Galactic MRI with Braginskii viscosity

Global MRI with Braginskii viscosity in a galactic profile


We present a global-in-radius linear analysis of the axisymmetric magnetorotational instability (MRI) in a collisional magnetized plasma with Braginskii viscosity. For a galactic angular velocity profile we obtain analytic solutions for three magnetic field orientations: purely azimuthal, purely vertical and slightly pitched (almost azimuthal). In the first two cases the Braginskii viscosity damps otherwise neutrally stable modes, and reduces the growth rate of the MRI respectively. In the final case the Braginskii viscosity makes the MRI up to times faster than its inviscid counterpart, even for asymptotically small pitch angles. We investigate the transition between the Lorentz-force-dominated and the Braginskii viscosity-dominated regimes in terms of a parameter where is the viscous coefficient and the Alfvén speed. In the limit where the parameter is small and large respectively we recover the inviscid MRI and the magnetoviscous instability (MVI). We obtain asymptotic expressions for the approach to these limits, and find the Braginskii viscosity can magnify the effects of azimuthal hoop tension (the growth rate becomes complex) by over an order of magnitude. We discuss the relevance of our results to the local approximation, galaxies and other magnetized astrophysical plasmas. Our results should prove useful for benchmarking codes in global geometries.

instabilities – accretion, accretion discs – Galaxy: disc – MHD – magnetic fields – plasmas.

1 Introduction

In the formation of compact objects (stars, planets and black holes) from accretion discs, turbulence driven by the MRI, and possibly the MVI, offers a promising mechanism for the necessary angular momentum transport (Velikhov, 1959; Chandrasekhar, 1960; Balbus & Hawley, 1991; Balbus, 2004). It has also been suggested that the observed velocity fluctuations in parts of the interstellar medium (ISM) with low star formation rates may, in part, arise from this process (Sellwood & Balbus, 1999; Tamburro et al., 2009). The evidence for this comes primarily from numerical simulations and a wide range of studies agree that weak magnetic fields and outwardly decreasing angular velocity profiles are an unstable combination (Balbus & Hawley, 1998; Balbus, 2003)

The most illuminating explanation for this comes from a shearing sheet analysis in which the mean azimuthal flow of the differentially rotating disc is locally approximated by a constant angular velocity rotation plus a linear velocity shear (Goldreich & Lynden-Bell, 1965; Umurhan & Regev, 2004). In the simplest possible setup, incompressible, isothermal, dissipationless, axisymmetric linear perturbations to a magnetic field with a weak vertical component, i.e. parallel to the rotation axis of the disc, are unstable when the angular velocity decreases away from the disc’s centre. Azimuthal velocity perturbations to fluid elements at different heights, tethered to each other by the magnetic field, increase (decrease) their angular momentum. This causes them to move to larger (smaller) radii as dictated by the gravitational field which sets the mean flow. In the frame rotating at constant angular velocity, this motion deforms the tethering magnetic field, provided it is not too strong, and this induces a prograde (retrograde) Lorentz force on the outer (inner) element thus destabilising the system as it moves to yet larger (smaller) radii. This mechanism is at the heart of the MRI.

However, although this model and description captures much of the essential physics, to fully understand the MRI, or at least the framework within which the shearing sheet should exist, a more nuanced approach is needed. In part, this is because the shearing sheet is formulated in a Cartesian coordinate system where curvature terms that arise naturally from the cylindrical geometry of the accretion disc are neglected. Indeed, in the local approximation, the over-stabilising effects of hoop tension (a curvature effect) associated with the azimuthal magnetic field, are totally ignored. Furthermore, the model predicts that the fastest growing, and therefore most physically relevant, linear MRI modes have a homogeneous radial structure on the scale of the shearing sheet in which the local approximation is made ( Guan et al. (2009) has shown the MRI is well localized in the non-linear regime). This means that the global disc structure, including boundary conditions, not captured by the local approximation may have a significant effect on these large scale modes in a way that cannot be determined locally. Other limitations exist too (Knobloch, 1992; Regev & Umurhan, 2008).

