Viscous, Resistive Magnetorotational Modes
We carry out a comprehensive analysis of the behavior of the magnetorotational instability (MRI) in viscous, resistive plasmas. We find exact, non-linear solutions of the non-ideal magnetohydrodynamic (MHD) equations describing the local dynamics of an incompressible, differentially rotating background threaded by a vertical magnetic field when disturbances with wavenumbers perpendicular to the shear are considered. We provide a geometrical description of these viscous, resistive MRI modes and show how their physical structure is modified as a function of the Reynolds and magnetic Reynolds numbers. We demonstrate that when finite dissipative effects are considered, velocity and magnetic field disturbances are no longer orthogonal (as it is the case in the ideal MHD limit) unless the magnetic Prandtl number is unity. We generalize previous results found in the ideal limit and show that a series of key properties of the mean Reynolds and Maxwell stresses also hold for the viscous, resistive MRI. In particular, we show that the Reynolds stress is always positive and the Maxwell stress is always negative. Therefore, even in the presence of viscosity and resistivity, the total mean angular momentum transport is always directed outwards. We also find that, for any combination of the Reynolds and magnetic Reynolds numbers, magnetic disturbances dominate both the energetics and the transport of angular momentum and that the total mean energy density is an upper bound for the total mean stress responsible for angular momentum transport. The ratios between the Maxwell and Reynolds stresses and between magnetic and kinetic energy densities increase with decreasing Reynolds numbers for any magnetic Reynolds number; the lowest limit of both ratios is reached in the ideal MHD regime.
Subject headings:black hole physics — accretion, accretion disks — MHD — instability — turbulence
The magnetorotational instability (MRI, Balbus & Hawley, 1991, 1998) has been widely studied in the inviscid and perfectly conducting, magnetohydrodynamic (MHD) limit. The departures from this idealized situation are usually parametrized according to the Reynolds and magnetic Reynolds numbers, where and stand for the relevant characteristic velocity and lengthscale and and stand for the kinematic viscosity and resistivity. The ideal MHD regime is then formally identified with the limit . There are many situations of interest in which the effects of dissipation need to be considered.
From the astrophysical point of view, accretion disks around young stellar objects constitute one of the most compelling reasons for investigating the MRI beyond the ideal limit. In particular, there is great interest in understanding to what extent can MHD turbulence driven by the MRI enable efficient angular momentum transport in cool, poorly conducting, protoplanetary disks (see, e.g. Blaes & Balbus, 1994; Jin, 1996; Gammie, 1996; Sano & Miyama, 1999; Salmeron & Wardle, 2005). Most of the studies addressing the effects of dissipation in non-ideal MRI have usually focused in inviscid, resistive plasmas. However, accretion disks are characterized by a wide range of magnetic Prandtl numbers, with varying by several orders of magnitude across the entire disk (see, e.g., Balbus & Henri, 2007). In order to understand the behavior of the MRI under these conditions it is necessary to relax the assumption of an inviscid plasma.
A large fraction of the shearing box simulations addressing the non-linear regime of the MRI have been carried out in the ideal MHD limit, i.e., without including explicit dissipation in the codes (see, e.g., Hawley, Gammie, & Balbus, 1995; Brandenburg et al., 1995; Sano et al., 2004). However, even in the absence of explicit viscosity and resistivity, finite difference discretization leads to numerical diffusion/dissipation. Therefore, even in this type of simulations, it is necessary to understand the impact of these numerical artifacts that lead to departures from the ideal MHD regime and how similar they are when compared with physical (resolved) dissipation.
A handful of numerical studies with explicit resistivity but zero physical viscosity have been carried out in order to understand the effects of ohmic dissipation in the saturation of MRI-driven turbulence (see, e.g., Sano, Inutsuka, & Miyama, 1998; Sano & Inutsuka, 2001; Fleming, Stone, & Hawley, 2000; Sano et al., 2004; Turner, Sano, & Dziourkevitch, 2007). In particular, Sano & Stone (2003) have shown that the saturation level of the stresses increases with increasing magnetic Reynolds number and seem to converge to an asymptotic value for magnetic Reynolds numbers larger than unity.
Recent work has pointed out problems with convergence in zero-net-flux numerical simulations of ideal MHD driven by the MRI (Pessah, Chan, & Psaltis, 2007; Fromang & Papaloizou, 2007) implying the necessity of incorporating explicit dissipation in the codes. Numerical studies with both resistivity and viscosity, in the presence of a mean vertical magnetic field (Lesur & Longaretti, 2007) and in the case of zero net flux (Fromang et al., 2007), have begun to uncover how the characteristics of fully developed MRI-driven turbulence depends on the Reynolds and magnetic Reynolds numbers. Even though the ranges in Reynolds and magnetic Reynolds numbers that can be currently addressed is still limited, the results obtained from the simulations suggest that the magnetic and kinetic energies contained in turbulent motions in the saturated regime depend on the values of the microphysical viscosity and resistivity. In particular, the mean angular momentum transport in the turbulent state increases with increasing magnetic Prandtl number.
From the experimental perspective, understanding the effects of non-vanishing resistivity and viscosity in the behavior of the MRI seems imperative, since the physical conditions achievable in the laboratory depart significantly from the ideal MHD regime (Ji, Goodman & Kageyama, 2001; Goodman & Ji, 2002; Sisan et al., 2004; Liu, Goodman, & Ji, 2006; Rüdiger, Schultz, Shalybkov, 2003). Liquid metals (such as sodium, gallium, and mercury) are often characterized by rather low magnetic Prandtl numbers (–). Although the regime of Reynolds numbers involved is still orders of magnitude smaller than any astrophysical system with similar magnetic Prandtl numbers, MRI experiments offer one of the few prospects of studying anything close to MHD astrophysical processes in the laboratory.
A number of analyses addressing some aspects of the impact of viscosity and resistivity on the MRI in various dissipative limits appear scattered throughout the literature on theoretical, numerical, and experimental MRI. More recently, Lesaffre & Balbus (2007) have found particular solutions of the viscous, resistive MHD equations (including even a cooling term) in the shearing box approximation. However, we are unaware of any comprehensive, systematic study addressing how the MRI behaves in viscous, resistive, differentially rotating magnetized plasmas for arbitrary combinations of the Reynolds and magnetic Reynolds numbers. The aim of this work is to carry out this analysis in detail.
The rest of paper is organized as follows. In § 2, we state our assumptions. In § 3, we solve the eigenvalue problem defined by the MRI for arbitrary Reynolds and magnetic Reynolds numbers. We provide closed analytical expressions for the eigenfrequencies and the associated eigenvectors. In § 4, we address the unexplored physical structure of MRI modes for finite Reynolds and magnetic Reynolds numbers and derive simple analytical expressions that describe these modes in various asymptotic regimes. In § 6, we calculate the correlations between magnetic and velocity MRI-driven perturbations that are related to angular momentum transport and energy densities. We find that some key results previously shown to hold in the ideal MHD limit (Pessah, Chan, & Psaltis, 2006) are also valid in the non-ideal regime. In particular, we show that even though the effectiveness with which the MRI disrupts the laminar flow depends on the Reynolds and magnetic Reynolds numbers, the instability always transports angular momentum outwards. We also find that magnetic perturbations dominate both the energetics and the transport of angular momentum for any combination of the Reynolds and magnetic Reynolds numbers. In § 7 we summarize our findings and discuss the implications of our study.
Let us consider a cylindrical, incompressible background characterized by an angular velocity profile , threaded by a vertical magnetic field . We work in the shearing box approximation, which consist of a first order expansion in the variable of all the quantities characterizing the flow at the fiducial radius . The goal of this expansion is to retain the most relevant terms governing the dynamics of the MHD fluid in a locally-Cartesian coordinate system co-orbiting and corrotating with the background flow with local (Eulerian) velocity . (For a more detailed discussion on this expansion see Goodman & Xu 1994 and references therein.)
The equations governing the dynamics of an incompressible MHD fluid with constant kinematic viscosity and resistivity in the shearing box limit are given by
where is the pressure, is the (constant) density, the factor
parametrizes the magnitude of the local shear, and we have defined the (locally-Cartesian) differential operator
where , , and are, coordinate-independent, orthonormal vectors corrotating with the background flow at . The continuity equation reduces to and there is no need for an equation of state since the pressure can be determined from this condition.
We focus our attention on the dynamics of perturbations that depend only on the vertical coordinate. Under the current set of assumptions, these types of perturbations are known to exhibit the fastest growth rates in the ideal MHD case (Balbus & Hawley, 1992, 1998; Pessah & Psaltis, 2005). The equations governing the dynamics of these perturbations can be obtained by noting that the velocity and magnetic fields given by
where the time dependence is implicit, constitute a family of exact, non-linear, solutions to the viscous, resistive MHD equations (1)-(2). As noted in Goodman & Xu (1994), even in the dissipative case, the only non-linear terms, which are present through the perturbed magnetic energy density, are irrelevant in the case under consideration.
We can further simplify equations (1) and (2) by removing the background shear flow111In the shearing box approximation, the dependence of the background flow on the radial coordinate is strictly linear and therefore viscous dissipation does not affect its dynamics. and by realizing that we can take without loss of generality. We then obtain
where the first term on the right hand side of equation (8) is related to the epicyclic frequency
at which the flow variables oscillate in a perturbed hydrodynamic disk. For Keplerian rotation the parameter is and thus the epicyclic frequency is .
It is convenient to define the new variables for , and introduce dimensionless quantities by considering the characteristic time- and length-scales set by and . The equations satisfied by the dimensionless perturbations, , , are then given by
where and denote the dimensionless time and vertical coordinate, respectively.
The dynamics of ideal MRI modes, with , is completely determined by the dimensionless shear (Pessah, Chan, & Psaltis, 2006). The effects of viscous and resistive dissipation introduce two new dimensionless quantities that alter the characteristics and evolution of the MRI. With our choice of characteristic scales, it is natural to define the Reynolds and magnetic Reynolds numbers characterizing the MHD flow as 222In a compressible fluid, the sound speed, , provides another natural characteristic speed to define the Reynolds and magnetic Reynolds numbers. These definitions, e.g., and , are related to those provided in equations (16) and (17) via the plasma beta parameter , with the Alfvén speed in the direction, simply by and for a polytropic equation of state , with and constants.
with associated magnetic Prandtl number
In order to simplify the notation, we drop hereafter the tilde denoting the dimensionless quantities. In the rest of the paper, all the variables are to be regarded as dimensionless, unless otherwise specified.
3. The Eigenvalue Problem for the Non-ideal MRI :
A Formal Analytical Solution
In this section we provide a complete analytical solution to the set of equations (12)–(15) as a function of the shear parameter , (or, equivalently, the epicyclic frequency, ) for any set of values defining the viscosity and resistivity.
It is convenient to work in Fourier space, as this provides the advantage of obtaining explicitly the basis of modes that is needed to construct the most general solution satisfying equations (12)–(15). Taking the Fourier transform of this set with respect to the -coordinate, we obtain the matrix equation
where the vector stands for
and represents the matrix
The functions denoted by correspond to the Fourier transform of the real functions, , and are defined via
where we have assumed periodic boundary conditions at , with being the (dimensionless) scale-height and the wavenumber in the -coordinate,
where is an integer number. In order to simplify the notation, hereafter we denote these wavenumbers simply by .
In order to solve the matrix equation (19), it is convenient to find the eigenvector basis, with , in which is diagonal. This basis exists for all values of the wavenumber (i.e., the rank of the matrix is equal to 4, the dimension of the complex space) with the possible exception of a finite number of values of . In this basis, the action of over the set is equivalent to a scalar multiplication, i.e.,
where are complex scalars.
In the eigenvector basis, the matrix has a diagonal representation = diag. The eigenvalues , with , are the roots of the characteristic polynomial associated with , i.e., the dispersion relation associated with the non-ideal MRI, which can be written in compact form as
where we have defined the quantities
The dispersion relation (25) is a fourth order polynomial with non-zero coefficients in and . In order to find its roots it is convenient to take this polynomial to its depressed form. This can be achieved by defining the new variables and such that333A physical interpretation of the variable is provided in § 3.3.
The resulting polynomial can then be written as
where the coefficients , , and are given by
The solutions to equation (30) are
where the subscripts and in the “” and “” signs label the four possible combinations of signs and we have defined the quantities444Defining the quantities and in this way allows us to show explicitly that in the limit the solutions to equation (25) converge smoothly to the solutions found in the ideal MHD case (see Appendix A).
and is any of the solutions to the cubic equation
which has closed analytic solutions
Note that the choice of either sign in is immaterial.
3.2. Four classes of solutions
All of the quantities , , and , depend on the viscosity and the resistivity only through . This has a series of important implications, in particular, there is always a range of wavenumbers for which the discriminant in equation (43) is positive, i.e., . It can also be seen that the last two terms between parentheses in equation (43) are always positive, i.e., , and thus they always produce damping. Because of this, we can classify the modes in four types: two (damped) growing and decaying “unstable” modes with eigenvalues
and two (damped) “oscillatory” modes with eigenvalues
We arbitrarily label these eigenvalues as
Figure 1 shows the growth rate as a function of the vertical wavenumber for different combinations of the Reynolds and magnetic Reynolds numbers for Keplerian rotation. These growth rates are more sensitive to changes in the resistivity than to changes in the viscosity. A qualitative understanding of this behavior can be obtained by realizing that viscosity tends to quench the instability, without altering the large scale magnetic field. Thus, as long as the resistivity is negligible, the range of unstable lenghtscales are the same in both ideal and viscous, perfectly conducting fluids. On the other hand, resistivity tends to destroy the magnetic field at small scales having a stronger impact on the stability of the perturbations at these scales.
Mathematically, the asymmetric response of the growth rate to changes in the viscosity or the resistivity originates in the different functional form of the terms that contribute to produce damping, i.e., and , in the exponential growth characterized by in equation (44). If the oscillatory modes, in equation (45), are considered instead, the roles of the plus and minus signs in these terms are interchanged. From this analysis we can infer that the “oscillatory” mode is affected (damped) more strongly by viscosity than by resistivity.
The simultaneous analysis of the various panels in Figure 1 leads to the conclusion that viscous, resistive unstable modes with magnetic Prandtl number equal to unity resemble more closely inviscid, resistive modes rather than viscous, conductive ones. In § 4 we provide analytical expressions to support this conclusion.
3.3. Normalized Eigenvectors: Geometrical Representation
The set of normalized eigenvectors, , associated with the eigenvalues (46) are given by
, , and the norms are given by
Here, is the -th component of the (unnormalized) eigenvector associated with the eigenvalue .
The set of four eigenvectors , together with the set of complex scalars in equation (46), constitute the full solution to the eigenvalue problem defined by the MRI for any combination of the Reynolds and magnetic Reynolds numbers.
A geometrical representation of the eigenvectors (47) can be brought to light by defining the angles and according to
It is important to remark that each of the four eigenvectors define, in principle, four sets of angles , for . We label the angles associated with the different types of modes discussed in § 3.2 according to
with similar definitions corresponding to for . Note that these angles are defined in spectral space and depend, in general, on the wavenumber , the epicyclic frequency, , the viscosity , and the resistivity . The angles associated with the modes labeled by and are always real while the ones associated with the modes and are in general complex. For the sake of brevity, in what follows we will refer to the set of angles describing unstable MRI modes simply as .
A normalized version of the MRI eigenvectors can be obtained by multiplying the set of vectors in equation (48) by the amplitudes
where we have defined
where, for the sake of simplicity, we have omitted the subscript on the left hand side. The expressions for the normalized eigenvectors , for , are then given by
3.4. Temporal Evolution
evolves in time according to
with and , for , given by equations (46) and (55). The initial conditions are related to the initial spectrum of perturbations, , via . Here, is the matrix for the change of coordinates from the standard basis to the normalized eigenvector basis555The eigenvectors (55) are not in general orthogonal, i.e., for . If desired, an orthogonal basis can be constructed using the Gram-Schmidt orthogonalization procedure (see, e.g., Hoffman & Kunze, 1971). and can be obtained by calculating the inverse of the matrix
The temporal evolution of a single MRI-unstable mode in physical space can be obtained from a linear combination of and as defined in equation (55). In particular, setting in equation (60) and substituting the result in equation (61) we obtain
These solutions are of particular importance for the linear late-time evolution of MRI modes. Note that any reasonable spectrum of initial perturbations of the type used in numerical simulations of shearing boxes will have a non-zero component along the unstable eigenvector . If the value of the magnetic field is such that the MRI can be excited for given values of the viscosity and resistivity then the exponentially growing perturbations in physical space will evolve towards a mode of the form (63) dominated by the lengthscale for which the growth rate reaches its maximum value .
Note that if a perturbation in physical space is composed by a single mode of the type described in § 3.2, no matter which class, then the angles defined in equations (50) and (51) are constant in time and are identical to the physical angles between the planes containing magnetic and velocity perturbations in physical space, see Figure 2, with
Finally, defining the angle such that
which implies that , and using the fact that
it is not difficult to show that
This means that provides a measure of how non-orthogonal velocity and magnetic field perturbations are.
It is evident that when the magnetic Prandtl number approaches unity viscous, resistive, MRI-driven magnetic and velocity perturbations tend to be orthogonal, i.e., , and therefore
for every wavenumber . This is illustrated in Figure 2 which shows the evolution of the angles and corresponding to the most unstable MRI mode as a function of the Reynolds/magnetic Reynolds number when the magnetic Prandtl number is equal to unity.
4. Physical Structure of MRI Modes
The evolution of the physical structure of a single growing MRI mode with wavenumber is characterized by its growth rate , the relative magnitude between the amplitudes of magnetic and velocity field perturbations, , and the two angles defining the planes containing them, and . For any reasonable spectrum of initial perturbations the mode that exhibits the fastest exponential growth, , which we refer to as , will dominate the dynamics of the late time evolution of the viscous, resistive MRI. It is therefore of particular interest to characterize the physical properties of this fastest growing mode in different dissipative regimes.
4.1. Marginal and Fastest Growing MRI-modes
Because the eigenvalue associated with the unstable growing mode, , is always real for any combination of the Reynolds and magnetic Reynolds numbers, it is possible to find the marginally stable mode such that . Setting in equation (25), we obtain a polynomial in
valid for any value of the viscosity and resistivity. Note that sets the minimum domain height for numerical simulations of viscous, resistive MRI-driven turbulence. Figure 3 shows the solutions of equation (70) in various dissipative regimes for Keplerian rotation. The analytic solutions of equation (70) are algebraically complicated but their asymptotic limits are rather simple. We find expressions for this critical wavenumber in several regimes of interest below.
In the ideal MHD limit, it is straightforward to find simple analytical expressions for the most unstable wavenumber, , and its associated growth rate, . However, the analytical expressions that we derived for the eigenfrequencies in the non-ideal case, equation (46), are not amenable to the usual extremization procedure. More precisely, it is very challenging to find the values of and that satisfy
Another possible path to find the values of the wavenumber , and the associated growth rate, is to use the fact that satisfies simultaneously the dispersion relation (25) and its derivative to eliminate between these two and obtain a polynomial in . The largest of the roots of this polynomial is the desired maximum growth rate. It is possible to find following a similar methodology, but eliminating between the two polynomials instead. However, for arbitrary values of the viscosity and resistivity, both procedures lead to a seventh degree polynomial whose roots must be found numerically, defeating altogether the attempt to find analytical expressions for and .
Using as a guide the results shown in Figures 4 and 5, we follow an alternative procedure. The goal is to find simple analytical expressions to describe the asymptotic behavior of the most unstable mode, , and the maximum growth rate, , in different dissipative regimes. It is evident from Figure 1 that and for all the non-ideal MRI modes. This information can be used to simplify the dispersion relation and its derivative so as to decrease their order without loosing vital information. This makes it possible to obtain manageable, but accurate, expressions for and in different limiting regimes.
Figure 6 shows contour plots for the critical wavenumber, , the most unstable wavenumber, , and the maximum growth rate, , as a function of the Reynolds and magnetic Reynolds numbers for Keplerian rotation. In all three panels, lighter gray areas correspond to larger values of , , and , respectively. Note that in all the cases, the functional form of the contours naturally divides the plane in three distinctive regions that we denote according to (ideal), (resistive), and (viscous). Note that when the most unstable wavenumber, , and the maximum growth rate, , are considered, these regions can be associated with the regions where , , and , respectively. The overlap between these regions is not as clear when the critical wavenumber is considered and some care is needed when deriving approximated expressions for it.
4.2. Ideal MRI Modes
Let us first demonstrate briefly how the formalism presented in § 3 reduces to previously known results in the ideal MHD limit. In the absence of dissipation, the eigenvalues , with , are the roots of the dispersion relation associated with the ideal MRI (Balbus & Hawley, 1991, 1998),
and are given by (Pessah, Chan, & Psaltis, 2006)
where we have defined the quantities and such that