Linear Vlasov Theory in the Shearing Sheet Approximation with Application to the Magneto-Rotational Instability
LLinear Vlasov Theory in the Shearing Sheet Approximationwith Application to the Magneto-Rotational Instability
Tobias Heinemann
1, 2, ∗ and Eliot Quataert † Department of Astronomy, University of California, Berkeley, CA 94720, USA Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA (Dated: January 6, 2015)We derive the conductivity tensor for axisymmetric perturbations of a hot, collisionless, andcharge-neutral plasma in the shearing sheet approximation. Our results generalize the well-knownlinear Vlasov theory for uniform plasmas to differentially rotating plasmas and can be used forwide range of kinetic stability calculations. We apply these results to the linear theory of themagneto-rotational instability (MRI) in collisionless plasmas. We show analytically and numericallyhow the general kinetic theory results derived here reduce in appropriate limits to previous results inthe literature, including the low frequency guiding center (or “kinetic MHD”) approximation, HallMHD, and the gyro-viscous approximation. We revisit the cold plasma model of the MRI and showthat, contrary to previous results, an initially unmagnetized collisionless plasma is linearly stable toaxisymmetric perturbations in the cold plasma approximation. In addition to their application toastrophysical plasmas, our results provide a useful framework for assessing the linear stability ofdifferentially rotating plasmas in laboratory experiments.
1. INTRODUCTION
The shearing sheet is a widely used model for studying the local dynamics of a differentially rotating disk (e.g. Hill1878, Goldreich and Lynden-Bell 1965, Hawley et al. local physics of rotating disks, independent of global boundary conditions. Studies of the shearing sheet have providedvaluable insight into gravitational instability in galactic and protostellar disks, excitation of shearing waves and spiralstructure, and angular momentum transport by turbulence driven by the magneto-rotational instability (MRI).In studies of self-gravitating stellar systems (in particular, galactic disks) the linear theory of the shearing sheethas been widely developed using both kinetic approaches as well as fluid approximations. However, despite theastrophysical importance of low collisionality ionized accretion flows onto compact objects (Rees et al. plasma inthe shearing sheet approximation. That is the aim of this paper. In particular, we derive the general conductivitytensor for such a plasma. We then apply this result to study the MRI in kinetic theory.The analysis in this paper can be viewed as the generalization to a plasma of the previous linear kinetic shearingsheet calculations for self-gravitating stellar systems (e.g. Julian and Toomre 1966). It is likewise a generalization to adifferentially rotating plasma of the standard linear Vlasov theory of uniform plasmas (e.g. Ichimaru 1973, Krall andTrivelpiece 1973, Stix 1992). Moreover, although we specialize to the MRI in the latter part of the paper, our resultsare likely to be useful for understanding other aspects of differentially rotating plasmas in kinetic theory.Our study of the MRI in kinetic theory generalizes and extends previous work. In particular, Quataert et al. (2002)studied the MRI in the “kinetic MHD” (or guiding center) approximation, which averages over the cyclotron motion ofparticles. They also focused solely on kinetic ions, neglecting the electron dynamics under the assumption that theelectrons are much colder than the protons. The formalism derived in this paper relaxes both of these restrictions.Our analysis contains as various limits many non-ideal MHD studies of the MRI, in particular those assessing theimportance of the Hall effect (Section 5.3.1, Wardle 1999) and gyro-viscous forces (Section 5.3.2, Ferraro 2007). Wedisagree, however, with the cold plasma analysis of Krolik and Zweibel (2006), see Section 4.4.This paper is organized as follows. We begin by providing a brief overview of kinetic theory in the shearing sheetapproximation (Section 2). We then describe in some detail the equilibrium solution of the shearing sheet equations,emphasizing in particular the properties of the equilibrium distribution function (Section 3). Section 4.1 containsthe key derivation of the conductivity tensor and the linear dispersion relation in the shearing sheet and shows howthose can be derived from standard uniform plasma results in the literature by an appropriate set of coordinatetransformations. Sections 4.4 and 5 derive analytically and determine numerically the MRI dispersion relation in various ∗ [email protected] † [email protected] a r X i v : . [ a s t r o - ph . H E ] J a n limits and discuss the connection between our results and previous derivations in the literature. Finally, Section 6summarizes our key results. The Appendices A through D contain various technical aspects of the calculations.
2. VLASOV DESCRIPTION OF CHARGED PARTICLE DYNAMICS IN THE SHEARING SHEET
The collisionless evolution of charged particles in a rotating frame of reference is described by the Vlasov equation ∂f s ∂t + v · ∇ f s + e s m s ( E + v × B ) · ∂f s ∂ v − (2 Ω × v + ∇ ψ ) · ∂f s ∂ v = 0 . (1)Here, f s ( r , v , t ) is the one-particle distribution function. This is normalized such that the number density n s and bulkvelocity u s are given by n s = Z d v f s and u s = Z d v v ˆ f s , (2)where ˆ f s = f s /n s is the reduced distribution function.The tidal potential ψ is the sum of the centrifugal potential and the gravitational potential. In the shearing sheetapproximation (Goldreich and Lynden-Bell 1965), this is given by ψ = − q Ω x + 12 ν z . (3)The angular velocity Ω is determined from radial force balance. We assume that in equilibrium, the magnetic field isfrozen into a bulk flow u common to all species ( E + u × B = 0). If the gravitational potential is due to a point massand radial pressure support is negligible, then q = 3 / ν = Ω.With the tidal potential as given in eq. (3), the dynamics are invariant with respect to translations along y . As wewill verify below, this implies conservation of canonical angular momentum. On the contrary, there is no invariancewith respect to mere translations along x . Instead, such translations need to be followed by a Galilean boost along y ,to wit x x + x , y y − q Ω tx , (4)cf. Wisdom and Tremaine (1988).
3. THE STEADY STATE
In light of the symmetries discussed in the previous section, we may seek an equilibrium in which the numberdensities are uniform in x and y and in which the bulk motion of each plasma species is given by u = − q Ω x e y . (5)This shear flow is the local representation of the global differential rotation. We take the magnetic field to be uniform.Ampère’s law µ J = ∇ × B then implies that the electric current density vanishes, which is consistent with thebulk flow (5) being common to all species. The frozen-in condition gives rise to an electric field E = q Ω x e y × B .Given this, the radial magnetic field B x must vanish according to Ferraro’s theorem (Ferraro 1937). The equilibriumelectromagnetic fields are thus given by E = q Ω xB z e x and B = B y e y + B z e z . (6)Let us now consider the motion of an individual particle with charge-to-mass ratio e/m . Given eq. (6), the equationof motion may be written as d v dt = em ( E + v × B ) − (2 Ω × v + ∇ ψ ) = ( v + q Ω x e y ) × S − ν z e z , (7)where we have introduced the vector S = em B + 2 Ω . (8)Note that for neutral particles we recover the equation of motion for test particles in Hill’s approximation (Hill 1878).The equation of motion (7) is generated by the Hamiltonian H = p x + ( p y − S z x ) + ( p z + S y x ) − q Ω S z x + ν z , (9)see Appendix A for a derivation. Expressing this in terms of the particle velocity v = ∂ H /∂ p yields2 E = v − q Ω S z x + ν z , (10)which we will refer to as the energy integral or simply the energy. According to Jeans’ theorem, any distribution functionthat depends only on the integrals of motion is in equilibrium. Thus, any f ( E ) is in equilibrium. However, any suchdistribution function describes an equilibrium in which the plasma is not rotationally supported but pressure-supportedagainst the tidal force: there is no mean flow and the density is non-uniform in x . An equilibrium distribution functionthat is compatible with the underlying assumptions of the shearing sheet is obtained as follows.Inspection of eq. (9) immediately reveals that y is a cyclic coordinate. A second integral of motion (in addition tothe energy integral) is thus given by the conjugate momentum p y . We will henceforth refer to p y = v y + S z x as the(canonical) angular momentum (cf. Wisdom and Tremaine 1988). For a given p y , the dynamics may be described inthe reduced phase space ( x, z, p x , p z ). A particle traveling in the mid-plane at the local shear velocity (5) correspondsto a stationary point in this reduced phase space. To see this, we note that the first variation of H is given by δ H = v x δp x + v z δp z + [ S y v z − S z ( v y + q Ω x )] δx + ν zδz. (11)A stationary point is reached if all coefficients in this expression vanish identically. This is the case if and only if z = 0 and v x = v y + q Ω x = 0 . (12)These orbits are the local representation of global circular orbits (see Appendix B). In the following this will be impliedand we will refer to them simply as circular orbits.For a given angular momentum p y , the energy of a particle on a circular orbit is2 E c = p y − S z /q Ω . (13)Note that E c depends only on the angular momentum and is thus an integral of motion. We now define the gyrationenergy K = E − E c as the difference between a particle’s energy and the energy of a hypothetical particle with thesame angular momentum but on a circular orbit. The gyration energy is thus2 K = v x + ( v y + q Ω x ) − ∆ + v z + ν z , (14)where we have introduced the tidal anisotropy ∆ = q Ω /S z . (15)The gyration energy K is the natural generalization of the epicyclic energy (Shu 1969) to charged particles. Inspectionof eq. (14) reveals that any distribution function of the form f ( K ) describes a rotationally supported plasma with anumber density uniform in x and y and a bulk flow equal to eq. (5). Note that K is manifestly positive definite as longas ∆ <
1. In this case the circular orbit thus minimizes the energy for a given angular momentum. We discuss thecase ∆ ≥ < B = B z e z ), the horizontal degrees of freedom are decoupled from the vertical ones. The corresponding integrals ofmotion are 2 K ⊥ = v x + ( v y + q Ω x ) − ∆ and 2 E z = v z + ν z . (16)The inhomogeneity in z can be dealt with using an action-angle formalism (e.g. Kaufman 1971, 1972). Here and in the following we will refer to both H and 2 H as the Hamiltonian. The same goes for any energy-like quantity. In the remainder of this paper we will, however, make no assumption about the inclination of B in the yz -planebut instead assume that ν = 0 so that the plasma is homogeneous in z . Neglecting stratification is a commonlymade simplification that greatly eases the analysis because we may assume that any disturbance of the equilibrium isperiodic in z . In order to integrate the particle equation of motion in this case, we introduce a new velocity ˚ v whosecomponents are given by ˚ v x = v x , ˚ v y = v y + q Ω x √ − ∆ , ˚ v z = v z . (17)Expressed in terms of ˚ v , the equation of motion (7) is d ˚ v /dt = ˚ v × ω g , (18)where the gyro-frequency vector ω g = S y e y + S z √ − ∆ e z . (19)The gyration energy defined in eq. (14) is simply 2 K = ˚ v . (20)The salient feature of having introduced ˚ v is that eq. (18) takes the same form as in a uniform plasma, with thegyro-frequency vector ω g defined in eq. (19) playing the role of the cyclotron frequency vector ω c = e B /m . Letus therefore adopt a cylindrical coordinate system (˚ v k , ˚ v ⊥ , ϑ ) in ˚ v -space that is aligned with ω g . In this and in thefollowing section, the labels ⊥ and k will always refer to the direction of ω g . The steady state Vlasov equation is[( v + q Ω x e y ) × S ] · ∂f∂ v = ( ˚ v × ω g ) · ∂f∂ ˚ v = − ω g ∂f∂ϑ = 0 , (21)where the square of the gyration frequency ω g is given by ω g = S y + S z ( S z − q Ω) . (22)The equilibrium distribution function is thus gyrotropic in ˚ v -space, i.e. of the form f (˚ v ⊥ , ˚ v k ). For this reason we willrefer to ˚ v as the gyrotropic velocity .The general form f (˚ v ⊥ , ˚ v k ) includes distribution functions of the form f ( K ). Such distribution functions are isotropicin ˚ v -space but anisotropic in v -space (unless q = 0 or e/m → ∞ ). This anisotropy is measured by ∆ as defined ineq. (15). It always lies within the orbital plane and is inevitable in the presence of the tidal force, hence the name tidal anisotropy . By contrast, the anisotropy due to different temperatures in the direction parallel and perpendicularto ω g is a free parameter and is distinct from the tidal anisotropy. Throughout much of the subsequent analysis we will make no assumptions about the specific form of f (˚ v ⊥ , ˚ v k ).Starting with section Section 4.3 we will, however, assume that the distribution function is given by a Maxwellian in ˚ v -space, to wit f ( K ) = n exp( −K /v t )(2 π ) / v t √ − ∆ , (23)where v t = const is the thermal velocity. In the guiding center limit ( e/m → ∞ ), this reduces to a drifting Maxwellianin v -space (as used in e.g. Quataert et al. e/m → ∞ , the tidally induced anisotropy ∆ ∼ m/e disappears. From the definitions given in eqs. (17) and (19) it then follows that ˚ v ≈ v + q Ω x e y and ω g ≈ ω c . In thelimit, any distribution function of the form f (˚ v ⊥ , ˚ v k ) thus becomes a drifting gyrotropic distribution in v -space, withthe anisotropy being due to the magnetic field only.It is also instructive to consider neutral particles ( e = 0). In this case the tidal anisotropy ∆ = q/ ω g = κ , where κ = p − q ) Ω is the epicyclic frequency. The distribution function (23) thus reduces tothe Schwarzschild distribution (see e.g. Julian and Toomre 1966). For this reason and for lack of a better name, wewill refer to f ( K ) as given in eq. (23) as the electromagnetic Schwarzschild distribution . Note that there is of course overlap between the two cases. If P = m R d v ( v − u )( v − u ) f ( K ) denotes the pressure tensor, then ( P xx − P yy ) /P xx = ∆ and P zz = P xx . One might expect that the tidal anisotropy alone can give rise to well known plasma instabilities such as the fire hose instability or themirror instability. We note, however, that the tidal anisotropy is unusual in the sense that it is, in general, not aligned with the magneticfield. Indeed, for a vertical field, the tidal anisotropy lies in the plane perpendicular to the magnetic field. ∆ > In the above derivation of the equilibrium distribution function we have assumed that the tidal anisotropy∆ = q Ω /S z <
1. We have done so because several quantities are singular if ∆ = 1. The validity of our analy-sis is also called into question when ∆ >
1. In this case the gyration energy is not sign definite, in which case a gyratingparticle may reach arbitrarily high velocities for any finite K . One consequence of this is that the electromagneticSchwarzschild distribution given in eq. (23) is not normalizable unless a seemingly arbitrary cutoff in velocity space isintroduced.What is the significance of ∆ >
1? To address this question, let us first consider neutral particles. In this case, thetidal anisotropy ∆ = q/ > κ = 2(2 − q )Ω < . (24)We thus recover the familiar Rayleigh criterion. As is well known, circular orbits are unstable if eq. (24) is satisfied.The local manifestation of this is that unperturbed orbits as described by eq. (7) are hyperbolic if κ <
0. It is doubtfulthat the shearing sheet approximation can be made sense of in this case. It certainly would not describe the localdynamics of a disk.The above discussion of neutral particle generalizes straightforwardly to charged particles orbiting in a vertical magnetic field, in which case the analogue of eq. (24) is ω g < , (25)where the gyration frequency is defined in eq. (22). By analogy with the neutral case, we are led to conclude thateq. (25) is the criterion for the motion of a charged particle on a circular orbit to be unstable. In Appendix B we shedlight on the stability of circular orbits from a global perspective and confirm that this conclusion is indeed correct forparticles orbiting in the mid-plane of the disk.The question of orbital stability becomes more involved if the magnetic field has a toroidal component, in whichcase the magnetic field is misaligned with the rotation axis. In this case, ∆ > ω g < S y (inclined magnetic field) , (26)where we remind the reader that S y = eB y /m . Equation (26) says that if the magnetic field has a toroidal component,then there is a region of parameter space, 1 < ∆ < S y /S z ) , where circular orbits are linearly stable eventhough they do not minimize the energy. In Appendix B we speculate that in spite of eq. (26), circular orbits aredestabilized by dissipation if the tidal anisotropy is in the range 1 < ∆ < S y /S z ) , and that ∆ < <
4. LINEAR THEORY
In this section we calculate the conductivity tensor of a collisionless plasma in the shearing sheet approximation. Theconductivity tensor is obtained from the solution of the linearized Vlasov equation and establishes a linear relationshipbetween the perturbed current density and the perturbed electric field. This together with an independent relationshipbetween current and electric field obtained from Maxwell’s equations leads immediately to the dispersion relation of acollisionless plasma.The calculation of the conductivity tensor of a uniform collisionless plasma is the subject of many textbooks (e.g.Ichimaru 1973, Krall and Trivelpiece 1973, Stix 1992). In generalizing this calculation to the shearing sheet, ourstrategy will be to manipulate the governing equations in a way that allows us to leverage as many results as possiblefrom the (already tedious) uniform plasma calculation.Let us first write the dynamical equations, i.e. the Vlasov equation (1) and Faraday’s law ∂ B /∂t + ∇ × E = 0, in aform that is manifestly invariant with respect to the Galilean transformation (4). For this we introduce the relative velocity and electric field ˜ v = v + q Ω x e y and ˜ E = E − q Ω x e y × B . (27)Note that ˜ v and ˜ E are the velocity and electric field as seen by an observer that is locally at rest with respect to thebackground shear flow. Expressed in terms of the relative velocity and electric field, the Vlasov equation is given by (cid:18) ∂∂t − q Ω x ∂∂y (cid:19) f s + ˜ v · ∇ f s + e s m s ( ˜ E + ˜ v × B ) · ∂f s ∂ ˜ v − (2 Ω × ˜ v − q Ω˜ v x e y + ν z e z ) · ∂f s ∂ ˜ v = 0 , (28)where it is understood that ∂/∂x is now to be evaluated at fixed ˜ v y . Faraday’s law expressed in terms of ˜ E may bewritten as (cid:18) ∂∂t − q Ω x ∂∂y (cid:19) B + q Ω B x e y + ∇ × ˜ E = 0 . (29)This form of Faraday’s law makes explicit that the magnetic field is stretched and advected by the background shear.We note that eqs. (28) and (29) are exact. No approximations have been made so far. We also note that the onlyexplicit x -dependence in eqs. (28) and (29) is multiplying ∂/∂y . Provided the dynamics are axisymmetric, we maythus take the distribution function f s , the magnetic field B , and the relative electric field ˜ E to be periodic in x . The conductivity tensor is obtained from the solution of the linearized Vlasov equation
DDt δf s + e s m s ( δ ˜ E + ˜ v × δ B ) · ∂f s ∂ ˜ v = 0 , (30)where D/Dt denotes the time derivative along unperturbed orbits discussed in Section 3. The formal solution is givenby δf s ( r , ˜ v , t ) = − e s m s Z t −∞ dt [ δ ˜ E ( r , t ) + ˜ v × δ B ( r , t )] · ∂f s ∂ ˜ v . (31)This is a so-called orbit integral: r ( t ) and ˜ v ( t ) describe the phase space trajectory of particle on an unperturbedorbit passing through ( r , ˜ v ) at time t .From now on we will assume no stratification ( ν = 0) and axisymmetric disturbances ( ∂/∂y = 0). The dynamicalequations then have constant coefficients and we may assume that all disturbances depend on position r and time t asexp( ik x x + ik z z − iωt ). For such disturbances, Faraday’s law (29) is easily solved for δ B , to wit δ B = k ω × (cid:20)(cid:18) − q Ω iω e x e y (cid:19) · δ ˜ E (cid:21) . (32)Here, the second term in parentheses arises from the stretching term in eq. (29). Substituting eq. (32) into eq. (31),taking the first moment, multiplying by e s n s , and summing over species yields the perturbed current density δ J ( r , t ) = − iω X s e s n s m s Z d v Z t −∞ dt ˜ v ∂ ˆ f s ∂ ˜ v · (cid:2) i ( ω − k · ˜ v ) + i k ˜ v + q Ω e x e y (cid:3) · (cid:18) − q Ω iω e x e y (cid:19) · δ ˜ E ( r , t ) , (33)where ˆ f s = f s /n s is the reduced distribution function.Equation (33) establishes a linear relationship between δ J and δ ˜ E . For this relation to be useful we need to carryout the orbit integral. Inspection of the unperturbed equations of motion might give the impression that in theshearing sheet approximation, this is a more complicated task than in a uniform plasma. However, we have seen in theprevious section that the equilibrium orbits closely resemble those in a uniform plasma when expressed in terms ofthe gyrotropic velocity ˚ v defined in eq. (17). This suggests that we try to express eq. (33) solely in terms of ˚ v . InAppendix C we show that this yields δ J = − iω X s e s n s m s ( Q s · Λ s · Q s + ∆ s e y e y ) · (cid:18) − q Ω iω e x e y (cid:19) · δ ˜ E , (34) More generally, if ∂/∂y = 0, then each field appearing in eqs. (28) and (29) may be written as a sum of so-called shearing waves (Goldreichand Lynden-Bell 1965, Thomson 1887). In the context of computer simulations, the relative electric field ˜ E as well as B and f s may betaken to be shearing periodic (Lees and Edwards 1972, Hawley et al. where the response tensor Λ s = Z d v Z ∞ dτ ˚ v ∂ ˆ f s ∂ ˚ v · [ i ( ω − k · ˚ v ) + i k ˚ v ] exp { i k · ( r − r ) + iωτ } (35)with τ = t − t . Here, ∆ s = q Ω / ( ω cs b z + 2Ω) is the tidal anisotropy, with both ω cs = e s B/m s and b z being signedquantities. The anisotropy tensor Q s in eq. (34) is given by Q s = e x e x + e y e y p − ∆ s + e z e z . (36)While everything to the left of δ ˜ E on the right hand side of eq. (34) could rightfully be referred to as the shearingsheet conductivity tensor, we choose to define this as σ = − iω X s e s n s m s Q s · Λ s · Q s . (37)Defined in this way, it most closely resembles the uniform plasma conductivity tensor (see e.g. Ichimaru 1973).We are now in a position to write down the general dispersion relation of a charge-neutral collisionless plasmain the shearing sheet approximation. Inserting the solution of Faraday’s law given in eq. (32) into Ampère’s law µ δ J = i k × δ B (neglecting the displacement current), and eliminating the perturbed current using eq. (34) yields D · (cid:18) − q Ω iω e x e y (cid:19) · δ ˜ E = 0 , (38)where the dispersion tensor D = ( k − kk − iωµ σ ) v a − q Ω X s n s m s ρ ω cs /b z ω cs b z + 2Ω e y e y . (39)Here, ρ = P s n s m s is the mass density, v a = B /µ ρ is the square of the Alfvén speed, and we have used thecharge-neutrality condition to rewrite the last term on the right hand side. From now on we will suppress the species index when there is no risk of confusion. As promised earlier, the orbitintegral in eq. (35) only involves the gyrotropic velocity. Let us adopt a right-handed Cartesian coordinate system( e , e , e ) with e = ω g /ω g . The unperturbed orbit in ˚ v -space is then simply ˚ v = ˚ v ⊥ [ e cos( ϑ + ω g τ ) + e sin( ϑ + ω g τ )] + ˚ v k e . (40)The gyro-frequency has no component along the x -direction. Thus we may write e = e z cos ϕ − e y sin ϕ. (41)The orientation of the coordinate axes in the plane perpendicular to ω g is arbitrary at this point. We can exploit thisfreedom by choosing a frame in which the wave vector k = k x e x + k z e z lies in the plane spanned by e and e . Thedesired orientation of the basis vectors is achieved if e and e are defined through k ⊥ e = k x e x + k z sin ϕ e × e x (42)and e = e × e . Indeed, projecting the wave vector onto the new basis shows that k = k ⊥ e + k k e , where theperpendicular and parallel wave number are given by k ⊥ = k x + k z sin ϕ and k k = k z cos ϕ, (43)respectively. In the following we will refer to the coordinate basis defined by eqs. (41) and (42) and illustrated in fig. 1as the Stix basis . x yz ϕχ ω g k Figure 1. Illustration of the coordinate transformations (41) and (42) used in calculating the plasma response tensor (45). The( x, y, z )-coordinate frame is the standard shearing sheet coordinate frame. The gyro-frequency vector ω g defines the 3-axis ofthe (1 , , x -axis by an angle χ about the 3-axis, which in turn is inclined with respect to the z -axis by an angle ϕ about the x -axis. The two Euler angles ϕ and χ are related through k x tan χ = k z sin ϕ . In order to determine the particle orbit in position space, we note that for axisymmetric disturbances it follows fromeq. (17) that k · d r /dt = k · v = k · ˚ v . Working in the Stix basis, it is easy to see that integrating this equationsubject to the condition r ( t ) = r yields k · ( r − r ) = k ⊥ ˚ v ⊥ ω g [sin( ϑ + ω g τ ) − sin ϑ ] + k k ˚ v k τ. (44)With this we can now carry out the integrals over τ and the gyro-phase ϑ in eq. (35). This calculation is entirelyanalogous to the corresponding calculation in a uniform plasma (see e.g. Ichimaru 1973) modulo the substitutions( v ⊥ , v k , ω cs ) → (˚ v ⊥ , ˚ v k , ω gs ). The result may be written as Λ = + ∞ X n = −∞ Z d v nω g ˚ v ⊥ ∂ ˆ f∂ ˚ v ⊥ + k k ∂ ˆ f∂ ˚ v k ! Π n (˚ v ⊥ , ˚ v k ) nω g + k k ˚ v k − ω . (45)In the Stix basis ( e , e , e ), the tensors Π n have the matrix representation Π n (˚ v ⊥ , ˚ v k ) . = n ω g k ⊥ J n i ˚ v ⊥ nω g k ⊥ J n J n ˚ v k nω g k ⊥ J n − i ˚ v ⊥ nω g k ⊥ J n J n ˚ v ⊥ J n − i ˚ v k ˚ v ⊥ J n J n ˚ v k nω g k ⊥ J n i ˚ v k ˚ v ⊥ J n J n ˚ v k J n , (46)where J n is the Bessel function of order n . Its argument is k ⊥ ˚ v ⊥ /ω g and J n denotes the derivative with respect tosaid argument.We stress that all we have assumed in deriving eqs. (45) and (46) is k y = 0. In order to carry out the remainingintegrations over ˚ v ⊥ and ˚ v k in eq. (45), we need to specify the particular form of the distribution function f (˚ v ⊥ , ˚ v k ). Inthe next section we will carry out this calculation for the electromagnetic Schwarzschild distribution defined in eq. (23). The integral over velocity space in eq. (45) is understood to be R d v = 2 π √ − ∆ R ∞ d ˚ v ⊥ ˚ v ⊥ R ∞−∞ d ˚ v k . The result given in eq. (45) holds for general equilibrium distribution functions of the form f (˚ v ⊥ , ˚ v k ). We nowspecialize to the electromagnetic Schwarzschild distribution defined in eq. (23) and given by f ( K ) = n exp( −K /v t )(2 π ) / v t √ − ∆ , (47)where v t = const is the thermal velocity. We remind the reader that the tidal anisotropy ∆ = q Ω / ( ω c b z + 2Ω) and thegyration energy K = 12 " v x + ( v y + q Ω x ) − ∆ + v z = ˚ v . (48)The second equality shows that we are simply considering a plasma that is Maxwellian in ˚ v -space.Given the equilibrium distribution function in eq. (47), the integrals over ˚ v ⊥ and ˚ v k that appear in the responsetensor (45) can now be performed. This calculation is again entirely analogous to the corresponding calculation in auniform plasma (see e.g. Ichimaru 1973). We may thus spare the reader the details and simply state the result Λ = ∞ X n = −∞ ζ ζ n { − W ( ζ n ) } T n − ζ e e , (49)where ζ n = ω − nω g | k k | v t . (50)The W -function used here and in Ichimaru (1973) is related to the standard plasma dispersion function Z defined byFried and Conte (1961) through W ( ζ ) = 1 + ξZ ( ξ ) with ζ = √ ξ . In the Stix basis, the tensors T n in eq. (49) havethe matrix representation T n . = n λ Γ n in Γ n k k | k k | n √ λ ζ n Γ n − in Γ n n λ Γ n − λ Γ n − i k k | k k | √ λ ζ n Γ n k k | k k | n √ λ ζ n Γ n i k k | k k | √ λ ζ n Γ n ζ n Γ n . (51)Here, the Γ n are related to the modified Bessel functions of the first kind I n throughΓ n ( λ ) = I n ( λ ) exp( − λ ) . (52)Their argument in eq. (51) is λ = k ⊥ v t ω g . (53) The cold plasma conductivity tensor may be obtained from the results of the previous section by letting v t →
0. Inthis limit we may drop all W -functions in eq. (49), since ζ n → ∞ , and most of the Γ n in eq. (51), since λ →
0. We arethen left with the response tensor Λ s = e e + 12 X ± ωω ± ω gs ( e ± i e )( e ∓ i e ) . (54)0The response tensor is singular for ω = ω gs . This can be understood as follows. The linearized cold plasma momentumequation in the presence of Coriolis and tidal forces is (cid:18) ∂∂t − q Ω x ∂∂y + 2 Ω × − q Ω e y e x (cid:19) · δ ˜ u s = e s m s ( δ ˜ E + δ ˜ u s × B ) . (55)Multiplying this with n s m s , summing over species, assuming axisymmetric perturbations, and using charge-neutralityyields X s (cid:20) − iω + (cid:18) Ω + e s m s B (cid:19) × − q Ω e y e x (cid:21) · n s m s δ ˜ u s = 0 . (56)This equation has a non-trivial solution (with at least one δ ˜ u s = 0) if the tensor in square brackets is singular, i.e. if iω ( ω − ω gs ) = 0 , (57)where the square of the gyration frequency ω gs is defined in eq. (22).The response tensor in eq. (54) is derived by expressing the perturbed current in terms of the perturbed electric field,which requires inverting the tensor in square brackets in eq. (56). The apparent singularity in the plasma response at ω = ω gs is a consequence of this inversion being possible only if said tensor is non-singular, i.e. if ω = ω gs . In whatfollows we will thus assume that ω = ω gs .Combining eqs. (37), (39) and (54) leads to the dispersion relation of a cold multicomponent plasma in the shearingsheet approximation. The dynamics of such a plasma were studied previously by Krolik and Zweibel (2006). Wedisagree with their analysis for the following reasons. First, they have an incorrect expression for the gyration frequency,which can be read off from the denominator of their eq. (3), which disagrees with our eq. (22). Second, their dispersionrelation does not reproduce the Hall MRI dispersion relation (Wardle 1999, Balbus and Terquem 2001) in the limit m e →
0, as it should (and as our dispersion relation does, see Section 5.3.1 below). Finally, Krolik and Zweibel (2006)claim that a differentially rotating, cold, and initially unmagnetized plasma is linearly unstable with a growth ratecomparable to the rotation frequency. We will now show that this is in fact not the case.
Here we consider the stability of an initially unmagnetized cold plasma in the shearing sheet. The dispersion relationfor this case is obtained from the general expression eq. (39) together with eq. (37) by substituting eq. (54) and thenletting the magnetic field strength tend to zero. The dispersion tensor in this limit reduces to v − a ‘ D = ( k − kk − iωµ σ ) ‘ + q/ e y e y , (58)where ‘ is the mean inertial length, defined through ‘ − = µ X s e s n s m s . (59)The limiting form of the cold plasma conductivity tensor is − iωµ σ = ‘ − Q · Λ · Q , (60)where the anisotropy tensor is given in eq. (36) with ∆ s = q/ B = 0 is Λ = e z e z + 12 X ± ωω ± κ ( e x ± i e y )( e x ∓ i e y ) . (61)With this, setting the determinant of eq. (58) equal to zero yields the dispersion relation ω (1 + k ‘ ) − κ k z ‘ ( q/ k ‘ ) ω − κ = 0 . (62) Multiplying the right hand side of eq. (55) by n s m s and summing over species yields of course the linearized Lorentz force δ J × B . Thisvanishes for k = 0 by Ampère’s law, in which case the solubility condition reduces to ω = κ , independent of any electromagnetic effects.This also follows in complete generality (assuming only charge-neutrality) from the Vlasov equation (28). ω = κ , the solution is ω = κ k z ‘ ( q/ k ‘ )(1 + k ‘ ) . (63)This is manifestly positive as long as q >
0. Thus there is no instability as long as the orbital frequency decreasesradially outward.
5. THE KINETIC MRI WITH FINITE ION CYCLOTRON FREQUENCY5.1. Cold and massless electrons
In most plasmas, electrons are much lighter than ions, and so the mass ratio m e /m i (cid:28)
1. If we assume that theelectrons are also much colder than the ions ( T e /T i (cid:28) σ e ≈ − iω e n e m e bb + en e B b × , (64)where we have dropped terms of order m e /m i and higher. The first term evidently diverges as the mass ratio goes tozero. Because this term is proportional to bb , the dispersion tensor becomes block-diagonal. In the cold and masslesselectron limit, the dispersion relation thus factorizes into δ ˜ E k = 0 and det D ⊥ = 0, where the perpendicular dispersiontensor D ⊥ = ( k − kk − iωµ σ i ) ⊥ v a − iω ¯ ω c b × − q Ω X s n s m s ρ ω cs b z ω cs b z + 2Ω ( b × e x )( b × e x ) . (65)We note that throughout this section, the labels ⊥ and k refer to the direction of b , not to the gyro-frequency vector ω g as in Sections 3 and 4. In eq. (65), σ i is the (multi-component) ion conductivity tensor and¯ ω c = en e Bρ = X s n s m s ρ ω cs (66)is the (mass density weighted) mean ion cyclotron frequency.In Appendix D we derive the dispersion relation of a so-called Vlasov-fluid in the shearing sheet. Vlasov-fluid theory(Freidberg 1972) assumes that the electrons are massless and isotropic, but not necessarily cold. The dispersion relationderived above agrees with the Vlasov-fluid dispersion relation for T e = 0, see eq. (D17).In the remainder of this section we will assume that there is only one ion species. The perpendicular dispersiontensor in this case is given by D ⊥ = v a ( k − kk ) ⊥ + ω c ( Q · Λ · Q ) ⊥ − iωω c b × − q Ω ω c b z /S z ( b × e x )( b × e x ) . (67)We focus on the Vlasov-fluid limit in this section to make contact with the existing kinetic theory calculations ofthe MRI, all of which have assumed cold electrons. We stress, however, that the results derived in the previoussection are fully general and can be used to quantify the effects of kinetic electrons on the MRI. The cold electronapproximation used here and in previous studies formally requires that ζ n → ∞ for all n (Section 4.4). However, ζ = ω/ | k k | v te ∼ v a /v te for the MRI (taking ω ∼ k k v a ). Thus ζ (cid:29) T e /T i (cid:28) ( m e /m i ) β − i . This is ingeneral not satisfied so that the cold electron approximation is not formally applicable when studying the MRI in hotaccretion flows. In future work, we will study the effect of (hot) kinetic electrons on the MRI in collisionless plasmas. In the limit e/m → ∞ we should be able to recover the results of Quataert et al. (2002) who analyzed the stabilityof a differentially rotating plasma in the shearing sheet approximation using the so-called kinetic MHD equations (see We note that since D = − iωµ σ e v a + . . . , it follows from the dispersion relation given in eq. (38) that the parallel component of the(relative) electric field must be at least first order in m e /m i . The parallel current density would otherwise diverge, a result that isconsistent with guiding center theory (e.g. Grad 1961). . . . . . . . . k x /k z . . . . . . . . γ / Ω − k k v a / Ω = p /
16 0 . . . . . k k v A / Ω0 . . . . . − − k x = 0 Figure 2. Dispersion relation of the MRI in the guiding center limit ( ω c / Ω → ∞ ) for an inclined magnetic field with b z = b y .Labels give the value of β z = β/b z . These results agree well with the analogous plots shown in figs. 3 and 4 of Quataert et al. (2002). e.g. Grad 1961, Kulsrud 1983, Hazeltine and Waelbroeck 2004). We start by expanding the response tensor in powersof m/e . Keeping mostly leading order terms, this yields Λ ≈ − ω ω c ( − e e ) + iωω c (cid:20) − (4 − q )Ω b z ω c (cid:21) e × − W − k ⊥ v t ω c e e + iζ W k ⊥ v t ω c e × − ζ W e e . (68)The argument of the W -function is ζ = ω/ | k k | v t . Expanding the anisotropy tensor Q yields Q ≈ − q Ω2 ω c b z (cid:20) − (8 − q )Ω4 ω c b z (cid:21) e y e y . (69)In order to derive the guiding center limit of the perpendicular dispersion tensor given in eq. (67), we need tocompute the tensor product Q · Λ · Q and project the result onto the plane orthogonal to the magnetic field. Doing sois complicated by the fact that for finite e/m and B y = 0, the magnetic field and the gyro-frequency vector ω g = ω g e are misaligned, see fig. 1. This makes taking the limit e/m → ∞ a tedious (albeit straightforward) task. At the end ofthe day, the matrix representation of the perpendicular dispersion tensor (67) in the basis ( e x , b × e x ) is, however,given by the relatively compact expressionlim e/m →∞ D ⊥ . = " k k v a − ω iω Ω b z (1 + W b y /b z ) − iω Ω b z (1 + W b y /b z ) k k v a − ω − q Ω − b y ζ W + v a { β (1 − W ) } " k z b y k x k z b y k x k z b y k x , (70)where β = 2 v t /v a . Figure 2 shows plots of the resulting dispersion relation det D ⊥ = 0 for different values of β z = β/b z .The results are in good agreement with Figs. 3 and 4 of Quataert et al. (2002). The physics of the MRI in the guidingcenter limit are discussed in detail in Quataert et al. (2002) (see also Balbus 2004). k k B k Ω) In this section we study modes whose wave vector is aligned with the magnetic field, which in turn is aligned withthe rotation axis. The ion response tensor in this case is given by Λ = 12 X ± Λ ± ( e ± i e )( e ∓ i e ) − ω k z v t W (cid:18) ω | k z | v t (cid:19) e e , (71)3with perpendicular components Λ ± = ωω ± ω g (cid:20) − W (cid:18) ω ± ω g | k z | v t (cid:19)(cid:21) . (72)At this point we note that even in the case of a purely vertical magnetic field ( B = B z e z ), the coordinate basis( e , e , e ) as defined in Section 4.2 does not in general coincide with ( e x , e y , e z ). This is because e = ω g /ω g , andwhether the gyro-frequency vector ω g is aligned or anti-aligned with the angular frequency vector Ω = Ω e z dependson the sign of S z = ω c b z + 2Ω. In this context we also note that the gyration frequency ω g > The cyclotron frequency ω c = eB/m can have either sign, but this depends only on the sign of the charge, not on thedirection of the magnetic field.In the absence of Coriolis and tidal forces, Λ + (Λ − ) would describe a right (left) handed circularly polarized wave.Insertion of eq. (71) with eq. (72) into eq. (67) yields the perpendicular dispersion tensor D ⊥ . = k z v a (cid:18) (cid:19) + ω c (cid:18) ω g /S z (cid:19) · (cid:18) < =−= < (cid:19) · (cid:18) ω g /S z (cid:19) − iωω c b z (cid:18) −
11 0 (cid:19) − q Ω ω c b z S z (cid:18) (cid:19) , (73)where the right hand side is the matrix representation of D ⊥ in the basis ( e x , e y ). We have also introduced the shorthands < = 12 (Λ + + Λ − ) and = = 12 i (Λ + − Λ − ) . (74)We note that Λ − = Λ ∗ + for purely imaginary ω . In this case, < and = are simply the real and imaginary parts of Λ + . The dispersion relation simplifies significantly in the cold ion limit v t →
0. In this case we may drop the W -functionin eq. (72). The dispersion relation det D ⊥ = 0 may then be written as k z v a ( ω − ω g ) − k z v a ω c b z [ q Ω( ω − κ ) − ω + q Ω ) ω c b z ] − ω ( ω − κ ) ω c ω − ω g = 0 . (75)The numerator reproduces the Hall limit of the dispersion relation given by Wardle (1999). We refer the reader to thatpaper as well as to Balbus and Terquem (2001) for discussions of the Hall MRI. Here we only summarize several keyfindings.Equation (75) is quadratic in both ω and k . Its discriminant with respect to ω is positive definite. Thus thereare no over-stable solutions. The maximum growth rate is γ = q Ω / k z = Ω v a p − α (cid:20) α ) 2Ω ω c b z (cid:21) − / with α = κ , (76)cf. eq. (29) of Wardle (1999). Remarkably, the maximum growth rate is independent of ω c . The system is thus unstableeven as ω c →
0. Below we will see that this is also the case if the ions are hot. This behavior might seem paradoxical inlight of the result of Section 4.4.1, where we have shown that an initially unmagnetized charge-neutral cold collisionlessplasma is linearly stable in the shearing sheet. The resolution of this apparent paradox is that here, it is only the ioncyclotron frequency ω c = ω ci that goes to zero whereas the electron cyclotron frequency ω ce is formally infinite sincewe have assumed the electrons to be massless. We note, however, that for realistic plasmas with finite mass ratios, it isat least questionable whether taking the limit ω ce → ∞ before letting ω ci → In Section 5.2 we have shown that we recover the kinetic MHD result obtained by Quataert et al. (2002) in theguiding center limit ω c → ∞ . For parallel modes it is easy to calculate the first order correction for finite cyclotronfrequency. In order to do so, let us assume the ordering ω ∼ k z v a ∼ Ω ∼ (cid:15)β Ω ∼ (cid:15) ω c , (77) Unless of course ω g <
0, in which case much of our analysis is invalidated, see the discussion in Section 3.2. Note that for parallel modes the MHD and kinetic MHD dispersion relations are identical (cid:15) (cid:28)
1. We thus assume that the plasma is hot but not too hot in the sense that1 (cid:28) β (cid:28) ω c / Ω. The perpendicular dispersion tensor to first order in (cid:15) is then given by D ⊥ = ( k z v a − ω )( − e z e z ) − iω (cid:18) − k z v t ω c b z (cid:19) e z × − q Ω e y e y . (78)Setting its determinant equal to zero yields the dispersion relation ω − ω [ κ + 2 k z v a (1 − β Ω /ω c b z )] + k z v a ( k z v a − q Ω ) = 0 . (79)The lowest order correction to the standard MHD dispersion relation is thus of order β Ω /ω c (cid:28)
1. This is the magnitudeof the gyro-viscous force relative to the magnetic tension force in the Braginskii equations, a fluid closure for weaklycollisional plasmas that includes FLR effects. Ferraro (2007) provides a numerical analysis of the MRI with thegyro-viscous force taken into account. Our results confirm that gyro-viscosity is indeed the dominant FLR effect in ahot plasma, not the Hall effect.Like the cold ion dispersion relation, eq. (78) is quadratic in both ω and k and may be analyzed in the same way.The maximum growth rate γ = q Ω2 (cid:20) − α ) β Ω ω c b z (cid:21) / (80)occurs at k z = Ω v a p − α " − α (1 + α ) β Ω ω c b z / , (81)where again α = κ/ Outside of particular limits such as those discussed in Sections 5.3.1 and 5.3.2 the dispersion relation even for thesimplest case of parallel modes with k k B k Ω , given in eq. (73), is transcendental in ω and k z and its roots can ingeneral only be found numerically. This is the goal of this section.Our numerical analysis has the caveat that we only look for purely growing modes. In support of this we note thatthe MRI is known to be purely growing in many fluid limits, including the cold ion/Hall MHD limit discussed inSection 5.3.1. In addition, the guiding center analysis carried out by Quataert et al. (2002) also yielded purely growingmodes only. For purely growing modes, the dispersion relation has a number of properties that make the root finding particularlyeasy. First, because the conductivity tensor vanishes as ω →
0, the dispersion relation at marginal stability is simply k z v a (cid:18) k z v a − q Ω ω c S z b z (cid:19) = 0 . (82)Second, the dispersion tensor is real for purely growing modes, see Section 5.3. Thus, in order to find the roots of thedispersion relation we may use Newton’s method for real valued functions. In addition, it is easy to get an overview ofhow many roots there are simply by plotting the dispersion relation as a function of γ = √− ω . Generally we findonly a single unstable root.Equation (82) determines the range of unstable wave numbers. This range is finite provided that ω c b z > ω c b z < − k z = 0 and k z = 1 v a s q Ω ω c S z b z . (83) In this context we remark that the gyro-viscous force in the Braginskii equations does in fact not depend on the collision frequency. We stress that it would be unreasonable to assume that all unstable modes are purely growing. We make this assumption only for theMRI. Our general dispersion relation described by eq. (39) encompasses and extends the uniform plasma dispersion relation, and auniform plasma is known to be subject to a host of instabilities that are not purely growing, particularly if said plasma is anisotropicwith respect to the magnetic field (see Schekochihin et al. et al. . β γ/ Ω γ/ Ω0 . ω c b z / . β k z v a / Ω − − − − ω c b z / k z v a / Ω0 . . . . . . . . . . . . . . . . . . . . . . . . Figure 3. Growth rate and wavenumber of the fastest growing MRI mode for k k B k Ω based on a numerical solution of thekinetic ion, cold, massless electron shearing sheet dispersion relation (73). The dashed black lines correspond to β Ω = ω c b z , theapproximate ratio of the gyro-viscous force to the tension force. We do not consider the region of parameter space where − ≤ ω c b z <
0. In support of this we note that for − ≤ ω c b z ≤ ( q − ≥
1, in which case the equilibrium orbits are unbound and the shearingsheet approximation rendered suspect, see the discussion in Section 3.2. For ( q − < ω c b z < q = 3 /
2) in our numerical solutions. The free parameters are then ω c b z and β .We sweep through this two-dimensional parameter space as follows. For each ω c b z we start with β (cid:28) β progressively and – for each k z in the unstablerange – use the previously obtained root as initial guess.Figure 3 shows pseudo color plots of the maximum growth rate and associated wave number. Note that the colorscales are not the same in the left and right panels, because the MRI physics depends on the sign of b z . The blackline in fig. 3 shows β Ω = ω c b z , suggesting that gyro-viscosity is the dominant FLR effect throughout a large portionof parameter space. Figure 4 shows cuts through fig. 3 at fixed β . The cusps in the maximum growth rate and theassociated discontinuities in wave number are a consequence of the dispersion relation having more than one localmaximum, as we discuss below.Figure 4 shows that for ω c b z (cid:38) Ω the growth rate of the fastest growing mode is always greater than the guidingcenter result γ = q Ω /
2. The growth rate increases with β and appears to asymptote for a certain ω c b z ∼ β Ω to γ = √ q Ω. This is the same as the maximum growth rate obtained by Quataert et al. (2002) and Balbus (2004)without FLR effects but with a 45 ◦ -inclined field; see also Section 5.2 and fig. 2. For 0 < ω c b z (cid:46) Ω the system is lessunstable. As ω c b z → γ → q Ω / k z → ω c b z (cid:46) −
4Ω the system is always lessunstable than in the guiding center limit. As ω c b z → −
2Ω from below, however, the growth rate of the fastest growingmode as well as its wave number diverge. While this might seem odd, we note that the tidal anisotropy also diverges60 . . . . . γ / Ω − ω c b z / . . . . . . . k z v a / Ω β = 10 β = 10 β = 10 β = 10 β = 10 − − − − − ω c b z / Figure 4. Growth rate and wavenumber of the fastest growing MRI mode as functions of ω c b z at fixed β . The cusps (upperpanels) and discontinuities (lower panels) correspond to one local minimum of γ ( k z ) taking over another, see fig. 5. in this limit. This anisotropy evidently constitutes a free energy reservoir that the system is able to tap into. We alsonote that the diverging growth rate as ω c b z → −
2Ω from below seems at odds with the cold ion result ( γ = q Ω / ω c b z < − β → γ ( k z ) for various values of ω c b z and β . The dispersion relation can have more than onelocal maximum. The discontinuities visible in the maximum growth rate and associated wavenumber in figs. 3 and 4correspond to one local maximum in fig. 5 overtaking another. The lower right panel demonstrates that for a fiducial ω c b z = − . β →
6. CONCLUSIONS
In this paper, we have studied the linear theory of a collisionless plasma in the shearing sheet approximation. Thisapproximation isolates the key local physics of differentially rotating disks and is free from the complications of globalboundary conditions. It is the ideal model for studying the local dynamics of accretion disks. Our study is motivatedby the application to the magneto-rotational instability (MRI) in such disks if the accreting plasma is collisionless.We have for the first time derived the general linear conductivity tensor for a collisionless plasma in the shearingsheet approximation. This accounts for the complete linear dynamics of kinetic ions and kinetic electrons with finiteLarmor radius/cyclotron frequency. Our calculation has leveraged the extensive literature on analogous calculations foruniform plasmas (Ichimaru 1973). In particular, we have shown that by a suitable choice of velocity space coordinates,the linear conductivity tensor for a collisionless plasma in the shearing sheet approximation can be expressed usingstandard results for a uniform plasma (Section 4.1 and Appendix C). In doing so, our primary assumptions arecharge-neutrality, axisymmetric perturbations, and no vertical stratification. The charge-neutral approximation ismade both for its applicability to realistic accretion disk dynamics and because the shearing sheet approximation in itsstandard form relies on Galilean invariance, which in turn necessitates charge-neutrality once the displacement currentis dropped (see e.g. Grad 1966).70 . . . . . . . . . . . . . . . . . . . . γ / Ω β = 1000 ω c b z / − . . . . . . . . . . . . . . . . . . . . ω c b z = 100 Ω β − − − k z v a / Ω0246810121416 γ / Ω ω c b z = − . β − k z v a / Ω0 . . . . . . . ω c b z = − . β − Cold
Figure 5. Plots of growth rate vs. wave number for select values of ω c b z and β . The lower panels illustrate the approach to thecold plasma dispersion relation as β → The resulting dispersion relation for a collisionless plasma in the shearing sheet approximation is given in eq. (39)with the conductivity tensor given in eq. (45). These results apply for any equilibrium distribution function. In general,the tidal force imposes a temperature anisotropy on the equilibrium distribution function, see eq. (15) and Section 3.This anisotropy vanishes only in the guiding center limit in which the ratio of the cyclotron frequency to the angularfrequency ω c / Ω → ∞ . We denote a “quasi-Maxwellian” distribution function with this minimal level of temperatureanisotropy as the “electromagnetic Schwarzschild distribution” because it is a generalization of the Schwarzschilddistribution of stellar dynamics to an ionized plasma. The conductivity tensor for this distribution function is given inequations eqs. (49) to (51). These results for the linear theory of the shearing sheet are likely to be applicable to awide range of linear stability calculations, beyond those we have carried out in this paper.We have focused on the application of our results to the linear theory of the MRI in collisionless plasmas. Inparticular, we have shown that the extensive existing literature on kinetic calculations of the MRI can all be understoodas approximations to the more general theory derived here. Specifically, the guiding center kinetic MRI results ofQuataert et al. (2002) follow from the high cyclotron frequency and cold, massless electron limit of our dispersionrelation (Section 5.2). The numerical gyro-viscous results of Ferraro (2007) – which represent finite cyclotron frequencycorrections to the guiding center results of Quataert et al. (2002) – follow analytically from our dispersion relation bykeeping the leading order cyclotron frequency corrections to the guiding center result (Section 5.3.2). And the HallMHD results of Wardle (1999) – which are particularly important in the context of protostellar accretion disks – followfrom the cold ion and cold, massless electron limit of our dispersion relation (Section 5.3.1).8In addition to providing a more general framework for understanding existing calculations of the MRI in kinetictheory, we have also derived a number of new results. In Section 5.1 and Appendix D we derive the dispersion relationfor the MRI in the limit of kinetic ions, but cold, massless electrons. We solve this dispersion relation for the simplestcase of k k Ω k B in Section 5.3. These numerical solutions of the full hot plasma dispersion relation confirm thatthe primary finite cyclotron frequency correction to the guiding-center MRI dispersion relation for β (cid:38) k k Ω k B . In this case, the kinetic corrections to the MRI (as well any other linear waves/modes thatmay be present) are due to finite cyclotron frequency effects. However, the full dispersion relation derived in Section 4also includes finite Larmor radius effects that become important when k ⊥ v t ∼ ω g . While differential rotation is unlikelyto be important on such scales if ω c (cid:29) Ω, the scale separation between ω c and Ω that can be achieved in numericalsimulations is always limited (unless of course the guiding center motion is averaged over, as in gyrokinetic simulations).Thus, even if the interplay between orbital motion and Larmor motion is insignificant in a given astrophysical system,our dispersion relation can be useful for benchmarking a code used to simulate said system.Possible future applications of our analysis include studies of the stability of linear shear profiles in collisionlessplasmas and the possibility of velocity space instabilities driven by the background temperature anisotropy presentin differentially rotating plasmas. In addition, all existing kinetic calculations of the MRI assume cold electrons,motivated by the fact that hot accretion flows onto compact objects are believed to have T i (cid:38) T e . However, thegeneral conductivity tensor given in Section 4.1 applies equally well to kinetic electrons and can be used to assess theapplicability of the cold electron approximation under realistic accretion disk conditions. Finally, our calculations canbe used to determine the properties of the MRI under the conditions achievable in laboratory experiments aimed atstudying low collisionality high- β plasmas (e.g. Collins et al. ACKNOWLEDGMENTS
We thank Omer Blaes, Martin Pessah, Mario Riquelme, Anatoly Spitkovsky, and Scott Tremaine for usefulconversations. In particular, conversations with Anatoly Spitkovsky inspired us to think in more detail about theunmagnetized limit (Section 4.4.1). We also thank Ellen Zweibel and Julian Krolik for a free exchange of ideasthat prompted us to sharpen our understanding of the cold plasma dispersion relation (Section 4.4). This workwas supported in part by NASA HTP grant NNX11AJ37G, NSF grants PHY11–25915 and AST-1333682, a SimonsInvestigator award from the Simons Foundation, and the Thomas Alison Schneider Chair in Physics at UC Berkeley.
Appendix A: The unperturbed Hamiltonian
In this section we derive the unperturbed Hamiltonian given in eq. (9). Let us start with the Hamiltonian governingthe dynamics of charged particles in a rotating frame. This is given by H = 12 (cid:16) p − Ω × r − em A (cid:17) + ψ + em Φ , (A1)where A is the magnetic vector potential and Φ is the electrostatic potential. The coordinate system is oriented suchthat Ω = Ω e z . In the shearing sheet approximation, the tidal potential ψ is given as in eq. (3) by ψ = − q Ω x + 12 ν z . (A2)It is straightforward to verify that ∂f /∂t + { f, H} = 0 is the same as the Vlasov equation (1) when the former isexpressed in terms of the particle velocity v = p + Ω × r + e A /m .We now consider the equilibrium discussed in Section 3, consisting of a linear shear flow and a uniform magneticfield. Without loss of generality we may (at first) work in the so-called symmetric gauge, in which the vector potential9is given by A = B × r /
2. The magnetic field is assumed frozen into the background shear flow. The correspondingelectric field derives from the electrostatic potential Φ = − q Ω xA y . With this, the Hamiltonian may be written as2 H = (cid:18) p − S × r (cid:19) − q Ω S z x + ν x . (A3)The Hamiltonian (9) is obtained after a canonical transformation to new coordinates ( r , p ) derived from the type-2generating function G ( r , p ) = xp x + yp y − S z xy/ . (A4) Appendix B: Global equilibrium
In this section we investigate the stability of charged particles on circular orbits in a global disk equilibrium. Wealso derive a generalization of the modified Schwarzschild distribution (Shu 1969) and show that for nearly circularorbits it reduces to the “electromagnetic Schwarzschild distribution” used in the main text.
1. Bulk flow and electromagnetic field configuration
We seek an axisymmetric equilibrium with a purely rotational bulk flow common to all plasma species. In acylindrical coordinate system ( r, ϕ, z ), this bulk flow may be written as u = r Ω ∇ ϕ, (B1)where Ω( r, z ) is the orbital frequency. Following Kaufman (1972), Cary and Littlejohn (1983) we write the magneticfield as B = ∇ w ( ψ ) × ∇ ϑ − ∇ ψ × ∇ ϕ, (B2)where ψ ( r, z ) is the poloidal flux function. The poloidal angle ϑ ( r, z ) may be chosen to increase clockwise from − π to π along poloidal field lines, which are isocontours of ψ . Note that the first term on the right hand side of eq. (B2)represents the toroidal field. The second term is the usual parametrization of an axisymmetric poloidal field. Themagnetic field given in eq. (B2) derives from the axisymmetric vector potential A = w ( ψ ) ∇ ϑ − ψ ∇ ϕ. (B3)The electric field is determined from the condition that the magnetic field be frozen into the equilibrium flow, to wit E = − u × B . (B4)This must be irrotational, thus ∇ × E = (cid:18) ∂ϑ∂r ∂ψ∂z − ∂ψ∂r ∂ϑ∂z (cid:19) ∂ Ω ∂ϑ r ∇ ϕ = 0 . (B5)The factor in parentheses is the Jacobian of the transformation from flux coordinates ( ψ, ϕ, ϑ ) to cylindrical coordinates( r, ϕ, z ). Assuming this is non-singular, it follows that the orbital frequency Ω( ψ ) is constant on poloidal field lines(Ferraro 1937), see also Papaloizou and Szuszkiewicz (1992). With this result, it easy to show that the electric fieldderives from the electrostatic potential Φ( ψ ) = − Z Ω( ψ ) dψ. (B6)
2. Stability of circular orbits
The equations of motion for charged particles orbiting in the equilibrium fields are generated by the Hamiltonian H = 12 (cid:20)(cid:16) p r − em A r (cid:17) + 1 r (cid:16) p ϕ + em ψ (cid:17) + (cid:16) p z − em A z (cid:17) (cid:21) + Ψ + em Φ , (B7)0where Ψ( r, z ) is the gravitational potential. The equilibrium Hamiltonian is independent of ϕ due to the assumed axialsymmetry. The canonical angular momentum p ϕ = r ˙ ϕ − em ψ (B8)is thus an integral of motion and for a given p ϕ , the dynamics may be described in the reduced phase space ( r, z, p r , p z ).We assume that the equilibrium also has equatorial symmetry. By this we mean that there is a plane of symmetryat z = 0 (the mid-plane ) and that the equilibrium fields are either even or odd functions of z . Specifically, we take thegravitational and electrostatic potentials as well as the poloidal flux ψ to be even, and the poloidal angle ϑ to be odd.This completely specifies the equatorial symmetry of the equilibrium.Using Hamilton’s equation, the first variation of eq. (B7) is δ H = ˙ rδp r + ˙ zδp z − ˙ p r δr − ˙ p z δz . All equilibrium pointsin phase space (for which δ H = 0) are thus circular orbits ( ˙ r = ˙ z = 0). For such orbits, the rates of change of thecanonical momenta are ˙ p r = em r ( ˙ ϕ − Ω) B z − ∂ Ψ ∂r + r ˙ ϕ (B9)˙ p z = em r (Ω − ˙ ϕ ) B r − ∂ Ψ ∂z , (B10)where we have made use of our assumption that the electric field is frozen into the bulk flow u = r Ω ∇ ϕ . Given theequatorial symmetry of Ψ, particles on circular orbits are (relative) equilibria if they travel in the mid-plane ( z = 0)and their angular velocity matches that of the bulk flow, i.e. if˙ ϕ = Ω( r, , (B11)provided the latter arises due to gravito-centrifugal balance, i.e. provided that r Ω ( r,
0) = ∂ Ψ( r, ∂r . (B12)From now on it will be understood that all quantities are evaluated at z = 0.In order to determine the stability of circular orbits we can invoke the so-called Lagrange-Dirichlet theorem (seee.g. Krechetnikov and Marsden 2007, and references therein). This states that an equilibrium point (in this case thecircular orbit) is stable to finite perturbations (nonlinear stability) if the second variation of H (a quadratic form) or,equivalently, its Hessian is positive definite. In the phase space coordinate system ( r, z, p r , p z ), the Hessian has thematrix representation D . = S z (cid:18) S z + r ∂ Ω ∂r (cid:19) + (cid:18) em ∂A z ∂r (cid:19) − em ∂A z ∂r ∂ Ψ ∂z + (cid:18) em ∂A r ∂z (cid:19) − em ∂A r ∂z − em ∂A r ∂z − em ∂A z ∂r , (B13)where S = em B + 2 Ω (B14)is defined in the same way as in the main text, see eq. (8). The Hessian matrix given in eq. (B13) is positive definite ifand only if S z (cid:18) S z + r ∂ Ω ∂r (cid:19) > ∂ Ψ ∂z > . (B15)Otherwise the Hessian is indefinite. If the above inequalities are satisfied, then circular orbits minimize the energy.Note that if we define ∆ = − rS z ∂ Ω ∂r , (B16)1in complete analogy to eq. (15), then the first inequality in eq. (B15) is equivalent to ∆ <
1. Note also that for circularorbits, the radial derivative of the canonical angular momentum is given by ∂p ϕ ∂r = r (cid:18) S z + r ∂ Ω ∂r (cid:19) = rS z (1 − ∆) . (B17)For neutral particles, the condition for stability ∆ < κ = 2Ω (cid:18)
2Ω + r ∂ Ω ∂r (cid:19) > . (B18)From eq. (B17) we recover the familiar result (e.g. Chandrasekhar 1961) that the angular momentum is a decreasingfunction of radius in this case. For charged particles, the (canonical) angular momentum in a stable configuration doesnot necessarily decrease with radius because S z is not sign definite. a. Linear stability Above we have shown that circular orbits are nonlinearly stable provided that ∆ < ∂ Ψ /∂z >
0. This doesnot necessarily imply that circular orbits are unstable if ∆ > To see this, we consider the linearized Hamiltoniansystem δ ˙ w = J · D · δ w , (B19)where w = ( r, z, p r , p z ), J = − − , (B20)and D is the Hessian given in eq. (B13). Assuming that δ w depends on time as exp( − iωt ), the linearized system hasnon-trivial solutions provided thatdet( iω + J · D ) = ω − ω (cid:20) S z (cid:18) S z + r ∂ Ω ∂r (cid:19) + S ϕ + ∂ Ψ ∂z (cid:21) + S z (cid:18) S z + r ∂ Ω ∂r (cid:19) ∂ Ψ ∂z = 0 , (B21)where denotes the identity matrix.Let us first consider the case of a vanishing toroidal magnetic field. For S ϕ = eB ϕ /m = 0, the solubility condition(B21) dramatically simplifies because the radial and vertical degrees of freedom decouple. The stability of the respectivemotions is determined by ω = S z (cid:18) S z + r ∂ Ω ∂r (cid:19) and ω = ∂ Ψ ∂z . (B22)Comparison with eq. (B15) shows that circular orbits are linearly stable if they are nonlinearly stable (i.e. if theminimize the energy), and linearly unstable otherwise.The situation is more complicated if there is a toroidal magnetic field. In this case there is a region of parameterspace where circular orbits are linearly stable even though they do not minimize the energy. To see this, let us assume(as in the bulk of the main text) that there is no stratification (Ψ = const). In this case, the solubility conditioneq. (B21) reduces to ω ( ω − ω g ) = 0 , (B23)where the square of the gyration frequency ω g is given by ω g = S z (cid:18) S z + r ∂ Ω ∂r (cid:19) + S ϕ . (B24) We use the terms linear and nonlinear stability rather loosely. For more precise definitions see e.g. Holm et al. (1985). < S ϕ /S z , (B25)which does not coincide with the nonlinear stability boundary ∆ <
1. In this context we remark that according to atheorem due to Bloch et al. (1994), a relative equilibrium that is linearly but not nonlinearly stable becomes linearlyunstable when a small amount of dissipation is added to the system. We are thus led to conclude that circular orbitsin realistic systems are always unstable if ∆ >
1, regardless of whether or not there is a toroidal magnetic field.
3. The modified electromagnetic Schwarzschild distribution
The electromagnetic Schwarzschild distribution (23) may be viewed as a generalization of the Schwarzschilddistribution (e.g. Julian and Toomre 1966). The latter is the local limit (nearly circular orbits) of the modified
Schwarzschild distribution as derived by Shu (1969). Here we generalize the modified Schwarzschild distribution andshow that it reduces to eq. (23) in the local limit.The integrals of motion are the particle energy E = ˙ r + r ˙ ϕ + ˙ z em Φ (B26)and the canonical angular momentum p ϕ defined in eq. (B8). By Jeans’ theorem, any distribution function of the form f ( E , p ϕ ) is a solution of the steady state Vlasov equation. Let us define the guiding center radius r c through p ϕ = r c Ω( r c , − em ψ ( r c , . (B27)Thus, r c is the radius of circular orbit with angular momentum p ϕ . Its energy is E c = r c Ω ( r c , r c ,
0) + em Φ( r c , . (B28)Note that E c depends only on the angular momentum and is thus an integral of motion. We define the modifiedelectromagnetic Schwarzschild distribution as f ∝ exp( −K /T ) , (B29)where K = E − E c (B30)is the gyration energy. The proportionality constant in eq. (B29) and the temperature T may in general depend on thecanonical angular momentum.We now show that for nearly circular orbits, eq. (B29) reduces to the electromagnetic Schwarzschild distributiongiven in eq. (23). Let us first introduce the peculiar velocity˜ v = ˙ r − r Ω( r, z ) e ϕ . (B31)Substituting r ˙ ϕ = 1 r h p ϕ + em ψ ( r, z ) i = 1 r h r c Ω( r c ,
0) + em ψ ( r, z ) − em ψ ( r c , i (B32)and expanding for small ˜ v ϕ /r Ω, the azimuthal component of the peculiar velocity is approximately given by˜ v ϕ = (cid:20) S z ( r,
0) + r ∂ Ω( r, ∂r (cid:21) ( r − r c ) , (B33) Circular orbits are relative equilibria. S is defined in eq. (B14). Corrections to eq. (B33) are quadratic in r − r c and z . Expanding the gyration energydefined in eq. (B30) to quadratic order in ˜ v ϕ /r Ω yields2 K = ˜ v r + S z ( r, (cid:20) S z ( r,
0) + r ∂ Ω( r, ∂r (cid:21) ( r − r c ) + ˜ v z + ν z , (B34)Here, the square of the vertical frequency ν is defined as ν = ∂ Ψ ∂z (B35)evaluated at z = 0. Combining eq. (B33) into eq. (B34) finally yields K = 12 " ˜ v r + ˜ v ϕ − ∆ + ˜ v z + ν z , (B36)which is just the electromagnetic Schwarzschild distribution defined in eq. (14). Appendix C: The conductivity tensor
Here we fill in some of the details of going from the formal solution to the linearized Vlasov equation, eq. (31), tothe “conductive” relationship between current and electric field, eq. (34). Note that we are assuming axisymmetricperturbations, so k y = 0.Using Faraday’s law in the form of eq. (32), the formal solution is given by δf s = − iω e s m s Z t −∞ dt ∂f s ∂ ˜ v · [ i ( ω − k · ˜ v ) + i k ˜ v + q Ω e x e y ] · δ ˜ E , (C1)where r ( t ) and ˜ v ( t ) satisfy d r /dt = ˜ v − q Ω x e y and d ˜ v dt = e s m s ˜ v × B − Ω × ˜ v + q Ω˜ v x e y (C2)subject to the “final conditions” r ( t ) = r and ˜ v ( t ) = ˜ v . Taking the first moment of eq. (C1), multiplying by e s n s ,and summing over species yields the perturbed current density δ J = ˜ σ · (cid:18) − q Ω iω e x e y (cid:19) · δ ˜ E . (C3)Here, we have introduced the tensor˜ σ = − iω X s e s n s m s Z d v Z ∞ dτ ˜ v ∂f s ∂ ˜ v · [ i ( ω − k · ˜ v ) + i k ˜ v + q Ω e x e y ] exp {− iφ ( τ ) } , (C4)where the new time variable τ = t − t and the short hand φ ( τ ) = k · ( r − r ) − ωτ. (C5)From now on we will suppress the species index when there is no risk of confusion. In order to express eq. (C4) interms of the “gyrotropic velocity” ˚ v = Q − · v , let us first note that ∂f∂ ˜ v = ∂f∂ ˚ v · Q + ∆ ∂f∂ ˜ v y e y , (C6)where the anisotropy tensor Q is defined in eq. (36). With this, the orbit integral in eq. (C4) may be rewritten as Z ∞ dτ ∂f∂ ˜ v · [ i ( ω − k · ˜ v ) + i k ˜ v + q Ω e x e y ] exp {− iφ ( τ ) } = − ∆ ∂f∂ ˜ v y e y + Z ∞ dτ ∂f∂ ˚ v · [ i ( ω − k · ˚ v ) + i k ˚ v ] · Q exp {− iφ ( τ ) } . (C7)4In obtaining this result we have used that q Ω ∂f∂ ˜ v x − ∆ ∂∂τ ∂f∂ ˜ v y = (cid:18) q Ω˜ v x − ∆1 − ∆ ∂ ˜ v y ∂τ (cid:19) v ⊥ ∂f∂ ˚ v ⊥ = 0 . (C8)The second equality in follows from eq. (C2). To get the first equality, notice that ∂f∂ ˜ v x = ˜ v x ˚ v ⊥ ∂f∂ ˚ v ⊥ and ∂f∂ ˜ v y = 11 − ∆ ˜ v y ˚ v ⊥ ∂f∂ ˚ v ⊥ + ˚ v k sin ϕ √ − ∆ (cid:18) v ⊥ ∂f∂ ˚ v ⊥ − v k ∂f∂ ˚ v k (cid:19) , (C9)where ϕ is defined through ω g = e z cos ϕ − e y sin ϕ . The second term in the expression for ∂f /∂ ˜ v y does not dependon time and therefore makes no contribution in eq. (C8). Using eq. (C7) in eq. (C4) and integrating the contributionfrom the first term on the right hand side of eq. (C7) by parts yields˜ σ = σ − iω X s e s n s m s ∆ s e y e y (C10)with the conductivity tensor σ as defined in eq. (37). Inserting eq. (C10) into eq. (C3) then finally yields eq. (34). Appendix D: Vlasov-fluid dispersion relation
The Vlasov-fluid equations (Freidberg 1972) describe the low-frequency electromagnetic interaction of collisionlessions with a massless and isotropic electron fluid. A detailed description is given in the recent work by Cerfon andFreidberg (2011). Here we show explicitly that the Vlasov fluid dispersion relation in the shearing sheet is equivalentto the cold, massless electron limit of the general hot plasma shearing sheet dispersion relation derived in Section 5.1(up to a term related to electron pressure that is included in the Vlasov fluid model as a crude representation of finiteelectron temperature effects). The Vlasov fluid calculation presented here is useful primarily in that it provides anindependent formalism for calculating the shearing sheet dispersion relation in the limit in which ion dynamics are themost important. In addition, the Vlasov fluid derivation is somewhat more similar to the MHD calculation in that itutilizes the Lagrangian displacement and the fact that field lines are frozen into the electron fluid.
1. The Vlasov-fluid equations
The electrons dynamics are described by the momentum equation en e ( E + u e × B ) + ∇ p e = 0 , (D1)which we will henceforth refer to as the generalized Ohm’s law. Here, we assume that the electron pressure p e is afunction of the electron number density n e alone, i.e. we assume a barotropic equation of state. Equation (D1) arisesout of the electron momentum equation in the limit of vanishing electron mass and assuming an isotropic pressure.The electrons must be collisional for this assumption to be reasonable. Indeed, Rosin et al. (2011) derive eq. (D1)rigorously through an expansion in the electron-to-ion mass ratio under the assumption of weakly collisional ions butstrongly collisional electrons.The crude electron model represents the most severe limitation of Vlasov-fluid theory. At the same time it may,however, also be viewed as the theory’s greatest merit. This is because the fast electron scales are not part of theproblem anymore: The plasma frequency and Debye length are not finite due to charge-neutrality and neither is theelectron Larmor motion due to the vanishing electron mass. The absence of these fast time scales simplifies the analysisand significantly reduces the computational burden in numerical simulations (Byers et al. et al. et al. µ J = ∇ × B thus lacks thedisplacement current. The system of equations is closed with the definition of the current density J = − en e u e + X s e s Z d v v f s (D2)and the charge-neutrality condition en e = X s e s Z d v f s , (D3)5where the summation is over ion species only.The electron velocity u e plays a more fundamental role in Vlasov-fluid theory than may be apparent at first sight.Differentiating eq. (D3) with respect to time yields after a little bit of algebra the continuity equation ∂n e ∂t + ∇ · ( n e u e ) = 0 . (D4)With the barotropic assumption, eliminating the electric field between Ohm’s law and Faraday’s law yields the inductionequation ∂ B ∂t = ∇ × ( u e × B ) . (D5)The magnetic field is thus frozen into the electron fluid. The form of eqs. (D4) and (D5) suggests that in linear theorywe introduce the Lagrangian displacement through δ u e + ξ · ∇ u e = (cid:18) ∂∂t + u e · ∇ (cid:19) ξ (D6)cf. Lynden-Bell and Ostriker (1967), Friedman and Schutz (1978). The linearized eqs. (D4) and (D5) are then easilyintegrated to yield δn e + ∇ · ( n e ξ ) = 0 (D7)and δ B = ∇ × ( ξ × B ) . (D8)We stress that eqs. (D6) to (D8) hold for general equilibria.
2. The dispersion relation in the shearing sheet
In this section we linearize the Vlasov-fluid equations around the equilibrium discussed in Section 3. The equilibriumbulk velocity common to all species is thus given by the linear shear flow (5). As in Section 4.1 we will work in termsof the invariant velocity ˜ v and electric field ˜ E as defined in eq. (27). We thus use the Vlasov equation and Faraday’slaw in the form of eqs. (28) and (29), respectively. Also introducing the invariant electron velocity˜ u e = u e + q Ω x e y , (D9)the generalized Ohm’s law reads en e ( ˜ E + ˜ u e × B ) + ∇ p e = 0 , (D10)and the electric current density is given by J = − en e ˜ u e + X s e s Z d v ˜ v f s . (D11)The dispersion relation is obtained by linearizing eq. (D11) and expressing it entirely in terms of the Lagrangiandisplacement ξ . As in Section 4.1 we consider an axisymmetric plane wave. Equation (D6) then reduces to δ ˜ u e = − iω ξ + q Ω ξ x e y . (D12)The ion current is given by the right hand side of eq. (34), to wit X s e s Z d v ˜ v δf s = σ − iω X s e s n s m s ∆ s e y e y ! · (cid:18) − q Ω iω e x e y (cid:19) · δ ˜ E , (D13)where the ion conductivity tensor σ is defined in eq. (37) with the understanding that the summation in this case isover ion species only. The electric field is obtained from the generalized Ohm’s law in the form of eq. (D10). Usingeqs. (D8) and (D12) to eliminate the electron density and velocity yields δ ˜ E = ( iω + q Ω e x e y ) · ( ξ × B ) + e − dp e /dn e ∇∇ · ξ . (D14)6Finally, we combine Ampère’s law and eq. (D8) to obtain µ δ J = ∇ × [ ∇ × ( ξ × B )] . (D15)Putting everything together yields after some algebra µ − ∇ × ∇ × ( ξ × B ) = en e iω ξ + σ · ( iω ξ × B + e − dp e /dn e ∇∇ · ξ ) − q Ω ξ x e y X s e s n s ω cs b z + 2Ω . (D16)Assuming a plane wave decomposition, this equation may be written as D · ξ = 0 with the dispersion tensor given by D = " ( k − kk − iωµ σ ) v a − iω ¯ ω c b × − q Ω X s n s m s ρ ω cs /b z ω cs b z + 2Ω e y e y × b + ¯ ω c (cid:18) iω bb − σ · kk e n e dp e dn e (cid:19) , (D17)where ¯ ω c = en e Bρ = X s n s m s ρ ω cs (D18)is the mass-weighted mean cyclotron frequency. The perpendicular components of eq. (D17) agree with eq. (65). G. W. Hill, American Journal of Mathematics , 129 (1878).P. Goldreich and D. Lynden-Bell, Monthly Notices of the Royal Astronomical Society , 125 (1965).J. F. Hawley, C. F. Gammie, and S. A. Balbus, The Astrophysical Journal , 690 (1996).M. J. Rees, M. C. Begelman, R. D. Blandford, and E. S. Phinney, Nature , 17 (1982).F. Yuan and R. Narayan, Annual Review of Astronomy and Astrophysics (2012).W. H. Julian and A. Toomre, The Astrophysical Journal , 810 (1966).S. Ichimaru, Basic principles of plasma physics: a statistical approach (Benjamin/Cummings Publishing Company, 1973).N. A. Krall and A. W. Trivelpiece,
Principles of plasma physics (McGraw-Hill, New York, 1973).T. H. Stix,
Waves in Plasmas (American Institute of Physics, New York, 1992).E. Quataert, W. Dorland, and G. W. Hammett, The Astrophysical Journal , 524 (2002).M. Wardle, Monthly Notices of the Royal Astronomical Society , 849 (1999).N. M. Ferraro, The Astrophysical Journal , 512 (2007).J. H. Krolik and E. G. Zweibel, The Astrophysical Journal , 651 (2006).J. Wisdom and S. Tremaine, The Astronomical Journal , 925 (1988).V. C. A. Ferraro, Monthly Notices of the Royal Astronomical Society , 458 (1937).F. H. Shu, The Astrophysical Journal , 505 (1969).A. N. Kaufman, Physics of Fluids , 387 (1971).A. N. Kaufman, Physics of Fluids , 1063 (1972).W. Thomson, Philosophical Magazine Series 5 , 188 (1887).A. W. Lees and S. F. Edwards, Journal of Physics C: Solid State Physics , 1921 (1972).J. F. Hawley, C. F. Gammie, and S. A. Balbus, The Astrophysical Journal , 742 (1995).B. D. Fried and S. D. Conte, The Plasma Dispersion Function (Academic Press, New York, 1961).S. A. Balbus and C. Terquem, The Astrophysical Journal , 235 (2001).H. Grad,
Microscopic and Macroscopic Models in Plasma Physics , Tech. Rep. (1961).J. P. Freidberg, Physics of Fluids , 1102 (1972).R. M. Kulsrud, in Handbook of Plasma Physics , Vol. 1, edited by M. N. Rosenbluth and R. Sagdeev (North Holland, New York,1983) Chap. 1.4, pp. 115–145.R. D. Hazeltine and F. L. Waelbroeck,
The Framework of Plasma Physics (Westview, Boulder, CO, 2004).S. A. Balbus, The Astrophysical Journal , 857 (2004).A. A. Schekochihin, S. C. Cowley, R. M. Kulsrud, G. W. Hammett, and P. Sharma, The Astrophysical Journal , 139 (2005).P. Sharma, G. W. Hammett, E. Quataert, and J. M. Stone, The Astrophysical Journal , 952 (2006).H. Grad,
The Guiding Center Plasma , Tech. Rep. (1966).C. Collins, N. Katz, J. Wallace, J. Jara-Almonte, I. Reese, E. Zweibel, and C. B. Forest, Physical Review Letters , 115001(2012).J. R. Cary and R. G. Littlejohn, Annals of Physics , 1 (1983).J. C. B. Papaloizou and E. Szuszkiewicz, Geophysical & Astrophysical Fluid Dynamics , 223 (1992).R. Krechetnikov and J. E. Marsden, Reviews of Modern Physics , 519 (2007).S. Chandrasekhar, International Series of Monographs on Physics, Oxford: Clarendon, 1961 (Dover Publications, Oxford, 1961).D. D. Holm, J. E. Marsden, T. Ratiu, and A. Weinstein, Physics Reports , 1 (1985). A. Bloch, P. S. Krishnaprasad, J. E. Marsden, and T. Ratiu, Annales de l’Institut Henri Poincaré (C) Analyse Non Linéaire ,37 (1994).A. J. Cerfon and J. P. Freidberg, Physics of Plasmas , 012505 (2011).M. S. Rosin, A. A. Schekochihin, F. Rincon, and S. C. Cowley, Monthly Notices of the Royal Astronomical Society , 7(2011).J. A. Byers, B. I. Cohen, W. C. Condit, and J. D. Hanson, Journal of Computational Physics , 363 (1978).D. Winske, L. Yin, N. Omidi, H. Karimabadi, and K. Quest, in Lecture Notes in Physics: Space Plasma Simulation , edited byJ. Büchner, M. Scholer, and C. Dum (Springer Berlin/Heidelberg, 2003) pp. 136–165.M. W. Kunz, J. M. Stone, and X.-N. Bai, Journal of Computational Physics , 154 (2014).D. Lynden-Bell and J. P. Ostriker, Monthly Notices of the Royal Astronomical Society , 293 (1967).J. L. Friedman and B. F. Schutz, The Astrophysical Journal221