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

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

and

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,

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}$),

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

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

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

where $\Omega_{\rm K}=\sqrt{GM_{\star}/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 $\Omega_{\rm 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 $\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.

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, $\Delta\Omega_{\rm eq}=\Omega_{\rm eq,-}-\Omega_{\rm eq,+}(>0)$, $\langle\cdots\rangle$ indicates the average over the $\phi-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 $\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. 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.

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 $\beta$, and we set

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)$.

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., 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 $\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).

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$:

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\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.

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:

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

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 (Brandenburg et al., 1995; Davis et al., 2010). 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).

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 ($\phi$) 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,

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 $[\alpha_{\rm 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 $[\alpha_{\rm R}]$ of the cylindrical cases are reliable. This is because it is not clear whether $\delta v_{\phi}$ defined as the deviation from the sub-Keprerian equilibrium rotation, $R\Omega_{\rm eq}$ (equation 7), is physically reasonable or not to estimate the Reynolds stress when the background rotation profile is deviated from $\Omega_{\rm eq}$ (Hawley, 2001, see also Section 4.4 for further discussion). For example, if we calculate $[\alpha_{\rm R}]$ by defining $\delta v_{\phi}$ as the deviation from the Keplerian rotation instead of $R\Omega_{\rm eq}$, then the Reynolds stresses obtained from both cylindrical cases are larger by $\sim 10^{-3}$ than those in Table 3.

Figure 2 presents the time evolution of the domain-averaged dimensionless Maxwell stress. 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.

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, 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 $H_{0}/R_{0}=0.1$ gives $\partial_{t}\ln[B_{R}^{2}]=0.028$, $\partial_{t}\ln[B_{\phi}^{2}]=0.21$, and $\partial_{t}\ln[B_{z}^{2}]=0.20$ (rotation^-1). These values are much larger than those expected from the small increase in $B_{i}^{2}$ during the time averaged period between $t=25$ and 300 rotations (Figure 2). Therefore, the net increase in $B_{i}^{2}$ by the positive $\partial_{t}\ln[B_{i}^{2}]$ 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 $\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., Brandenburg et al., 1995; Balbus & Hawley, 1998); $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 (blue arrows in Figure 6).

One can find that $\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,R]}$ is also positive. This is due to the transfer of angular momentum. The magnetic fields are distributed preferentially along the $-r\phi$ 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_{\phi}$ (red arrows in Figure 6), namely positive $\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,R]}$ is obtained. Moreover, $\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,R]}$ ($=+12.2$ and $+12.9$) in the cylindrical simulations are much larger than $\overline{[y\Rightarrow_{{}_{y}}~{}\,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 dominant component of $B_{\phi}$ ($B_{y}$) to $B_{z}$ also gives positive $\overline{[\phi\Rightarrow_{{}_{\phi}}~{}\,z]}$ ($\overline{[y\Rightarrow_{{}_{y}}~{}\,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_{{}_{x}}~{}\,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 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 $\propto v_{R}\delta v_{\phi}$ 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.

$\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 epicyclic frequency is a key for the positive $\overline{[z\Rightarrow_{{}_{{}_{R}}}~{}\,z]}$ in the cylindrical simulations, which will be discussed in Section 4.1.

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}]$ (top panel) and variation rates of the most dominate component of the toroidal field, $\partial\ln([B_{\phi}^{2}])/\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).

The middle panel of Figure 7 compares the variation rates of $B_{\phi}^{2}$ due to the radial compression, $[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]$ (red solid) and the winding by radial shear, $[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]$ (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 stays at a relatively low level.

The bottom panel of Figure 7 shows the sum of the compressible terms, $\left\{[\phi\Rightarrow_{{}_{z}}~{}\,\phi]+[\phi\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]\right\}$ (red solid) and the shearing terms, $\left\{[z\Rightarrow_{{}_{z}}~{}\,\phi]+[R\Rightarrow_{{}_{{}_{R}}}~{}\,\phi]\right\}$, (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 expansion.

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 maximum peak at $t\approx 91$ rotations.

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]$, 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 reaches more than 1.5 times of the time-averaged rate. In the descending phase of $90\lesssim t\lesssim 95$ rotations, although the transfer rate is reduced as the energy source, $B_{\phi}^{2}$, has already declined, it is maintained at a level comparable to the time-averaged rate, which further reduces $[B_{\phi}^{2}]$. A qualitatively similar tendency is seen before the later peak at $t\approx 105$ rotations. After this time, the transfer rate drops to a low level because the energy source, $B_{\phi}^{2}$ has sharply decreased earlier than $B_{z}^{2}$. The drop in $[\phi\Rightarrow_{{}_{\phi}}~{}\,z]$ reduces $B_{z}^{2}$, and the high magnetic activity phase is finally terminated. We can conclude that the enhanced leakage of the magnetic energy on the $R-\phi$ plane to the $z$ component is a trigger for the end of the high magnetic activity.

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 $f_{\rm K}=H^{2}/R^{2}=H_{0}^{2}/R_{0}^{2}=$ const. (equation 10), which is 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 epicyclic oscillations in the case with $H_{0}/R_{0}=0.1$, giving $\Omega_{\rm eq}=0.99\Omega_{\rm K}$, where an arbitrary amplitude is set to be $\xi_{0}=0.01R_{0}$. The figure clearly illustrates that the displacements in 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 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 $\phi$ and $z$ averaged (equation 16) and mass-weighted $v_{R}$,

where $\xi_{R}=0$ is set at the initial time, $t=70$ rotations, of the presented period. In the inactive period of $t\lesssim 83$ rotations, the gas accretes inward in a quasi-steady manner. In contrast, when the magnetic activity is enhanced after $t\gtrsim 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 $\xi_{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 triggering 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_{\phi,{\rm init}}$ (or $B_{z,{\rm 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_{\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 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_{\phi,{\rm init}}$, occupy 50% of the volume (the shaded regions in the right picture of Figure 11) 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 radially “neighboring” regions each other by $\pi$, the compressible amplification works most effectively as shown in Figure 11. 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

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 ($\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 (29) is rather conservative because the effect of the radial compression is regarded to be already effective well before the phase difference reaches $\pi$. 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\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 (29) 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. If the numerical resolution in the radial direction is insufficient and the difference in $\kappa$ between neighboring cells is too large, it appears that intermittent variability cannot be captured (Akatsuka & Suzuki 2023, in preparation).

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-\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 (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. 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 (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, $\alpha_{\nu}$, is assumed to depend separately on density and magnetic field:

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

where the detailed derivation is presented in Appendix B.

We note that in Riols & Lesur (2019) $\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 (Salvesen et al., 2016; Scepi et al., 2018) and global simulations (Mishra et al., 2020). There is a rough consensus of $q\approx 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 13 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 density¹¹1In 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 (30) 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 31).

Figure 14 shows scatter plots between $\alpha_{\rm 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_{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 31). 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 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.