The extent to which these limitations matter should, and under a variety of assumptions have, been investigated by global analyses that take into account the full radial structure of the disc and its boundaries (Dubrulle & Knobloch, 1993; Curry et al., 1994; Curry & Pudritz, 1995; Ogilvie & Pringle, 1996; Ogilvie, 1998), and specifically for the galaxy by Kitchatinov & Rüdiger (2004). The conclusions of these investigations largely confirms the local picture of a large radial scale instability driven by the differential rotation of the disc. This suggests that whilst a local MRI analysis is generally correct, its regime of validity must be checked globally. It is the purpose of this work to do just that for the MRI operating in a collisional, magnetized (the ion cyclotron frequency ion-ion collision frequency ) plasma (like the ISM).

In such a plasma Braginskii (1965) has shown that to lowest order in , the deviatoric stress tensor is diagonal and anisotropic. This leads to different parallel and perpendicular viscosities or, more fundamentally, pressures (and thermal conductivities3) with respect to the local direction of the magnetic field. Of the important physical consequences of this, it will be the effect of the Braginskii viscosity in the presence of a galactic shear flow that will concern us here.

This is not a new topic and in recent years the study of magnetized accretion discs has attracted attention in both the collisionless (Quataert et al., 2002; Sharma et al., 2003; Sharma et al., 2006, 2007) and collisional regimes (Balbus, 2004; Islam & Balbus, 2005; Ferraro, 2007; Devlen & Pekünlü, 2010), and a well developed code to simulate them now exists (Parrish & Stone, 2007; Stone et al., 2008). However, a number of fundamental questions remain unanswered. Primarily, what is the non-linear fate of the MRI in a collisional magnetized plasma? Does it transport angular momentum and if so, is the transporting stress primarily viscous, Maxwell or Reynolds4? On what scales do the most unstable modes emerge and how does this vary with the system parameters? In what regime will a local analysis become untenable and global effects (either radial or vertical) become important (Gammie & Balbus, 1994)? What effect does the presence, or absence, of a net vertical field have given Cowling’s anti-dynamo theorem and the dissipative properties of the Braginskii viscosity (Moffatt, 1978; Lyutikov, 2007)? Could viscous heating from the Braginskii viscosity lead to secondary magnetized or unmagnetized thermal instabilities (Balbus, 2001; Quataert, 2008; Kunz et al., 2011)? Do channel solutions, or something approaching them, emerge (Goodman & Xu, 1994)? If they do, the associated field growth will generate pressure anisotropies that could feed new parasitic instabilities such as the mirror. What would their effects be at this stage and in the inevitable turbulence where the mirror and firehose instabilities will both arise (Schekochihin et al., 2005)?

Addressing these questions will require a two-pronged approach involving both numerical and analytic studies. It may transpire that much of the existing work on the unmagnetized and collisionless MRI is directly applicable, but this needs determining. As such we conduct a global linear stability analysis for three separate background magnetic field orientations: purely azimuthal, purely vertical and pitched (magnetic field lines follow helical paths on cylinders of constant radius). We embed these in a galactic rotation profile. In agreement with earlier, local studies, we find that when the field has both a vertical and an azimuthal component, a linear instability with a real part up to times faster than the MRI emerges (Balbus, 2004; Islam & Balbus, 2005). In contrast to local studies we also find that it has a travelling wave component, and its growth rate depends on the viscous coefficient.

In the presence of a vertical field, we also recover the standard inviscid MRI and show that upon introducing the Braginskii viscosity, its growth rate is always reduced. A similar effect is found for a purely azimuthal field where we find that the Braginskii viscosity damps modes that are, inviscidly, neutrally stable.

The layout of this paper is as follows. In Section 2, we introduce and perturb a series of equilibrium solutions to the Braginksii-MHD equations and this forms the basis of our global stability analysis. Relegating the manipulation of the perturbed equations to Appendix A, we obtain a single ODE governing the perturbed modes. We proceed by solving this for azimuthal, vertical and pitched field orientations in Sections 3, 4 and 5 respectively. For each case we contrast the behaviour with and without Braginskii viscosity. In Section 6 we discuss the physical mechanism of the instability, where it may occur astrophysically, and the relation of our results to the local approximation. Finally, we conclude in Section 7 with some thoughts on open questions relating to magnetized astrophysical plasmas.

