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

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 marginally unstable condition regarding ``viscous-type'' instability.

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 ±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 (Ebrahimi & Blackman 2016).Although these problems are not inherent in global simulations (e.g.,

Suzuki
Armitage 1998; Machida et al. 2000;Flock et al. 2013;Suzuki & Inutsuka 2014;Béthune et al. 2017;Takasao et al. 2018;Zhu et al. 2020;Jacquemin-Ide et al. 2021), 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 ±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 (Brandenburg et al. 1996;Klahr & Bodenheimer 2003;Obergaulinger et al. 2009).Along the same lines, Suzuki et al. (2019, S19 hereafter) 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 epicyclic 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.

Basic setup
We solve ideal MHD equations in cylindrical coordinates, (R, ϕ, z), and where d dt and ∂ ∂t denote Lagrangian and Eulerian time derivatives, respectively; ρ, p, v and B are density, gas pressure, velocity, and magnetic field, G is the gravitational constant, M ⋆ 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 = ρc 2 s , where c s is isothermal sound speed that depends only on R and does not change with time.We employ the following radial dependence of temperature (∝ c 2 s ), 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, δv = v − RΩ eq φ, or (δv R , δv ϕ , δv z ) = (v R , v ϕ − RΩ eq , v z ), (7) if expressed in components, evaluated in the frame corotating with the equilibrium rotation frequency, Ω eq .If we ignore magnetic field and assume v R = v z = 0, the radial component of equation ( 2) yields under the equilibrium condition.Defining Ω eq ≡ v ϕ,eq /R, from the radial force balance equation ( 8) and assuming the radial profile of density, ρ ∝ R −µ , we obtain where Ω K = GM ⋆ /R 3 is the Keplerian rotation frequency and is a sub-Keplerian parameter (Adachi et al. 1976, S19).
Equations ( 1) -( 5) are solved with the second-order Godunov + CMoCCT method (Sano et al. 1999;Evans & Hawley 1988;Clarke 1996).The azimuthal advection due to Ω eq (equation 7) is handled with the "FARGO (Fast Advection in Rotating Gaseous Objects)"-type algorithm (Masset 2000;Benítez-Llambay & Masset 2016), which is an improved treatment from S19 to loosen the Courant-Friedrichs-Lewy condition (Courant et al. 1928) 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 ϕ boundaries.
We use the same simulation units as in S19: the physical variables are normalized by the units of R 0 = 1, ρ 0 = 1, and Ω K,0 = 1; the magnetic field is normalized by R 0 Ω K,0 √ 4πρ 0 , which cancels the √ 4π factor in the cgs-Gauss units.In the following sections, we conventionally define 2π/Ω K,0 as "one rotation time", which is slightly shorter than the actual time, 2π/Ω 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.

Shearing radial boundary condition
The radial boundary condition we are adopting is an extension of the shearing periodic condition in a local Cartesian box (Hawley et al. 1995).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, where the subscripts, − and +, stand for the inner and outer radial boundaries, ∆Ω eq = Ω eq,− − Ω eq,+ (> 0), ⟨• • • ⟩ indicates the average over the ϕ − z plane, and 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 δ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, δS is rescaled in order that the root-mean-squared (rms) amplitude, √ δS 2 , in the ghost cells averaged over the ϕ − z surface is matched to the surface averaged √ δS 2 in the adjacent cells inside the simulation domain.We have free parameters regarding this amplitude matching, which are determined to reproduce steady accretion structure during magnetically inactive periods (see Sections 4.4 and 4.5).The detailed numerical implementation for the boundary treatment is described in Appendix A.

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: and Equations ( 6), (13), and ( 14) guarantee a constant initial plasma β, and we set In order to trigger MRI, we add random velocity perturbations in v R and v ϕ to the equilibrium velocity distribution, (v R , v ϕ , v z ) = (0, RΩ eq , 0).

Simulation domain & resolution
The simulation domain is a local cylindrical region that covers (R − ∼ R + , ϕ − ∼ ϕ + , z − ∼ z + ).The azimuthal and vertical spacings, ∆ϕ and ∆z, of each grid cell are constant.The radial grid size, ∆R, is proportional to R.
We consider two cases with thin-disk conditions of c s,0 = 0.1R 0 Ω K,0 and 0.01R 0 Ω K,0 at R = R 0 .We define the scale height, H 0 = c s,0 /Ω K,0 , at R = R 0 .The same box size per H 0 is adopted in these two cases, 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., Bodo et al. 2008;Johansen et al. 2009;Shi et al. 2016, for discussion on the numerical domain size in Cartesian simulations).
The numerical resolution = 64/H 0 for the R and z directions is the same as that in S19.It is ≈ 61/H 0 in the ϕ direction, which is slightly higher than ≈ 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 ϕ , L z ), per H 0 and the same numerical resolution, (N x , N y , N z ) = (N R , N ϕ , N z ) (Table 1).

