Revisiting linear dynamics of non-axisymmetric perturbations in weakly magnetized accretion discs

We investigate linear dynamics of non-axisymmetric perturbations in incompressible, vertically stratified Keplerian discs with a weak vertical magnetic field in the shearing box approximation. Perturbations are decomposed into shearing waves whose evolution is followed via numerical integration of the linearized ideal MHD equations. There are two basic modes in the system -- inertia-gravity waves and magnetic mode, which displays the magnetorotational instability (MRI). As distinct from previous studies, we introduce `eigenvariables' characterizing each (counter-propagating) component of the inertia-gravity and magnetic modes, which are governed by a set of four first order coupled ordinary differential equations. This allowed us to identify a new process of linear coupling of the two above non-axisymmetric modes due to the disc's differential rotation. We did a comparative analysis of the dynamics of non-axisymmetric and axisymmetric magnetic mode perturbations. It is shown that the growth of optimal and close-to-optimal non-axisymmetric harmonics of this mode, having transient nature, can prevail over the exponential growth of axisymmetric ones (i.e., over the axisymmetric MRI) during dynamical time. A possible implication of this result for axisymmetric channel solutions is discussed. Specifically, the formation of the channel may be affected/impeded by non-axisymmetric modes already at the early linear stage leading to its untimely disruption -- the outcome strongly depends on the amplitude and spectrum of initial perturbation. So, this competition may result in an uncertainty in the magnetic mode's non-linear dynamics. It is also shown that a maximum non-axisymmetric growth is at vertical wavelengths close to the scale-height for which compressibility effects are important. This indirectly suggests that compressibility plays a role in the dynamics of the non-axisymmetric MRI.


INTRODUCTION
The interplay of magnetic field and differential rotation in accretion discs offers a realistic hope of explaining energetic processesturbulence and angular momentum transport -in this kind of astrophysical flows. Despite more than twenty years of research in this field, after the seminal paper by Balbus & Hawley (1991) in which the magnetorotational instability (MRI) was proposed as the key element for angular momentum transport in accretion discs, many questions are still open. A fully consistent dynamical picture is still ⋆ E-mail: george.mamatsashvili@tsu.ge missing, mainly because of the complex nature of the ingredients entering into dynamical phenomena in magnetized discs.
In MRI studies, it is customary to consider two different kinds of configurations: one in which there is a non-zero net magnetic flux threading the disc and the second in which there is no net flux. In these two cases open problems and dynamical processes are somewhat different, in fact, in the first case we have a well defined linear instability, while in the second case the MRI must set in as a subcritical instability and dynamo action must be effective in order to sustain the field needed for the MRI to operate. In this paper, we will focus our attention on the net flux case: in this case in order to understand the dynamics one has to consider the relative importance and the interplay of several different processes: (I) exponential instability of axisymmetric perturbations (axisymmetric magnetorotational instability, Balbus & Hawley 1991;Goodman & Xu 1994); (II) transient growth of non-axisymmetric perturbations caused by the shear/non-normality of the disc's differential rotation (Balbus & Hawley 1992;Tagger et al. 1992;Brandenburg & Dintrans 2006;Johnson 2007;Pessah & Chan 2012, transient dynamics includes: Coriolis forces, the basic disc flow shear, magnetic tension forces); (III) non-linear interaction of axisymmetric perturbations (channel flow) with non-axisymmetric (parasitic) ones (Goodman & Xu 1994;Pessah & Chan 2008;Pessah & Goodman 2009;Longaretti & Lesur 2010;Pessah 2010); (IV) non-linear interaction between non-axisymmetric perturbations and the onset and character of turbulence (e.g., Hawley et al. 1995;Brandenburg et al. 1995;Lesur & Longaretti 2011;Bodo et al. 2008;Guan et al. 2009;Davis et al. 2010;Bodo et al. 2011); (V) linear coupling of different non-axisymmetric modes/branches (like the transient growth, this coupling is also induced by the non-normality of the disc shear flow, see e.g., Bodo et al. 2005); (VI) viscous and diffusive effects (e.g., Sano et al. 1998;Sano & Stone 2002;Fleming & Stone 2003;Pessah & Chan 2008;Longaretti & Lesur 2010;Pessah & Goodman 2009).
Certainly, consistent numerical simulations are able to cover the whole set of the above processes, however in simulations all the processes are intermixed and, in order to understand the dynamics, one has to find a way to disentangle them. In analytical studies, one is compelled to restrict attention only to some of the above processes. For instance, the basic scenario in the above-mentioned several recent papers focusing basically on I, III and VI processes is the following: one considers the exponential growth of axisymmetric perturbations, the fastest growing of which gives rise, in the non-linear regime, to the so-called channel solution. In essence, this is process I. Afterwards, on top of this channel flow, a secondary instability develops and non-axisymmetric, so-called parasitic, modes grow. This analysis should mimic process III, albeit with some simplifications, like neglecting the basic flow shear and Coriolis force, due to the complexity of the secondary instability. Dissipative effects (process VI) have been also considered, but processes II, IV and V have been in general neglected, however, there are cases where they may be appreciable and even decisive. The goal of this paper is in fact the analysis of the importance of processes II and V, through a study of the linear dynamics of perturbation in an incompressible stratified disc flow with a vertical magnetic field. A detailed numerical investigation of the non-linear process IV is out of the scope of the present study and will be covered in a subsequent paper. We would like to draw particular attention to process Va linear coupling of non-axisymmetric modes. This phenomenon, caused by shear flow non-normality, was revealed in the 1990's and is widely accepted by the hydrodynamical and geophysical/meteorological communities (e.g., Chagelishvili et al. 1997 Kalashnik et al. 2006). It is even more significant in astrophysical disc flows, where background differential rotation (i.e., flow non-normality) is strong and, hence, the linear mode coupling is an inevitable occurrence. Examples of the studies taking into account the coupling in astrophysical discs are Tevzadze et al. (2003); Bodo et al. (2005); Tevzadze et al. (2008); Mamatsashvili & Chagelishvili (2007); Heinemann & Papaloizou (2009a); Tevzadze et al. (2010); Paardekooper et al. (2010); Mamatsashvili & Rice (2011). This paper also demonstrates yet another manifestation of the coupling phenomenon in magnetized discs.
In the present paper, we primarily perform a comparative analysis of the exponential instability of axisymmetric perturbations and transient growth of non-axisymmetric ones, i.e., a quantitative comparison of the processes I and II. We also focus on the linear coupling of various non-axisymmetric physical modes of perturbations existing in the disc flow, i.e., on the process V. As a matter of fact, in complex/unstable systems, the dynamics of perturbation modes is driven by various physical phenomena. For instance, in the considered here magnetized incompressible disc flow, there are two types of non-axisymmetric modes: 1. magnetic mode that directly taps into the energy store of disc flow shear and experiences transient growth and 2. inertia-gravity waves that only weakly directly exchange energy with the disc flow (Balbus 2003;Ogilvie 1998). However, as we demonstrate here, the latter mode is coupled linearly with the magnetic mode and draws energy from it. As a result, inertia-gravity waves can also become energetically intensive and affect dynamical processes, particularly, angular momentum transport.
The disc model we consider is similar to that of Balbus & Hawley (1992); Goodman & Xu (1994); , i.e., it represents a local patch of incompressible, vertically stratified ionized Keplerian disc penetrated by spatially uniform background magnetic field. We work in terms of the shearing waves (Kelvin modes) that are a natural basis for the shearing box models of discs for describing the dynamics of waves (e.g., Goldreich & Tremaine 1978) and vortices (e.g., Lominadze et al. 1988;Chagelishvili et al. 2003;Johnson & Gammie 2005) as well as their linear coupling (e.g., Bodo et al. 2005;Mamatsashvili & Chagelishvili 2007;Heinemann & Papaloizou 2009a). As mentioned above, the considered system contains two distinct physical types/modes of perturbations: magnetic mode and inertia-gravity waves. In turn, each of these modes consists of two components. The magnetic mode has two counter propagating components that exhibit growth due to the MRI, but during a limited time interval. The second, inertia-gravity mode also consists of two counter propagating waves that can grow due to the coupling with the magnetic mode. To analyse the dynamics of these four components, we reduce the linearized ideal magnetohydrodynamic (MHD) equations in the shearing box to a set of four first order differential equations for eigen-variables characterizing each component. Then, analyzing the dynamics using the eigen-variables, we grasp the transient growth as well as the mode linear coupling processes. We also contrast the exponential growth of axisymmetric magnetic mode perturbations to transient growth of non-axisymmetric ones.
The paper is organized as follows. In Section 2, we describe the mathematical formalism -linearize basic equations and cast them into a system of four first order differential equations for mode eigen-variables. Section 3 is devoted to the dynamics of magnetic and inertia-gravity modes and to the comparison of axisymmetric and non-axisymmetric growths of the dominant magnetic mode. Summary and discussion are given in Section 4.
In the shearing box model, disc dynamics is studied in a local Cartesian reference frame co-rotating with the disc's angular velocity at some fiducial radius r0 from the central star, so that curvature effects due to cylindrical geometry of the disc are ignored. In this rotating coordinate frame, the unperturbed differential rotation of the disc manifests itself as a parallel azimuthal flow with a constant velocity shear in the radial direction. The Coriolis force is included to take into account the effects of the coordinate frame rotation. The vertical component of the gravity force due to the central star is also present and responsible for disc's vertical stratification. The background magnetic field B0 is assumed weak (subthermal) in the sense that the Alfvén speed is much less than the sound speed, so that we can work in the Boussinesq approximation, as in Balbus & Hawley (1991. As a result, the basic dynamical equations of incompressible ideal MHD in the shearing box are (e.g., Goodman & Xu 1994): where u is the velocity in the local frame, ρ is the density, B is the magnetic field, P = p + B 2 /8π is the total pressure, Ω is the angular velocity of the local rotating reference frame, equal to the disc's angular velocity at r0: Ω(r0). x, y, z are, respectively, the radial, azimuthal and vertical coordinates with the corresponding unit vectorsx,ŷ,ẑ in these directions. The shear parameter q = 1.5 for the Keplerian differential rotation considered in this paper. The set of equations (1)-(4) has an equilibrium solution that is stationary and axisymmetric. In this unperturbed state, the velocity field of Keplerian rotation represents, as noted above, a parallel azimuthal flow, u0, with a linear shear in the radial direction u0x = u0z = 0, u0y = −qΩx.
In the shearing box model, equilibrium density ρ0 and total pressure P0 depend only on the vertical z−coordinate and satisfy the hydrostatic balance equation while the background magnetic field along the z-axis, B0 = B0zẑ, is assumed to be spatially uniform, as is typical for configurations with non-zero net flux.

Perturbation equations
Consider now perturbations u = u0 + u ′ , ρ = ρ0 + ρ ′ , P = P0 + P ′ , B = B0 + B ′ to the equilibrium described above. Assuming the perturbation of density to be small, ρ ′ ≪ ρ0, from equations (1)-(4) we obtain the following system governing the perturbation dynamics where The form of equations (5)-(8) permits decomposition of the perturbations into shearing plane waves, or spatial Fourier harmonics (SFHs) with time-dependent amplitudes and phases where F ≡ (u ′ , ρ ′ , P ′ , B ′ ) denotes the perturbed quantities and F ≡ (ū,ρ,P ,B) are the amplitudes of corresponding SFHs. The azimuthal, ky, and vertical, kz, wavenumbers remain unchanged, whereas the radial wavenumber kx(t) varies with time at a constant rate qΩky if ky = 0 (i.e., for non-axisymmetric perturbations) due to sweeping of wave crests by the background shear flow. In other words, an initially leading SFH (with kx(0)/ky < 0) eventually becomes trailing (with kx(t)/ky > 0) as time goes by. This change of SFH's orientation from leading to trailing at kx(t) = 0 is called 'swing'. In the local disc model, shearing plane waves are a natural platform of decomposition of perturbations compatible with the shearing sheet boundary conditions (see e.g., Hawley et al. 1995;Heinemann & Papaloizou 2009b) and, in fact, represent the basic/simplest 'elements' of dynamical processes at linear shear (Yoshida 2005). Substituting (9) into equations (5)-(8), all the non-linear terms vanish identically in these equations and we arrive at the following system of linear first order ordinary differential equations that govern the dynamics of SFHs of perturbations We note that vanishing of non-linear terms for perturbations of a single Fourier harmonic form (9), that reduces full non-linear equations exactly to linear ones, is a general property of incompressible MHD discovered by Goodman & Xu (1994). Thus, equations (10)-(16) and the subsequent analysis are, in fact, valid for any values of the amplitudesū andB.
Using the continuity equation, the density amplitudeρ can be found through the z-velocity amplitude dρ dt = −ūz dρ0 dz and then from equations (13)-(16), be expressed through the vertical magnetic field, Before proceeding further, we make a non-dimensionalization of quantities in terms of the orbital time Ω −1 , sound speed cs and scale-height H = cs/Ω specific to the disc's local dynamics: Despite working in the incompressible limit, we still prefer to use the sound speed instead of the Alfvén speed vAz = B0z/(4πρ0) 1/2 in the normalizations. A more natural normalization length scale to be used in the present incompressible problem would obviously be a characteristic length Hm = vAz/Ω along the z−axis associated with the background magnetic field, but the reason for using the scale-height, H, instead is that this will be useful below in making a connection with the compressible case by means of the plasma β = 8πρ0c 2 s /B 2 0z = 2H 2 /H 2 m parameter. Typically, this parameter is large for weak fields (i.e., Hm ≪ H) and we take a value β = 400 often adopted in simulations of MRI-driven turbulence in discs (e.g., Guan et al. 2009).
In the Boussinesq approximation employed here, all the equilibrium quantities and their vertical derivatives entering equations (10)-(17) are approximately uniform in the vertical z−coordinate. This holds if the vertical length-scale of perturbations is much shorter than the disc scale-height, implying that our incompressible analysis is, strictly speaking, self-consistent for kz ≫ 1 with the above scaling. Still, as will be shown below, it can give an idea about the trend of mode dynamics (amplification) at kz > ∼ 1, near the compressibility regime.
Next, substituting the velocities from equations (13)-(15) and density from equation (17) into equations (10)-(12), we arrive at the set of second order differential equations only for the magnetic field components, which after the above non-dimensionalizations take the form (henceforth bars over the amplitudes will be omitted) supplemented with the divergence-free constraint where N 2 = −d ln ρ0/d ln z is the squared Brunt-Väisälä frequency normalized by Ω 2 ; it is positive here and therefore the disc is convectively stable. Everywhere below we mostly use a value N 2 = 0.8, although show that the results do not change much with different stratifications. Equations (18)-(21) form the basis for our subsequent analysis. They contain all the information on various perturbation modes existing in the shearing box model of a stratified, incompressible disc threaded by a constant background vertical magnetic field. We will classify these modes below by making use of the WKBJ approach.

Total perturbation energy and α parameter
For further analysis it is useful to introduce the total spectral energy density of SFHs and derive its evolution equation. From equations (10)-(17), we obtain the following equation governing the energy evolution (see e.g., Balbus 2003) where asterisks denote complex conjugate and the non-dimensional total spectral energy density E = E k + Eg + Em is composed of three basic terms: the kinetic energy density E k = |ux| 2 + |uy| 2 + |uz| 2 , the potential, or buoyancy energy density and the magnetic energy density The right hand side of equation (22) is the shear parameter multiplied by the sum of Reynolds and Maxwell stresses, that is, by the well-known Shakura-Sunyaev α parameter. These stresses ensure energy exchange process between the disc mean flow and perturbations; in the absence of the shear (i.e., in a rigidly rotating disc), the total energy of perturbations is conserved. So, the growth of small perturbations, and therefore the MRI, is entirely fed with the energy drawn from the disc's differential rotation by means of the stresses.

Eigen-variables and classification of modes -WKBJ approximation
The shear parameter q enters equations (18)-(21) implicitly through the time-dependent radial wavenumber kx(t) and explicitly in the right hand side of equation (18). Therefore, we can distinguish between two kinds of effects induced by the disc flow shear on the perturbation dynamics. The implicit dependence on the shear leads to finite-time phenomena -transient growth and coupling of different modes -in the dynamics of non-axisymmetric perturbations that arise from the drift (time-variation) of kx(t) in the wavenumber space (see next subsection); they vanish for axisymmetric perturbations with ky = 0. The explicit dependence on the shear is not related to non-axisymmetry and equally affects both axisymmetric and non-axisymmetric perturbations. This term sets the frequency of epicyclic oscillations and, in particular, is responsible for the existence of the MRI in differentially rotating discs.
To classify and characterize all the modes involved in equations (18)-(21), first we consider the regime where the radial wavenumber is much larger than the azimuthal one, |ky| ≪ |kx(t)|. This condition implies that the time-dependent radial wavenumber kx(t) varies a little during the shear/dynamical time, (qΩ) −1 |dkx/dt| = |ky| ≪ |kx(t)|, so that the WKBJ, or adiabatic approximation holds with respect to the small parameter ky/kx and hence the effect of disc flow shear stemming from non-axisymmetry of perturbations does not play a role. We will see later that in the opposite regime |ky| > ∼ |kx(t)| (which ensues eventually as kx(t) drifts along the kx−axis), the dynamics is non-adiabatic (non-WKBJ) and shear-induced transient effects become important. So, in the present case, we ignore for the moment the variation of kx with time in equations (18)-(21), retain only the explicit shear term and look for solutions in the standard form ∝ exp −i t ω(t ′ )dt ′ . Substituting this into equations (18)-(21) and assuming that the time-dependent frequency, ω(t), also satisfies the same WKBJ condition |dω(t)/dt| ≪ ω 2 (t) for |ky| ≪ |kx|, we get a linear system of equations Associated eigenfrequencies are obtained by equating the determinant to zero, which yields the formal WKBJ dispersion relation: In the limit ky = 0, equation (27) reduces to the well-known dispersion relation of Balbus & Hawley (1991) describing the axisymmetric MRI. So, we can view it as a generalized non-axisymmetric WKBJ dispersion relation at small ky.
Dispersion relation (27) is a quadratic equation with respect to ω 2 and has two distinct solutions corresponding to two physically different types of perturbation modes: 1. An inertia-gravity wave with the frequency the restoring force for which is mainly provided by rotation and vertical gravity (buoyancy), but is modified by magnetic field tension. For very large |kx| ≫ |ky|, |kz| the effect of rotation is negligible and this frequency asymptotically tends to ω 2 g ≈ N 2 +2k 2 z /β. Generally, inertia-gravity waves play an important role in overall disc dynamics and they were studied extensively in the past mainly in the context of unmagnetized discs (Lubow & Pringle 1993;Korycansky & Pringle 1995;Ogilvie 1998;Balbus 2003;Tevzadze et al. 2003Tevzadze et al. , 2008). Obviously, these waves have ω 2 g > 0, since the disc is stratified stably. 2. A magnetic mode with the frequency/growth rate 1 the restoring force for which is mainly provided by magnetic tension, but modified by rotation and vertical buoyancy. Due to the disc's differential rotation, this mode becomes unstable (ω 2 m < 0), that is, exhibits the MRI when k 2 z < qβ, at radial wavenumbers smaller than a certain critical value, as it follows from expression (29). At larger kx, the magnetic mode is stable and oscillatory; it is essentially of the same nature as Alfvén waves in classical MHD. Again, in the limit of very large radial wavenumbers |kx| ≫ |ky|, |kz| the effect of rotation is negligible in the dynamics of the magnetic mode and its oscillation frequency tends to Alfvénic, ω 2 m ≈ 2k 2 z /β.
Thus, in a stratified incompressible magnetized disc, perturbations can be classified into two basic types -the inertia-gravity wave and magnetic modes and any general perturbation can be decomposed into the sum of these two modes. Eigenvalue equations (23)-(26) can be transformed into a slightly different form. We eliminate the total pressure P through Bx and By from equations (25)-(26) and switch to new independent variables h1 = Bx, h2 = dBx/dt, h3 = By, h4 = dBy/dt, which will be our basic variables in the subsequent analysis. After substituting these into equations (23)-(24), we obtain where h = [h1, h2, h3, h4] T is the state column vector and the components of a new matrix A are Eigenfrequencies associated with equation (31) are found from Calculating the determinant obviously gives the same dispersion relation (27). So, the eigenfrequencies ±ωm and ±ωg, given by (28) and (29), multiplied by −i represent the eigenvalues of the matrix A. Since the modes' frequencies and therefore the eigenvalues are generally distinct and corresponding eigenvectors are linearly independent, we can use the spectral decomposition, or diagonalization of A (Golub & Van Loan 1996), where Λ is the diagonal matrix formed from the eigenvalues iωm, −iωm, iωg, −iωg of the A matrix, and the columns of a matrix C are the eigenvectors of A, such that the first column corresponds to the eigenvalue iωm, second column to −iωm, third to iωg and fourth to −iωg, where the auxiliary functions G1(ω) and G2(ω) are defined as An advantage of the spectral decomposition is that it allows us to reduce equation (31) to a set of decoupled first order differential equations, i.e., to corresponding canonical equations. Substituting (32) into equation (31) and defining a new vector variable g ] T as ψ = C −1 · h, we obtain the canonical form of equation (31) for ψ, are dynamically decoupled and evolve independently from each other, respectively, with the frequencies ωm, −ωm, ωg, −ωg. Physically,  (28) and (29), respectively, versus the time-dependent ratio kx(t)/kz at ky/kz = 2, kz = 6 and N 2 = 0.8. ω 2 m is negative during a finite time, in the interval |kx(t)/kz | < 8.1 (i.e., |kx(t)/ky | < 4.05 as given by condition 30), indicating the MRI of the non-axisymmetric magnetic mode. At large describe inertiagravity waves -and in this respect are very convenient and revealing. Below we will refer to them as mode eigen-variables. The absence of coupling among modal equations (33) in the linear analysis means that in the WKBJ limit, the perturbation modes evolve without exchanging energy among each other, i.e., initially exciting one either mode with a specific characteristic time-scale does not lead to the excitation of other modes with different time-scales. In the following subsection, we will see that in the non-WKBJ regime, these equations are no longer independent from each other, i.e. become coupled that, in turn, results in the mutual linear coupling between inertia-gravity waves and the magnetic mode.

The dynamical equations in the non-WKBJ regime
In the WKBJ/adiabatic regime considered above, the radial wavenumber kx(t) and hence the mode eigenfrequencies have been assumed approximately independent of time by virtue of the limit |ky| ≪ |kx(t)|. However, after some time, kx(t), drifting along the kx−axis, enters the range |ky| > ∼ |kx|, where the WKBJ treatment breaks down and one must take into account, neglected so far, the dependence of kx(t) on time due to the disc flow shear in equations (18)-(21). In other words, in this non-WKBJ/non-adiabatic regime, transient effects arising from the implicit presence of the shear in equations (18)-(21) manifest themselves in the mode dynamics.
Combining equations (18)-(21), we eliminate the total pressure Figure 2. Real (solid lines) and imaginary (dashed lines) parts of the effective, or modified frequencies/growth rates, ωm −iM 11 , ωm +iM 22 , ωg − iM 33 , ωg + iM 44 , versus kx(t)/kz at the same ky/kz, kz and N 2 as in Fig. 1. For comparison, we also plot Re(ωm) (dashed-dotted line) and Im(ωm) (dotted line) in the upper two panels and purely real ωg (dashed-dotted line) in the lower two panels. Sharp peaks near kx/kz = ±8.1(kx/ky = ±4.05) in the upper two panels correspond to the points ω 2 m = 0, where the WKBJ approximation breaks down and the mode dynamics becomes non-adiabatic. Consequently, the difference between the effective frequencies/growth rates and ωm, ωg found in the WKBJ limit is most appreciable in the non-adiabatic region, between these two peaks. from equation (21), back into equations (18) and (19), we finally arrive at the system of two second order differential equations Equations (34) and (35) Balbus & Hawley (1992), but cast instead in terms of the components Bx and By. They permit to more easily and conveniently carry out parallel analysis of non-axisymmetric and axisymmetric dynamics (in the latter case there is no need for an involved procedure for recovering the axisymmetric limit as done in Balbus & Hawley (1992)). In a vertically stratified incompressible disc, there is no characteristic horizontal length-scale and therefore kx and ky enter equations (34) and (35) in combination kx/k and ky/k, whereas kz enters explicitly as multiplied by Hm. So, these equations in fact contain three free parameters, kx/kz, ky/kz and kz by means of which we will describe the dynamics of the modes. Without loss of generality, we can take ky and kz to be positive, ky, kz > 0.
Equations (34) and (35) rewritten in terms of the variables h1, h2, h3, h4 introduced above take the form where A is the former matrix as in equation (31) except now timedependent through kx(t). A1 is a new matrix originating from the time-dependence of kx and is therefore proportional to the shear parameter q, We see that in the non-WKBJ regime, where the variation of kx with time is allowed for, new shear-related terms appear in basic equation (36) compared to its counterpart WKBJ equation (31); these terms are proportional to ky/k and vanish in the WKBJ limit and/or for axisymmetric perturbations.
To extend canonical equations (33) to the non-WKBJ case, we should take into account that the matrix C in the spectral decomposition of A(t) also varies with time through kx(t). As a consequence, substituting representation (32) into equation (36), we obtain a more general set of canonical equations governing the linear dynamics of eigen-variables in both adiabatic and non-adiabatic regimes where the generalized eigen-variables are given, as before, by They are a natural extension of the eigenvariables defined above in the WKBJ regime. Λ(t) is the former eigenvalue matrix which has also become time-dependent through kx(t). The matrix is a new term, compared to canonical WKBJ equations (33), introduced by the shear. Like A1, it is also proportional to ky/k and thus was absent in the WKBJ regime considered above. In order to better understand its role in the mode dynamics, it is convenient to write equation (37) componentwise It is clear that the matrix M modifies the mode frequencies as well as couples the canonical equations for each eigen-variable, so we appropriately refer to it as a coupling matrix. It involves the timederivative of C(t) and the matrix A1, which are both proportional to the shear parameter q and to the azimuthal wavenumber ky. As a result, the coupling matrix identically vanishes in the shearless limit q = 0 (rigidly rotating disc) and/or for axisymmetric perturbations with ky = 0, reducing equations (38)  (33). In other words, the axisymmetric magnetic mode and inertiagravity waves are not mutually coupled. So, the linear dynamics of non-axisymmetric perturbations, involving possible couplings and interactions among various modes, is richer and much more complex than that of axisymmetric ones. The homogeneous (left-hand side) parts of equations (38) describe the individual dynamical behaviour of the magnetic and inertia-gravity modes. The inhomogeneous (right-hand side) parts, originating from the non-diagonal components of the M matrix, describe mutual couplings between the magnetic mode bearing the MRI and inertia-gravity waves as well as interaction between counter-propagating components contained in each of these modes. It should be stressed that these shear-induced inhomogeneous terms in the governing equations (responsible for the linear mode coupling phenomenon) are relatively new in comparison with related studies on the linear dynamics of perturbations in magnetized discs (e.g., Balbus & Hawley 1992;Brandenburg & Dintrans 2006;Johnson 2007). It is seen that compared to the WKBJ limit, the corresponding instantaneous frequencies/growth rates, ωg and ωm, of both modes are modified by the diagonal components M11, M22, M33, M44. Now the quantities ωm − iM11, ωm + iM22, ωg − iM33, ωg + iM44 appear as effective, or modified instantaneous frequencies/growth rates and determine amplification properties (instability) of the magnetic and inertia-gravity modes. Figure 1 shows the squared frequencies/growth rates ω 2 m and ω 2 g as a function of time during the swing of a perturbation SFH from leading to trailing. ω 2 g remains always positive, whereas ω 2 m becomes negative during a limited time, indicating the transient non-axisymmetric MRI of the magnetic mode SFH. Figure 2 shows a similar time history of the modified frequencies/growth rates, together with ωm and ωg for comparison. It is seen that the modified frequencies/growth rates deviate noticeably from the WKBJ values ωm and ωg mostly in the non-WKBJ range |kx(t)/kz| < ∼ 10 (i.e., |kx(t)/ky| < ∼ 5), implying that transient effects induced by the shear due to non-axisymmetry play an important role in the dynamics (instability) of individual modes at moderate and small radial wavenumbers; at large |kx(t)/kz| → ∞ the dynamics is mostly adiabatic and the modified frequencies go to ωm and ωg. Notice in Fig. 1 that the non-adiabatic interval |kx(t)/kz| < ∼ 10 also covers the range where the MRI would occur (i.e., where ω 2 m < 0) if it were the WKBJ limit, so the growth rate of the non-axisymmetric MRI is actually determined not by ωm, as in the WKBJ and/or axisymmetric cases, but rather by the modified frequencies ωm ∓ iM11,22. Therefore, calculating the growth rate of the non-axisymmetric MRI using only ωm, that is, concentrating on the WKBJ limit during an entire evolution, as done in several papers (e.g., Balbus & Hawley 1992; Goodman & Xu 1994; Johnson 2007), may not give a full information on the non-axisymmetric instability of the magnetic mode. This implies that when studying the dynamics of non-axisymmetric MRI in discs, one must necessarily take into account the effects of disc flow shear in full -first as being responsible for the existence of the MRI itself (the explicit shear term in linear equations 18-21) and second as bringing about (implicitly through time-variation of kx) transient amplification and possible mutual couplings of different perturbation modes, including the magnetic mode present in a magnetized disc. Figure 3 shows the non-diagonal elements of M as a function of time. M12, M21 and M34, M43 describe the coupling, respectively, between the counter-propagating components 2 , ψ of inertia-gravity waves. The elements M13, M14, M23, M24 and M31, M32, M41, M42 describe the coupling between the magnetic and inertia-gravity modes. In other words, initially imposed only one either of these two modes can act as a source for another mode and excite it in the course of evolution. As expected, all these non-diagonal components reach their maximum values in the non-adiabatic range |kx(t)/kz| < ∼ 10, where the effects of the shear due to non-axisymmetry of the modes play a role. So, the coupling of the magnetic and inertiagravity modes, being induced by the disc flow shear, is expected to occur in this wavenumber range; at large |kx(t)| → ∞ the non-diagonal elements vanish -there is no energy exchange in the WKBJ regime. Because of the existence of such a coupling, the magnetic and inertia-gravity modes actually have strictly separate identities only in the WKBJ regime, when the characteristic time-scales (frequencies/growth rates) of these modes are different (Fig. 1) and the right-hand side terms of equations (38), responsible for the mode interaction, disappear. By contrast, in the non-WKBJ regime, the mode time-scales are comparable and the coupling terms are appreciable (Figs. 2 and 3), so that the modes cannot be quite disentangled from each other and we have some mixture of both modes at such times. Finally, we note the following relations between the various components of the M matrix: 43 . Thus, the shear associated with the disc's differential rotation plays a twofold role in the dynamics of perturbation modes: 1. modifies the evolution (stability) of each mode. In particular, causes the magnetic mode to become magnetorotationally unstable in a certain range of radial wavenumbers (i.e., we have transient MRI), 2. introduces the phenomenon of linear mode coupling. In the numerical analysis performed below, we will see the consequences of this mode coupling in the evolution of the perturbed magnetic field components.

DYNAMICS OF NON-AXISYMMETRIC MAGNETIC AND INERTIA-GRAVITY MODES
In this section, we separately investigate the evolution and amplification properties of non-axisymmetric SFHs (with ky = 0) of basic modes in the system -inertia-gravity wave and magnetic mode -by numerically solving equation (36) (or, equivalently equations 34 and 35). Then, we concentrate mainly on the magnetic mode because it determines the MRI in the disc and perform a comparative analysis of the transient growth of non-axisymmetric and the exponential amplification of axisymmetric magnetic modes SFHs during a dynamical time for a range of wavenumbers. For the numerical integration of equation (36) we used a standard Runge-Kutta scheme (MATLAB ode45 RK implementation). We start integration by initially imposing a tightly leading (i.e., with kx(0) < 0, |kx(0)| ≫ |ky|, |kz|) SFH of purely magnetic mode or purely inertia-gravity waves on the disc flow and trace the subsequent evolution of the perturbed quantities corresponding to these initial conditions until kx(t) ≫ |ky|, |kz|, when the SFH becomes tightly trailing. At |kx(0)| ≫ |ky|, |kz|, the SFH starts out from the WKBJ regime, where the coupling coefficients are small, so the modes are not coupled with each other and therefore can be cleanly separated. This, in turn, permits us to isolate either of these two modes and to insert it as an initial condition into equation (36). In this adiabatic regime, the eigen-variables corresponding to counter-propagating components of the magnetic and inertia-gravity modes evolve independently from each other and are given by the asymptotic solutions of the homogeneous parts of equations (38), To prepare initial conditions to equation (36) corresponding to a purely magnetic mode imposed at t = 0, we select out only the eigen-variable ψ m as an initial condition is dictated by the fact that it relates to an exponentially growing branch, i.e., −iωm = γ > 0 when ω 2 m = −γ 2 < 0, where γ is the growth rate in the WKBJ regime and/or in the axisymmetric case. These initial conditions ensure that there is no contribution from inertia-gravity waves at the beginning. Similarly, to impose purely inertia-gravity waves in the equations, we select out only ψ  Figures 4 and 5 show the subsequent time-development of the magnetic field components, Bx, By, Bz, and normalized by their initial values kinetic and magnetic energies E k /E k (0), Em/Em(0) during the evolution of SFHs corresponding to the above initial conditions from being tightly leading to tightly trailing. For the initially imposed magnetic mode SFH (Fig.  4), being still in the WKBJ regime at kx(t)/kz ≪ −1, the azimuthal magnetic field By oscillates with the Alfvén frequency ωm ≈ kz 2/β and constant amplitude, whereas the radial Bx and vertical Bz components remain much smaller than By, vanishing  Figure 4. Evolution of the magnetic field components pertaining to an initially imposed magnetic mode SFH with kx(0)/kz = −240, ky/kz = 2, kz = 6, N 2 = 0.8 (solid lines correspond to real parts, dot-dashed to imaginary parts). The azimuthal field By is dominant over the other two components in the magnetic mode. At the early adiabatic stage of evolution, By oscillates with a constant amplitude, while the amplitudes of Bx and Bz are small. Then, in the non-adiabatic region 77 < t < 83 (i.e., |kx(t)/kz | < 8.1), all the three components undergo transient growth and the corresponding kinetic and magnetic energies amplify by a large factor of 100, indicating the transient non-axisymmetric MRI exhibited by the magnetic mode. Afterwards, in the next adiabatic region on the tightly trailing side at t > 100, Bx falls off, while By and Bz continue to oscillate with the frequencies ωm and ωg, respectively, and with the amplitude of the former being much larger than that of the latter. In other words, at the non-adiabatic stage, the dominant azimuthal component of the magnetic field generates vertical one, which was asymptotically zero at the initial tightly leading stage. This is a consequence of the linear mode coupling phenomenon -generation of the inertia-gravity mode SFH by the magnetic mode SFH.
in the limit kx(0) → −∞, that is, the magnetic mode is characterized by the dominant azimuthal magnetic field perturbation. This is not surprising if we recall that in a classical Alfvén wave, magnetic field perturbation lies perpendicular to the background magnetic field and wavevector. Since, in the magnetic mode, the main restoring force is magnetic tension, the kinetic and magnetic energies are comparable and larger than the potential (buoyancy) energy. By contrast, for the imposed inertia-gravity mode SFH (Fig.  5), at the initial WKBJ stage, Bz oscillates with the frequency ωg ≈ N 2 + 2k 2 z /β and constant amplitude, whereas Bx and By remain much smaller than Bz, vanishing at kx(0) → −∞, that is, in inertia-gravity waves, dominant is the vertical magnetic field perturbation. This is due to the fact that the restoring force for these waves at large kx is provided mainly by vertical buoyancy (but is modified by magnetic tension) that induces vertical magnetic field oscillations via induction equation. As a result, for the inertia-gravity wave SFH, the kinetic energy is of the same order as the buoyancy energy and both are much larger the magnetic energy.
As time passes, the radial wavenumber kx(t) of the SFHs, drifting along the kx−axis, approaches the non-WKBJ region, where the elements of the M matrix are appreciable and therefore the shear-induced effects become important. The components   Figure 5. Evolution of the magnetic field components pertaining to an initially imposed inertia-gravity mode SFH with kx(0)/kz = −240, ky/kz = 2, kz = 6, N 2 = 0.8 (solid lines correspond to real parts, dot-dashed to imaginary parts). In the inertia-gravity wave SFH, the vertical field Bz dominates over Bx and By. Initially, in the adiabatic region, Bz oscillates with a constat amplitude, while Bx and By are asymptotically zero. Then, during crossing the non-adiabatic region 77 < t < 83 (i.e., |kx(t)/kz | < 8.1), Bx and By components undergo small growth, so that the corresponding kinetic energy is almost unchanged and the magnetic energy exhibits a small transient amplification by a factor of 1.25. Afterwards, on the tightly trailing side at t > 100, Bx falls off, while By and Bz continue to oscillate with the frequencies ωm and ωg, respectively. In other words, at the non-adiabatic stage, the dominant vertical component of the magnetic field generates azimuthal one, which was asymptotically zero at the initial tightly leading stage. This is again a consequence of the linear mode coupling phenomenon -generation of the magnetic mode SFH by the inertia-gravity mode SFH.
Bx, By and Bz for the magnetic mode SFH (Fig. 4) and Bx and By for the inertia-gravity wave SFH (Fig. 5) start to grow transiently. Then, at about t1 = 77, kx(t1)/kz = −8.1 the SFHs enter the non-WKBJ region, cross it and leave at about t2 = 83, kx(t2)/kz = 8.1 (see Fig.2). For the magnetic mode, the modified frequency ωm + iM22 associated with the eigen-variable ψ (−) m has a positive imaginary part in the non-WKBJ region (Fig. 2), causing this mode to undergo transient-exponential amplification in this interval. As a result, at around t = 80, all the three components, Bx, By and Bz, reach maximum values about 10 times larger than those on entering the non-adiabatic region, but the azimuthal component, By, still remains dominant. Accordingly, the magnetic Em and kinetic E k energies of the magnetic mode SFH grow about 100 times during the transient amplification event. This is essentially the transient MRI experienced by the non-axisymmetric magnetic mode SFH. On leaving the non-adiabatic region, in the next adiabatic region, where the SFH is now tightly trailing (with kx(t)/kz >> 1), Bx gradually decays and asymptotically, at kx(t) → ∞, tends to zero; By continues to oscillate with the amplified amplitude at an Alfvén frequency and the amplitude of Bz levels off at a constant value, though its oscillation frequency clearly differs from that of By and is equal to the frequency of inertia-gravity waves, ωg ≈ N 2 + 2k 2 z /β. This is a direct consequence of the coupling between the magnetic and inertia-gravity modes. As noted above, in the inertia-gravity mode SFH dominant is Bz, whereas in the magnetic mode SFH By component. In the initially imposed magnetic mode SFH, the Bz component has been negligible at the tightly leading stage, because of the absence of inertia-gravity waves in the initial conditions (see Fig. 4). Then, in the non-adiabatic region, the non-diagonal coupling elements of M matrix become appreciable, resulting in the generation of the inertia-gravity mode SFH by the original magnetic mode SFH. 3 So, the Bz component that keeps oscillating with constant amplitude at large times, in the tightly trailing adiabatic regime in Fig. 4, corresponds to this newly generated inertia-gravity wave SFH and hence its oscillation frequency coincides with ωg. It can be said that during the evolution of the magnetic mode, azimuthal magnetic field generates the vertical one due to shear of the disc's differential rotation. However, the intensity of the excited inertia-gravity waves is still small compared with that of the original magnetic mode, whose dominant azimuthal field By is about 50 times larger than Bz of the excited inertia-gravity wave.
For the initially imposed inertia-gravity wave SFH (Fig. 5), the effective frequency ωg + iM44 acquires a positive imaginary part during a short time in the non-WKBJ region because of the shear-induced component M44, however, Bz does not exhibit any noticeable increase in amplitude and oscillates with nearly ωg, during an entire course of evolution of the SFH, from tightly leading to tightly trailing. Both Bx and By undergo moderate transient growth in the vicinity of t = 80, but remain smaller than Bz. The kinetic energy, E k , shows almost no growth and the magnetic energy, Em, grows only about 1.5 times. After leaving the non-adiabatic region, when the SFH wraps up tightly trailing, Bx gradually decays, as in the above case, but By continues to oscillate with a constant amplitude. As seen from Fig. 5, the oscillation frequency of By clearly differs from that of Bz and is equal to the Alfvén frequency, ωm ≈ kz 2/β. This is again a result of the coupling between the inertia-gravity and magnetic modes. The initially imposed inertiagravity mode SFH with dominant Bz component generates, during crossing the non-WKBJ region, the magnetic mode SFH with dominant By component. However, the amplitude of the Bz component related to the waves is about 5 times more than that of By belonging to the magnetic mode -the strength of the generated magnetic mode SFH is smaller than that of the original inertia-gravity wave SFH. We obtained a similar result for other values of ky and kz.
Having separately analysed the SFH evolutions of the magnetic and inertia-gravity modes, we have found that generally the magnetic mode tends to grow by orders of magnitude larger factors than the inertia-gravity waves do, i.e., it is the magnetic mode that displays the non-axisymmetric MRI and thus remains a dominant driver of dynamical processes (e.g., angular momentum transport) in magnetized discs. It should be noted that the amplification of the magnetic mode is not just transient -the energy drawn from the differential rotation during transient amplification phase is not returned back to the disc, but remains stored in the magnetic mode which then continues to oscillate with constant amplitude (Fig. 4). So, below we focus only on the magnetic mode and explore the ef-  (0)), corresponding to an initially imposed magnetic mode SFH, as a function of ky/kz and kz at N 2 = 0.8. The initial energy, Em (0), is taken at the tightly leading stage kx(0)/kz ≪ −1. Em,max is the maximum value achieved by the magnetic energy during the swing of the SFH from tightly leading to tightly trailing. The values of the growth factor below kz = 2 (dashed lines) are not strictly precise, because the incompressibility requirement is not valid at such kz. ficiency of its growth for various values of the free parameters of the problem.
To quantify the transient growth of the magnetic mode SFH, we define the amplification factor -a ratio of the maximum value achieved by the magnetic energy over an entire evolution of this SFH from beginning (tightly leading) to end (tightly trailing) to its initial value, Em,max/Em(0). The way we characterize the amount of amplification, and hence the non-axisymmetric MRI, differs to that of Goodman & Xu (1994), who used an ex-  Figure 8. Normalized α/α(0) parameters vs. kx(t)/kz corresponding to initially imposed magnetic mode (m) and inertia-gravity wave (g) SFHs at ky/kz = 2, kz = 6, N 2 = 0.8. The total stress related to the magnetic mode is positive, while that of inertia-gravity waves is negative in the dynamically interesting non-adiabatic region |kx(t)/kz | < ∼ 10 and both decay at late times. ponential growth factor integrated over an entire evolution time, exp( ∞ −∞ Im[ωm(t ′ )]dt ′ (in their notations time changes from t = −∞ to t = ∞) to characterize non-axisymmetric parasitic instabilities afflicting axisymmetric channel solution. Figure 6 shows the amplification factor calculated for various ky and kz. A maximum value at ky/kz < ∼ 0.1 occurs near kz = 10, which approximately coincides with the critical wavenumber of the most rapidly growing axisymmetric MRI equal to kz,cr = (15β/32) 1/2 = 13.7 ≫ 1 with our normalization. This is easy to understand: at small ky/kz, the WKBJ approximation holds during almost an entire evolution of SFH and one can find the amplification factor by a simple integration of the WKBJ growth rate given by relation (29) over the time interval [t1, t2], where t1 and t2 are the moments when drifting kx(t), respectively, enters and leaves the unstable region (30), exp( In this case, the extremely large growth factors (∼ exp(10)) result from the fact that SFHs with small ky/kz ≪ 1 take long (hundreds of orbits) to cross the unstable region. In this respect, the growth of non-axisymmetric SFHs with such small ky, despite being huge, do not actually play a role because it occurs over times much larger than the dynamical/orbital time. At ky/kz > ∼ 1, the growth factor is independent of ky/kz and increases with decreasing kz (see also Fig. 10). This indicates that as opposed to the (nearly) axisymmetric mode, which has its maximum exponential growth rate at larger vertical wavenumbers, most of growth of the non-axisymmetric magnetic mode occurs preferably at smaller kz < ∼ 5. However, the very large (∼ exp (8)) values of the growth factor below kz < ∼ 2 (dashed lines in Fig. 6) are not really precise, since they are related to SFHs with vertical wavelengths comparable to the disc scale-height for which compressibility effects, neglected here because of the Boussinesq approach, are important (see also next subsection, Fig. 10).
We also analysed the role of vertical stratification in the dynamics of non-axisymmetric modes. It is well-known that (con-vectively stable) stratification does not affect the fastest growing axisymmetric perturbations (Balbus & Hawley 1991). In Fig. 7, we show the time-development of the normalized magnetic energy, Em/Em(0), for the magnetic mode SFH at various stratifications at fixed ky and kz. It is seen that the growth factor of energy weakly depends on the stratification, especially during the exponential-transient amplification phase around kx(t) = 0. This is easily explained: the magnetic mode's amplification is mainly provided by the combined action of rotation and magnetic tension, with buoyancy playing only a minor role. At later times, at the tightly trailing stage (kx(t)/kz = 2kx(t)/ky ≫ 1), the energies oscillate with amplitudes that are a bit larger at smaller stratifications. Thus, stratification has a weak influence on the growth of the non-axisymmetric MRI in the linear regime, however, it plays an important role in the non-linear development and saturation of the MRI (Davis et al. 2010).
In Fig. 8 we also plot the time-development of the total stress, a sum of Reynolds and Maxwell stresses, given by equation (22), corresponding to the initially imposed leading harmonics (with kx(0)/ky < 0) of inertia-gravity waves and magnetic mode. Being still at the leading stage (kx(t)/ky < 0), both are positive and comparable to each other, but at the trailing stage, at 0 < kx(t)/kz = 2kx(t)/ky < ∼ 10, transport due to the magnetic mode is positive and larger by absolute value than the negative transport due to inertiagravity waves. So, inertia-gravity wave harmonics in this range of wavenumbers, like convection, tend to transport angular momentum inwards (see also Balbus 2003). However, the harmonics of both modes oscillate at large times, in the tightly trailing regime, and their respective α, being also oscillatory, gradually decay. So, in the linear theory, the sign of transport cannot be determined with certainty. In the non-linear regime, however, inertia-gravity waves may amount to moderate fraction of or become even comparable to the magnetic mode, despite the prevalence of the latter in the linear regime. In this case, both these dynamically intertwined modes (not only the magnetic mode) are to be taken into account when studying a nonlinear (turbulent) state in incompressible, stratified and magnetized discs. Hence, understanding the true nature of turbulent angular momentum transfer in such discs requires a scrupulous analysis of the contributions of magnetic mode and inertia-gravity waves in the transport.
Recently, Pessah & Chan (2012) have carried out an analogous study of the linear dynamics of non-axisymmetric shearing MHD waves, but with kz = 0, in an incompressible unstratified shearing box threaded by a weak magnetic field. Although they focused on accretion disc boundary layers with negative shear, q < 0, and therefore not subject to the standard MRI, their results regarding the evolution of energy and total stress for a magnetic mode SFH resemble those described above. So, we would like to draw a parallel between the findings presented in that paper and ours. In the shearing box without vertical stratification, for 2D perturbations with kz = 0, only the magnetic mode -Alfvén wave -survives, whereas inertia-gravity wave disappears. The transient amplification of these Alfvén waves, which was investigated by Pessah & Chan, has non-exponential character and can be attributed to the shear/non-normality of the boundary layer flow. 4 As a result, the corresponding amplification factors in the non-adiabatic region, in the vicinity of kx(t) = 0, are smaller than those of the magnetic mode SFHs in our case with the MRI, as shown in Figs. 4 and 6; outside the non-adiabatic region, the magnetic mode SFH behaves like that of 2D Alfvén wave does. Specifically, the evolution of the total stress for the magnetic mode SFH shown in Fig.  8 and for 2D Alfvén wave SFH are qualitatively similar: both experience transient growth near kx(t) = 0, however, the one related to the 2D Alfvén wave is negative there due to the negative shear, and oscillate around zero with decreasing amplitude at large |kx(t)/ky| ≫ 1. Also, the SFH energies, starting with a constant value, after a transient amplification phase, do not decay and keep on oscillating about a time-independent value at late times. In our 3D case with stratification, inertia-gravity waves also come into play. Although their transient amplification is small compared with that of the magnetic mode, as demonstrated above, they couple with the magnetic mode in the non-adiabatic regime near kx(t) = 0 and must be taken into account for a correct/fuller description of the mode dynamics.

comparative analysis of non-axisymmetric and axisymmetric SFHs' growth
An analysis of the competition between the non-axisymmetric and axisymmetric modes' linear dynamics is of primary interest, because this competition will determine the main cause/ingredient of the transition to turbulence and the resulting non-linear turbulent state. If the initial amplitude of perturbations is very small, axisymmetric magnetic mode, characterized by permanent exponential growth with constant rate, eventually wins in this competition and forms an axisymmetric channel flow (Hawley et al. 1995;Bodo et al. 2008;Longaretti & Lesur 2010). However, if the initial amplitudes of both non-axisymmetric and axisymmetric modes are sufficiently large, then outcome is not so obvious. Below we examine possible dynamics in the latter case. It is well-known that the most unstable axisymmetric magnetic mode SFH has the growth rate γmax = 0.75Ω, i.e., the dynamical (e-folding) timescale tc = 1/γmax = 1.33Ω −1 (Balbus & Hawley 1991;Goodman & Xu 1994). This is the shortest characteristic timescale of the axisymmetric MRI and effectively a dynamical time too, as it is of the same order as the shear time (qΩ) −1 = 0.67Ω −1 . So, a natural choice for the time interval during which the non-axisymmetric and axisymmetric growths are to be compared is therefore tc.
We have seen above that the magnetic energy of the nonaxisymmetric magnetic mode SFH, evolving from far tightly leading (kx(0)/kz ≪ −1) to tightly trailing phases (kx(t)/kz ≫ 1), increases by orders of magnitude (see Fig. 6). The amplification factor of the energy increases with decreasing the azimuthal wavenumber, ky. At small ky, as noted above, the amplification factor is determined to a good approximation by the WKBJ growth rate given by expression (29) and the corresponding duration of this amplification is the time the radial wavenumber, kx(t), takes to cross the unstable interval (30) in the wavenumber space. However, this swing time can be quite long in comparison with the dynamical/e-folding timescale, tc, and therefore the most powerful growth occurring for such small ky does not really play a decisive role as far as disc's dynamical instability is concerned. In some sense, "truly" non-axisymmetric dynamical growth, that occurs over times of order tc, starts from about ky/kz > ∼ 1 and kz > ∼ 1. Thus, to analyze the relative importance of the transient growth of non-axisymmetric magnetic mode SFHs compared with the fastest growing axisymmetric one during tc and to understand their role in the system dynamics, we adopt the concept of "optimal" perturbations (Butler & Farrell 1992;Farrell & Ioannou 1996, 2000Bakas et al. 2001), which is widely employed by geophysical/meteorological community. As applied to our problem, the strategy of the optimal perturbation approach is to consider the evolution of a single magnetic mode SFH with some initial kx(0), ky and kz within the time interval [0, tc], to calculate the resulting transient amplification and then to explore how it changes with varying these wavenumbers. Specifically, at each point in the (kx(0), ky, kz) space, at the initial moment of time, t = 0, we prepare initial conditions corresponding to a purely magnetic mode SFH (i.e., initially impose ψ m (0) = 0) and trace subsequent evolution during the time tc. From Figs. 4 and 7 we see that the energy of a magnetic mode SFH undergoes most of its growth on the trailing side at kx(t)/kz < ∼ 10 and, afterwards, turns into oscillatory with constant amplitude. Obviously, optimal SFHs are expected to be in this range (see below). Each non-axisymmetric SFH drifts in the (kx, ky)-plane along the kx-axis with speed qky, so that the optimal SFH with the initial radial wavenumber kx(0) will have the wavenumber kx(tc) = kx(0) + qkytc at tc. Hence, we can regard the energy of the optimal perturbation as a function of kx(tc), ky and kz.
We define the transient amplification factor of the magnetic mode SFH in terms of the associated magnetic energy, Em(t), as the ratio of its final, at t = tc, to the initial, at t = 0, values. Since this is to be compared with the largest exponential growth factor of the magnetic energy of the most unstable axisymetric magnetic mode during the same time interval, that is, with exp (2γmaxtc) = exp(2), we introduce a new parameter -a relative amplification factor fe(kx(tc), ky, kz) = exp(−2) Em(kx(tc), ky, kz) Em(kx(0), ky, kz) .
This perameter conveniently measures the non-axisymmetric growth of the magnetic energy against axisymmetric one during the dynamical time and thus characterizes relative significance of non-axisymmetric magnetic mode SFHs with respect to axisymmetric ones. SFHs with those wavenumbers that yield the maximum (or close to maximum) value of fe represent optimal (or close-to-optimal) perturbations. Figure 9 shows isocontours of the relative amplification factor, fe, in the (kx/kz, ky/kz)-plane at kz = 7 and N 2 = 0.8. In this map, the values range from 0.17 to the optimal (largest) value fe,opt = 1.65 at kx/kz = 1.8, ky/kz = 1.03. The values above unity fe reaches in the vicinity of the line ky = 0.57kx, which can be labeled as the line of extremum (maximum) of fe (obviously, the optimal SFH also lies on this line). Farther up along the line of extremum, it quickly levels off to an asymptotically constant value fe,max = 1.56. All the magnetic mode SFHs with fe > ∼ 1 are referred to as "close-to-optimal" harmonics. These harmonics are situated in the vicinity of the extremum line and located on the trailing side of the (kx, ky)-plane. Thus, it is evident that there is a fairly broad range of kx(tc), ky for which the energy of nonaxisymmetric magnetic mode SFHs can transiently grow by factors larger than a purely exponential amplification of energy for the most unstable axisymmetric magnetic mode during the dynamical time tc. This does clearly demonstrate that non-axisymmetric magnetic mode perturbations can prevail over, or compete with axisymmetric ones during tc at those wavenumbers. Contrasting Figs. 6 and 9, both describing non-axisymmetric amplification of the magnetic mode, we see different trends with respect to ky at fixed kz. In Fig. 6, where amplification factor, being calculated over much larger time interval for the magnetic mode SFH evo-  Figure 9. Isocontours of the relative amplification factor fe(kx(tc), ky, kz) in the (kx/kz, ky/kz)-plane at kz = 7 and N 2 = 0.8. It achieves an extremum (maximum), fe,max = 1.56, on the trailing side (kx/kz > 0), on the line ky = 0.57kx. This extremum is almost constant all along this line, reaching a slightly larger peak, or optimal value fe,opt = 1.65 at the lower end kx/kz = 1.8, ky/kz = 1.03. This peak corresponds to the optimal, while the values fe > 1 near the line of extremum, to close-to-optimal SFHs. As is clear, the amplification factors of the non-axisymmetric magnetic mode SFHs can be comparable to or even larger than that of the maximally growing axisymmetric one during the characteristic time of the latter. lution, increases with decreasing ky, because SFHs with smaller ky stay much longer in the unstable region (30). By contrast, in the second case in Fig. 9, the evolution time is fixed to tc, which is not enough for nearly axisymmetric (ky ≪ 1) harmonics with various initial kx(0) to undergo considerable growth, whereas nonaxisymmetric ones with moderate ky/kz > ∼ 1 manage to get amplified more within the same time.
These findings should be important for further development in the non-linear stage and, ultimately, for turbulence (and possibly coherent structures) that may emerge in this system. Of course, if the amplitudes of initial perturbations are very small and the linear regime holds for a long time, the axisymmetric/exponential MRI will eventually dominate. Optimal non-axisymmetric perturbations undergo most of growth only during the dynamical time interval. However, if the optimal and close-to-optimal SFHs have sufficiently large initial amplitudes, the non-linearity will come into play during the dynamical time and the further exponential growth of the axisymmetric MRI (e.g., axisymmetric channel flow) may not persist. In other words, in such a course of events, the optimal and close-to-optimal non-axisymmetric SFHs can prevail over axisymmetric SFHs and determine the onset of a non-linear regime and its dynamics, including angular momentum transport parameter α.
As is well-known, the exponential growth rate of an axisymmetric magnetic mode SFH reaches a maximum at the critical wavenumbers kx,cr = 0, kz,cr = (15β/32) 1/2 , decreases at kz < kz,cr and kz > kz,cr and reaches zero at kz = 0 and kz = √ qβ (Balbus & Hawley 1991). Let us now examine how the non-axisymmetric relative amplification factor (found above at given kz) varies with kz. Figure 10 shows the maximum, fe,max, along the above-mentioned line of extremum as a function of kz. Note that fe,opt is close to fe,max, so we have chosen to plot only the latter, because as evident from Fig. 9, those close-to-optimal SFHs with values of fe near to fe,max generally occupy larger area in the wavenumber space than the optimal SFH. From this figure it is seen that the relative amplification factor of non-axisymmetric magnetic mode SFH monotonically increases with decreasing kz. For kz < 12, fe,max > 1 indicating that at such vertical wavenumbers the transient growth of non-axisymmetric magnetic mode perturbations always exceeds that of the most unstable axisymmetric ones during the dynamical time. Therefore, non-axisymmetric MRI becomes dominant over, or at least competitive with the axisymmetric one in the disc dynamics. Note that at kz < ∼ 2, when vertical wavelength of perturbations becomes comparable to the disc scaleheight, compressibility is at work and the part of this curve at such kz, calculated here in the incompressible limit, does not give quite reliable values of the non-axisymmetric amplification factor. Nevertheless, Fig. 10 correctly predicts a trend -the amplification of non-axisymmetric magnetic mode SFH is the larger, the smaller kz is. By contrast, as mentioned above, the axisymmetric MRI attains its largest growth rate at critical vertical wavelengths a few times smaller than the scale-height, λz,cr/H = 2π/Hkz,cr = 2π(32/15β) 1/2 = 0.46. Extrapolating this to the compressible case, non-axisymmetric magnetic mode perturbations, which transiently grow the most, are expected to be more influenced by compressibility than axisymmetric ones. So, in order to accurately characterize the non-axisymmetric MRI in discs with net vertical magnetic field, one should take into account compressibility. This may be related to the saturation of magnetic energy and stresses in MRI-driven turbulence that exhibits dependence on compressibility and hence the latter plays an important role in the non-linear devel-opment of the non-axisymmetric MRI with non-zero net magnetic flux (e.g., Sano et al. 2004;Pessah et al. 2007;Guan et al. 2009).

SUMMARY AND DISCUSSION
In this paper, we have investigated the linear dynamics of nonaxisymmetric perturbations in a vertically stratified Keplerian disc with a weak uniform vertical magnetic field. We mainly considered incompressible perturbations, but, in the final part, speculated on the dynamics in the compressible case. An analysis was performed in the framework of the local shearing box approximation. As is usually done in this case, perturbations were decomposed into spatial fourier harmonics (SFHs), or shearing plane waves, and substituted into linearized ideal MHD equations of the shearing box. The temporal evolution of the amplitudes of SFHs was followed by numerical integration in time of the resulting set of linear differential equations. In contrast to previous linear studies on the non-axisymmetric MRI in discs (e.g., Balbus & Hawley 1992;Tagger et al. 1992;Brandenburg & Dintrans 2006;Johnson 2007), we took a different approach and rewrote these equations in the canonical form for characteristic variables, or eigen-variables (equations 38) that allowed us to identify and describe all the modes/types of perturbations contained in this disc model. In this respect, our analysis offers a deeper insight into the specific dynamics and amplification properties of individual modes, allows us to easily capture a new effect of their mutual couplings brought about by the disc's differential rotation. There are two basic modes in the considered incompressible disc model: inertia-gravity waves and magnetic mode. The magnetic mode, which is driven mainly by magnetic tension force but modified by rotation, displays the MRI. Inertia-gravity waves are driven by (convectively stable) vertical buoyancy and rotation and are spectrally stable, though nonaxisymmetric ones can still undergo transient growth at small radial wavenumbers due to disc flow shear. For this reason, we focused primarily on non-axisymmetric modes, since the dynamics of axisymmetric modes in magnetized discs is quite well understood (e.g., Balbus & Hawley 1991;Ogilvie 1998). Besides, shear effects manifest themselves only for non-axisymmetric perturbations and hence one can say that the dynamics of non-axisymmetric modes is more diverse and richer than that of axisymmetric ones. We compared the transient growth of initially imposed non-axisymmetric inertia-gravity waves and magnetic mode SFHs during evolution from tightly leading to tightly trailing and found that the latter exhibits much larger transient growth factors. Also, the stresses associated with it, dominate those due to inertia-gravity waves. In other words, the magnetic mode is a primary mode taking part in energetic processes and determines the MRI in discs. Another important and relatively new effect we observed is the coupling of magnetic mode and inertia-gravity waves. The evolution of an initially imposed magnetic mode SFH has shown that, for a certain range of vertical and azimuthal wavenumbers, the transient amplification of the magnetic mode SFH is accompanied by the generation of inertia-gravity wave SFHs, though with smaller amplitude -the non-axisymmetric MRI acts as a new source of inertiagravity waves. This is a consequence of the linear mode coupling phenomenon resulting from disc flow shear/non-normality (Chagelishvili et al. 1997). It demonstrates that, although inertiagravity waves are not a main participant in energetic processes defining disc instability, still their presence should be taken into account, together with the more important magnetic mode, in order to form a complete dynamical picture of the non-axisymmetric MRI in stratified discs. In magnetized discs, inertia-gravity waves carry predominantly vertical magnetic field perturbation, whereas in the dominant magnetic mode its mostly azimuthal. The linear coupling between these two modes implies that vertical magnetic field can be generated by azimuthal one and vice versa due to the disc flow shear.
Another goal was a parallel analysis of the dynamics of non-axisymmetric and axisymmetric magnetic modes. The accepted course of events in the presence of non-zero net vertical magnetic field is the following: initially the linear axisymmetric MRI arises whose subsequent non-linear development leads to the formation of a channel solution. Afterwards, on this channel flow, the secondary instabilities of (parasitic) non-axisymmetric perturbations develop (Goodman & Xu 1994;Pessah & Goodman 2009). Typically, in the stability analysis of a channel flow, one considers small non-axisymmetric perturbations to the channel solution once it has been well formed, i.e., reached amplitudes much larger than the background magnetic field (Goodman & Xu 1994;). In this case, the growth rate of parasitic instabilities turns out to be much larger than that of the channel itself. This fact considerably simplifies the analysis, because the background shear and rotation can be neglected and the channel can be assumed steady. In the above papers, a case of not very large channel amplitudes (compared to the background magnetic field) is also studied. Stability of this (intermediate amplitude) channel solution against non-axisymmetric perturbations is complicated by shear, rotation and the growth (nonsteadiness) of the channel itself. This problem is attacked using the technique of shearing waves (Goodman & Xu 1994). But what the above papers have in common is that they all consider the linear stability of already well-developed channel flows being subject to non-axisymmetric (parasitic) perturbations with amplitudes smaller than that of the channel itself.
The above dynamical picture can be dramatically altered if there is a competition between axisymmetric and transient nonaxisymmetric growths from the early linear stage of evolution during the characteristic e-folding time, tc, of the most unstable axisymmetric MRI (since a typical time for the development of the non-linear channel is of the order of tc). If the initial amplitude of non-axisymmetric perturbations is not small, their transient growth comes into play and these perturbations may become a key ingredient of the overall dynamics. We compared the amplification factor of the magnetic energy of non-axisymmetric magnetic mode SFHs with various radial, azimuthal and vertical wavenumbers with that of the most unstable axisymmetric magnetic mode SFH during tc. From this comparative analysis it turned out that at vertical wavelength larger than about the critical wavelength of the axisymmetric MRI (i.e., for kz < kz,cr, see Fig. 10), the growth factor of nonaxisymmetric magnetic mode SFHs is larger than that of the most unstable axisymmetric ones. In other words, the transient growth of these optimal and close-to-optimal non-axisymmetric SFHs exceeds the maximal exponential growth of axisymmetric SFHs over the dynamical time. When forming a complete dynamical picture, one should also take into account the fact that the set of the closeto-optimal SFHs is substantially larger than the set of the exponentially growing axisymmetric ones, as evident from Fig. 9. As a result, axisymmetric SFHs can be overwhelmed by the close-tooptimal SFHs if they have sufficiently large initial amplitudes. In this case, non-linear interaction between non-axisymmetric and axisymmetric SFHs will start just after transient amplification, leading directly to turbulent state, without/before the formation of the channel flow. In other words, non-axisymmetric magnetic pertur-bations can prevent the development of the channel in its early growth stage. This course of events may be another reason for the absence of channels in some simulations of the MRI-driven turbulence with non-zero net vertical magnetic flux configurations (e.g., Hawley et al. 1995;Bodo et al. 2008).
Our investigation also demonstrated that amplification factor of non-axisymetric magnetic mode SFHs increases with decreasing vertical wavenumber and reaches its maximal value at vertical and azimuthal wavelengths comparable to the disc scale-height (see Figs. 6 and 10). In weakly magnetized discs, the axisymmetric MRI is essentially incompressible, as the largest growth rate comes at wavelength smaller than the disc scale-height. By contrast, according to our incompressible calculations, the largest values of the non-axisymmetric amplification factor are achieved at vertical wavelengths of perturbations comparable to the disc scale-height for which compressible wave modes become inevitably involved (via mode linear coupling) in the dynamics . Consequently, for a realistic analysis of the non-axisymmetric MRI in discs, one should take into account compressibility. In support of this, we would like to mention that the saturation level of MRI-driven turbulence, involving non-axisymmetric Fourier harmonics, in numerical simulations with box sizes of the order of scale-height does depend on compressibility (Sano et al. 2004;Pessah et al. 2007;Guan et al. 2009). In addition, our preliminary analysis indicates that the linear coupling of spiral density waves and non-axisymmetric MRI occurs in the compressible case. Due to this coupling, non-axisymmetric magnetic mode with azimuthal and vertical wavelength comparable to the disc scale-height can effectively excite spiral density waves. This is of the same nature as the generation of high-frequency density waves by aperiodic (non-oscillatory) modes in non-magnetized discs (Bodo et al. 2005;Mamatsashvili & Chagelishvili 2007;Mamatsashvili & Rice 2011) and in MRI-driven turbulence (Heinemann & Papaloizou 2009b). However, in the latter paper, a set up with zero net magnetic flux was considered, where the linear MRI was apparently absent (or at least was weak) and density waves were primarily generated only by the vortical mode. The situation can be different in the case of non-zero net vertical magnetic flux we considered, where the main mode generating density waves is not the vortical but the magnetic mode.