2 Global Stability Analysis

2.1 Governing equations

The simplest set of equations that capture the physics of the collisional magnetized MRI are those those of isothermal ideal MHD with the Braginskii viscosity (Lifshitz et al., 1984). Explicitly these equations are the momentum equation, including the Braginskii stress tensor; the induction equation that describes the evolution of the magnetic field; the incompressibility condition (because the perturbations are linear, the Mach number can always be made small enough to ensure this); and the stress tensor itself:

T (4)

where the constant mass density has been scaled out of the problem, is the velocity field, is the mass-density scaled magnetic field, i.e. the Alfvén velocity, is the direction of the magnetic field, is the gravitational potential, is the total (gas plus magnetic) pressure, ‘:’ is the full inner product, I the identity, and T is the full Braginskii stress tensor whose form we now explain.

In a magnetized plasma the total pressure tensor (whose divergence appears in the momentum equation) is given by

where are the perpendicular and parallel scalar pressures with respect to the magnetic field direction. When the plasma is also collisional, the pressure anisotropy can be related to the rate-of-strain by a Chapman-Enskog style perturbation theory (Braginskii, 1965; Chapman & Cowling, 1970).

Microphysically, the anisotropy arises from the conservation of the magnetic moment (first adiabatic invariant) of a gyrating particle in a magnetic field. However, collisions break this conservation and relax the anisotropy by pitch-angle scattering particles in velocity space (this dissipative process will turn out to be important in isolated field configurations). The competition between these two processes is governed by

where we have used the BGK operator to approximate the full collision operator.

In the presence of time-varying magnetic fields, the pressure anisotropy tends to a steady state that tracks the fields’ rate of change. It follows

where is the coefficient of the Braginskii viscosity, and we have used equation (2) in the final equality. Equation (4) follows directly.

2.2 Equilibrium solutions

Working in cylindrical polar coordinates , we introduce equilibrium solutions to equations (1)-(4) that describe a differentially rotating global shear flow constrained by gravity and threaded by a magnetic field that lies on cylinders of constant radius. Our equilibrium solutions are uniform in , and the plasma motion is restricted at an inner boundary of finite radius . (The validity of these assumptions with be discussed in Section 6.2.)

We allow gravity to dictate the rotation profile of the equilibrium flow where is the rotation frequency at the inner boundary and is a dimensionless measure of the shear. Of interest to us here is the case of that is both analytically treatable and physically corresponds to a galactic disc where, unlike the Keplerian case of , the gravitational potential of the dark matter halo sets the rotation profile (Rubin & Ford, 1970; Sofue & Rubin, 2001). In modelling this we set the dark matter mass distribution (whose sole purpose is to ultimately set the rotation profile) to a Mestel (1963) profile so, via Poisson’s equation, . In this case, the flow does not vary with radius.

We decompose the magnetic field into vertical and azimuthal components, (so as to construct a time independent equilibrium, we neglect radial magnetic fields5). That is, where is the pitch angle of the magnetic field. We do not specify yet, however, for mathematical simplicity, we demand that both and are independent of radius so remains constant. Whilst this implies a vertical current is associated with , a simple super-galactic magnetic field can account for . In all, our equilibrium fields take the form


and are constant in space and time. As is physically relevant, we restrict .

It is a mathematically convenient feature of our equilibrium solutions that there is no evolution of the magnetic field strength. From equations (4) and (5), T is absent from the unperturbed state and the system will be stable to pressure anisotropy driven microscale instabilities, e.g. firehose and mirror (Schekochihin et al., 2005). In contrast to the case where Laplacian viscosity (or indeed resistivity) is present, we can construct an ideal MHD solution independent of any radial flows (Kersalé et al., 2004).

2.3 Perturbed equations

Figure 1: Radial mode structure of the fastest growing branch MRI modes with and . The (purely real) growth rates are and associated wavenumber .

To determine the stability of this system we linearise equations (1)–(4) about (5) with axisymmetric velocity perturbations and similarly for the magnetic and pressure fields. Here is the wavenumber in the z direction and the growth rate. We retain curvilinear terms from the cylindrical geometry but neglect self-gravity. We non-dimensionalise with respect to time-scales and length-scales .