Spatial and temporal average
To analyze numerical data, we take various averages of simulated quantities.We express ⟨A⟩ as the ϕ-and z-integrated average of a variable, A, at R: We define as the volume-integrated average.For the average over the entire simulation domain, we simply express We write 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 ≳ 20 rotations.Thus, we take the temporal average from t 1 = 25 rotations to the end of the simulation at t 2 = t final = 300 rotations.

Evolution of magnetic energy
Taking the inner product of the induction equation (3) with B/4π, we have the equation for the evolution of magnetic energy: For numerical investigation, we rewrite equation ( 20) for B 2 R , B 2 ϕ , and B 2 z separately in a logarithmic derivative form: To examine in detail the evolution of magnetic field, we analyze the (i ⇒ 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 ⇒ k j], in Table 2. Similar numerical investigation on the induction equation has been done to understand dynamo-like magnetic evolution (Brandenburg et al. 1995;Davis et al. 2010).An essential extension from these works is that we separate the electromotive force, v × B, into the shearing and compressible parts (Section 3.2).

Time Evolution
We perform numerical simulations of the three cases in Table 1 until t = 300 rotation time.Figure 1 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 (ϕ) component with a moderate level of density fluctuations.Table 3 shows various quantities regarding magneto-turbulence averaged over the time and simulation domain.The 2nd column presents the time-and domain-averaged dimensionless Maxwell stress,

