# Global MRI with Braginskii viscosity in a galactic profile

## Abstract

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.

###### keywords:

instabilities – accretion, accretion discs – Galaxy: disc – MHD – magnetic fields – plasmas.^{1}

^{2}

## 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
conductivities^{3}

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 Reynolds^{4}

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:

(1) | ||||

(2) | ||||

(3) | ||||

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.

### 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 fields^{5}

(5) |

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

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

(6) |

where

(7) | ||||

(8) |

and

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 force^{6}

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

(9) |

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

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

(10) |

where

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

(11) |

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) simplify^{7}

(12) |

with

(13) |

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

(14) |

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

(15) |

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

(16) |

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

(17) |

Demanding that is real, this has solutions

(18) |

and to ensure at , it must satisfy the boundary condition

(19) |

which determines the spectrum of solutions.

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

(21) |

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 .

### 4.2 Viscous MRI,

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

(22) | ||||

(23) |

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

(24) | |||||

(25) |

where is the Euler constant.

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

(26) |

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

(27) |

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

(28) |

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

(29) |

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

(30) |

and at

(31) |

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

(32) |

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 .

In summary, our orderings are

(33) |

### 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 quartic^{8}

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

(37) |

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

(38) |

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

(39) |

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

(40) | ||||

(41) |

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

(42) |

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

(43) |

Now taking , we write

(44) | ||||

(45) |

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

(46) |

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

(47) |

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

(48) |

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

(49) |

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

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