As detailed in Appendix A, the linearised equations combine into a single complex second order ordinary differential equation for , the Modified Bessel equation





Here is the vertical Alfvén frequency and ( in dimensional form) is a measure of the relative sizes of the anisotropic pressure force and Lorenz force6. Here is the number density, the ion mass and the plasma-beta.

To obtain solutions to equation (6), we choose our boundary conditions to be an impenetrable wall at so and the requirement decays at infinity (see Furukawa et al. (2007) for an inviscid treatment that includes the shear singularity at the origin).

In this case, the eigenfunctions of equation (6) are modified Bessel functions of the second kind with argument and order (Watson, 1944; Abramowitz & Stegun, 1964). The spectrum of solutions is discrete and we index with , . In the special case when (or from equation (7)) is real (in general it is complex) the problem is Sturm-Liouville and is an infinite ordered set of eigenvalues . To determine and therefore , it is necessary to solve


From solutions to this and equations (7) and (8), the full set of flow, magnetic and pressure fields can be constructed, Appendix A. In Fig. 1 we show the radial structure of when .

In general, determining must be done numerically. However, for the most physically relevant magnetic field configurations, the problem becomes, in part, analytically tractable. These cases, a purely azimuthal field (in the galactic plane), a purely vertical field (a super-galactic field) and a slightly pitched field (very slightly out of the galactic plane), exhibit categorically different behaviour such that are singular limits.

In the first case, linear perturbations are damped; in the second, the system is unstable to the MRI but the Braginskii viscosity reduces the maximum growth rate below the Oort-A value maximum , (Balbus & Hawley, 1992); in the third, the system is unstable with , even for asymptotically small (Balbus, 2004). We present the details of these calculations now.

3 Azimuthal Field,

The stability of inviscid axisymmetric perturbations to a purely azimuthal magnetic field in the presence of a shear flow are well known. When and the magnetic field has the system is always stable (as ours is). When only one criterion is met, depending on the form of the fields, the system may still remain stable by the modified Rayleigh criterion (Rayleigh, 1916; Michael, 1954; Chandrasekhar, 1961). (It is worth noting however that global non-axisymmetric perturbations are unstable (Ogilvie & Pringle, 1996).)

Figure 2: Left panel: Growth rate of the n=0 (fastest) branch of the inviscid global MRI as a function of the vertical Alfvén frequency for several values of (time and length are non-dimensionalised with respect to and ). In the weak field limit, we recover the behaviour of the local MRI (dotted line). For all figures, the top-to-bottom order of the varied parameter corresponds to the top-to-bottom order of the main-panel curves. Right panel: Maximum growth rate as a function of . As asymptotes to the local Oort-A maximum (dotted line). Inset: The wavenumber at which occurs as a function of along with the local maximum (dotted line). Note: In the inviscid weak-field limit, two configurations with the same vertical magnetic field will exhibit identical behaviour – Section 5.3

Do these results for our inviscidly stable system, , persist in the presence of the Braginskii viscosity? The answer is no.

Setting in equations (6) to (8) and making a change of variables, , we obtain the simple expression



We multiply equation (10) by the complex conjugate of , and integrate between and infinity. Boundary terms vanish, so


where is non-negative, as is .

Equation (11) is a quadratic in whose roots depend crucially on . When (the inviscid limit), is purely imaginary (neutrally stable, travelling waves), whereas when the result is quite different. In this case and the only question is whether perturbations are purely damped, or damped and travelling.

It follows that if the system is stable in the absence of the Braginskii viscosity, it remains so in its presence.

4 Vertical Field,

4.1 Inviscid MRI,

For simplicity, we start by considering the inviscid limit of a purely vertical field , i.e. . In this case equations (6) to (8) simplify7:




If , has no nontrivial zeros and it follows that for an instability (hence our sign convention in equation (6)). This implies is bounded from above by .

Numerical solutions to equations (12) and (13) are obtained using Newton’s method (implemented in Mathematica). There is an unstable solution (and three stable ones) which is shown in Fig. 2. We find is real, positive, and of order the shear rate. The instability is the global MRI, and the branch has the largest growth rate.