Label
Unabbreviated expression Table 2. Summary of the volume-integrated averages for the evolution of magnetic energy.See equations ( 21) -( 23) for the detail.
and the 3rd -7th columns represent different components of the magnetic energies and their ratios.While the three cases give similar magnetic properties, 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).
The 8th column of Table 3 shows the time-and volume-averaged Reynolds stress: We find that the cylindrical cases yield much lower [α R ] than the Cartesial case by nearly an order of magnitutde.We infer that this large difference between the cylindrical and Cartesian cases is due to the substantially different channels for the magnetic field amplification as will be discussed in Section 3.2.However, we do not suppose that the quantitative values of [α R ] of the cylindrical cases are reliable.This is because it is not clear whether δv ϕ defined as the deviation from the sub-Keprerian equilibrium rotation, RΩ eq (equation 7), is physically reasonable or not to estimate the Reynolds t (rotation) stress when the background rotation profile is deviated from Ω eq (Hawley 2001, see also Section 4.4 for further discussion).For example, if we calculate [α R ] by defining δv ϕ as the deviation from the Keplerian rotation instead of RΩ eq , then the Reynolds stresses obtained from both cylindrical cases are larger by ∼ 10 −3 than those in Table 3.
Figure 2 presents the time evolution of the domainaveraged dimensionless Maxwell stress.The cylindrical case with H 0 /R 0 = 0.1 (red solid) exhibits very large temporal variability; [α M ] reaches nearly 0.3 in the active period between 80 ≲ t ≲ 105 rotations, while [α 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 ≲ [α M ] ≲ 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 [α M ]'s of these three cases are similar (Table 3) because larger [α M ] in the active phases and smaller [α M ] in the inactive phases are obtained in higher-variability cases.

Evolution of magnetic field
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. Al-though some terms, e.g., [z ⇒ z ϕ] and [R ⇒ R ϕ] of the cylindrical case and [x ⇒ x z] and [y ⇒ 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 ⇒ 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, when the steady-state condition is satisfied, the sum of numerical values from all four arrows entering each B i is 0. Figure 4 indicates that, while this is approximately satisfied in these three cases, the sums yield tiny positive values.For example, the cylindrical case with ).These values are much larger than those expected from the small increase in B 2 i during the time averaged period between t = 25 and 300 rotations (Figure 2).Therefore, the net increase in B 2 i by the positive ∂ t ln[B 2 i ] is considered to be compensated by the numerical dissipation of the magnetic fields in a sub-grid scale; in other words, this analysis can quantify the numerical dissipation of magnetic fields in ideal MHD simulations (see also Section 4.2 for discussion on the magnetic diffusion in the ideal MHD condition).
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 epicyclic 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).
We examine each [i ⇒ k j] term in more detail, particularly focusing on the similarities and differences between the cylindrical and Cartesian simulations.Let us begin with the [i ⇒ 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) 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.The 8th column presents Reynolds stress. -10 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 H0/R0 = 0.1 (left) and the Cartesian case (right).What are shown from top to bottom are ), and ∂ ln [Bz ] 2 ∂t in units of (rotation time) −1 .Smoothing with spline interpolation is applied to the plotted lines to display the long-time evolution.By this treatment, the original numerical data are averaged over ≈ 3 rotation times.
of the magnetic fields are similar at least in a qualitative sense.Positive [z ⇒ z R] and [R ⇒ R ϕ] reflect the standard pathway triggered by MRI (e.g., Brandenburg et al. 1995;Balbus & Hawley 1998); B R is generated from B z by the MRI, and eventually, B ϕ is amplified from the fluctuating B R by the winding due to the radial differential rotation (blue arrows in Figure 6).
One can find that [ϕ ⇒ ϕ R] is also positive.This is due to the transfer of angular momentum.The magnetic fields are distributed preferentially along the −rϕ direction by radial differential rotation as shown in Figure 1.In stronger-field regions, the angular momentum is more effectively transported outwardly along the field lines.This causes the inner (outer) fluid parcel to move inward (outward), whereas the volume averaged bulk flow keeps the steady accretion structure (see Section 4.4).Therefore, B R is generated from B ϕ (red arrows in Figure 6), namely positive [ϕ ⇒ ϕ R] is obtained.Moreover, [ϕ ⇒ ϕ R] (= +12.2 and +12.9) in the cylindrical simulations are much larger than [y ⇒ y x](= +0.468) in the Cartesian simulation.The difference is probably be-cause the angular momentum cannot be defined in the Cartesian setup owing to the assumed symmetry with respect to the ±x directions (See Section 1).The leakage from the dominant component of B ϕ (B y ) to B z also gives positive [ϕ ⇒ ϕ z] ([y ⇒ y z]) with the similar tendency concerning the difference between the cylindrical and Cartesian cases.
Next, we inspect the compressible terms, [i ⇒ j i], which correspond to the round arrows around each B i in Fig- ure 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 [i ⇒ j i] terms are different in the cylindrical and Cartesian simulations.Most remarkable difference is found in the compressible amplification of B ϕ along R; [ϕ ⇒ R ϕ] 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 [y ⇒ x y](= −1.05) in the Cartesian case.Surprisingly, the amplification rate by [ϕ ⇒ R ϕ] is slightly larger than that of [R ⇒ R ϕ], the winding term by radial differential rotation.In other words, the con-  , per rotation time averaged over 25 -300 rotations (See Table 2).The positive (negative) rates with ≥ 1(≤ −1) are shown by red solid (blue dashed) arrows.We note that the sum of the values on the arrows entering each Bi is = 0 when the time-steady condition is achieved.
tribution from the compressible flows dominates that from the shear flows in the amplification of B ϕ 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 epicyclic frequency, which cannot be considered in the Cartesian setup.The relative importance of the compression against the shear is a possible and probably plausible reason for the small Reynolds stress in the cylindrical simulations (Section 3.1).From geometrical considerations, the Reynolds stress ∝ v R δv ϕ is expected to be generated by shearing motions.If the compressible effect dominates, the sheared stress will be perturbed or interrupted, which will reduce the Reynolds stress.
[R ⇒ ϕ R] is also different between the cylindrical and Cartesian simulations.Although [x ⇒ y x] is nearly 0 in the Cartesian simulation, [R ⇒ ϕ R] takes a relatively large negative value in the cylindrical simulations, which partially cancels positive [ϕ ⇒ ϕ 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 epicyclic frequency is a key for the positive [z ⇒ R z] in the cylindrical simulations, which will be discussed in Section 4.1.

Onset of enhanced magnetic activity
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 [i ⇒ k j] terms.In this subsection, we scrutinize the rising phase of magnetic activity Figure 7 shows the time evolution of [α M ] (top panel) and variation rates of the most dominate component of the toroidal field, ∂ ln([B 2 ϕ ])/∂t, (middle and bottom panels) during t = 65 − 97 rotation time, which covers the rising and saturated phases of the largest peak in [α M ] (Figure 2).
The middle panel of Figure 7 compares the variation rates of B 2 ϕ due to the radial compression, . An example for the amplification of a magnetic field by shearing terms.The magnetic field is distributed with −BRB ϕ > 0, which is the usual situation in the system with inner-fast differential rotation.B ϕ increases by radial differential rotation of v ϕ , which is represented by (R ⇒ R ϕ) (blue).BR 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 (ϕ ⇒ ϕ R) (red).

solid) and the winding by radial shear, [R ⇒
R ϕ] (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 2 ϕ in the rising phase of the magnetic activity when the radial shear term stays at a relatively low level.
The bottom panel of Figure 7 shows the sum of the compressible terms, [ϕ R ϕ] , (black dashed).As shown in Figures 3 & 4, the vertical compression term, [ϕ ⇒ z ϕ], 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 ≲ 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 ⇒ x y] and [y ⇒ z y] are always negative, namely the compressible terms work as the decay of B 2 y by divergent expansion.
After t ≳ 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 maximum peak at t ≈ 91 rotations.