The fastest growing modes have , so in the strong field regime the mode is large scale (physically, smaller scale modes are suppressed by magnetic tension) and is peaked away from the inner boundary. The result is that is reduced below the Oort-A maximum that occurs when and the mode is localised at (Curry et al., 1994).

In this weak field regime the problem is amenable to asymptotic analysis. It has been shown by Cochran (1965) and Ferreira & Sesma (1970, 2008) that in this limit, the zeros of are given by


where , is the modulus of the real negative zero of the Airy function (Abramowitz & Stegun, 1964) and omitted terms are of the form with . This result can be used to find the asymptotic-in- form of by inverting equation (13) to form a bi-quadratic


into which we substitute equation (14). Retaining the first two terms in , we solve exactly for


in which the positive root corresponds to the instability. Numerically we find the fastest growing mode (for a given ), and the wavenumber at which it occurs , obey .

To draw an analogy with the local approximation (see Section 6.3), that indexes the number of zeros in the domain is like a radial wavenumber. For large enough , the solutions are plane waves. To see this, we apply a WKB analysis to equation (6) using the small parameter (Dubrulle & Knobloch, 1993; Ogilvie, 1998). We have


Demanding that is real, this has solutions


and to ensure at , it must satisfy the boundary condition


which determines the spectrum of solutions.

Combining equations (13) and (19) we have


where and we have neglected factors of .

In the large limit, on small enough scales, the mode structure given by equation (18) is especially simple. Reverting to as our radial coordinate, we expand equation (18) about . To leading order we find


where and .

These solutions describes rapidly oscillating modes with frequency and a slowly varying amplitude . For sufficiently large , the solutions are plane waves whose growth rate decreases with .

Figure 3: Growth rate of the n=0 (fastest) branch of the global MRI with Braginskii viscosity for a two different weak, , field configurations and a range of . Left panel: Vertical Field: For fixed , decreases as increases and is bounded above by the inviscid growth rate (dotted line), but the critical wavenumber above which is independent of . Right panel: Slightly pitched field: In contrast to , for a given , the viscous case inviscid case (dotted), but is less than the local Lorentz-force free, or MVI, limit (dashed). As , but, unlike the vertical case,

4.2 Viscous MRI,

Allowing for viscosity whilst retaining a vertical field, equations (7) and (8) again reduce to a simple form:


Unlike the inviscid case, we can no longer guarantee the reality of as can, and indeed sometimes does, take complex values. The problem is generally not of Sturm-Liouville form and so the roots of are complex.

Numerically we again find four branches of which three are stable, and one unstable. The unstable mode (the only one of interest) is real, so if then . In this case, the problem is Sturm-Liouville. The unstable branch can be traced from the inviscid MRI (and is maximized for ), Figs. 1, 3 and 5.

Asymptotic expressions for and can be found using the results of Cochran (1965) and Ferreira & Sesma (2008) for the (complex) roots of equation (9):


where is the Euler constant.

For we combine equations (22) and (23) to get


where is the inverse of the Reynolds number. We expand equation (26) in and solve order by order. To lowest order we find

and to

and so


Numerically we find which, combined with equations (22) and (24), implies . It follows that to lowest order


and, differentiating equation (27) with respect to k and neglecting leading order logarithmic variations,


For , equation (25) applies and so is governed by

where we have used . We expand in so , and again solve order by order. To lowest order we find satisfies the inviscid equation


and at


To see the effect of variations in on we perform a WKB analysis in which we order so is a small parameter – see Section 4.1. The boundary condition for the viscous modes is then

Combining this with equation (23), we expand in so . To lowest order is given by the inviscid equation (20), and to


In the various limits, our asymptotics confirm the numerical results that , and . Contrasting the results of this section with the inviscid case when , we have

5 Pitched Field,

5.1 Ordering assumptions

When the magnetic field has both a vertical and an azimuthal component, the perturbed Braginskii stress tensor exerts an azimuthal ‘tension’ force on separating plasma elements. For arbitrary , the system is neither Sturm-Liouville (equation (6) is complex) nor is it amenable to the kind of polynomial inversion used in Section 4.