Termination of magnetic activity
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 2  i ] is normalized by the time-averaged [B 2 i ].One can recognize that [B 2 z ]/[B 2 z ] (red dashed line) reaches the highest peaks at t ≈ 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 ϕ components.Moreover, in the subsequent declining period, ϕ ] (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 ϕ components rapidly leaks to the z component.
In order to verify this speculation, we inspect the transfer rate, [ϕ ⇒ ϕ z], from B 2 ϕ to B 2 z in the bottom panel of Figure 8.It monotonically increases in the rising phase of t ≤ 90 rotations as [B 2  ϕ ] increases, and it stays at a high level, which reaches more than 1.5 times of the time-averaged rate.In the descending phase of 90 ≲ t ≲ 95 rotations, although the transfer rate is reduced as the energy source, B 2 ϕ , has already declined, it is maintained at a level comparable to the time-averaged rate, which further reduces [B 2 ϕ ].A qualitatively similar tendency is seen before the later peak at t ≈ 105 rotations.After this time, the transfer rate drops to a low level because the energy source, B 2 ϕ has sharply decreased earlier than B 2 z .The drop in [ϕ ⇒ ϕ z] reduces B 2 z , and the high magnetic activity phase is finally terminated.We can conclude that the enhanced leakage of the magnetic energy on the R − ϕ plane to the z component is a trigger for the end of the high magnetic activity.

Radial dependence of epicyclic frequency
A significant difference of the cylindrical approach from the Cartesian one is the presence of the radial variation in epicyclic frequency.For the equilibrium rotation profile, equation ( 9), the epicyclic frequency is written as where we used (equation 10), which is practically satisfied in our simulations.While in the Cartesian box obviously κ is spatially constant, in the cylindrical box κ is radially dependent, κ ∝ R −3/2 .Figure 9 demonstrates radial displacements, ξ R = ξ 0 sin(κt), arising from epicyclic oscillations in the case with H 0 /R 0 = 0.1, giving Ω eq = 0.99Ω K , where an arbitrary amplitude is set to be ξ 0 = 0.01R 0 .The figure clearly illustrates that the displacements in neighboring radial locations gradually become out of phase with time because ∂κ ∂R ̸ = 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 fields.
Figure 10 presents the radial displacements of fluid in the cylindrical case with H 0 /R 0 = 0.1 in a time-distance diagram between 70 and 90 rotations.Here the radial displacement of each radial cell is calculated from the ϕ and z averaged (equation 16) and mass-weighted v R , where ξ R = 0 is set at the initial time, t = 70 rotations, of the presented period.In the inactive period of t ≲ 83 rotations, the gas accretes inward in a quasisteady manner.In contrast, when the magnetic activity is enhanced after t ≳ 83 rotations, the outer gas moves outward while the inner gas continues to accrete; the gas in the domain radially expands, which will be discussed in Section 4.4.
Compared with the ideal setting of ordered epicyclic oscillations in Figure 9, radial oscillations in the numerical simulation (Figure 10) are more randomly excited by turbulence in a stochastic fashion.Accordingly, the picture of the phase mixing of initially ordered epicyclic oscillations based on Figure 9 has to be generalized.Indeed, one can recognize randomly and ubiquitously formed concentrations of ξ R trajectories, which are a characteristic feature of converging regions, as seen in the ideal setup (Figure 9).Furthermore, these converging regions are more frequently seen in the elevating phase of the magnetic activity during t = 80 − 90 rotation time, which is consistent with the argument on the importance of the compressible amplification in trig- gering the high magnetic activity (Figure 7 and Section 3.3).
Let us suppose a simple situation in which toroidal (or vertical) magnetic field with constant strength, B ϕ,init (or B z,init ), is initially distributed (left of Figure 11).If we pick out the radial compression term from equation ( 22), the corresponding term of the induction equation is This is essentially the continuity equation for B ϕ , whereas the geometrical curvature term, 1 compared to the mass continuity equation (1).Therefore, B ϕ is compressed (rarefied) in the converging (diverging) regions (right of Figure 11).If we compare the left and right pictures of Figure 11, 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 ϕ,init , occupy 50% of the volume (the shaded regions in the right picture of Figure 11) so that [B ].The consideration based on this simple model indicates that the radial compression term should systematically increase B 2 ϕ , which does not occur in the Cartesian setup with the spatially constant κ.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, [ϕ ⇒ R ϕ], takes the relatively large positive values in the cylindrical simulations, which is in clear contrast to the negative [y ⇒ 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 , [z ⇒ R z], is positive in the cylindrical cases, which is also in contrast to the negative corresponding term, [z ⇒ x z], in the Cartesian case.
The abovementioned argument indicates that the presence of ∂κ ∂R ̸ = 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 ⇒ 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 radially "neighboring" regions each other by π, the compressible amplification works most effectively as shown in Figure 11.We can estimate the time, t = τ comp , to reach the phase difference of ∆ϕ κ = π from the uniform oscillation at t = 0 via ∆κ 0 t = ∆ϕ κ , where the deviation of κ 0 is derived from ∆κ = | ∂κ ∂R 0 |∆R 0 by using the radial spacing between the neighboring regions, ∆R 0 .Since | ∂κ ∂R 0 | = 3Ωeq,0 2R , we have where the radial spacing is normalized by H 0 for a typical scale of the system.This equation shows that the epicyclic oscillation becomes completely out of phase (∆ϕ κ = π) at radial spacing of ∆R = H 0 when τ comp ≈ 3 rotations for the cylindrical case with H 0 /R 0 = 0.1.For H 0 /R 0 = 0.01, longer τ comp ≈ 30 rotations are required.However, we would like to note that the estimate based on equation ( 29) is rather conservative because the effect of the radial compression is regarded to be already effective well before the phase difference reaches π.Additionally, oscillations are expected to be excited rather randomly by turbulence, as discussed with Figure 10.Thus, even in this nearly Cartesian case, the compressible amplification is already working at t ≈ 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 ( 29) is that, as time goes on, the compressible amplification is going to work for smaller ∆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 ∂κ ∂R ̸ = 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.If the numerical resolution in the radial direction is insufficient and the difference in κ between neighboring cells is too large, it appears that intermittent variability cannot be captured (Akatsuka & Suzuki 2023, in preparation).

Substructures and viscous-type instability
Figure 11 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., Suriano et al. 2019;Cui & Bai 2021, see Lesur et al. (2022) for recent review).
Figure 12 shows a two-dimensional (2D) r − ϕ slice of density, ρ/ρ 0 (left), "magnetization", B 2 z /4πp (middle), defined as the inverse of β z (c.f., equation 15), and Maxwell stress, α M (right) on z = 0 at t = 82.4rotation 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 weakerfield 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 (Steinacker & Papaloizou 2002) and discussed from a viewpoint of viscous-type instability (Lightman & Eardley 1974) generalized by including the physical properties of the MRI (Hawley 2001).
One of the key ingredients of the viscous-type instability in the MHD framework is magnetic diffusion.Al-Suzuki c s,0 0.8 0.9 .1 at z=0, t=82.4(rot.)though 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 (Lazarian & Vishniac 1999;Lazarian et al. 2020).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 Riols & Lesur (2019), 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 Riols & Lesur (2019) is that dimensionless viscosity, α ν , is assumed to depend separately on density and magnetic field: where we can practically assume α ν ≈ α 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 where the detailed derivation is presented in Appendix B.
We note that in Riols & Lesur (2019) α ν is assumed to depend on magnetization, B 2 z /4πρc 2 s , namely q = q B = q ρ is imposed, and they derived the instability condition as q > 1.The relation between time- and volume-averaged α M and B 2 z /4πρc 2 s have been extensively examined with Cartesian shearing box simulations (Salvesen et al. 2016;Scepi et al. 2018) and global simulations (Mishra et al. 2020).There is a rough consensus of q ≈ 0.5 − 0.7, which does not satisfy the instability condition in Riols & Lesur (2019).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 14 shows scatter plots between α M (vertical axis) and various quantities (horizontal axis) of each grid point displayed in Figure 12.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 2 z (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-α M and lower-density (upper right) side.Moreover, the slope is slightly steeper than q ρ = 1 and is in the unstable regime (equation 31).We examined scatter plots for ρ-α M in other time frames and found that q ρ ≳ 1 is kept in most of the time.Therefore, we can interpret that the substructures seen in our simulations (Figure 12) 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.

Intermittency
Cartesian shearing box simulations for MRI often show quasi periodicities with t ∼ 10 rotation time (Davis et al. 2010;Gressel 2010;Guan & Gammie 2011), which is related to the recurrent growth and breakup of largescale channel-mode flows (Gogichaishvili et al. 2018).While quasi-periodic magnetic activity is more clearly seen in vertically stratified simulations, being frequently associated with mass outflows (e.g., Suzuki & Inutsuka Suzuki 2009;Suzuki et al. 2010;Wissing et al. 2022), similar periodicity is also observed in unstratified simulations (Sano & Inutsuka 2001;Shi et al. 2016).Our Cartesian case is also showing mild time variation in [α M ] (black dashed line in Figure 2).
On the other hand, the duration ∼ 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 ∼ 10 rotations described above.A speculative mechanism for this weak and long-term periodicity is related to the compressible effect arising from the radial variation of κ (Section 4.1).If we utilize equation ( 29), the time to reach the phase difference, ∆ϕ κ = 2π, between neighboring radial regions can be estimated as a possible source of the periodicity.If we adopt ∆R ≈ 0.2R 0 , referring to the radial width of a ring (or a gas) in Figure 12, we obtain t ≈ 33 rotation time, which is roughly consistent with the observed cycle of the enhanced magnetic activitye.On the other hand, this estimate based on the initially ordered epicyclic oscillations (Figure 9 Section 4.1) may not be relevant, because oscillations driven constantly and randomly by turbulence interact each other as illustrated in Figure 10.

Accretion structure in active and inactive phases
We compare the radial dependences of various physical quantities during magnetically active and inactive phases in Figure 15.In order to see the difference between the active and inactive periods, we take the temporal averages for [α M ] < 0.13 (left) and [α 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 15 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 ϕ-z plane is 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 appears to be bro-ken.This is mainly because the increased magnetic pressure pushes gas inward (outward) from the inner (outer) boundary as discussed below.However, the averaged radial velocity, ⟨v R ⟩, near the inner and outer boundaries is only a few % of the sound velocity2 .This gradual expanding velocity is much smaller than the fluctuating velocity, which is nearly an order of the sound speed in the active phase.
The middle panels of Figure 15 show each component of the rms magnetic field strength, ⟨B 2 i ⟩, in comparison to the initial B z,init .The left panel indicates that the vertical magnetic field strength, ⟨B 2 z ⟩, averaged over the inactive periods increases about 10 times from the initial value with almost keeping the initial profile, ∝ 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), although the field strength is larger as expected, the ratios between different components is not so different from those in the inactive phases.The toroidal component exhibits a weakly concave down profile.Then, the direction of the magnetic pressure gradient is inward near the inner boundary, while it is directed outward in the intermediate and outer regions.These cause the slow expanding flows seen in the top right panel.
In the bottom panel of Figure 15, density and the inverse of plasma β, 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, ⟨β⟩ −1 exhibits a weak radial dependence, which can be more easily seen in the active phase.The radial distribution of the gas pressure is also modified according to the altered density profile.This also leads to the alternation of the time-averaged rotation profile from Ω eq , which affects the estimate of the Reynolds stress in equation 25 as discussed in Section 3.1.
We present triangle diagrams for the active and inactive periods in Figure 16.One can find that these two phases exhibit quite similar tendencies.Careful comparison shows that most of the variation rates, , are slightly smaller in the active phase.This is because The middle panels compare rms BR (magenta dash-dotted), B ϕ (blue dashed), and Bz (black solid).The initial vertical field strength, Bz,init, is multiplied by a factor of 10 (red dotted), in order to compare with 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 β = 8π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.
these numerical values are normalized by larger [B 2  i ]; if we compare the change of the magnetic energies, ∂t , the active phase yields larger values.

Boundary effect
There is a freedom to set the fluctuation amplitudes at the radial boundaries (Section 2.2).We determined the parameters to control these boundary amplitudes in order that the steady accretion structure is reproduced when we take the time average over the inactive periods (right panels of Figure 15), which is described in detail in Appendix A. On the other hand, as discussed above, the steady accretion is not achieved for the radial structure averaged over the active periods.This raises doubts whether the obtained results may be affected by the boundary treatment.Hence, we perform cylindrical simulations with H 0 /R 0 = 0.1 for different values of the amplitude parameters in Appendix A.
Although these different cases show different radial flow structures and yield stochastically enhanced magnetic activities at different times, we obtain similar properties of the intermittency; magnetic bursts occasionally appear in low-activity phases that occupy most of the simulation time.Consequently, these different cases give very similar time averaged Maxwell stresses and magnetic field strengths (Figures 19 & 20 and Table 4 in Appendix A).Furthermore, the obtained variation rates of the magnetic energies, [i ⇒ k j], in the triangle diagrams are also very similar each other (Figure 21).
So far, we have examined the [i ⇒ k j] terms integrated over the whole simulation domain.However, these terms may be spatially dependent particularly near the the radial boundaries.To inspect the boundary effect, we compare all [i ⇒ 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 ) to those averaged over Suzuki H 0 /R 0 = 0.1, Active periods  the entire simulation domain shown in Figure 4. Two winding terms, [R ⇒ R ϕ] and [z ⇒ z ϕ], in the case with H 0 /R 0 = 0.1 are slightly affected; specifically, the averages in the narrow region are [R ⇒ R ϕ] = 4.12 and [z ⇒ z ϕ] = −4.41,which are both reduced by ≈ 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 contribution from the active periods causing the mismatched radial boundaries 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.

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 epicyclic frequency, κ, 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 κ (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 κ (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, α M ∝ ρ −qρ , with q ρ ≳ 1 is obtained (Section 4.2).We revisited the viscous-type instability by considering the dependence of α M separately on density and vertical magnetic field, and found that the instability condition is q ρ > 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 κ; 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.
The author thanks the anonymous referee for valuable comments to the original draft.Numerical computations were carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan, and Yukawa-21 (Dell Pow-erEdge 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.(famp,+, famp,−) (1.   Therefore, larger f amp,+ (f amp,− ) tends to induce inward (outward) flows.We are adopting f amp,+ = 1.02 and f amp,− = 0.93 for the cylindrical case with H 0 /R 0 = 0.1 and f amp,+ = 0.94 and f amp,− = 0.93 for H 0 /R 0 = 0.01.We demonstrate how different choices of f amp,± affect the results of the case with H 0 /R 0 = 0.1 in Figures 19 -21 and Table 4.We perform simulations with (f amp,+ , f amp,− ) = (1.06,0.96) and (1, 1) in addition to the original case with (f amp,+ , f amp,− ) = (1.02,0.93).From Figure 19 we find that, although the timings of magnetic enhancements are different, these three cases exhibit similar intermittent properties.
The top panels of Figure 20 indicates that the radial accretion profile is affected by the values of f amp,± .In the original case (left), although the steady accretion profile is achieved in the inactive periods, the radial flow shows a gentle expanding trend by the contribution from the active periods (Figure 15).On the other hand, the case with (f amp,+ , f amp,− ) = (1.06,0.96) yields a steady accretion feature4 , because the expanding flow is confined by a small increase in both f amp,+ and f amp,− as discussed above.At the same time, one can also find from the middle and bottom panels of Figure 20 and Table 4 that these two cases give very similar time-averaged magnetic fields.
The case with f amp,± = 1 (right panels of Figure 20) exhibits an expanding trend with positive v R in the outer region, which is expected from the smaller f amp,+ and larger f amp,− than those in the original case.The magnetic properties of this case are similar to those of the other two cases, although the time-averaged field strength is slightly larger (Table 4) by the larger contribution of the high magnetic activity periods (Figure 19).
The triangle diagrams of these three cases show very similar shearing and compressible amplification rates (Figure 4); in particular, the critical importance of the radial compression in the amplification of B ϕ is universally attained.In summary, although f amp,± controls the accretion profile, the intermittency of the magnetic activity and the timeaveraged field strength are not significantly affected, provided that f amp,± ≈ 1 is employed.

B. LINEAR PERTURBATION ANALYSES FOR VISCOUS-TYPE INSTABILITY
We conduct linear perturbation analyses on the 2D r − ϕ plane under the axisymmetric approximation.In the momentum equation (2), we do not explicitly include magnetic fields but consider them through α prescription (Shakura & Sunyaev 1973).Then, the radial and angular momentum equations can be written as in Eulerian forms, where (w R , w ϕ ) = (v R , v ϕ − RΩ K ) is the velocity deviation from Keplerian rotation and α ν is dimensionless viscosity, equation (30).We note that, in equation (B5), the inward gravity, −ρ GM R 2 , and the centrifugal force, ρRΩ 2 K , should be present but they are exactly canceled out each other.We adopt the radial profiles of c 2 s ∝ R −1 (equation 6), ρ ∝ R −1 , (equation 13) and α ν ∝ R 1/2 (S19) to match our simulation setting.Then, from the unperturbed state (∂ t = 0) of equations ( B5) and (B6), we obtain which is consistent with equation ( 9), and α ν is assumed to depend separately on ρ and B z (equation 30), and thus, we have Here, density is calculated from the continuity equation ( 1) and vertical magnetic field from the induction equation (3) with resistivity, η, explicitly included: where we assumed v z = 0. We also adopt the α prescription for η with the same dependence on ρ and B z , We note that η ∝ R ensures that equation (B10) has an unperturbed-state solution of B z ∝ R −1 , which is consistent with equation ( 14) (see also S19).
We expand ρ, w R , w ϕ , and B z with the radial dependences that match our simulation setting: We apply these expansions to equations (1), (B5), (B6), and (B10) and further apply plane-wave expansion, δ ∝ exp(iωt − ikR), (B16) to the first-order variables.Then, the corresponding four equations are

Suzuki
From equations (B17) -(B20), we can derive a fourth-order equation with respect to ω.Here, we ignore the solutions of ω 2 ≈ Ω 2 K + k 2 c 2 s describing sound waves in shearing systems, and focus on secular modes by assuming ω 2 ≪ Ω 2 K .In addition, we also assume H/R ≪ 1 and ignore the terms ∝ 1/R, 1/R 2 , • • • arising from the curvature of cylindrical coordinates.Then, we finally obtain where we omitted the subscripts, 0, and the superscripts, ′, because we ignored the R-dependent terms.Substituting equations (B8) and (B11) into equation (B21), we can derive is the magnetic Prandtl number and we once again disregarded the terms ∝ 1/R that originally involved w R .As shown in Figure 14, q ρ ≈ 1, and then, in the square-root of equation (B22) 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 B16), one can find a growing mode, if it is present, only for the negative sign of equation ( B22).The growth rate can be calculated as ) 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 Ω K α ν for Pm ≪ 1 and Ω K α η for Pm ≫ 1.In the former case, an unstable mode is present if q ρ > 1 (equation 31), 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 ρ > 1 because (1 − q ρ + 2q B ) > 0 is probably satisfied in usual situations (Figure 14).

Figure 2 .
Figure 2. Comparison of the time evolutions of [αM] averaged over the simulation box for the cylindrical cases with H0/R0 = 0.1 (red solid) and 0.01 (blue dotted) and the Cartesian case (black dashed).

Figure 4 .
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 H0/R0 = 0.1 and 0.01, and the Cartesian case, respectively.The numerical values indicate the rate of the change of magnetic energy, ∂ ln [B 2 i ] ∂t

Figure 5 .
Figure 5. Left: The shearing term, [i ⇒ i j], stands for the generation (or decay) of Bj from Bi by the i derivative of vj.Right: The compressible term, [i ⇒ j i], represents the convergent amplification (or divergent disperse) of Bi by the j derivative of vj.

Figure 7 .
Figure 7. Time evolution of Maxwell stress (top) and variation rates of [B 2ϕ ] (middle and bottom) averaged over the domain in the rising and saturated phases of the most active period in the cylindrical simulation with H0/R0 = 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.

Figure 8 .Figure 9 .
Figure 8.Time evolution of magnetic energy (top) and the variation rate of B 2 z by shearing from ϕ to z (bottom) averaged over the domain in the declining phase of the same active period shown in Figure 7.What are shown in the top panel are the three components of [B 2 i ] with i = R (black dotted), ϕ (blue solid), and z (red dashed), normalized by the respective time-averaged values, [B 2 i ], during t = 25 − 300 rotations.[B 2 i ]/[B 2 i ] = 1 is plotted by the black thin solid line.In the bottom panel, the time average, [ϕ ⇒ ϕ z] = +9.71 is also plotted by the black dotted line.

Figure 10 .Figure 11 .
Figure 10.Radial distance (horizontal) vs. time (vertical) diagram for radial displacements taken from the numerical data of the cylindrical case with H0/R0 = 0.1.

Figure 14 .
Figure 14.Scatter plots of dimensionless Maxwell stress, αM, against B 2z /4πρc 2 s (left), (ρ/ρ0) −1 (middle), and B 2 z /4π (right) at t = 82.4rotations.The data are averaged over 2 radial and 16 azimuthal grid cells of the 2D plots in Figure12.The slope of y ∝ x is overlaid in the left and middle panels by the red solid line, and y ∝ x 0.5 is in the left panel by the blue dashed line.

Figure 15 .
Figure 15.Comparison of time-averaged radial profiles of various physical quantities of the case with H0/R0 = 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 BR (magenta dash-dotted), B ϕ (blue dashed), and Bz (black solid).The initial vertical field strength, Bz,init, is multiplied by a factor of 10 (red dotted), in order to compare with 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 β = 8π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.

Figure 17 .
Figure 17.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 −,d (R +,d ) is the radial location of the innermost (outermost) cells in the simulation domain, which are illustrated by gray shade.R−,g (R+,g) is the radial location of the inner (outer) ghost cells.The physical quantities of the ghost cell at R = R±,g (red) are determined from those of the corresponding sheared cells at R = R ∓,d (blue) and the radially adjacent cell(s) at R = R ±,d (green).See text for the detail.

Figure 21 .
Figure 21.Same as Figure 4 but for cases with different famp,± for H0/R0 = 0.1.The left, middle, and right triangles correspond to those in Figure 20.Note that the left triangle is the same as the top right triangle in Figure 15.

Table 1 .
Summary of the simulation domains and resolutions.