However, assuming matters simplify considerably. Physically, this choice of pitch angle represents the most realistic non-isolated galactic magnetic field configuration (Beck et al., 1996). Mathematically, as we now show, it is a singular limit that constitutes the stability threshold between actively damped modes, Section 3, and an unstable configurations that grows even faster than the inviscid global MRI (to which the following stability threshold also applies)

Here is a small parameter with respect to which we order the pitch angle and the remaining quantities in equation (6): and .

Letting we can retain only the first few terms in a series expansion of our trigonometric functions

To include their effects, we assume a strong magnetic field and Braginskii viscosity . To retain vertical magnetic tensions to order 1 (failing to do so removes the high wavenumber cutoff and leaves the equations ill-posed), we set . Then, in anticipation of unstable modes that grow at the shear rate, we order too. Balancing terms in equation (6), we find .

Figure 4: (if it exists), and the corresponding vertical Alfvén frequency for two field configurations. Left panel: Vertical field: and vs. for a range of . These are bounded above by the local limits of and (dotted lines), and well matched by the asymptotic results (dashed lines) given by equations (28) and (29). Right Panel: Slightly pitched, : and are bounded from above by the local limit of (main panel, top dotted line) and the inviscid limit (top inset, dotted line). The imaginary part of the growth rate (bottom inset) is well described in the small and large limits by equations (36) and (48). Note that (bottom curve, main panel) demonstrates the combined effects of hoop tension and the Braginskii viscosity.

In summary, our orderings are


We apply these scalings to equations (6) to (8) and retain the lowest order terms. We find is still governed by Bessel’s equation and, since , is determined by equation (12) and



We start by solving the inviscid problem.

5.2 Pitched inviscid MRI.

In the inviscid weak-field regime, we recover the vertical weak-field instability of Section 4.1 and Fig. 2. appears via the vertical Alfvén frequency only so so the unstable mode will be confined to a boundary layer of width .

If is small but finite (implicitly all orderings are subsidiary to equation (33)) the governing dispersion relation is the complex quartic8

Writing as a series in we find


with a maximum value


which occurs at , Fig. 4.

5.3 Pitched MRI with Braginskii viscosity,

In the presence of the Braginskii viscosity we find a complex instability whose real part exceeds that of its inviscid counterpart. Numerical solutions described by equations (12) and (34) are shown in Figs. 3, 4 and 5. When the growth rate tends to the shear rate and, as before, the fastest growing modes occur for .

Asymptotic expressions for and can be found by substituting equation (14) into equation (34). For , is governed by

Differentiating this with respect to ( and appear together only in one combination and so is treated as a single independent variable), and setting , the stationary points of obey:


By considering the branch that maximizes the stationary value of , and introducing , we find


and substituting this into equation (37) we obtain a polynomial whose solutions describe :


Because it is a sixth order polynomial, there is no general formulae for its roots. However, the presence of the asymptotic parameter recommends a series solution. The power of the expansion parameter depends on whether it is large or small.

Taking first, we write


Substituting equation (40) into equation (39), and solving order by order, we find, to lowest order . To next order and so


Substituting this into equations (38) and (41), we find , so the most unstable mode is


Now taking , we write


Following the same procedure, to lowest order and we take the positive root corresponding to the instability. At first order we obtain no information, but at second order we find that . Taking the negative root (see Fig. 5) we have


Substituting this and equation (45) into equation (38), to lowest order , and at next order so , and


Now considering the effect of finite magnetic fields, (and which is now present) increase with . Guided by the numerical results, we assume and take the leading order balance of terms in equation (34). We find

and the resultant growth rate is


which agrees well with the numerics, Fig. 4.

To see the effects of on we adopt the same WKB approach and orderings as Section 4.2 where is now given by the appropriate limit of equation (34). We assume and expand in so . To lowest order is given by equation (20), and to


so the Braginskii viscosity increases the growth rate beyond the inviscid limit.

Figure 5: Direct comparison of the (main panel) and (insets) behaviour of the unstable mode for (top solid line and inset) and (bottom solid line and inset) with for a weak magnetic field, . The dotted lines are the local maxima, and the dashed lines are the asymptotic scalings given by equations (28) and (46) for , and equations (29) and (47) for .

Contrasting the results of this section with the inviscid case, for we have