Magnetic Reconnection in Plasmas (Cambridge Monographs on Plasma Physics) - PDF Free Download (2024)

the first term coming from interlinkage, the two others from the individual twists. Thus Hf = 0 2 = H, the helicity of the original twisted ribbon 0. [By the way, cutting a Moebius (T = | ) ribbon lengthwise results in a single ribbon twisted twice (T = 2) with helicity

2(0/2)2 = ±02].

This is a simple example showing that the helicity of a closed configuration can be expressed by those of the constituent closed flux tubes. If there are N constituents of flux 0,-, J2^=i 0i = 0? w e h a v e i,

(2.49)

where L^ gives the number of linkages of tubes / and j and Hi are the internal helicities of the constituents given by their twists, kinks and knots. Relation (2.49) also shows that for large N the contribution from the internal helicities becomes negligible, since if,- oc 0 2 = O(N~2), hence YltHi = OiN'1), while the first term, the interlinkage sum, is finite. We now return to the question raised previously, when we computed the helicity of a trefoil knot: Is the helicity conserved in reconnection? At first sight this does not seem to be the case. If we naively reconnect two interlinked field lines, we seem to be able to reach different configurations like those in fig. 2.10(a) or (b) where, in neither process, is H conserved. This would, however, imply an oversimplification of the concept of field lines, which in a conducting medium are to be viewed as thin flux tubes with an internal structure, the twist. One might think that in the limit of vanishing tube radius, field lines in the tube become parallel, such that the

27

2.4 Conservation of magnetic helicity in reconnection

(a)

(b)

Fig. 2.10. Different ways of reconnecting an unstructured field line.

Fig. 2.11. Ribbon model to demonstrate helicity conservation in reconnection (after Pfister & Gekelman, 1991). twist loses its meaning. The helicity of a tube does, however, not depend explicitly on the tube radius, as seen in (2.44). The helicity of a twisted tube is of the same order as that of two interlinked tubes of similar fluxes. We may say that field lines in a plasma are not "naked" but "dressed with twist".

28

2 Basic kinematic concepts

This property has to be accounted for in the reconnection process. Reconsidering the reconnection of two interlinked flux tubes, fig. 2.9, into a single tube (which is only possible if fa = 02 = ), we have to distinguish between the field lines which reconnect first and those reconnecting subsequently, i.e., treat the flux tubes as extended flux bundles. This process introduces a twist, which is, however, difficult to guess just by looking at the configuration. I therefore suggest that you perform the corresponding experiment using two interlinked closed untwisted paper ribbons as done by Pfister & Gekelman (1991) (fig. 2.11). To make the reconnection process definite, one should draw B vectors along the ribbons on both sides. The resulting ribbon has a twist T = 2, hence Hr = 22, see (2.44), and hence H! = if, the helicity of the original interlinked tubes (2.47). The reconnection implied in the transition from fig. 2.8(b) to 2.8(c) occurs between flux tube parts which can be considered plane. Reconnection proceeds in a natural peeling manner, which does not introduce additional field line twists. Though it seems to be difficult to give a mathematical proof we have at least made very plausible that reconnection does not change the helicity. The term on the right side of (2.16) can be made arbitrarily small by localizing the volume where En =£0.

3 Current sheets

In an electrically conducting magnetized plasma even slow motions do not, in general, preserve a smooth magnetic field, but give rise to sheetlike tangential field discontinuities, called current sheets. These are the natural loci of magnetic reconnection. In this chapter we consider the properties of current sheets and reconnection via current sheets in the traditional framework of resistive magnetohydrodynamics (MHD). This restriction is justified, since macroscopic current sheets mostly occur, when the reconnection process R in Ohm's law (2.1) is dissipative. Collisionless reconnection processes, which are treated in chapters 6 and 7, usually give rise to microcurrent sheets. The chapter starts with a brief introduction to MHD theory, discussing the basic equations, magnetostatic equilibria, and linear MHD waves. In low-/? plasmas it is often convenient to eliminate the fast compressional MHD wave by using a reduced set of equations. In section 3.2 we first consider the conditions under which a current sheet arises and how it is formed by rapid thinning, which continues until finite resistivity leads to a stationary sheet configuration. Then the structure of a resistive current sheet, called a Sweet-Parker sheet (Sweet, 1958; Parker, 1963), is analyzed. While the global properties of such sheets follow from the basic conservation laws, the detailed structure requires a more specific analysis. Section 3.3 deals with the role of current sheets as centers of reconnection in a magnetic configuration. Syrovatskii (1971) has developed a simple and elegant theory of current sheets which captures many features of the fully dynamic resistive theory, in particular the complicated structure of the sheet edges, the so-called Y-points. We then discuss the scaling laws of stationary current-sheet reconnection. Since in a macroscopic sheet this is a rather slow process, Petschek's reconnection model (Petschek, 1964), predicting a much higher reconnection rate, had attracted strong interest, being considered for many years as the fundamental reconnection 29

30

3 Current sheets

configuration. Though it has now become clear that Petschek's model is not valid in the most interesting limit of small resistivity, it may nevertheless be valid as a phenomenological model in the presence of a more efficient reconnection process. 3.1 Resistive MHD 3.1.1 The MHD equations The MHD equations can be obtained in a systematic way from Braginskii's two-fluid theory (see chapter 6) at sufficiently large spatial scales. Here, we introduce these equations in a simple heuristic way. The acceleration of a fluid element is due to the sum of the pressure gradient and the magnetic, gravitational and viscous forces, p-r = ~VP + - j x B + pg + p/iV2v, (3.1) at c where d/dt = dt + v • V, and p is the mass density, which obeys the continuity equation dtp + V • py = - / + p V • v = 0. (3.2) at On macroscopic scales the plasma is electrically neutral, there are just as many positive as negative charges in a fluid element, hence the charge density vanishes and there is no direct electric force in (3.1). Since the gravitational force pg is only important in certain astrophysical systems, we will usually ignore it in (3.1). The viscous term contains only a scalar kinematic viscosity /a, ignoring the anisotropy of the viscosity tensor induced by the magnetic field. In this form the viscosity is not meant to describe real momentum transfer but only to prevent singularities of the flow. The pressure tensor is assumed to be isotropic, corresponding to conditions close to local thermodynamic equilibrium. In most cases of interest, plasmas are sufficiently dilute to follow the ideal gas equation p = (m + ne)kBT, where n^e are the particle densities of ions and electrons, ks is the Boltzmann constant, and p ~ m^ is the mass density. The Lorentz force in (3.1) can also be written in the following way:

- j x B = - v f - + bb (B VB) + (I - bb) • (B • VB), C

o7T

(3.3)

b = B/B, where we have used Ampere's law

J = fVxB,

(3.4)

3.1 Resistive MHD

31

i.e., the magnetic force acts as an isotropic pressure added to the thermal pressure p, a tension which strives to shorten the field lines and a stress which tends to straighten the field lines. On the large scales over which the MHD approximation is valid, heat conduction effects are negligible, hence changes of state are adiabatic,

or, using the continuity equation (3.2), dtp + v • Vp + ypV • v = 0,

(3.5)

where y = Cp/Cy, Cp,Cy = specific heats, y = 5/3 in the usual case. The change of the magnetic field is given by Faraday's law (2.3) with E from Ohm's law (2.1), dtB = V x (v x B) - cV x R

(3.6)

The basic assumption of the resistive MHD approximation is that collisions are sufficiently frequent, such that the collisional effects in R dominate, of which resistivity is usually the most important one, R = r,j.

(3.7)

Here, too, a scalar resistivity r\ is assumed, which is, however, a much less drastic approximation than the assumption of a scalar viscosity, since f/jL and j/|| differ only by a factor of order one, t]± ~ 2ri\\ (see Braginskii, 1965). Inserting Ampere's law, (3.4), (3.6) becomes dtB = V x (v x B) + V • (vlVB).

(3.8)

The coefficient X = rjc2/4n has the dimensions of a diffusion coefficient and is called the magnetic diffusivity. As an alternative to adiabatic changes of state, (3.5), one often assumes incompressible motions V • v = 0. The conditions for the validity of this assumption are discussed in section 3.1.2. For V • v = 0, the continuity equation reduces to the simple convection equation dp/dt = 0, such that (at least in the absence of a gravitational force) we can take p = po = const for simplicity. In the incompressible case the total pressure P = p + B2/Sn is no longer an independent dynamic variable, but is determined from Poisson's equation obtained by taking the divergence of (3.1). Taking the curl of this equation eliminates the pressure. The incompressible MHD equations can thus be written in terms of the vorticity w = V x v and B: + v • Vw - w • Vv = — ( B • Vj - j • VB) + juV2w, cpo dtB + v VB - B Vv = AV2B.

(3.9) (3.10)

32

3 Current sheets

These equations assume a more symmetric form if written in terms of the Elsasser fields (Elsasser, 1950): z ± = v + B/v/4^po,

(3.11)

which gives Stz* + z+ • Vz± = - — V P + ^+V2z± + /i_V2z+, Po

(3.12)

3.1.2 Equilibrium and linear modes Magnetic plasma configurations, both in laboratory and in space, vary in time because of changes of the external circuit and because of internal transport processes. These time-scales are usually long compared to the dynamics described by the MHD equation (3.1), such that the system changes quasi-statically evolving through a sequence of states which satisfy the equilibrium equation V p - - j x B = 0,

(3.13)

neglecting the gravitational force. Equation (3.13) can also be written in terms of the total pressure P = p + B2/8n: VP--B

4TT

VB = 0.

(3.14)

From (3.13) it follows immediately that p is constant along magnetic field lines and along current density lines, B • Vp = 0,

j • Vp = 0.

Since, for a finite pressure gradient, B and j are not parallel, we can locally construct a surface \p(x9y,z) = const, where p is constant, p = p{\p). In an equilibrium configuration this property must, however, hold globally, i.e., there must be surfaces of constant pressure on which a field line runs on for ever, a condition which is rather stringent for general nonsymmetric configurations. That infinitely long field lines span smooth magnetic surfaces can be shown directly for configurations with a continuous symmetry, where a coordinate system £,*/,£ exists such that all physical quantities including the elements of the metric tensor depend only on two coordinates £,?/. If this condition holds, (3.13) can be reduced to an elliptic differential equation for a scalar function

3.1 Resistive MHD

33

such that surfaces xp = const have the desired properties of magnetic surfaces (Edenstrasser, 1980a). It turns out that the requirements on the coordinates £,*7,£ are rather restrictive. The most general case is that of helical symmetry (Edenstrasser, 1980b), of which plane and axisymmetry are special cases. One-dimensional equilibria can easily be obtained analytically. In plane geometry, with p(x),B = (0,By(x\Bz(x)), the configuration is a sheet pinch, which follows the equation

since B • VB = 0. In cylindrical geometry, where we have p{r),B = (0,Be(r),Bz(r)), the most general configuration is a screw pinch, dp

1

Tr = c A quantity often used in laboratory plasma physics is the safety factor q measuring the field-line twist T,

where 2nR is the length of the plasma column. A screw pinch with q ~ 1 serves as a simple model of a tokamak plasma, while a reversed-field pinch has q [y/(y — l)]po/gPo- Only in the special case y = 1, i.e., for isothermal conditions, is the extent not limited, p(z) = p 0 exp{-gp o z/p o }.

(3.21)

Given a certain equilibrium it is important to know how it responds to small perturbations. If the equilibrium is stable, the perturbation generates a propagating wave, which is most conveniently discussed for a hom*ogeneous plasma embedded in a uniform magnetic field. By linearizing the MHD equations writing B = Bo + B, p = po 4- p, p = po + P, v = v, the perturbations B,v,p,p can be Fourier transformed in space and time, B oc exp{ik • x — icot} etc., such that the differential equations reduce to a hom*ogeneous algebraic form, the solubility condition of which gives the dispersion relation for three types of modes a>i, 002,003: (a) The shear Alfven wave

(oj = kyA,

(3.22)

where vA = Bo/\/4npo

(3.23)

is the Alfven velocity and k\\ = k • b = k cos 6. The plasma motion v associated with this mode is incompressible k-v = 0, i.e., the mode is transverse, v is also perpendicular to Bo, and so is the perturbation of the magnetic field B = +^/4npov, which causes a bending of the field lines. Hence we see that the Elsasser fields, (3.11), describe Alfven waves propagating in the direction of or opposite to the mean magnetic field Bo. (b) The (fast) magnetosonic wave, also called compressional Alfven wave,

[v2A + c] + yl(v\ + <W-4o\ ((02/k)2 > max{i;,4,cs}. Here, the left-hand equality refers to perpendicular

3.1 Resistive MHD

35

Fig. 3.1. Phase velocities v\^ = (o\^(B)/k of the MHD modes (3.22), (3.24), (3.25). propagation k _L Bo, where the mode is longitudinal v || k leading to pure field compression. Since both B2 and p are simultaneously compressed, the restoring force and hence the frequency are particularly high. For parallel propagation, the mode degenerates to the shear Alfven wave, (3.22), for the usual case of VA > cs in a magnetized plasma. The phase velocity of the slow mode satisfies 0 < (co^/k)2 < min {VA,CS}. In the case of perpendicular propagation (left-hand equality) it is longitudinal, but since the variations of B2 and p are now just opposite in phase, SB2/$n = —dp, the restoring force and hence the frequency vanishes, corresponding to a quasi-static change of equilibrium. For parallel propagation (right-hand equality) the slow mode is again longitudinal, the nonmagnetic sound wave co2 = k2c2. (In the case cs > VA the parallel modes of (b) and (c) just interchange roles.) MHD waves are nondispersive, the phase velocity co/k equalling the group velocity dco/dk which depends only on the angle 9 between k and B. This is to be expected, since ideal MHD theory contains no intrinsic spatial scale. The frequencies co/(k) of the three modes are illustrated in fig. 3.1.

36

3 Current sheets

At finite amplitude, compressional waves steepen to form shocks. Magnetosonic shocks are generated as blast waves in eruptive events, for instance in a solar flare or a supernova explosion, or as stationary standoff shocks in front of an obstacle in a super-Alvenic plasma stream such as the Earth's bow shock. In reconnection theory, slow-mode shocks, i.e., finite amplitude slow modes, are more important. They play a crucial role in Petschek's model (section 3.3.4). In an inhom*ogeneous plasma the three types of modes are in general coupled, but they remain nevertheless useful for a qualitative classification of small perturbations. While a hom*ogeneous system is stable, co being real, as follows from the dispersion relations (3.22), (3.24), (3.25), sufficiently strong pressure and current density gradients give rise to instabilities. (For a quick overview of MHD stability theory see Biskamp, 1993a, chapter 4.) Finally, we briefly discuss the validity of the incompressibility assumption. Though this assumption is often made for linear modes, it is basically a consequence of the nonlinear equations since, for instance, any linear perturbation in a neutral fluid corresponds to a (compressible) sound wave. As a general rule, the dynamics of the fluid can be assumed incompressible, V • v = 0, if v is small compared with the (fastest) compressional mode in the direction of v. In addition, the time variation of v must be weak, i.e., the frequency must be small compared with the frequency of the compressible wave, such that the dt terms are small compared with the advective terms v • V. Consider first a neutral fluid or a high-/?-plasma. From the pressure equation (3.5) we have ypV • v ~ —v • Vp. Substituting Vp from the equation of motion (3.1) for high /?, — Vp ~ pv-Vv, we find

-Vt,~4f

c2s L where L is a typical gradient scale-length. In a low-/? plasma one has to distinguish between perpendicular (to B) and parallel motions. While in the latter case (3.27) applies, in the former case Faraday's law (3.6) gives B2V - v ~ - v • VJS2/2. From the equation of motion (3.1) for low /J one has pv • Vi;2/2 ~ - v • VB 2 /8TT, hence

3.1 Resistive MHD

37

As discussed in more detail in section 3.1.4, the perpendicular motion is usually incompressible or quasi-incompressible, while the parallel motion is often compressible. In MHD turbulence theory, one often assumes full incompressibility V • v = 0 (see section 5.4). While this assumption is valid in 2D, if one can rely on the presence of a strong axial field, in the general 3D case the assumption is less well founded and more a matter of convenience. 3.1.3 Conservation laws The resistive MHD equations (3.1), (3.2), (3.5), (3.8) satisfy a number of conservation laws, which can easily be derived. The mass M in a volume V(t) moving with the fluid is conserved dM dt

d

= i. /

ai jv(t) at Jv(t)

pdV = 0.

(3.29)

The change of momentum P is due to the total stress acting on the surface F of the volume, ^ = *- [ pvdV = - / ST - dF. dt dt Jv(t) JF Here 2T = {Tij} is the stress tensor Tij =\P + ^n) hj ~ -^BiBJ

+ \vWvj

+ dPil

(3.30)

( 3 - 31 )

consisting of the thermal, the magnetic, and the viscous stresses. The energy balance is dW _ d r dt dt Jv(

= -£(v#~)'dF-ri

f j2dV - \i f p Y^i^jfdV.

(3.32)

In an isolated system, where p, v, B vanish at the boundary, the momentum is conserved and the energy decays due to ohmic and viscous dissipation.* In most cases of interest, however, B is not zero at the boundary and there is an exchange of momentum and energy with the external system. Also, the magnetic helicity H = J A • BdV is conserved for r\ = 0 as discussed in section 2.1. Note that this conservation law is a property of * The boundary terms also vanish for periodic boundary conditions, which are usually assumed in studies of hom*ogeneous turbulence.

38

3 Current sheets

the frozen-in magnetic field, which is valid for any velocity v, not only for motions satisfying the MHD equations. In the incompressible case V v = 0, there is an additional ideal invariant of the MHD equations, the cross-helicity K = f \-BdV. IV(t) Jv(t)

(3.33)

The cross-helicity of a flux tube moving with the fluid is conserved up to resistive and viscous effects, 1 _ = -(ji + X) / > (diV,)(diB,)dV.

dt

(3.34)

J jf

Note that the kinetic helicity HK — / v • (V x vW3x H 35) which is an invariant of the Euler equation, is not conserved in MHD. Finally, there is an important linear invariant, if the dynamics is restricted to helical geometry, i.e., if all physical variables depend only on r and 8 + az. Such motions occur in the development of a kink instability in a cylindrical plasma column, which plays an important role for instance in the sawtooth phenomenon in a tokamak plasma or the instability of a twisted coronal loop. Integrating Faraday's law (2.4) for r\ = 0, dtA = v x (V x A) - V0,

(3.36)

one can easily derive that the helical flux function ifH = ocrAg — Az

(3.37)

is conserved along the orbit of a fluid element: •+V-VV>H = 0.

(3.38)

From (3.38) it follows that for incompressible motions the total helical flux Q>H is conserved, dt

a

In the limiting case of plane geometry, a = 0, xpn reduces to — Az.

(3.39)

3.1 Resistive MHD

39

3.1.4 Reduced MHD models In many cases the resistive MHD equations are not considered in full generality. Instead, approximations are made concerning both physics and geometry. We often encounter the situation that the system is embedded in a strong magnetic field Bo generated by external currents. Such a field is called a potential field, V x Bo = 0. By "strong" we mean that thermal and kinetic plasma energies are small compared with the magnetic field energy:

p~ pv2 < Bl/%n,

i.e., the ratio of thermal to magnetic energy, called /?, is small,

l < 1.

(3.40)

Bo is, in general, not hom*ogeneous. The gradient of the field intensity is related to the field line curvature K = b • Vb by the equation V ± B 0 = S O K,

(3.41)

which is easily derived from Bo x (V x Bo) = 0. Since for fl < 1 the frequencies associated with the plasma motions are much slower than the frequency of the magnetosonic mode, dt < t^Vj_, the latter decouples from the plasma dynamics and can be eliminated (Drake & Antonsen, 1984). This implies that the inertia and the viscous terms in the perpendicular part of (3.1) can be neglected, such that there is instantaneous equilibrium perpendicular to the field, | j

- ^B2K

=0

(3.42)

or decomposing B into a large potential field Bo and a component Bi generated by plasma currents, B = Bo + Bi, and using (3.41). Hence we find that the variation of the longitudinal part Bi\\ = Bi • Bo/Bo is small, 5Bn~~5p.

(3.44)

The perpendicular velocity V_L is no longer determined directly from the force balance. Instead, the incompressible (or solenoidal) and the irrotational parts of vj_ are computed separately. The incompressible part is obtained from the parallel component of the curl of (3.1)

" ( v x pJt)

b

=

Tv" I +

2(b x K) Vp + b (v x

'

"

(3>45)

40

3 Current sheets

using the force balance, (3.42), to express VB2 in terms of Vp. To derive (3.45) use j = b/'n + jj_ and the quasi-neutrality condition V • j = 0. To compute the compressible part of V_L, multiply Faraday's law (3.8) by B to obtain (neglecting resistive effects) \{dt + v ± • V)B2 = -B2(V •

VL

+ K • v ± ),

where we have used v±

= yE = c - ^ - ,

(3.46)

obtained from Ohm's law (2.2). Substitution of ¥j_B2 from (3.42) gives \dtB2 - AnyL • Vp = -B 2 (V • v ± + 2K • v ± ).

(3.47)

In a low-/? plasma the terms on the l.h.s. are proportional to p by use of (3.44), hence we find to lowest order in /? V ••]_ = - 2 K • v ± .

(3.48)

We call this property quasi-incompressibility. The result can easily be interpreted. For low /? the perpendicular field E^ in v# is essentially electrostatic, 1 (3.49) c since 3tA_L oc 8B\\ gives a correction O(/3) to the electrostatic field, hence

n-c^.

(3.50)

Taking the divergence and using (3.41) it is easily seen that v# satisfies (3.48). The parallel dynamics is determined by the parallel component of (3.1) h

To further simplify matters we make use of the fact that B± —oo. While the full solution ^{t),r\(t\ in particular the value of to, must be computed numerically, the most interesting behavior in the vicinity of the singularity can easily be obtained analytically. Since £ —• 0 for t —• to, while r\ reaches a finite value rj(to), as can be checked a posteriori by inserting the solution £(i)9 integration of the £ equation in (3.78), I ~ -ri(to)/e, gives

{^Qij(to)) and hence y^ — pec

2 1 , 3 to — t 1

fa>-02/3, 5 ^ const,

(3.79) (3.80)

Note that, even if both 70 and 0 implies unlimited expansion, while for yo = .

3.2 Resistive current sheets

47

This system has the similarity solution y2

\

(183)

W))' = A(t)xy,

(3.84)

corresponding to j = j(t) and w = 0. (Note also that the compressible flow, (3.73), is irrotational.) From the xp equation follows t = -A£,

ii=Ati,

(3.85)

where A(t) is a free function. Since (3.82) is satisfied identically, the solution is not completely determined, in contrast to the compressible case, (3.75)-(3.77). A special choice is A = const, Z=

fo-A\

fl = meA\

(3.86)

corresponding to exponential flattening of the X-point configuration and exponential growth of the current density j oc e2At. Numerical simulations (Sulem et al.9 1985; Friedel et ah, 1997) show that this behavior is indeed typical in 2D MHD. The result can be understood in terms of a continuous weakening of the flow nonlinearity, as the configuration becomes more and more one-dimensional, as can be seen directly from the 2D MHD equations. Taking the Laplacian of (3.70) gives

8tj + v ' V; - B • Vw = 2 ] £ ez • (V^y x V5^) = 2 ]T emd& • dkv. (3.87) i,k

i

Combining this equation with (3.82) and introducing the Elsasser variables z^, (3.11), and their derivatives w± = w + j9 we obtain dtw* + z T • Vw± = ^

= ±2^2 eiikdiB • dk\.

(3.88)

Uk

For a sheet-like structure only the cross-sheet derivative becomes large. Hence, because of the cross product, the source term in (3.88) contains only one large derivative, such that the equation can be written symbolically in the form w = Cw, whence exponential growth. The presence or absence of a finite-time singularity (FTS) in a fluid model has been the subject of considerable debate, since this property may shed some light on the process of turbulence generation from a smooth initial state and even elucidate the character of the turbulence itself. The simplest example is the compressible flow in the ID Burgers' equation dtu + udxu = 0, where the presence of FTS can easily be shown. Taking the x-derivative yields dtw + udxw = w = —w2,

w = dxu,

(3.89)

48

3 Current sheets

which has the solution w[xo(£)] = (* ~~ ^o)"1- The simplest case of incompressible flow, the 2D plane Euler flow, does not exhibit FTS, as is to be expected from the equation dtw + v • Vw = w = 0. This is, however, no longer true for more general 2D Euler flows, since, for instance, a rotationally symmetric flow can lead to FTS (Grauer & Sideris, 1991). For nonsymmetric, i.e., three-dimensional Euler flows, FTS seems to be generic (except for high-symmetry initial conditions such as the TaylorGreen vortex, see, e.g., Brachet et a/., 1992), which has recently been studied numerically using the method of adaptive mesh refinement to reach the necessary high local resolution (Grauer et a/., 1998). By contrast, 3D MHD flows show only exponential growth, much as in the 2D case, since narrow structures always become sheet-like eventually, even if initially the magnetic field is weak (Grauer & Marliani, 2000). The absence of a finite-time singularity seems to weaken the relevance of the similarity solution (3.86) since, for long times, spatial regions far away from the X-point, where the behavior of the solution is not universal depending on the particular configuration, may affect the behavior at this point. In practice, however, exponential growth up to the resistive saturation level is sufficiently fast that non-universal effects remain small. Though the similarity solutions given above provide strong evidence for rapid current-sheet formation, their validity is restricted to the center region of the sheet, the immediate vicinity of the X-point, where j ~ const and w ~ 0. As seen in the following sections, current sheets in their full extent are also the location of strong vorticity. In contrast to the current distribution in the sheet, which has a flat elliptical monopole character, the vorticity distribution has the structure of a flat quadrupole (hence w = 0 in the center). This behavior of j and w can be interpreted in terms of the variables w± = w + j , such that w = | ( w + + w~), j = \(w+ — w~~). In a current sheet, w + and w~ both form elliptical structures, slightly tilted against each other. (Since the ideal MHD equations are symmetric in the Elsasser fields, w + and w~ must be similar in contrast to j and w.) Hence one of the two fields, w or j , must have a quadrupole character. Numerical simulations reveal that current sheets indeed develop in a selfsimilar way (Biskamp, 1993b). Figure 3.3 illustrates the thinning process. The system starts from the smooth state at t = 0 (configuration Ai from Biskamp & Welter, 1989) {x, y) = cos(x + 1.4) + cos(y + 0.5), xp(x9 y) = cos(2x + 2.3) + cos(y + 4.1),

(3.90) (3.91)

3.2 Resistive current sheets

49

(a)

(b)

(c)

V

J

Fig. 3.3. Development of current sheets from the initial configuration, (3.90) and (3.91): (a) t = O;(b)t= 1.0; (c) t = 1.3 (from Biskamp, 1993b). which is a generalization of the Orszag-Tang vortex (Orszag & Tang, 1979). Four current sheets develop at the X-points of this configuration. Since relaxation to equilibrium inside the narrow sheets occurs quasiinstantaneously, the magnetic field obeys the equilibrium relation B • V/ = 0, such that j = j(xp). In general, j(ip) is not identical on both sides of the sheet because of the asymmetry of the configuration, but has two different branches j±(xp). Figure 3.4 gives scatter plots In j vs xp for one of the sheets taken at four different times plotted on top of each other, where 1 = j/jmax(t),

xp = oc(t)(xp - xpo)9 and ; m a x = j(y>o)is

the

maximum current

50

3 Current sheets

u

i

A

¥

*

Fig. 3.4. Normalized In j versus xp scatter plots taken at t = 1.0,1.1,1.2 and 1.3, plotted on top of each other, for the upper-right current sheet of fig. 3.3. The straight full line is j(\p) for the exponential profile, (3.92), the dashed line for a Gaussian profile, (3.94) (from Biskamp, 1993b). value. In good approximation the points fall on two straight lines, i.e., the two branches of j±(\p) are j± ^ jmax{t)exp[-2cc±(t)(\p

-

(3.92)

The maximum current density jmaxM is found to increase exponentially in time in agreement with the similarity solution (3.86), and the thickness S±(t) oc a+(£)~1//2 decreases exponentially. In configuration space the profile j = 7maxexp(—2<xxp) corresponds to jo(x) =

jmaxsech2(x/8).

(3.93)

As seen from fig. 3.4, the numerically observed current profiles are distinctly different from a Gaussian profile, the dashed line in fig. 3.4, which is obtained from the parametric representation = ax erf(x/a)

(3.94)

3.2 Resistive current sheets

51

where a is chosen such that j(\p) has the same derivative at the origin as the corresponding branch of the scatter plot. The current profile j = exp(—2xp) can be "derived" theoretically (Montgomery et al, 1979; Taylor, 1993). Here the idea is to use statistical arguments to determine the most probable profile, which should be set up approximately under typical conditions. Using a thermodynamic approach one would calculate this state as the minimum of the free energy 3F = W — ST under certain constraints. To avoid functional integrals we consider a discrete model (Montgomery et a/., 1979): divide the area F containing the current sheet into M cells AF, F = MAF, and decompose the current density into N filaments of current k such that / = Nk is the total current flowing across F with N >• M. For a given current distribution there are Nt filaments in cell i. Assuming equal a priori probability of states the probability P{Nt} to find a certain set of occupation numbers {Nt} is proportional to the number of states with these occupation numbers, i.e.,

P{Nt} = nWO" 1 -

(3-95)

i

The entropy of the system is S = —lnP ~ £),• JV,-In iV,-, using Stirling's formula. The energy W is

W = J(Vxp)2d2x = -Jxpjd2x = -J

Jd2xd2x'j(x)g(x,x')j(x'\

where g(x,x r ) is the Green's function, xp(x) = J g(x, x!)j(x!)d2xf. In discretized form we have xpt = k Y,j gijNt and therefore jNj.

(3.96)

We now minimize the free energy under the constraint of constant total current by varying the occupation numbers Nu =0, j

(3.97)

J

where a is a Lagrange multiplier. Variation gives TlnNt + cc = -2k2 ] T SijNj = -l^xpu j

or, passing to the continuum limit, j = CQ-2V/T

(3.98)

52

3 Current sheets

with C = (/l/AF)e~ a/r . Hence xp obeys the differential equation V 2 V = Ce~2xp/T.

(3.99)

For the simplest case of a symmetric current sheet the solution of (3.99) is given by the expression (3.93) with 5 = (T/C) 1 / 2 , C = j m a x . Hence the "temperature" T is a measure of the sheet width. It is interesting to consider the Fourier transform of the profile (3.93), Z100 coskx 1 nk / j - dx = -—n . (3.100) 7-oo cosh x sinh(^nk) Hence the magnetic energy spectrum is asymptotically exponential, B\k) = j^~>=

.

*

„ * 4e"*,

\k\ > 1,

(3.101)

which agrees with the numerically observed energy spectrum (Frisch et a/., 1983). By contrast, a Gaussian profile would lead to a Gaussian Fourier spectrum. 3.2.2 Basic properties of resistive current sheets The growth of the current density during sheet thinning saturates because finite resistivity leads to field diffusion, which balances the convective magnetic flux transport into the sheet. The result is a quasi-stationary resistive current sheet, called a Sweet-Parker sheet (Sweet, 1958; Parker, 1963). Resistivity introduces reconnection - magnetic flux is transported into the sheet, reconnected and swept out of the sheet. The configuration is illustrated in fig. 3.5. Assuming incompressible motions, the stationary quasi-one-dimensional state is characterized by six quantities, three describing the dynamics: (i) the (poloidal) magnetic field Bo immediately outside the sheet, called the upstream field (the downstream field is negligible), (ii) the upstream flow uo perpendicular to the field, and (iii) the downstream flow VQ along the field taken at the sheet edge; two geometric quantities: (iv) the sheet length A and (v) the width S (in the literature one often finds a different terminology S = thickness, A = width, implying a three-dimensional sheet structure with the length measured along the third direction); and, finally, (vi) the diffusion coefficient r\. We will see below that the effect of the viscosity is usually weak. These six quantities are connected by three relations derived from the continuity equation, Ohm's law, and the equation of motion, assuming stationarity. Furthermore, the analysis implies that the profiles of the dynamic quantities along and across the sheet have certain self-similar

53

3.2 Resistive current sheets

Fig. 3.5. Schematic drawing of a Sweet-Parker sheet. properties. Integrating the continuity equation V • p\ of the sheet and assuming constant density gives = v0S.

0 over a quadrant (3.102)

Consider Ohm's law (2.1) along the x-axis. Stationarity requires Ez to be uniform in space. In the upstream region outside the sheet, where the current density is small, the resistive term is negligible, while in the sheet center, where j = j m a x , and the velocity vanishes, the resistive term dominates, which gives the relation u0B0 = /77'max ^

rjB0/d.

(3.103)

Since, usually, uo < Bo as can be seen a posteriori, the inertia term is negligible in the force balance across the sheet. Hence dx(p + \B2) = 0, which gives ^ 2 p 0 . ( 3 - 104 ) Here po is the upstream fluid pressure and p m a x the maximum pressure in the sheet center, where the (poloidal) field vanishes. Now consider the force balance along the mid-plane of the sheet. Since Bx is negligible, so is the magnetic force, such that only the pressure force accelerates the fluid along the sheet, vydyvy = —dyp. Integration between center and edge yields \»l =Pm. For stationary conditions, (3.65) and (3.66) read m

n

dx(f>dyW - dy(f)dxw - dxxpdyj + dyxpdxj = n(dxw + d2yw), - dy(j)dx\p = t](d2xy) + d2yxp) - E.

(3.113) (3.114)

56

3 Current sheets

Here E = Ez = dtxp is the axial electric field, which is uniform and constant in the stationary case. Since at the neutral point, which in our symmetric case is also a flow stagnation point, B = v = 0, (3.114) gives */(V20 + V02) = E.

(3.115)

Differentiating (3.114) twice with respect to x and twice with respect to y at the origin one obtains, respectively, + */(V40 + W22) = 0, rj(ip22 + t/>04) = 0.

(3.116) (3.117)

Differentiating (3.113) once with respect to both x and y at the origin gives -V20(^22 + V04) + V02(V40 + V22) = M051 + 2033 + 015),

(3.118)

which by use of (3.116) and (3.117) becomes 4 --011V>2OV>O2=/J(051 + 2033 + 015)-

(3.119)

Cowley (1975) considered the inviscid case ju = 0, where the r.h.s. of (3.119) seems to vanish. Assuming 0 n ^ 0, i.e., stream lines to form hyperbolae, either ip2o or tpo2 (not both because of (3.115)) must vanish. This would imply that field lines do not form hyperbolae, in particular the separatrix branches do not intersect at a finite angle but osculate, meaning that the normal magnetic field component does not increase linearly but cubically Bx oc y3. The argument is, however, not strictly valid (Uzdensky & Kulsrud, 1998), since the solution, in particular the velocity field, does not remain regular in the inviscid limit, such that the Taylor expansion (3.112) breaks down. Uzdensky & Kulsrud perform a boundary layer analysis for small viscosity / * < * / , where the velocity exhibits a viscous sublayer. Asymptotic matching of the solutions in the outer inviscid region and the inner viscous region yields that l i m ^ o ^051 is finite, a result which they could corroborate and quantify by a fully numerical solution. It follows that, in general, both tpo2 =/= 0 and \p20 ^ 0, their values depending on the solution of the entire current sheet with the upstream field Boy(y) as boundary condition. Hence in the limit \i —• 0 the separatrix branches do not osculate, but intersect at a finite angle, which is of course very small in a long current sheet. For finite viscosity \i ~ r\ both 0 and \p are regular with converging Taylor series, (3.112). But now one cannot draw definite conclusions from relation (3.119). Numerical simulations show that, although ip2o (or T/>02) does not vanish, it is small, hence the separatrix branches nearly osculate.

3.2 Resistive current sheets

57

To consider the global sheet configuration (apart from the edge regions, where a complicated behavior occurs, which is discussed separately in section 3.3.2), we can use a Taylor expansion in y along the sheet, which is valid for finite /n and r\: xp(x9 y) = ipo(x) + y2xp2(x)/2 + ..., (x9y) = y4>i{x) + y*fc(x)/l

! + •...

(3.120) (3.121)

Insertion into (3.113) and (3.114) gives a set of differential equations for the coefficients \pn(x)9(j)n(x). If xpo(x)9 or the current profile jo(x) = \PQ(X)9 is given, the remaining functions ip2(x),...,4>i(x),... can be determined successively in terms of y)Q. Let us assume the profile x/(5),

(3.122)

whence integration gives Voto = ^ax<S2ln[cosh (*/. As an approximation we could try to assume one-dimensionality for xp = xp(x). This would, however, imply that the Lorentz force term B • Vj in (3.65) vanishes and hence w = 0, which contradicts the fact that a current sheet is also the location of high vorticity. Instead, we may consider (3.66) only along the symmetry axis y = 0, where vy = 0, assuming a linear velocity profile vx = —ux/L with some length-scale L. In this case Ohm's law,

E - —xB = r\B\ JLJ

B = \p\x\

3 Current sheets

58

Fig. 3.6. Theoretical current-sheet profiles: (3.125) (full), seen2* (dashed), Gaussian (dot-dashed).

has the solution (Sonnerup & Priest, 1975):

B

=

n

i:

S2=EL/u.

(3.125)

Since here B —> 0 for x2 —» oo, the total current vanishes and j(x) has negative ("return") current parts adjacent to the central positive one, as seen in fig. 3.6. In typical resistive current sheets, however, the net current is finite and the negative contributions are small or absent, hence the solution (3.125) can, at most, model the central positive part of a current sheet. Alternative current profiles are a Gaussian e~~* or the sech2x profile, (3.93), which are also shown in fig. 3.6. To compare these with the current profiles obtained in numerical simulations, consideration must be restricted to the central part |x| ;$, 2(5, since only here simulations exhibit a self-similar profile shape. It appears that the numerical results agree best with a Gaussian, though the difference between Gaussian and sech2x is small.* More importantly, the profile does not depend on the type of resistive dissipation, a hyperresistivity r\2, obtained by the substitution r\V2\p —> — 7/2VV, gives the same profile as a normal resistivity r\. This shows, that the internal structure of the current sheet is primarily determined by the flow dynamics, and not by the character of the dissipation. * This behavior seems to be at variance with fig. 3.4, where the profile (3.92), i.e., sech2x, agrees significantly better with the simulation result than a Gaussian. However, fig. 3.4 refers to the self-similar behavior in the ideal phase of current-sheet formation, while the present results deal with stationary resistive sheets. We could have chosen a Gaussian in (3.122), but the sech2x profile is algebraically more convenient.

3.3 Driven current-sheet reconnection

59

3.3 Driven current-sheet reconnection In reconnection theory, one traditionally distinguishes between driven and spontaneous processes, where, however, the definition of these terms varies considerably. Driven reconnection refers to open, externally forced systems, while we call reconnection spontaneous if it arises by some internal instability or loss of equilibrium, the dynamic evolution of which depends only weakly on the external coupling. Hence in the former case the dynamics is determined by the boundary conditions, while it is rather insensitive to their choice in the latter. The concept of driven reconnection can, however, be applied more generally. In many systems, reconnection occurs in well-defined localized regions of space. If we are primarily interested in the physics of the reconnection process, we may restrict consideration to a small region of linear dimensions L around such a location instead of the global system of scale A > L. On the other hand L should not be too small, for instance L > A in the case of a single current sheet. The main advantage of the restriction to the subsystem L is that it allows us to simplify the geometry. In addition we can assume stationarity, even though the global system A is usually not stationary. Since the coupling of the L-system to the A-system occurs through the boundary conditions imposed on the L-system, which change on the global time-scale ~ A/VA, while the subsystem adjusts to these changes on the much shorter timescale ~ L/VA, we may consider the latter as stationary (if a stationary state exists). Boundary conditions are as indicated, for example, in fig. 3.2. Fluid and magnetic flux are injected from above and below, while the fluid leaves the system on both sides horizontally carrying along the reconnected field. The small region where resistivity is important is called the diffusion region, which is surrounded by the ideal external region, where resistive effects can be neglected. What can be learned from considering a driven reconnection configuration? The interpretation of such a model has led to some confusion. Since a fundamental issue of reconnection theory is to account for rapid processes in systems where the resistivity is very small, a figure of merit of a reconnection model is a weak or possibly no dependence of the reconnection rate on the resistivity. In stationary driven reconnection this point needs further specification, since the reconnection rate is determined by the boundary conditions for the inflow velocity and magnetic field intensity and hence is independent of r\ by definition. For the reconnection process to be independent of r\ one has to require instead that at fixed boundary conditions the configuration remains unchanged if Y\ is varied, at least for sufficiently small values of r\. A consequence of such behavior would be that the ratio of outflow energy flux to input energy

60

3 Current sheets

flux should be independent of r\ (essentially unity) and so should be the ratio of the energy dissipation rate to the input energy flux (essentially zero). 3.3.1 Syrovatskii's theory of current sheets Before we deal with the full resistive dynamics, which in general requires a numerical approach, we discuss a particularly simple theory due to Syrovatskii (1971), which deals with the ideal equations but allows sheetlike singularities. Though the theory does not describe real reconnection dynamics, it accounts for many important features thereof. The basic equations use a vanishing pressure gradient Vp = 0 and so are somewhat different from those of two-dimensional incompressible MHD, to which the major part of this chapter is confined. Since, however, the dynamics is only computed a posteriori from the change of the magnetic field, the difference is of no importance. The main assumption is that all currents in the system are localized in isolated current sheets. Hence \p satisfies Laplace's equation V2tp = 0. (3.126) The solution is determined by the boundary conditions. If these vary in time, xp obtains a parametric time dependence xp{x,y, £), which then determines the perpendicular component \± of the velocity (3.127) from the frozen-in condition dtxp + v • Vxp = 0,

while the parallel component v\\ is calculated from the equation As I At xVxp=0,

(3.128)

which follows from the equation of motion (3.71) for Vp = 0. (Equation (3.128) implies that the current density and hence the Lorentz force does not vanish identically. Equation (3.126) has to be regarded as an approximation in the sense that the effect of the distributed currents is small compared to that of the sheet currents). Hence the evolution of the system is described by (3.126), i.e., a simple boundary value problem. The crucial point is that regular solutions \p(x9y,t) do not always exist. If the change of the boundary conditions leads to a change of \p, dtip =£ 0, at a neutral point V\p = 0, the frozen-in condition cannot be satisfied at this point, which is therefore called a singular neutral point. In the presence of such points current singularities are expected to appear. To discuss these

3.3 Driven current-sheet reconnection

61

singularities in more detail it is convenient to introduce the potential F(z) in the complex plane z = x + iy, F(z, t) = ip(x, y, t) + ix(x, y, t\

(3.129)

which is analytic in the region considered except for isolated singular points and branch cuts. The conjugate harmonic function % is determined by the Cauchy-Riemann relations SxX = -dyV>,

Syx = dxxp.

(3.130)

From (3.129) and the definition B = ez x Vip we obtain the magnetic field in the form = By + iBx, (3.131) as can easily be verified by choosing a special direction of the derivative, e.g., dF/dx, since owing to (3.130) the complex derivative is independent of this choice. If z = zo(O is the position of a neutral point, dF/dz\z=ZQ = 0, the complex potential in the vicinity of this point has the form F(z, t)=°^-[z-

zo(t)}2 + P(t).

(3.132)

Here we restrict consideration to the only practically relevant case of second-order or X-type neutral points. (Note that there is no nonsingular O-type neutral point, since this would require a nonvanishing current density V2ip =/= 0.) If the change of the boundary conditions is such that jS ^ 0, the neutral point z = zo is singular and F(z9t) cannot remain analytic at this point. The natural conclusion seems to be that the electric field dt\p = ft induces a line current / corresponding to adding a singular contribution in (3.132), F(z, t)=°^[z-

zo(t)]2 + P(t) + ^

In [z - zo(f)].

(3.133)

Here I(t) is the total current induced at the neutral point, which is determined by the boundary conditions. It is, however, easy to see that this generates closed \p lines, i.e., an O -point topology, around the neutral point z = zo and two additional X-points, as illustrated in fig. 3.7. The change of topology occurs at these secondary X-points, which is not allowed since it would require dtxp ^ 0 at these points, too, and hence additional singular contributions of the form (3.133) and so forth. Thus it is impossible to construct a solution with isolated singular points. The alternative is to introduce a branch cut corresponding to a current

62

3 Current sheets

(a)

(b)

Fig. 3.7. Generation of a singular current at an X-point. (a) Initial nonsingular configuration; (b) effect of an induced singular current line in the original Xpoint, leading to a fictitious O-point and two adjacent X-points. The heavy line indicates the current sheet actually arising. sheet carrying the induced current I(i). The location of this branch cut is determined by the two fictitious X-points arising by adding the line current to the original X-point, (3.133). As indicated in fig. 3.7(b), the cut passes through this point connecting the new X-points, which drift further apart as I(t) increases. Addition of the current sheet splits the original X-point into two Y-points. Let us discuss the structure of the field in the vicinity of the branch cut. At distances large compared with the length of the cut the potential F(z9t) has approximately the form (3.133). The condition at the cut is that B does not intersect the cut, i.e., xp = const at the cut. To simplify notation, assume zo = 0 and a straight cut extending along the y-axis between the points y = +A. The solution with the asymptotic form (3.133) for \z\ > A and xp = 0 at the cut is given by F(z) =l-az

(3.134)

with the derivative dF

Iz

— By -\- iBx

(aA 2 /2)

=

otz2

(3.135)

While the magnetic potential xp(x,y) = Re{F(z)} is continuous, the magnetic field By has a jump across the cut, the line density of the current carried by the sheet,

= By{0+9y)-By(0-,y)

=.

(3.136)

with

fA J(y)dy = I.

(3.137)

3.3 Driven current-sheet

(a)

63

reconnection

(b)

Fig. 3.8. Contours of xp = Re{F}, where F is given by (3.134); heavy line is the current sheet, dashed line is the separatrix. (a) General case yo < A exhibiting a singularity at the current-sheet endpoints; (b) limiting regular case yo = A (from Syrovatskii, 1971).

The current distribution, (3.136), shows an interesting feature. It is positive, i.e., in the direction of the total current /, in the center part \y\ < yo, and negative in the outer parts \y\ > yo, where y$ = (I/2TKX)+A2/2. The points z+ = (0, ±yo) are neutral points of magnetic field, dF/dz = 0, where the separatrix branches off the y-axis, while the current sheet continues along the y-axis, fig. 3.8(a). At the end points |y| = A the current density J(y) becomes singular giving rise to infinite magnetic field intensity. Only in the special case / = nocA2, where the neutral points coincide with the current sheet endpoints, the singularity vanishes, fig. 3.8(b). This value is in fact an upper limit of the current / (or a lower limit of the sheet width A), since for larger values the points z+ would form isolated neutral points connected by a separatrix similar to the dashed line drawn in fig. 3.7(b), which is topologically not possible because of the frozen-in property. The characteristic features of the Syrovatskii current sheet, (3.136), agree well with those of dynamic current sheets in a fully resistive theory, as discussed in section 3.3.2. The velocity field v corresponding to the magnetic configuration (3.135), (3.127), (3.128), cannot be given in simple analytical form, but must be computed numerically, even in the stationary case dt\p = const. The qualitative behavior close to the current sheet is, however, easily understood. It follows from (3.127) that there is a plasma flux into the sheet, vx(0+9y) = —vx(0-9y) =fc 0. Mass conservation then requires the plasma to flow along the sheet, leaving the sheet at high speed close to the neutral points z±, while the plasma flow vanishes at the singular sheet endpoints. The downstream velocity in the cone formed by the two branches of the separatrix is of the order of the upstream flow and is not related to the cone angle. The technique of conformal mapping used by Syrovatskii is restricted

3 Current sheets

64

(b)

1.5

Fig. 3.9. Construction of a current sheet in an axisymmetric configuration, (a) Initial flux distribution (3.138); (b) modified flux distribution resulting from a change of the parameters ipo and a. Since this increases the flux inside the separatrix, a current must be induced at the X-point to undo this change. This is achieved in (c) by adding a single current ring, while the superposition of 16 current rings in (d) already gives a good approximation to a continuous current sheet (from Longcope & Cowley, 1996). to 2D plane geometry. In axisymmetry for instance, d^ = 0 in cylindrical coordinates r, z, (/>, an X-line is a closed line with coordinates r = ro, z = ZQ. If an electric field E^ acts along this line, a current sheet will be generated, say in the plane z = zo (the orientation of the sheet depends on the boundary conditions), which extends between the closed lines through the points r',zo and r",zo forming a toroidal ribbon. Longcope & Cowley (1996) devised an algorithm to compute the current density in the sheet, approximating the latter by a number n of linear ring currents, which produce a chain of magnetic islands, the width w/ of which shrinks to zero rapidly, w/ oc n~2 as n —• oo. As a particular example the authors consider the purely poloidal field V(/> x Vtp, Wo [4r2z2 - (r2 - 2a2)r (3.138) 2na4 Changing the parameters xpo and a corresponds to a change of the external field configuration, which gives rise to change of xp at the X-line and hence to a current sheet, an example being given in fig. 3.9. Let us briefly discuss the case of a fully three-dimensional magnetic field. Imagine a configuration characterized by two magnetic nulls, one xp =

3.3 Driven current-sheet

reconnection

(a)

65

(b)

Fig. 3.10. (a) Schematic plot of the field, (3.139), including the current ribbon, inside which the current flows in a clockwise direction indicated by the wide gray arrows, (b) The magnetic field lines at the top and the bottom of the current sheet, shown separately (only half of the ribbon is plotted). The current in the sheet arises from the difference in angle between the field lines in the two surfaces (from Longcope & Cowley, 1996).

^4-type and one 5-type (see fig. 2.5). If there is a nonvanishing electric field E along the null-null line AB, a current sheet is expected to form. It appears all we have to do is to split the null-null line thus creating pairs of nulls A\ A" and B\ Bn', which now span a ribbon of finite length, consisting of pieces of S^ and Z# glued together. However, because of the boundaries of the sheet, it cannot carry a net current in the direction of the applied field £, while in reality there will be such a current, which continues to flow in a diffuse way outside the sheet, a feature not allowed in Syrovatskii's theory. Hence it does not seem possible in general to extend the theory to three-dimensional configurations. Longcope & Cowley considered the special case of a closed sheet forming an annular ribbon, generalizing the axisymmetric case, (3.138), by superimposing a weak hom*ogeneous field in the y-direction, (3.139)

= V(/> x

This field has two nulls A, B located at rA,B = (1 + €2/%)a, (\)A$ = +7C/2,

= +ea/2,

where e = na2b/xpo. Now, a current sheet may form in a way similar to the axisymmetric case, leading to a splitting of the nulls into pairs

66

3 Current sheets

of nulls, A —• A',A" and B —> B' ,B". The current sheet consists of two half-ribbons, one extending from n/2 to 3n/2, the other from TC/2 to — n/2 joining smoothly at the lines A'A" and B'B", forming an annular ribbon, see fig. 3.10. The magnetic topology on the top and on the bottom of the sheet is different, corresponding to the properties of A'9 A11 being ,4-type nulls and B\ B" being B-type nulls. Thus the top of the sheet represents part of l^B' a n d equivalently 2^/, the bottom Z#" and equivalently Z^. This field topology modulates the predominantly toroidal sheet current. (In the presence of a sufficiently strong toroidal field B^ the addition of the hom*ogeneous field By no longer gives rise to magnetic nulls and the field topology at the sheet is qualitatively as in the axisymmetric case.) 3.3.2 Dynamic structure of Y-points We have seen in the preceding section that a singular X-point, where the induced electric field dtxp ^ 0, is transformed into a current sheet. Topologically, the separatrix is pulled apart, forming two Y -type neutral points connected by the sheet. Syrovatskii's theory, however, shows that only in a special case are these regular Y -points coinciding with the endpoints of the current sheet. In general, the sheet continues beyond these neutral points with reversed current direction, extending up to the singular endpoints, where the magnetic field becomes infinite. It is interesting to compare these predictions with the structure of the edge regions in a dynamic resistive current sheet. Since the Taylor expansion along the sheet, (3.120) and (3.121), breaks down at the sheet edge y ~ A, the only reliable information is obtained from numerical simulations. Previous analytical treatments of the diffusion region (see, e.g., Vasyliunas, 1975) assume a smooth transition to the ideal external region in the downstream cone, the plasma continuing to flow at the upstream Alfven velocity, to which it has been accelerated along the sheet. Such highly super-Alfvenic flow (the local Alfven velocity in the downstream cone is much smaller than that in the upstream region) should be prone to shock formation, which would increase the field intensity and slow down the plasma to sub-Alfvenic speed. In fact, simulations show that the flow in the downstream region is clearly sub-Alfvenic and apparently unrelated to the high flow speed inside the diffusion region. Closer inspection of the edge of the diffusion region reveals a complicated structure, as is illustrated in figs. 3.11 and 3.12. Figure 3.11 contains stereographic plots of the current distribution viewed from (a) the upstream side and (b) the downstream side, which show the main features of the configuration, the diffusion layer represented by the central current sheet along the y-axis and a weaker sheet current along the separatrix both joining in a region of rather complex behavior. As seen in fig. 3.1 l(b),

3.3 Driven current-sheet reconnection

67

(a)

ib)

Fig. 3.11. Stereographic plots of the current distribution of a stationary reconnection configuration, viewed from (a) the upstream side, and (b) the downstream side (from Biskamp, 1986).

68

3 Current sheets

-0.8

Fig. 3.12. Contours of j , (/>, \p in the edge region of the diffusion layer, (a) r\ = (b) rj = rio/y/2; (c) rj = rjo/2 (from Biskamp, 1986).

the current density in the diffusion layer changes sign, the positive central part is followed by a negative part, which ends in a quasi-singular spike (numerically well resolved, though). Figure 3.12 contains contour plots of j , (j>, xp in the edge region of three stationary simulation states with decreasing values of r\, showing the rapid increase of complexity. The point where the current density of the diffusion layer changes sign (marked by the arrows in the j-plots) coincides with the location where the separatrix (dashed in the t/;-plots) branches off, as in Syrovatskii's current sheet configuration, fig. 3.8(a). The dynamics can readily be interpreted when considering the ^-contours, high streamline density indicating high velocity. The flow, which is accelerated in the central current sheet up to the upstream Alfven speed, is decelerated in the following part of reversed current density and finally completely blocked and turned backwards at the current-sheet endpoint singularity, the spike in fig. 3.1 l(b). The flow is subsequently again accelerated, forming two secondary current sheets adjacent to the primary one, with the same characteristics, j > 0 and j < 0 parts and a flow-blocking shock-like structure. In the limit *;—•() a hierarchy of higher-order current sheets is generated with a self-similar scaling behavior, as drawn schematically in fig. 3.13.

3.3 Driven current-sheet reconnection

69

Fig. 3.13. Schematic drawing of the self-similar hierarchy of current sheets in the diffusion layer: thick line = positive j ; thin line = negative j ; open circle = quasi-singularity at sheet endpoint. These properties of the diffusion layer are consistent with Syrovatskii's current-sheet model, the multi-current-sheet edge behavior representing the dynamically resolved singularity predicted in Syrovatskii's quasi-static theory. The picture indicated in fig. 3.13, however, rests on two assumptions, viz., symmetry and stationarity. As will be discussed in section 4.7, a dynamic current sheet, though more stable than a static one, becomes unstable if the aspect ratio A = A/d is sufficiently large, which gives rise to nonstationary behavior. In a less symmetric configuration such non-stationary behavior will be very complex and finally lead to fully developed turbulence consisting of a statistical distribution of microcurrent sheets, see section 5.4.

3.3.3 Scaling laws of stationary current-sheet reconnection In Syrovatskii's theory, discussed in section 3.3.1, the system is completely determined by the boundary conditions on xp. The corresponding plasma flows are passive, i.e., they do not affect the magnetic configuration, and the physical reconnection mechanism, in particular the value of the resistivity, does not enter the theory. By contrast, numerical solutions of the resistive MHD equations exhibit a strong influence of the magnitude of r\ on the entire system. Hence, to determine the dependence of reconnection on r\ the dynamic equations must be solved. More specifically we consider an open system, where reconnection is driven by imposing suitable inflow and outflow boundary conditions. Unfortunately the problem appears to be too complicated to permit an analytical treatment without severe approximations. In particular, the usual boundary layer approach of matching the solution of the resistive region to that of the ideal external region is practically impossible, since the complicated shape of the resistive layer (section 3.3.2) makes the problem non-separable, i.e., fully two-dimensional.

70

3 Current sheets

Fig. 3.14. Computational system as used in numerical simulations of driven reconnection. Hence the problem has to be treated numerically. The simplest and also most reliable method is that of dynamical simulation, i.e., to follow the system from some initial state with stationary boundary conditions until a steady configuration is reached, which automatically eliminates unstable solutions. Accurate numerical solutions in the relevant parameter regime are state of the art, and since there are only two essential parameters their scaling laws can be determined from a limited number of computer runs. (In practice, we do not have to restart the system from the initial state for each parameter value. Instead, the run is simply continued, changing, for instance, the resistivity in not too large steps, which requires only a short time for the system to relax to the new steady state, provided such state exists.) Since primary interest is in understanding the qualitative behavior, we can choose the simplest possible geometry with up-down and right-left symmetry, such that only one quadrant must actually be computed. The computational system is indicated in fig. 3.14. The basic equations to be solved are the 2D incompressible MHD equations (3.65) and (3.66). Boundary conditions have to be assigned to ip,,j and w. While at the internal boundaries, the x-axis and the y-axis, these conditions follow from the imposed symmetry, \p and j being symmetric, (j) and w antisymmetric, at the upper (x = Lx) and the left-hand (y = Ly) boundaries they correspond to the inflow and outflow conditions, respectively. These should be chosen in a way conforming with the concept of an open system, such that

3.3 Driven current-sheet reconnection

71

boundaries, in particular the outflow boundary, do not obstruct the flow. Open boundaries are well defined for linear waves, requiring that waves are not reflected at the boundary, or, in mathematical terms, that for a hyperbolic system of differential equations all characteristics are outgoing at the boundary, which guarantees that perturbations arising due to the presence of the boundary do not propagate into the system. However, (3.65) and (3.66) of incompressible MHD are of mixed hyperbolic-elliptic type. Moreover, the equations are nonlinear and two-dimensional, and to date there is no rigorous method to determine whether a particular set of boundary conditions is admissible. (A discussion of the boundary conditions for different MHD systems has been given by Forbes & Priest, 1987.) Let us therefore apply a more practical procedure. While it seems to be sufficient to require continuity of w and j at the boundaries in (3.65), dnw = dnj = 0, the system is more sensitive to the boundary conditions for the potentials \p in (3.66) and 1, which means that inertia effects now give rise to a current sheet of length A, which modifies the configuration on a scale A >• 5. The properties of this current sheet Bo,uo, (5, A indicated in fig. 3.14, are determined by the asymptotic conditions on the configuration, in particular the imposed reconnection rate £, and by the resistivity. We can rewrite (3.107) and (3.108) such as to express Bo,uo,S in terms of E,r\ and the sheet length A, Bo ~ {E2&/n)ll\

(3.148)

uo - (£>//A) 1/3 ,

(3.149)

5 ~ (r)2A/E)1/3. (3.150) The length A, which depends on the specific dynamics at the sheet, cannot be determined in a simple way. A series of simulation runs has been performed (Biskamp, 1986) for identical boundary profile functions varying only the parameters E(= Woo) and r\, which yield the scaling relation A~E4/rj2. (3.151) Inserted into relations (3.148)—(3.150) we obtain the current-sheet scaling laws in a stationary driven reconnection process in terms of the two free parameters E, r\: Bo ~ E2/n, (3.152) (3.153) S~E.

(3.154) n

Hence, increasing the driving E increases the field J?o i front of the sheet (flux pile-up) and correspondingly lowers the velocity wo because of UQBO = E. The most conspicuous feature is the increase of the sheet length A, when r\ is decreased. This is illustrated in fig. 3.16, which gives the flow patterns (j)(x,y) and the magnetic configurations tp(x,y) for three cases differing only in the value off/. The figure clearly shows the increase of A, the flux pile-up, and the slowing down of the inflow velocity. The scaling laws (3.151)—(3.154) can be related to the physics in the current sheet, which is dominated by inertia, i.e., the acceleration along

74

3 Current sheets

Fig. 3.16. Contours of 0 and xp for three stationary configurations of currentsheet reconnection. (a) r\ = rjo', (b) rj = rjo/2; (c) Y\ = rjo/4 (from Biskamp, 1986). the sheet. As shown in section 3.2.3 the velocity vy increases linearly, vy ~ Boy /A, such that the average force along the sheet is ~ COUSt

(3.155)

using (3.151) and (3.152). Hence the scaling relations imply that the force along the sheet is independent of E and r\, in particular remains finite for In Syrovatskii's theory, the current-sheet length A is independent of the reconnection dynamics inside the sheet, in contrast to the result, (3.151), obtained for a fully dynamic system, which is due to the fact that in Syrovatskii's approach the reaction of the flow dynamics on the magnetic field is ignored, the flow being passive, determined only a posteriori. The question of a finite reconnection rate in the limit r\ —> 0 has also been considered by Grad et a\. (1975). Instead of the full MHD equations these authors treat a system of equations suitable for slow diffusive processes, in which plasma inertia is omitted, which reduces the dynamics to a sequence of MHD equilibria changing quasi-statically because of changes in the boundary conditions and magnetic diffusion. By solving these equations numerically it is found that reconnection always occurs at a finite-angle X-point and that the reconnection rate seems to remain finite for Y\ —> 0. No indication of current-sheet formation is observed. The crucial objection, however, concerns the omission of the inertia term, which is only justified if fluid velocities remain much smaller than the Alfven velocity. As we have seen, this is not true for the flow along the current sheet, which limits the validity of Grad's reduced equations to

3.3 Driven current-sheet reconnection

75

reconnection rates E = 0.* Solving these equations, subject to boundary conditions and the jump conditions at the shock, determines the configuration, an example being given in fig. 3.19(a). Petschek derives an approximate analytical solution for a oo. 4.1.1 Linear tearing instability Writing ipi(x,y) = xpi(x)&ky, (4.8) and (4.9) in the resistive layer become ixkBr^l+nxpl = ixkB'oxp'{-ikfoxpu

(4.11) (4.12)

where we make the approximations V 2 0i ~ (///, ¥2xpi ~ 1///, Byo(x) ~ xB'o, the constant-t/; assumption xp\ ^ Vi(0)> and the coordinate system such that xs = 0. The individual terms in each equation are a priori of the same order. The f0 term in (4.12) can, however, be neglected, since it would only lead to a small correction of the stability threshold replacing A' by

86

4 Resistive instabilities

A! + O(r]2/5) (Bertin, 1982). The functions 4>i(x),ip\(x) can then be chosen with definite parities, 4>\ odd and \p\ even (or the reverse, but since we are only interested in the reconnecting mode with ipi(O) ^ 0, the former parity is the relevant one). This parity choice makes the following simple analysis semi-quantitative. With the approximations $[ ~ —cj)\/d2,\pf[ ~ A'xpi/S, (4.12) and the two relations obtained from (4.11) by equating the l.h.s. with either the first or the second term on the r.h.s. give y and 3 as functions of A', y ~ f/3/5(A')4/5(kBf0)2/5, (4.13) S ~ rjA'/y ~ r}2/5(Af)1/5(kBf0)-2/5.

(4.14)

(From a rigorous treatment one obtains a numerical factor of 0.55 in the expression for y9 see section 4.2.) In the analysis we have tacitly assumed that A! > 0. It can in fact be shown (Furth et a/., 1963) that, in the framework of incompressible MHD, the tearing mode is unstable if, and only if, Af > 0. Relations (4.13), (4.14) are consistent with the constant-t/; assumption, which requires that the skin time of the layer T$ ~ S2/f] is shorter than the growth time y" 1 . We find from (4.14)

rs ~ nh2 = o(T 1 / 5 ) < y"1 = 0(>T3/5).

(4.15)

The quantity A' is obtained from the solution in the external region. Since the instability is weak, y ~ x^ T:~A < t j 1 , where xn = a2/r\ is the global resistive diffusion time, the inertia term in (4.9) is negligible, such that the equation reduces to the perturbed equilibrium equation B o - V j i + B i - V ; o = O.

(4.16)

In the vicinity of the resonant surface, (4.16) assumes the form t// 1 '-^ 1 =0, where K = fo/Bfo.

The solution has the form y i ( x ) M ( 0 ) = /i(x) + Cfi(x),

with

(4.17)

CK*2),

f2 = x + O(x2\

(4.18) (4.19)

where C may have different values C+ on both sides of the singularity, which are determined by the outer boundary conditions. To obtain an overall smooth behavior the resistive layer solution must be matched to the external solution by identifying the logarithmic derivative of \p\ of the former for |x| -> oo, A\ with that of the latter for x —• 0, A' = c+ - C_.

(4.20)

87

4.1 The resistive tearing mode

(a)

(b)

Fig. 4.4. Schematic plots of the eigenfuction xpi for (a) fo(xs) = 0 as in the symmetric sheet pinch; (b) fo(xs) ^ 0 as in the cylindrical pinch. Note that the logarithmic singularity of xp[, which is proportional to f0, cancels in A'. In the case of a plane sheet pinch, (4.13) becomes (4.21) For the standard symmetric sheet pinch BQ(X) = Botanh(x/a),jo = (Bo/a)sech2(x/a), and boundary conditions \p\ = 0 for |x| —• oo, (4.21) can be solved analytically (see fig. 4.4(a)) (4.22)

ka hence A'= 2

(4.23)

— ka

a ka We thus find that the tearing mode is unstable, A' > 0, for long wavelength ka < 1 and stable, Ar < 0, for short wavelength ka > 1. These properties remain qualitatively valid for rather general jo(x) profiles, in particular the behavior A' oc k~l for small k. The constant-y; property breaks down for k ~ 3/a2, where A' ~ 3~{. A more general analysis without the constant-tp assumption shows that the growth rate, which increases with decreasing /c, y ~ (A')4/5fc2/5 ~ /c~2//5, reaches a maximum .,

.,//, \ _,

/ m a x — yy^O)

M l/2

i\

^ _

9 K>0

'I

•

OTl/4

M 9/i\ V^"*^ v

For still smaller /c, y decreases again, as can be seen also from (4.13) and (4.14) by inserting Ar ~ 3~x. In the general nonsymmetric case, where JQ(XS) ^ 0, for instance in the presence of a mean field B, BQ(X) = B + BQ tanh (x/a), the eigenfunction

88

4 Resistive instabilities

exhibits a logarithmic singularity in the derivative, (4.19), see fig. 4.4(b). The growth rate is, however, not affected by this singularity. In a cylindrical screw pinch, (3.16), with periodicity length Lz = 2nR the resonant surface of a helical mode with poloidal and axial (or toroidal) mode numbers (m, n) is given by the relation k B = -B0(r) - ^Bz = -Bo (l - -q(r)) = 0, * (4.25) r K r \ m J where q = rBz/RBe is the safety factor, the inverse of the field-line twist T or the rotational transform i. The reduced equations can be written in terms of the helical flux function \p* = (n/m)(r/R)Ao — AZ9 (3.37),*

since, in the framework of reduced MHD, Bz = const, hence AQ = rBz/2 and Az = —xp. Equations (4.6), (4.7) thus become (4-27)

= nh (dt + v V)w = B* • V/,

(4.28)

where B* = ez x Vtp* is called the helical magnetic field. The equation for the external solution (4.16) becomes

dV dr2

lftyi _ /^m2 ! r dr \r2

7o(r) \ Be(r)[l-nq(r)/m]J

= 0

The short-wavelength stability of the tearing mode, A' < 0 for ka >, 1, has important consequences in a cylindrical plasma. For equilibrium profiles with gradient scale a, tearing instability arises only for ma/r | ^ |x|vo = y ^ ^ v - t p i c o s y ,

(4.38)

and dl ^ dy, so that, for instance, cos ydl

^4 f

Jo

™y*y Jxp -xpicosy

{439)

with u>1

Inserting j(xp) from (4.34) into expression (4.36) and using (4.35) gives 4

(4.40)

91

4.1 The resistive tearing mode

\ \ \\«

ii)

Fig. 4.5. Flow pattern around the X-point in the tearing mode evolution. where

Ju — cos y

= 1.827.

(4.41)

Hence we obtain (4.42)

Rutherford's theory is exact for small island size. This does, however, not imply that inertia effects are completely negligible, since on the separatrix j need not be a flux function, B vanishing at the X-point. But the detailed behavior in this region does not influence the global evolution, which is determined by (4.34). From (4.42) one finds the reconnection rate

E=

for

t 1. This behavior is related to the vacuum bubbles in a shearless plasma column (Kadomtsev & Pogutse, 1974). The stability of the tearing mode depends rather sensitively on the current-density gradient in the vicinity of the resonant surface. A local steepening may destabilize the mode, while a flattening leads to stabilization. The latter effect is indeed the dominant self-stabilization process of an unstable tearing mode when growing to finite amplitude. The process may, however, destabilize other modes owing to the current-profile steepening in the region adjacent to the island, which appears to play an important role in the major disruption in tokamak plasmas (see, for instance, Biskamp, 1993a, chapter 8). A particular shape of the ^-profile close to q = 1 can give rise to cascades of high-rc tearing modes with n up to 20 (Hallatschek et al, 1998). By current-profile tailoring around the major resonant surfaces completely tearing mode-stable profiles can be constructed (Glasser et al, 1977). 4.1.4 Effect of dynamic resistivity In the linear theory of the tearing mode, Y\ is usually assumed hom*ogeneous and constant in time (resistivity gradient effects give rise to a different type of instability, the rippling mode, see Furth et al., 1963). A hom*ogeneous resistivity distribution violates, strictly speaking, the equilibrium condition rjj = const, but in linear theory this effect is negligible, since the instability time-scale is faster than the resistive diffusion time of the equilibrium. In the nonlinear regime, however, the behavior of the resistivity becomes important, in particular for large island width. Let us therefore assume that the resistivity adjusts dynamically such that f]{x,y) = f](Te) = f](xp) owing to fast parallel electron heat transport, where the xp -dependence is determined by slow perpendicular heat conduction and local heat sources and sinks. Consider the difference of (4.34) taken at the X-point and the O-point, where V\p = 0, so that the advection term vanishes: —(xpx-xp0) = rjxjx - rjojo-

(4.51)

96

4 Resistive instabilities

Since \px-xp0 = 2xpu jx - jo - 2xpi A'/wj and r\x-r\oobtains approximately w/ ~ 77 A' + j-r-wi.

2xpidr\/dxp, one (4.52)

The drj/dxp-term is related to the a-term in (4.44). We first discuss the tearing mode in a sheet pinch. The physical meaning of (4.51) can easily be understood. Because of the choice \p$ > 0 the magnetic flux in the island is positive, xpx — xpo > 0. The variation of this flux is given by the difference between the injection rate fjxjx due to reconnection at the X-point and the dissipation rate rjojo due to resistive dissipation in the island. Lowering r\o with respect to rjx, drj/dxp > 0 ("island heating"), leads to enhanced growth and larger saturation width, while increasing rjo, drj/dxp < 0 ("island cooling"), slows down the growth and reduces the final width. In a cylindrical tokamak-like configuration, the effect of a dynamic resistivity on helical tearing modes is just the opposite of that in a sheet pinch. Since xp is now replaced by the helical flux function xp* = xp - (n/m)r2Bz/2R, (4.26), xp'm(r) = Be(r)[l - nq{r)/m], we have xp'l{rs) < 0 for the usual case of radially increasing q(r) and jo > 0. This implies that xp* has its maximum at the 0-point, in particular xp*x—xp*o < 0, such that (d/dt)(xp*x — xp*o) > 0 corresponds to island shrinking. Hence it follows from (4.52) that island heating r\o < r\x, meaning drj/dxp* < 0, leads to smaller islands, while island cooling leads to larger ones. Numerical simulations illustrating this effect are shown in fig. 4.7. The initial current profile is the rounded profile (v = 2) with ro = 0.5, rs = 0.65, and the resistivity rj(xp*) is defined such that the profile along the radius through the X-point remains constant, rj(r,9x) = rjo/jo{r,t = 0). The resistivity in the island is still free. For t < 10 3 T^ we assume a flat distribution rji(xp*) = r\x, the value at the X-point. For t > 103TA the island resistivity is shaped, f]i(xp*) = rix[l + e(xp* - xp*x)}•

Choosing e > 0, i.e., r\o > r\x corresponding to island cooling, gives rise to larger islands (the upper branch in fig. 4.7(a)), whereas e < 0, i.e., r\o < ^x corresponding to island heating, leads to smaller islands (the lower, dashed branch). The corresponding */-profiles are plotted in fig. 4.7(b). One can see that a relatively small change of r\ (~ 10 percent) has a dramatic effect on the island width. 4.1.5 Neoclassical tearing mode The main influence of the island resistivity profile rji(xp) on tearing mode evolution occurs through its effect on the current profile in the island. Since

4.1 The resistive tearing mode

97

(a) Fig. 4.7. Effect of dynamic resistivity, (a) Island width W21 for three island resistivity profiles; (b) radial 77-profiles across the O-point (full line) and the X-point (dashed line) for the final states in (a).

in equilibrium r]I(xp)jI(xp) = const, a decrease of r\\ implies an increase of jj and vice versa. Hence if ji is reduced compared to the ambient value, the tearing mode is destabilized leading to larger saturated islands. In hot tokamak plasmas the parallel (~ toroidal) current density is increased compared with the resistive value £||/*/|| by the bootstrap current, which is proportional to the pressure gradient (see, e.g., Kikuchi & Azumi, 1995),

J\\-

c dp

:

Ye~dr'

(4.53)

The last term arises from the friction between electrons circulating around the torus and those trapped on the low-field outer side. The rigorous calculation of these collision processes is the subject of neoclassical transport theory (see, e.g., Hinton & Hazeltine, 1976). Since in the island the pressure profile p ~ p(\p) is essentially flat, there is no bootstrap-current contribution from the island. Hence the current density in the island is reduced compared to the value outside the island, which has a destabilizing effect on the tearing mode. Extending Rutherford's theory, (4.42), to

98

4 Resistive instabilities

Fig. 4.8. Schematic drawing of the phase-space diagram of the neoclassical tearing mode. include the neoclassical effect gives (Carrera et a/., 1986),

* / = 1.22, (A'+ ^ ) , where Lq = (dlnq/dr)'1, Lp = —(dlnp/dr)'1, fipe = %npe/Bl is the electron poloidal /J, and a is a numerical factor of order unity. (An expression similar to (4.54), replacing ^Je —> icLq9 is obtained for the nonlinear tearing mode under the influence of magnetic field-line curvature JC, see Kotschenreuther et al., 1985). It follows from (4.54) that for small island width the growth of the tearing mode is dominated by the neoclassical term, w/ ~ t1^2. The largeamplitude behavior and saturation depend on the sign of A'. For A' > 0 these are determined by the nonlinear modification of A', A'(wj) —> 0, as discussed before. For Af < 0, however, (4.54) indicates saturation at w/ = A/\A'\9 which seems to imply that the neoclassical tearing mode grows to finite amplitude for any negative A'. Under real conditions, however, the form of the neoclassical term in (4.54) is only correct for not too small island width, since the flattening of the pressure profile in the island is counteracted by cross-field diffusion. One can introduce this effect phenomenologically by the following substitution in (4.54) (Fitzpatrick, 1995):

with Wd oc (Ka/Kj)1/4, where K±, K\\ are the perpendicular and parallel heat diffusivities. Hence there is, in general, a critical island width W/,Crit =

\A!\\V2d/A,

which must be exceeded, before the mode starts growing. (A detailed discussion of the threshold conditions is given by Wilson et a/., 1996.) Growth may be triggered by toroidal coupling to modes of different

4.2 The double tearing mode

99

helicities, which are excited in a tokamak discharge, for instance by the m = 1 mode associated with the sawtooth oscillation (see section 6.5). The phase-space diagram of the neoclassical tearing mode is shown infig.4.8. 4.2 The double tearing mode The basic assumption in the theory of the tearing instability, given in the preceding section, is the constant-tp approximation, the validity of which requires that the skin time of the resistive layer is shorter than the growth time, T$ < y~l. There are, however, circ*mstances when this condition is not satisfied. Formally this occurs when A' —> d~x9 which we have already encountered in the long-wavelength range k ~ t]1^4 at the maximum tearing mode growth rate, where T £, k}, —> /c, Boy —> l?o- In the external region, resistivity and inertia can again be neglected. Hence from (4.56) one has xpi = — #o£, which on substitution in

W / V

V

(4.58)

gives — (BQ-^-)

=k2B^.

(4.59)

100

4 Resistive instabilities

(b)

I/A

I

\

)

V

(0 Fig. 4.9. Schematic plots of the equilibrium field Bo and the eigenfunction £ for (a) the standard tearing mode, (b) the double tearing mode, and (c) the resistive kink mode.

To be definite, consider a symmetric equilibrium field Bo(x) with just two resonant surfaces B0(±xs) = 0 (Pritchett et al, 1980). For long wavelength kxs < 1 the factor on the r.h.s. is small, such that we can solve (4.59) in terms of a power series expansion in the parameter k2x2, £ = £(°)+£(1)+---. To lowest order, d^/dx = 0 except at the resonant surfaces x = +xs. Since for an isolated system £ must vanish for |x| —> oo, we have £o = const, \x\ < xS9 0, |x| > xs.

(4.60)

While in the case of a single resonant surface, boundary conditions enforce fo = 0, such that £ is localized at the resonant surface, fig. 4.9(a), in the case of two neighboring resonant surfaces, ^ may be finite in between,

4.2 The double tearing mode

101

fig. 4.9(b). The first-order term is easily obtained: ( k2 fx 2 -^ / Bkdx, dx

Xs

Bldx

\x\ < x s , \x\ > xs.

In the vicinity of x = xS9 one has BQ(X) ~ (x — xs)Bf0, such that (4.61) can be written in the form 1

^

dx

kH

n{xxr n{x-xsr

[

}

where r B2Qdx (4.63) Bo Jo is a measure of the free magnetic energy of the mode. In the resistive layer (4.11) and (4.12) have to be solved which, when written in terms of £, read 2 * (4.64) AH = - T T - ^ 2

yxpi = -(x - xs)yBfQZ + ij\p'{.

(4.65)

Asymptotic matching of the logarithmic derivatives of the solution to the external solution, (4.62), yields the growth rate. The problem can be treated analytically (Coppi et a/., 1976; Ara et al, 1978), which gives the following expression for the normalized growth rate y: (466)

Here y = y/{nk2B^ = y(x^A)1^ and 1H = ^ 1 1 f knk{xnjXA) ^, IA~ = B 0, xn~l = rjk2. From the general expression (4.66), y can be obtained explicitly in three special cases: (a) For strong ideal instability XE > 1 use of the asymptotic properties of the F-function, F(z + a)/T(z + b) ~ za~b> gives the ideal growth rate y = 2H, i.e., y = *HkBf0; (4.67) (b) Marginal instability Xu — 0 corresponds to the first pole l)/4], which is y = 1 or y = (t1kIB'oy/i;

(4.68)

102

4 Resistive instabilities

(c) For deeply ideally stable MHD conditions JH < 0, \XH\ > 1> corresponding to small growth rate y < 1, (4.66) gives (4-69) Hence in the MHD-stable regime the mode reduces, as expected, to the standard tearing mode. In fact XH is directly related to A'. From (4.62) one obtains (

S

) ;

n(x - xs)

(470)

x-x s -+0 +

or, using the ideal relation \pi = —(x — xs)B'o^,

x -xs + (XH/n)],

x- x

whence A' = xp[/xpi\o+ -tpi/tyilo. = -K/AH

(4.72)

(This analysis is only valid for A' > 0, since for A' < 0 the displacement £ differs from expression (4.70), see below.) Evaluating the coefficient in (4.69) the growth rate thus becomes y = 0.55rj3/5A'4/5(kB'0)2/5.

(4.73)

Equation (4.63) indicates that 1# is always negative for the double tearing mode. (Because of the assumed symmetry it is sufficient to consider only one resonant surface x = xs.) In fact, in a plane configuration all modes are ideally stable in the framework of reduced MHD. For long wavelength kxs 9o, where g(#o) = fu By contrast, the last non-reconnected flux surface has the value fe and extends all the way around the axis, since fe > gmax- For — 6Q < 9 < 9Q, where the surfaces labeled x = ±xs touch, the flux is continuous since xp*(—x) = tp*(x), but the derivative dr\p*9 the helical field 5*, exhibits a symmetric ^-dependent discontinuity from — ^Jfe — g(6) on the inside (x = —xs) to +sjfe — g(6) on the outside (x = xs). Note that the current singularity is not limited to the strip between the Y -points, —6o < 9 < 9o, the current sheet proper, but extends around the island along the separatrix, since y/fe — g(6) — yjfi — g(0) 7^ 0, though its value is smaller. This property is not captured by Syrovatskii's theory, but clearly appears in the driven reconnection simulations, fig. 3.11. On the other hand, the present theory misses Syrovatskii's singular behavior at the sheet edge, (3.136). Equation (4.94), resulting from the equilibrium condition, and the conservation of the helical flux, (4.95), have been solved numerically for g(6) and /(x), which determines the configuration in terms of the core displacement £o (Zakharov et a/., 1993). One finds in particular 90 = 60°, w7 = 1.84ft, *s = 0.4ft.

(4.96)

(xs = 0.5ft would correspond to 9Q = 0, i.e., osculating circular surfaces.)

110

4 Resistive instabilities

Though the instantaneous configuration can be calculated without considering the physics of the reconnection process, the time evolution, of course, depends thereon. Since reconnection occurs in a resistive current sheet, we adjust the Sweet-Parker model discussed in section 3.2.2 to the geometry of the kink mode. In the laboratory frame the flow into the sheet is not symmetric, but enters only from the core side. Since reconnection of the helical field is symmetric, involving equal amounts of internal and external helical flux, the sheet itself is moving outward radially at half the plasma speed u ~ to in the laboratory frame, such that in the sheet frame the inflow velocity is uo = u/2. Since the plasma entering all along the sheet leaves at the sheet edge, the continuity equation gives the relation

I

uOnd6 = uorx sin 0O = t>o.\ = lK|x s ^ lvolft/2.

(4.98)

(Note that in our normalization B* = VA 1, Lp = po/pfo = pressure scale-height, one has Bi • Vjo < Bo • V/i, hence the kink-mode term Bi • V/o has been omitted in (4.115). In the local approximation one neglects the variation of the coefficients, which corresponds to the Boussinesq approximation in fluid dynamics (see, e.g., Drazin & Reid, 1981). The perturbation can be Fourier transformed, V/ —• ikf, such that (4.115)—(4.117) become algebraic. Choosing K = 7cex, Vp0 = Po^x, b = ez, and neglecting r\ for the time being, we obtain the dispersion relation ! where v\ = B^/po andfcj_,fc||refer to the direction of the magnetic field. For K pointing along the pressure gradient, Kp^ > 0, convective instability arises, if the line-bending stabilizing term kh)\ is sufficiently small, icjff'>fcf

(4.119)

with P' = 2pf0/BQ. (If KPQ > 0, the curvature is called unfavorable or bad, since the resulting convective instability is obviously detrimental to confinement in a plasma device. The opposite case Kpf0 < 0 is called favorable curvature.) In real magnetic configurations, however, k\\ cannot, in general, become arbitrarily small. A lower limit of k\\ arises because of two different effects: (a) the extent along a field line, where K is destabilizing, is finite; (b) neighboring field lines are not strictly parallel, since there is, in general, magnetic shear. For a mode localized in the region of destabilizing curvature we have fe|| ~ L" 1 , where LK is the extent of this region along the field line. In a tokamak plasma, field-line curvature is mainly due to toroidal curvature, hence K ~ — £R/R, where e^ is in the direction of the major torus radius R. A field line winding around the plasma column with a twist T = q~l thus passes alternately through regions of unfavorable curvature on the torus outside and favorable curvature on the torus inside, the distance * We use the same normalization as before, but since in a stratified medium the mean density may depend on position, p 0 is written explicitly.

120

4 Resistive instabilities

between both being LK = nRq. A perturbation will grow, if it is confined to the outside. Hence, from (4.119), a necessary condition for this so-called ballooning instability is, qualitatively,

-f

> j±,.

(4.120)

Modes not localized in regions of destabilizing curvature feel only the average curvature of the field line K. These modes are called interchange modes, since they correspond simply to an interchange of the positions of neighboring flux tubes together with their plasma content. Interchange modes may be unstable in a cylindrical pinch, where the curvature is constant along the field lines, K = Kcer, KC ~ —Bg/rB2 = —r/(Rq)2. In a tokamak plasma K consists of two parts, the average toroidal curvature, which is stabilizing, since the favorable contribution on the torus inside is slightly larger than the unfavorable on the outside, K1 ~ r/R2, and the cylindrical curvature, KC. For a tokamak of small inverse-aspectratio e = r/R < 1 and circular cross-section one obtains (Shafranov & Yurchenko, 1968) K = KC{\ - q2) = - ^ j ( 1 " q2).

(4.121)

Hence interchange modes are not excited in a tokamak plasma if q2 > 1, which is the case for the major part of the plasma cross-section. In a reversed-field pinch, however, where the toroidal field is much weaker and hence the field-line twist much larger than in a tokamak, q2 < 1, the cylindrical curvature dominates, such that interchange modes may be unstable. The second effect limiting k\\ from below, the magnetic shear, is rather universally present. If the perturbation is aligned with a particular field line k • B = 0, it is in general no longer so on the neighboring field lines. Let us first discuss this effect qualitatively for a cylindrical plasma. In the vicinity x = r — rs of the resonant radius rs we have -q(r)) ~ -™Be(rsJsx, (4.122) m ) rj where ? = dlnq/dlnr is the shear parameter or simply the shear, see (3.18). Since the radial mode width is typically of the order of the poloidal wavelength, (m/rs)x ~ 1, we find from (4.122) fc|| - We/Bzrs = 3/Rq,

(4.123)

which on substitution in (4.119) gives the condition for interchange instability in a cylindrical plasma KCP' > a - ^ ,

i.e., -pf > a y ,

(4.124)

4.5 Pressure-driven instabilities

121

where a is a numerical factor of order unity. If ? = 0(1), as in a tokamak, we see that interchange modes require much higher /?, /? ~ 1, than ballooning modes, (4.120). For a rigorous treatment of localized modes in a general equilibrium configuration it is convenient to use flux coordinates chosen in such a way that one of the two coordinates spanning a flux surface is in the direction of the magnetic field. The problem then reduces to the solution of an ordinary differential equation along the field line. Consider a flux tube of small but finite diameter encompassing the field line, about which the mode is localized. In a local Cartesian coordinate system the sheared magnetic field has the form Bo = Bo[ez + (x/L 5 )e y ],

(4.125)

where Ls is the shear length, L~{ = H/Rq in a cylindrical plasma. We now transform to a twisted coordinate system (Beer et a/., 1995): £ = *, r , = y + (x/Ls)z9

f=z,

(4.126)

such that £ is the coordinate along the field line, which gives B V = Bodc,

Vi = d2x + d) = [de + (£/L s )^] 2 + 4

(4.127)

Neglecting r\, we can substitute \p\ and p\ from (4.116) and (4.117), 2

i

= Bo • VVi(B 0 • V0) + 2(b x K) • V[(b x Vp0) • V #

(4.128)

In the twisted coordinates, (4.126), the coefficients depend only on the parallel coordinate £, the dependence on x in the operator Bo * V being eliminated. Hence we can Fourier transform in the perpendicular directions d% —• ik%9 dn —> ikn, such that (4.128) becomes an ordinary differential equation. k^kn are the covariant components of kj_,

fe,]ex + V , , Vi = -k\ = -[fc{ + (C/Ls)/c,]2 - k\ Introducing normal and geodesic curvature components,

the operator in the last term in (4.128) becomes (b x K) - V ± = i(b x K) • k ± = i{Knkn - Kg[k{ + (ULs)k,]}. Hence (4.128) assumes the form,

(4.129) (4.130)

122

4 Resistive instabilities

where we have used the notation 2

k\

(C - Co)2

F = p = l+ —-j—,

._ h

Co--rrLs.

Equation (4.131) is called the ballooning equation. The curvature components Kn, Kg and the shear length Ls are functions of £. Equation (4.131) must, in general, be solved numerically, but some special cases can be discussed analytically. The condition for instability of interchange modes, Mercier's criterion (Mercier, 1960), is obtained by considering the asymptotic solution for £ —• oo (see, e.g., Freidberg, 1987). In the cylindrical case this reduces to Suydam's instability criterion (Suydam, 1958),

(4132)

iw-

which specifies the numerical factor in the qualitative relation (4.124), or (4.133)

P'KCL2S=DR>X-.

In the case of a large-aspect-ratio tokamak of circular cross-section, Mercier's criterion becomes (Shafranov & Yurchenko, 1968) K/r^l-gV^^r

(4-134)

Ejfect of finite resistivity Let us now include finite resistivity, which permits the magnetic field to slip across the plasma and thus reduce the stabilizing line-bending effect. The dispersion relation in the local approximation, (4.118), now reads

(4135)

^-Tinh-po 1 + r\kyy

1 + isr\kyy hence for rjk]_/y > 1 the stabilizing term reduced. We now show that in an inhom*ogeneous system, where k\\ arises because of magnetic shear, the instability threshold is, in fact, independent of the shear effect, only the growth rate depends on rj. From (4.116) we see, by use of (4.130), that the resistive term does not introduce additional differential operators, ^

(4-136)

such that the ballooning equation (4.131) becomes > = 0.

(4.137)

4.5 Pressure-driven instabilities

123

Equation (4.137) can be solved in a way methodically similar to that used in tearing mode theory, section 4.1.1. For rjk]_/y < 1 the problem exhibits two disparate scale-lengths along the field line, the equilibrium scale L s , and the resistive scale defined by rjk]_F2/y ~ 1, Ln = Ls(y/r]k])1/2 > Ls. In analogy to the tearing mode boundary layer problem, one solves (4.137) separately on both scales using appropriate approximations. Since the general theory, which also includes finite compressibility (Correa-Restrepo, 1982), would lead beyond the scope of this presentation, we only outline the method and sketch the results. The theory applies to toroidal configurations, where field lines wind around the plasma indefinitely, making the coefficients in the ballooning equation (4.137) periodic functions of C, apart from the explicit secular £ dependence. In the resistive range one formally introduces two scales £ and z = e~1^ en = (t]k\/y)1/2, expanding the solution in the form 0 ( 0 = 0 o (z) + ^ i ( C , z ) + • • •,

(4.138)

where the dependence of the /s on £ is assumes to be periodic. (j)o(z) is the envelope of the eigenfunction. Consider first a cylindrical plasma column, where Kg = 0 and Kn = KC = const. In this case we obtain the equation for the envelope o(z) directly from (4.137), since there is no explicit dependence on the equilibrium scale £. Let us write (4.137) in the following non-dimensional form (4.139)

dz 1 + (y\ly)zl dz where

z = f , y = y-, n = i/fci^, DR = pKcLi Ls

VA

(4.140)

VA

and Ls = Rq/!$ in cylindrical geometry. Moreover, we restrict ourselves to conditions, where y is sufficiently small such that (r\/y)z2 > 1. In this case (4.139) reduces to the Weber equation ^ ° + ( A - u 2 ) ( / > o = O,

(4.141)

with u = z(Jzy)1/4 and A=

£1/2

n

wDR.

The harmonic oscillator solutions lead to the eigenvalues A = \ n = 1 + 2n, n = 0,1,....

(4.142)

124

4 Resistive instabilities

The largest growth rate corresponds to n = 0, where 0o = exp(—u2/2), y = >71/3Dj/3 or

(4.143)

2/3

y= Since the eigenvalues Aw are positive, it follows from (4.142) that the condition for instability is DR > 0. (4.144) Compared with Suydam's criterion, (4.132) or (4.133), the stabilizing linebending effect induced by the magnetic shear is eliminated due to reconnection. Equation (4.141) is only valid for sufficiently small y, rjF2/y > 1, or DR along the field line and the growth rate, (4.143). For larger DR the unity in the denominator of F2/(l + rjF2/y) cannot be neglected. In this case DR in the growth rate, (4.143), has to be replaced, DR->\-

y/\-DR,

(4.145)

which also indicates that the present analysis of the resistive interchange mode is valid only for DR < \, i.e., for ideally stable conditions, see (4.133). Ideally unstable modes are practically not affected by a small resistivity. In a general magnetic configuration, where the coefficients Kn,Kg,Ls in (4.137) depend on the equilibrium scale £, the envelope 4>o(z) is determined by the solubility condition for 4>2 in the expansion (4.138). The resulting equation, containing only averages over the equilibrium scale, can be solved analytically, but the algebra is rather involved (Correa-Restrepo, 1982). On the equilibrium scale £, inertia and resistivity can be neglected, since we are only interested in ideally stable configurations where the resistive growth time is much longer than the Alfven time, so that the solution for marginal ideal stability suffices, 0(£) = id(0- Since, however, the coefficients depend on the details of the equilibrium, a numerical treatment is, in general, required. To obtain a smooth overall solution, the separate solutions must be matched asymptotically by identifying the logarithmic derivatives of (froiz) for z —> 0 with that of id(0 for |£| —• oo, which gives the dispersion relation. Because of finite compressibility the growth rate may become complex (see also section 4.5.2 below). The main results are: (a) If the average curvature is destabilizing pfKn > 0 or DR > 0, there are always unstable modes, the resistive interchange modes with a

4.5 Pressure-driven instabilities

125

growth rate scaling y ~ r\1^ as in the case of cylindrical symmetry. The instability threshold is DR = 0, hence shear stabilization is eliminated by resistivity. Contrary to the long-wavelength tearing modes, small-scale pressure-driven modes are dominantly electrostatic, Bo • V0 + Y\j\ ~ 0 in (4.116), which means that the magnetic perturbation \p\ produced by the current density j \ can be neglected. In a tokamak, the average curvature Kn = KC(\ — q2) is destabilizing only for q < 1, i.e., close to the magnetic axis, while in the reversed field pinch with \q\ < 1, Kn = KC is always destabilizing. (b) Also in the case of stabilizing average curvature fj'lcn < 0 or DR < 0, unstable modes exist. These are no longer driven by the basic electrostatic Rayleigh-Taylor mechanism, but by magnetic reconnection. Growth rates are slower than resistive interchange modes, scaling with larger exponents of r\. 4.5.2 Nonlinear evolution of pressure-driven modes The nonlinear dynamics of pressure-driven modes is quite different from that of current-driven kink modes. While the latter are dominated by reconnection, which slows down the nonlinear evolution compared with the linear growth rate as discussed in section 4.3.2,* pressure-driven instabilities tend to accelerate nonlinearly by a self-focusing process, which avoids reconnection. This behavior has been analyzed by Cowley & Artun (1997) for the simple model of a plane stratified plasma with a RayleighTaylor-unstable density profile g • Vp < 0. Gravity is chosen in the negative x-direction g = — gex. The system is stabilized by a horizontal magnetic field B = Bo(x)ez, which is line-tied at two boundary plates z = 0,L z , hence fc|| ~ L~l in (4.118). (A linear stability analysis has been performed by Zweibel & Bruhwiler, 1992.) Cowley & Artun use the Lagrangian form of ideal MHD to derive an approximate theory for the early nonlinear evolution. After some rather tedious algebra the authors arrive at an equation for the amplitude £(x9y,i) of the vertical plasma displacement y

^

?d2yZ + c3d2ye^

(4.146)

The function f(z) describes the sinusoidal variation along the field accounting for the line-tying boundary condition /(0) = f(Lz) = 0. The parallel mode structure is not expected to change strongly during the early nonlinear phase, since the nonlinear dynamics affects primarily the * This is, however, only true for resistivity-governed reconnection considered in this chapter. Noncollisional kink modes may exhibit an explosive behavior, see chapter 6.4.3. t Here we give Cowley's equation in a slightly simplified form.

126

4 Resistive instabilities

A Ay

o

-)

^

—

i \

y

V2_ Fig. 4.21. Schematic plot of the finger-like structure of the asymptotic solution of (4.147). perpendicular scales. Hence /(z) is assumed constant in time, which allows us to average over z thus reducing the problem to the two-dimensional equation (4.146). The first term on the r.h.s. of this equation is the linear RayleighTaylor drive (4.118), y2 > 0 corresponding to instability, y2 < 0 to stable oscillation. The second term comprises the quasi-linear effect, which is stabilizing, while the last term acts as a nonlinear drive leading to explosive growth. Ignoring the quasi-linear term, the x-dependence becomes parametric. Integration over y and further simplification leads us to the following basic equation (4.147) where the overline indicates the y-average. In the nonlinear regime (4.147) has the similarity solution -2

= 6(t0 - 1 )

(4.148)

The time to = to(y) depends on the initial conditions. If the maximum of £ is at y = 0, then near the maximum to(y) ^ to + \tf^y2, such that the width of £ shrinks to zero as t approaches the finite-time singularity to, Ay ~ J(to — t)/tfQ, as can easily be seen. Hence an initially smooth function £(y9t = 0) develops a finger-like structure shown schematically in fig. 4.21. In the solution (4.148) the_ayerage £2 in (4.147) has been neglected, which is justified for t —• to, where £2/£max ~ ty/Ly < 1. Cowley & Artun have also solved numerically the full nonlinear equation (4.146). While the scaling observed for the maximum amplitude is similar to (4.148), the

4.5 Pressure-driven instabilities

127

Flow

(a)

Expansion

Flow

(b) Contraction

Fig. 4.22. Schematic drawings of the essential processes in (a) a rising flux tube, and (b) a falling flux tube. width of the "finger" shrinks much more rapidly, Ay ~ (to — t)2, than predicted by the approximate equation (4.147), Ay ~ (to — t)1/2. Moreover, the vertical mode width increases Ax ~ (to—t)~0A, hence the mode expands into the linearly stable layers of the system. The physics leading to this explosive growth is illustrated in fig. 4.22. The stratification is such that Vp points upward, and V(p + \B2) and g point downward. The rising flux tube moves into a region of reduced ambient pressure, hence the cross-section broadens giving rise to reduced density at the tube apex, which makes the tube even more buoyant. In addition, the decreased magnetic field reduces the restoring field-line bending effect. In the adjacent (in y) falling flux tube the nonlinear effects are stabilizing, as the field is compressed. Hence the upward motion is strongly favored, leading to the finger-like structures just discussed. The dynamics found in the Rayleigh-Taylor-unstable magnetized slab studied by Cowley & Artun is qualitatively very similar to the nonlinear evolution of ballooning modes in a high-/? tokamak plasma, as demonstrated in three-dimensional MHD simulations of a toroidal plasma column by Park et a\. (1995). Here, a slowly growing n = 1 kink mode leads to a helically distorted three-dimensional quasi-equilibrium. Where the helical shift is toward the torus outside, a steep pressure gradient is generated, on which high-n ballooning modes become unstable, exhibiting explosive nonlinear growth with a finger-like structure much in the same way as in Cowley's model. The presence of magnetic shear, not included in the simple slab model, also induces a self-focusing of the pressure bulge along field lines. Because of this localization, shear stabilization is not ef-

128

4 Resistive instabilities

fective in saturating mode growth, which finally leads to global destruction of pressure confinement, contrary to the traditional idea that small-scale modes saturate at low amplitudes. Park et al show, by comparison with experimental data, that nonlinear ballooning modes can account for the characteristic features of high-/? disruptions in tokamaks. 4.53 Finite pressure effects on tearing modes In section 4.1.1 the tearing mode was treated in the simplest possible framework, low-/? reduced MHD, (4.6) and (4.7), where the instability threshold is A! = 0. For finite /? and more complicated toroidal geometry, threshold conditions and growth rates are considerably modified due to the effects of average magnetic curvature and parallel compressibility. A general theory has been developed by Glasser et al (1975). Consideration is restricted to the solution in the resistive layer, which must be matched to the ideal external solution. The quantity A', which depends on the equilibrium configuration, must, in general, be obtained numerically and is considered as a given parameter in the analytic treatment. The situation is further complicated by the presence of several resonant surfaces, the usual case in a tokamak plasma because of toroidal coupling (Connor et al, 1991). Here, we limit consideration to a low-/? circular-cross-section, large-aspect-ratio tokamak, where the dispersion relation reads (Glasser et al, 1975)

r(-)v 5/4 /

A' = 2TC-—^—-^

1-DR^J2

nl/2\

(4.149)

with the same definition of y and fj as in (4.140), and DR = f$rKnL2s, Ls = Rq/% Kn = KC(1 — q2), Af = Af/k, k = m/r. In particular kvA/Ls = (nqf/q)vA/R. If DR > 0, there is always one unstable root for any value of A', the resistive interchange mode. Since A! is independent of r\, this root is given by the vanishing of the factor in brackets on the r.h.s., as the factor in front, y 5 / 4 /^ 3 / 4 oc ?/~1//3, becomes large for small r\, hence 7 = which is identical with the dispersion relation for resistive interchange modes (4.143). For DR = 0 we find the standard tearing mode in a cylindrical plasma y = 0.55J/3/5(A')2/5, which is the usual expression (4.73) inserting the definitions of y,fj9 A'. The condition for instability is Af > 0.

4.6 Shear flow instability

o

129

-

-2x10

2x10 CO

Fig. 4.23. Solution of the dispersion relation (4.149): y and co as functions of S. The numbers are the corresponding S values (from Hender et al, 1987). The case DR < 0 is the one usually encountered in a tokamak plasma. Here the instability requires A' to exceed a finite positive threshold A' > Ac. At the threshold y is purely imaginary, y = +icoc.

The quantities coc and Ac are easily obtained by separating (4.149) into real and imaginary parts i2/3

(4.150) ,\DR\VI

A, =

(4.151)

We see that Ac —• oo for rj —> 0. Since A' is independent of ?/, there are no unstable modes for negative DR in this limit. However, the r\ -dependence of Ac is relatively weak, such that, in practice, for >/-values relevant in tokamak research, Ac is still reasonably small. Figure 4.23 gives the complex growth rate as function of S = r\~l. The standard tearing mode scaling is approximately valid for y > fj1/3\DR\2/3 or \DR\ < (rjAf3)2/5. 4.6 Shear flow instability In conventional MHD stability theory, consideration is restricted to magnetostatic configurations. By contrast, magnetic reconnection involves strong plasma flows, in particular along current sheets and separatrices. An important question therefore concerns the stability of such systems. In this section we first give a brief introduction to shear flow instabilities in neutral fluids, which is a central theme in fluid mechanics (see, e.g.,

130

4 Resistive instabilities

Drazin & Reid, 1981), and then consider the stabilizing effect of a parallel magnetic field. 4.6.1 Shear flow instability in neutral fluids Neglecting viscosity, the stability of an incompressible flow is determined by the linearized vorticity equation dtwi + v0 • Vwi + vi • Vw0 = 0

(4.152)

with wi = V2 0. Equation (4.152) now becomes

(c - V)(f - /c20) + F"0 = 0,

(4.153)

which is called the Rayleigh equation. Contrary to small-scale pressuredriven modes discussed in section 4.5.1, (4.153) does not give rise to localized modes depending only on the local value of F", but the instability involves the entire shear-flow profile V(x). A necessary condition for instability is that V(x) should have at least one inflexion point V" = 0 (Rayleigh, 1880). This condition can be derived from the energy integral obtained by dividing (4.153) by (c— F), multiplying by 0* and integrating over the entire system. The criterion excludes, for instance, linear or parabolic velocity profiles from instability. In magnetic reconnection processes, regions of high velocity shear or high vorticity are confined to narrow channels or sheets with no nearby boundaries. Here inflection points of the velocity, or extrema of the vorticity always arise, so that Rayleigh's criterion is satisfied. Consider first a simple vortex sheet corresponding to the velocity profile - 0 .

(4 154)

'

The perturbation of the configuration is conveniently described by the perpendicular displacement 0;x = £, which is related to the velocity perturbation, (4.155) Vx = -ik = dt£ + ikV£, * In this section c is the phase velocity of a fluid perturbation, not the velocity of light. This notation is chosen to conform with the literature.

4.6 Shear flow instability

131

hence u

eh = \ )

\

(4.156)

I (c - V2)£

x > 0.

v

7

Integration of (4.153) across the sheet gives f [(c - V)(f - k2cj)) + V"\ dx = cc//|i e + f (7" 0.46, where the kink mode becomes stable. The perturbations

4.7 Instability of a resistive current sheet

137

0.020 I—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r

0.015 -

0.010 -

0.005

0.5

1.0

1.5

2.0

Fig. 4.26. Growth rate y of the pinch-tearing mode for Bo = 0.8, rj = 3 x 10 4 (full), the ideal pinch mode (dashed), and the tearing mode in a static sheet pinch (dash-dotted) (from Biskamp et ah, 1998a). subsequently decay due to viscous and resistive dissipation. By contrast, the nonmagnetic jet is completely disrupted. Figure 4.27 illustrates this difference. It shows gray-scale plots of the vorticity: (a) the (slightly perturbed) initial jet, (b) the nonlinear disruption of the jet for Bo = 0, and (c) the nonlinearly broadened jet for finite field Bo = 0.3. While, in a hom*ogeneous or weakly inhom*ogeneous field, the jet is completely stable for Bo/Vm2LX^OA, in the strongly sheared field, (4.173), characteristic of a current sheet, only the Kelvin-Helmholtz mode is stable. Here, the pinch-tearing mode remains unstable leading to magnetic islands which, due to their advective motion, assume a drop-like shape with a blunt leading edge. Since small islands moving with the central jet velocity are faster than bigger ones, they catch up and coalesce with the latter, until only the largest islands remain. Such islands are called plasmoids, to which we return in section 4.7.2. 4.7 Instability of a resistive current sheet 4.7.1 Threshold condition for tearing instability The results of the preceding section indicate that a resistive current sheet is stable against Kelvin-Helmholtz instability of the sheared flow along

138

4 Resistive instabilities

(a)

(b)

Fig. 4.27. For caption see facing page.

4.7 Instability of a resistive current sheet

139

Fig. 4.27. Gray-scale plots of the vorticity illustrating the nonlinear evolution of the Bickley jet. (a) Initial state; (b) nonlinear state of the non-magnetized jet; and (c) nonlinear broadening of a magnetized jet for Bo = 0.3 (from Biskamp et al, 1998a).

the sheet, since v < VA, as explained in section 3.2.2. However, the most important candidate for instability in a sheet current, the tearing mode, is not stabilized by the velocity shear. In a static sheet pinch, for which the tearing mode is usually considered (section 4.1), the instability condition ka < 1 implies that the configuration becomes unstable if the sheet length A exceeds one or two wavelengths of the marginally stable tearing mode, corresponding to an aspect ratio A/a ^ 10. To account for the apparent stability of current sheets of much larger aspect ratio, as observed in numerical simulations, a mechanism other than velocity shear must be active. The problem has been studied by Bulanov et al (1979), who realized that the only stabilizing effect is due to the non-uniformity of the flow along the sheet. The authors assume a parallel velocity increasing linearly along the sheet

My) =

(4.174)

140

4 Resistive instabilities

which agrees with the behavior of the outflow in a Sweet-Parker sheet, see (3.124). The velocity perturbation is again assumed incompressible. Inserting (4.174) in the linearized vorticity equation d*wi + v0 • Vwi + wiV • v0 = Bo • V/i + Bi • V;o* and in the linearized induction equation, (4.6), we obtain (dt + Yydy + r)wi = Bo • V/i + Bi • V;o,

(4.175)

(dt + Tydy)xp1 = B o • V0i + rjV2xph

(4.176)

When riding on the equilibrium flow, (4.174), the wavenumber of a sinusoidal perturbation changes in time, the wavelength, i.e., the distance between points of equal phase, increases

A= n implying that k = 2n/X decreases k = -Tk.

(4.177)

This effect cancels the advective term vo • V, i.e., the explicit y-dependence in (4.175) and (4.176). With the ansatz

we obtain dtipi = (y + iky)xp\, hence (8t + vo • V)xpi = y\pi, using relation (4.177). We also have (dt + Yydy Thus (4.175) and (4.176) become (y + T ) ^ - (y - T)k2$ = ikB0(xp" - k2xp) - ikB'^xp, yxp - ikB^cf) = r\{\p" - k2xp).

(4.179) (4.180)

Bulanov et al. write the equations in terms of the perturbed magnetic field Bx = —ikxp and the perpendicular plasma displacement £ = fvxdt, where vx = —ikcj). While, in a hom*ogeneous system, this would only mean multiplying the equations by zfc, in the presence of inhom*ogeneous flow, where k = fcoe~n, the growth rate of Bx and vx is reduced by F, y1 = y—F, and £ = (y'+T)vx. In these variables, (4.179) and (4.180) assume * Because of the compressibility of the equilibrium flow, (4.174), the V • v0 term has to be included in the vorticity equation.

4.7 Instability of a resistive current sheet

141

the form

(/ + F)(/ + 2F)£ - / ( / + F)£ = ik(B'; - k2Bx) - ikBZBx, (/ + r)(Bx - ikBot) = nW -

klB

*)-

(4.181) ( 4 - 182 )

Applying the same procedure as in section 4.1.1 or 4.2.1, we find the dispersion relation

(/ + r) 4 (/ + 2r) = y05,

(4-183)

where yo is the growth rate in the static sheet pinch. Solving (4.183) one finds that the tearing mode becomes stable for T > O.87yo- By contrast, the dispersion relation following from the \p9(j> system, (4.179) and (4.180), y4(y + F) = yo> predicts a significantly higher threshold value of P. Hence we are faced with an ambiguity in the stability problem, which requires us to specify whether we wish to consider stabilization of B x , i.e., the magnetic energy associated with the tearing mode, or of xp, i.e., the magnetic flux. We prefer the former definition, since MHD stability is usually discussed in the energy picture. The stabilizing effect of the inhom*ogeneous outflow can easily be understood. The tearing instability corresponds to a local condensation or bunching of the current density, which is counteracted by the stretching due to the flow. One may assume that the mode is effectively stabilized if the change of the wavelength during one growth time exceeds a substantial fraction of the wavelength, say one-half, F/yo > 0.5. Hence the sheet is tearing-mode stable for T > 0.5ymax. In the asymptotic limit of large S one has ymax — 0.5(1^1^)~1^2 for k§5 ^ S" 1 / 4 , (4.24). As discussed in sections 3.2.2 and 3.2.3, the following relations are valid in a resistive current sheet:

T = vA/A ~ u/8 ~ nib1 = i" 1

and the local Lundquist number is S o = ?n/?A = vAS/rj = A/S.

Hence stability requires

or 8/A > 0.05. Since, however, for such low values of So the maximum growth rate is substantially smaller than predicted by the asymptotic expression, we may estimate that a resistive current sheet is susceptible to tearing instability only for aspect ratio So = A/8 > 102. (4.184)

142

4 Resistive instabilities

Fig. 4.28. Plasmoid generation in a tearing-unstable Sweet-Parker sheet. The figure shows tp-contours at times t\ < ti < h < U (from Biskamp, 1986). The tearing instability problem in a current sheet has also been investigated by Phan & Sonnerup (1991), who assumed that the flow is diverted in the sheet along the z-direction instead of the y-direction. Since there is no stabilizing flow in this case, the stability limit obtained, So ^ 12.5, is close to that of a static sheet. 4.7.2 Plasmoids From the scaling laws of driven reconnection, (3.151) and (3.154), A/d ~ r\~2 (for variable sheet length A), or, (3.150), A/S ~ Y\~2^ (for fixed A), it is clear that the tearing mode will arise in a Sweet-Parker sheet for sufficiently small r\. Let us therefore consider its nonlinear behavior. The

4.7 Instability of a resistive current sheet

143

Fig. 4.29. (a) Magnetic structure and (b) flow pattern of a plasmoid (from Biskamp, 1986). usual appearance of the tearing mode in a static sheet pinch is that of a chain of magnetic islands. In a Sweet-Parker sheet the more typical case is that of one or a few isolated islands, which are convected along the sheet, growing in a self-similar way both in width and in length being stretched by the inhom*ogeneous flow. Such isolated islands are called plasmoids. Figure 4.28 illustrates the repetitive generation and convection of plasmoids. Shown is a configuration of driven reconnection such as in fig. 3.16 (for clarity the lower symmetric quadrant is added in the plot). Changing the boundary conditions at y = 0 slightly, from dy\p = 0 to \p(x, t) = \p(x)+Et9 where E is the externally imposed average reconnection rate and \p(x) is the xp-profile taken from case (c) of fig. 3.16, is sufficient to destroy stationarity and give rise to repetitive generation of plasmoids. The dynamics of the plasmoid becomes clear from fig. 4.29(b), where the density of streamlines demonstrates the high speed of the plasmoid along the original sheet compared with the velocity in the ambient medium. Plasmoid formation and ejection is a general phenomenon, which occurs if the sheet aspect ratio A/d becomes large enough. Figure 4.30 shows such an event during the nonlinear evolution of the resistive kink mode. It also occurs in the extended noncoUisional current sheets dominated by electron inertia. Fully developed plasmoids assume an asymmetric droplike shape with a blunt front and a long tail corresponding to a current sheet, which may be thinner than the original stationary one and which is hence particularly sensitive to further plasmoid generation. Since such secondary slim plasmoids propagate faster than the primary bulky one, the former will catch up and finally coalesce with the latter. Thus the tearing

144

4 Resistive instabilities

Fig. 4.30. Plasmoid dynamics in the current sheet of the nonlinear resistive kink mode (from Biskamp, 1986).

instability can strongly increase the reconnection efficiency compared with that in a stationary current sheet. Generation of plasmoids is a widespread phenomenon. Plasmoids are an essential feature in many rapid MHD processes, such as an erupting magnetic arcade in the solar corona, the standard model of a large flare (see, e.g., Biskamp, 1993a, chapter 10), or the dynamics in the geomagnetic tail connected with the substorm phenomon, which we discuss in chapter 8. Tearing instability (and the resulting plasmoid dynamics) are also an important mechanisms for the generation of small-scale fluctuations in MHD turbulence (see, for instance, Biskamp, 1993a, chapter 7).

Dynamo theory

Dynamo theory deals with the generation or, more precisely, the amplification and sustaining of magnetic fields by motions in an electrically conducting fluid. The topic is of eminent practical importance, since, as is generally believed, the magnetic fields in astrophysical objects as diverse as planets, stars, accretion discs and galaxies are generated in this way. The property of the conducting medium to be a fluid is important, since dynamo action by the motion of solid ducts is commonplace as the general means of electric power generation. Here the simplest model is the disc dynamo proposed by Bullard (1955) sketched in fig. 5.1, which shows the essential features of dynamo action. The basic equations are (5.1)

LI +RI = Cl = F-

c2l\

(5.2)

where L and R are the inductance and the resistance of the system, F is the external torque on the disc, / is the current generating the magnetic

Fig. 5.1. Bullard's disc dynamo. 145

146

5 Dynamo theory

(b) Fig. 5.2. Disc dynamo models exhibiting spontaneous field reversals: (a) shunted disc dynamo; (b) coupled disc dynamos. field penetrating the disc, and c\, ci are geometrical factors. For small / the angular frequency Q of the disc is constant, such that / =/ o exp(ciQ — R)t. I is growing exponentially, if Q is sufficiently large to overcome the resistive damping. In the nonlinear regime there is a steady-state solution

It is easy to see from fig. 5.1 that a seed field in the opposite direction will also be amplified. If, however, the rotation is reversed Q —> —Q, field amplification does not occur, hence dynamo action depends on the relative sense of the current loop winding and the disc rotation. The general nonlinear solution is oscillatory, but does not reverse the sign of / or B. Field reversal does, however, occur in a shunted disc dynamo (Malkus, 1972; Robbins, 1977), illustrated in fig. 5.2(a). Here, the current / flowing in the disc is branched into the current Ic in the coil around the the axis and the current Is in the shunt. Written in terms of /, / c , and the angular velocity Q, the generalized disc dynamo equations are equivalent to the well-known Lorenz model, a truncation of the equations of two-dimensional convection (Lorenz, 1963). Hence for certain values of the external parameters, there are solutions oscillating in a nonperiodic form between two states of opposite field orientation. Field reversals also occur in a system of two coupled Bullard disc dynamos (Rikitake, 1958), illustrated in fig. 5.2(b). In a disc dynamo the current is flowing in a multiply connected region consisting of the wire and the disc. Dynamo action in a conducting fluid

5 Dynamo theory

147

occupying a singly connected region is more difficult to achieve, even the very existence of a dynamo effect in such a fluid has been questioned for a considerable period. We start this chapter with an outline of kinematic theory, section 5.1, discussing the topological properties of the fluid motions necessary for dynamo action. The essential requirement is finite (kinetic) helicity - the flow must, loosely speaking, be sufficiently twisting - since plane or axisymmetric flows cannot lead to field amplification, which is the essence of the anti-dynamo theorems. A major step from the proof-of-principle models to a theoretical tool for practical applications has been the introduction of mean-field electrodynamics (MFE). Here the fluid motions are assumed small-scale and turbulent, as they are expected to be in most systems of interest. The details of these motions have no effect on the behavior of the large-scale magnetic fields, only two average quantities being important, the turbulent kinetic energy leading to an anomalous magnetic diffusion and dissipation, and the helicity which through the a-effect may give rise to field amplification. MFE is the framework for most theoretical studies of stellar or galactic dynamos. While the basic mechanisms of small-scale field amplification, namely flux tube stretching and twisting, do not require reconnection, the generation of large-scale observable fields from small-scale ones do, hence reconnection is a basic ingredient of dynamo theory. Though not introduced explicitly, MFE of course involves reconnection, as will easily be seen. In its best-known and commonly used form MFE is a kinematic theory, where the fluid turbulence is determined by nonmagnetic mechanisms, for instance thermal convection in stellar dynamos. Since the magnetic field is primarily amplified at small scales, the small-scale fluid motion will be modified by the magnetic force even at still low mean magnetic energy, which reduces very efficiently their contribution to the a-effect, such that progressively larger turbulent scales become involved. This process requires a self-consistent MHD description. Section 5.2 introduces self-consistency in the framework of mean-field theory by modeling the reaction on the fluid motion in a phenomenological way. Then an overview is given of the present theory of the solar dynamo, which is essentially restricted to MFE. Parker's original idea of an aQdynamo working in the solar convection zone (Parker, 1955a), which was considered as a solid theoretical basis for many years, is again put into question by recent helioseismological evidence, such that solar dynamo theory remains a challenging problem. In section 5.3, direct numerical simulations of convection in a spherical shell are presented. Though the properties of the fluid are highly idealized and the Rayleigh numbers are still rather low, these simulations give a reliable picture of the basic dynamo mechanism which is not affected by ad

148

5 Dynamo theory

hoc assumptions. Of all the astrophysical dynamo systems, the properties of the Earth's liquid core come closest to the conditions simulated in these computations. I therefore briefly review the present understanding of the geodynamo. High-Rayleigh number convection typical, for instance, for stellar convection zones, involve fully developed MHD turbulence, which is discussed in section 5.4. This is essentially a sequel of chapter 7 of my previous book (Biskamp, 1993a). After a brief introduction to the basic properties of MHD turbulence, the section is concentrated on recent developments in 3D turbulence. Dynamo theory is a vast field of its own, which cannot fully be covered in this brief presentation. Many aspects are touched only cursorily or not at all, in particular the more mathematical aspects of dynamo theory, which have developed into a topic in their own right rather independent of the more practically oriented fluid mechanical approach (see, e.g., Bayly & Childress, 1987, 1988; Finn & Ott, 1988). For a more comprehensive treatment of dynamo theory I refer the reader to several special treatises, e.g., Moffatt (1978), Parker (1980), Zeldovich et al, (1983), Roberts & Soward (1992).

5.1 Kinematic dynamo theory In kinematic theory the fluid velocity v(x, t) is given, and the problem consists in determining the development of the magnetic field, i.e., we have to consider only the induction equation dtB = V x (v x B) + >/V2B,

(5.3)

which is a linear hom*ogeneous equation for B. Kinematic dynamo theory answers the question about the stability of a nonmagnetic flow to infinitesimal magnetic disturbances, which is formally similar to the stability theory discussed in chapter 4. We will learn about the conditions that must be satisfied for magnetic field amplification and to a certain extent also about the spatial structure of this field. Nonlinear saturation by reaction of the magnetic field on the fluid motion is considered in section 5.2. The basic physics of the terms in (5.3) is similar to that in (5.1) for the simple disc dynamo, which suggests that for sufficiently large v dynamo action should be possible. Since, on the other hand, the topology of the setup in the disc dynamo plays an important role in its functioning, we expect that in the fluid problem the structure of the flow must satisfy certain requirements.

5.1 Kinematic dynamo theory

149

Fig. 5.3. Magnetic configuration illustrating Cowling's theorem. 5.1.1 Nonexistence theorems and special dynamo solutions In fluid dynamics, in particular in MHD, problems are often simplified by assuming two-dimensional geometry. In dynamo theory, however, this proves to be too restrictive. Consider the simplest case of plane geometry dz = 0, where B has the form = e7 x such that (5.3) reduces to two equations for xp and Bz dtxp + v • Vxp = rjV2xp,

(5.4)

dtBz + v • VBZ = Bp • Vvz + rjV2Bz,

(5.5)

assuming, as commonly done, incompressibility V • v = 0. While the axial field Bz can be amplified by the poloidal field B p , via the first term on the r.h.s. of (5.5), there is no reaction of Bz on Bp. In fact the poloidal magnetic energy decays, as can be shown directly by multiplying (5.4) by xp and integrating over space and time < - / xp2dV

t=Q

= const,

(5.6)

hence the driving term in (5.5) tends to zero and so does the axial field energy J B2dV. A similar argument can be used in the axisymmetric case. It is, however, instructive to recapitulate briefly the proof presented by Cowling (1934) to demonstrate his famous theorem "A steady-state axisymmetric field cannot be maintained". An axisymmetric field has at least one magnetic axis C, where the poloidal field vanishes, see fig. 5.3. Integration of Ohm's law along this closed line gives

J)E'dl+(l) \xB-d\=

(5.7)

150

5 Dynamo theory

Since

J>E'dl= hxE-dF

= - fdtB-d¥ = 0

because of the assumed stationarity, and

I v x B • d\ = I v x Bp • d\ = 0, J

J

because of B p = 0 on C, also the current density on the r.h.s. of (5.7) must vanish on the magnetic axis, which implies that Bp = 0 not only on C but in a full surrounding thereof. Extending the argument to this region, it follows that Bp must vanish throughout. (There are some subtleties, which must be considered in a rigorous proof, but these do not change the basic result.) The preceding discussion might suggest that the anti-dynamo theorem applies to all symmetric field configurations, i.e., those exhibiting an ignorable coordinate. This is, however, not true for the most general case, that of helical symmetry r,u = Icj) + fez, where B(r, u) has the form B = h x WxpH + h/, h

= /TTfev V r x V u >

h v

' / = °>

— IAZ being the helical flux function, (3.37), which follows the equation

+ v • Vxpn = n \LxpH -

/2 + k2r2

fj .

(5.8)

Here L is the generalized Laplacian

while / , the field component in the helical direction, obeys the equation (5.9) F is a rather complicated linear functional of / and \pn generalizing the r.h.s. of (5.5), which we do not need to specify. The important point to note is that the xpn equation (5.8) now contains a reaction of / on \pn absent both in plane geometry (k = 0) and axisymmetry (/ = 0), which allows dynamo action. An explicit steady-state helical dynamo solution has been given by Lortz (1968). Such helical fields have a certain relevance for the dynamo effect in the reversed-field pinch (RFP) (for a review of RFP theory see Biskamp, 1993a; Ortolani & Schnack, 1993).

5.1 Kinematic dynamo theory

151

Steady spatially periodic flows have received particular attention as efficient dynamos, that most commonly studied being the so-called ABC flow v = (A sin z + C cos y, B sin x + A cos z, C sin y + B cos x), which has the property of maximal helicity since v = w, w = V x v. (A flow with w = h is called Beltrami flow, which is the hydrodynamic analog of a force-free magnetic field j = XR). A special case of a two-dimensional ABC dynamo was first analyzed by Roberts (1972). Helicity is the essential requirement for dynamo action, as is obvious from the a-effect introduced in section 5.1.3. It should, however, be noted that the presence of helicity does not guarantee field amplification or sustaining, since an axisymmetric flow, which does not give rise to dynamo action, also has finite helicity for nonzero azimuthal velocity v^ HK = IJdVv^W x vp)^, as can easily be checked. The first dynamo solution in a conducting sphere with a nonvanishing dipole field outside the sphere has been developed by Herzenberg (1958). Here the velocity consists of two spherical rotors of constant angular velocity embedded in the sphere, with the fluid otherwise at rest. Assuming the rotor radii to be small compared with the distance between their centers makes an analytic calculation of the field possible. (The analysis can be generalized to n rotors, see for instance Moffatt (1978), where references to the original papers are found.) While the actual calculations have been restricted to special rotor configurations and orientations, dynamo action seems to occur for rather general conditions, though the magnetic field may not be stationary but oscillating. Models like Herzenberg's rotor configuration are of course very artificial, serving only as a proof-of-principle of the existence of motions that give rise to dynamo action. Before moving on to consider systems of more practical relevance, it is useful to illustrate the basic process of field amplification in a simple mechanistic model. 5.1.2 The rope dynamo Vainshtein & Zeldovich (1972) proposed the following elementary process called a rope dynamo, which demonstrates the possibility of exponential growth of the magnetic field, see fig. 5.4. Consider a closed flux tube of cross-section F carrying the flux (/>. Stretching the tube to twice its original length reduces the cross-section by one-half because of mass conservation (assuming incompressibility), while the flux remains unchanged. Twisting and folding the tube as indicated in the figure leads to a tube with the original cross-section but twice the flux. Repeating the process results in

152

5 Dynamo theory

2D Fig. 5.4.

3D

Basic mechanism of the rope dynamo in 2D and in 3D.

exponential growth of the flux 711112

It is interesting to note that this stretch-twist-fold process does not involve reconnection, it occurs on the fastest possible time-scale - the advective time - and appears to be typical for the initial amplification of a seed field in real systems. Though of course under realistic conditions the process of folding back cannot be perfect, which will make the field configuration more and more filamentary, the VainshteinZeldovich process is, at least in principle, possible. In two-dimensional geometry, by contrast, only stretching can occur. Whereas a flux tube of the original cross-section can be generated, as indicated in fig. 5.4, not only does this require reconnection but, more importantly, the field

5.1 Kinematic dynamo theory

153

inside the combined tube is now anti-parallel, hence only the magnetic energy not the flux has increased. By repeating the process, the cross-sections of the unidirectional-field filaments shrink, until resistive diffusion leads to field annihilation, prohibiting further growth of the energy, which ultimately decays, consistent with the anti-dynamo theorem. The rope dynamo is the prototype of a so-called fast dynamo,* as contrasted with a slow dynamo wherein the growth rate y depends on the resistivity, or more precisely on the magnetic Reynolds number ReM = vL/ri, y ~ Re]J,0 < v < 1. This point has attracted considerable attention in the astrophysical community since, because of the typically huge value of ReM slow dynamos may need more time than is available to produce, for instance, the galactic magnetic field. However, the strict distinction between fast and slow dynamos results from the use of strongly simplified models and should not be overemphasized. It will become clear subsequently that the generation of observable large-scale fields, such as the dipole field of the Earth, requires reconnection and hence cannot be fast in the strict sense, but that, on the other hand, turbulence may give rise to reconnection rates independent of the Reynolds number. 5.1.3 Mean-field electrodynamics In the dynamo models considered so far the velocity field v(x, t) was treated as a known function, which has been chosen in a rather ad hoc way. We now take a first step toward self-consistency. In most systems of practical relevance, both in astrophysical and laboratory plasmas, the fluid motions responsible for dynamo action are caused by some instability buoyancy-driven in stars and planets or current-driven in the reversed-field pinch - and are therefore expected to become turbulent for the typical high Reynolds numbers. Moreover, these fluid motions have spatial scales which are small compared with the observed large-scale magnetic fields, these motions generate. This suggests a two-scale statistical approach called mean-field electrodynamics (MFE), which was developed in a pioneering paper by Steenbeck et a/., (1966) (for a more recent review see Krause & Radler, 1980). One assumes that v and B can be decomposed into an average part * The Vainshtein-Zeldovich picture of the rope dynamo is oversimplified, as noticed by Moffatt & Proctor (1985), since by the twisting (= kinking in the language of section 2.4) of the rope the field inside the rope becomes more and more twisted because of magnetic helicity conservation. At some point this internal twist will relax, which involves reconnection. Hence the rope dynamo can only be fast, i.e., independent of rj, for a finite time.

154

5 Dynamo theory

varying only on the large scale and a small-scale fluctuating part: v = vo+?,

B = B 0 + B,

(5.10)

(v) = ( B ) = 0 . Insertion into the induction equation yields equations for the mean and the fluctuating parts Bo and B: dtB0 = V x (v0 x Bo) + V x e + ^/V2B0, dtB = V x (v0 x B + v x Bo) + V x [(? x B) - (v x B)] + ^V2B,

(5.11) (5.12)

where e = (vxB)

(5.13)

is the mean electromotive force (the negative electric field, note the sign convention) induced by the turbulent fluctuations. By solving (5.12), e can be expressed in terms of v and Bo. An explicit solution is obtained in the quasi-linear approximation, where the nonlinear term in the fluctuation equation (5.12) is neglected, which is formally justified when the fluctuation amplitudes are sufficiently small. Moreover, one usually assumes, for convenience, that the dominant scales of the turbulent fluctuations are still large compared to the dissipative scales, so that the diffusion terms in (5.12) can be omitted. We thus find

B(x,t) = f J— 00

di!V x [v(x',0 x B0],

(5.14)

where the time-integral is along the Lagrangian orbit x' = x(tf) and the weak space-time dependence of the mean field is neglected. Insertion into (5.13) gives

e= f

*'(?(x,0 x {V x [?(x',0 x Bo]}}.

(5.15)

J — 00

On the assumption that v is statistically independent of Bo, e is a linear functional of Bo. Making use of the separation of scales we can conclude that e depends only on the local properties of Bo, i.e., on Bo and its first-order spatial derivatives, ei = *ijBoj + PijkdjBok.

(5.16)

A further simplifying assumption usually made for convenience is that the turbulent field v is isotropic, whence ay = a> Ca the a-effect dominates and the differential rotation can be ignored in (5.23). Such a system is called an a 2 -dynamo. In steady state, (5.22) gives the estimate Bp ~ CJi^ and (5.23) gives BAJ ~ CaBp, hence the value of Ca for field sustainment is Ca ~ 1, and the poloidal and toroidal magnetic fields are of the same order of magnitude, BA

~ B

(5 27)

In the opposite limit Ca > Ca one can neglect the a-term in the I n equation, which now gives the steady-state condition B^ ~ CaBp. Together with Bp from (5.22) we thus have D = 1, where D is called the dynamo number, and BA

Bp.

(5.28)

Such a system is called an aQ-dynamo, which carries a toroidal field much larger than the poloidal field. An aQ-dynamo is therefore also

5.1 Kinematic dynamo theory

159

characterized as "strong-field dynamo", compared with an a 2 -dynamo, where both fields have similar magnitude, which is called "weak-field dynamo". In many astrophysical objects, such as the Sun or the Earth, the toroidal field is unobservable, concealed in the interior, since only the poloidal field extends outside the conductive region. In this case, the toroidal field can only be calculated or estimated theoretically using the observable properties of the object. From these it is concluded that both the geodynamo and the solar dynamo are of afi-type. 5.1.5 Free dynamo modes The kinematic dynamo equations (5.22), (5.23) are linear in the magnetic field and thus for stationary mean velocity fields v and dynamo parameters a, rjr and given boundary conditions have eigenmodes proportional to eyt. In general, v and a are functions of space, which would require to solve the equations numerically. Some important qualitative properties can, however, already be disclosed in the local or WKB approximation, where vp,VQ, a, t]T and the geometrical factors r can be considered constant. The concept of dynamo waves has been introduced by Parker (1955a) and Yoshimura (1975) to explain the period and migration direction of the solar magnetic field. With the usual Fourier ansatz A, B^ oc exp{ik-x} one obtains (y + r]Tk2 + ik - yv)A - ocB^ = 0, (5.29) (y + rjTk2 + ik • yp)B^ + [(ik x rVQ)^ - ock2]A = 0, 2

(5.30)

2

using V x Bp = —k A. The simplest case is that of the a -dynamo in the absence of a mean flow, where the solvability condition of (5.29) and (5.30) gives the dispersion relation (y + r]Tk2)2 = a2k2.

(5.31)

Hence the mode is purely growing or decaying, Im{y} = 0. Dynamo action y > 0 occurs, if |a| > r\jk, which is independent of the sign of a. In the opposite case of an aQ-dynamo, where the a-term in (5.30) can be neglected, we find the dispersion relation p

(5.32)

where T = ±a(k x rVQ)^

(5.33)

The character of the solution depends on the sign of F, i.e., the direction of k for a given sign of a. First take F > 0. With yj—2i = +(1 — i) the solution of (5.32) is y

= -r,Tk2 ± F 1 / 2 - i(k • \p ± F 1 / 2 ),

(5.34)

160

5 Dynamo theory

i.e., Re{y} > 0 for the upper sign and F 1 / 2 > v\jk2. Contrary to the a 2 -dynamo the magnetic field is in general oscillating with frequency co = — Im{y} = k • \p + F 1 / 2 corresponding to so-called dynamo waves. For sufficiently small poloidal velocity \p, the wave propagates in the +k-direction, but may reverse its direction if \p is large, k • yp = —F 1 / 2 leading to a stationary wave pattern. In the case F < 0 the solution of the dispersion relation reads, using

y = -nrk2 ± IFI1/2 - i(k • v, + IFI1/2),

(5.35)

such that the frequency of the dynamo wave is co = k\p—|F|1//2. For \p = 0 the dynamo wave (upper sign in (5.35)) now travels in the —k-direction. Besides the frequency and the propagation direction of the dynamo wave, the eigenmode equations (5.29), (5.30) also predict the phase relation between the poloidal and the toroidal fields. For the most relevant case of marginal stability and \p = 0, (5.29) gives A = i((x/a>

and hence for Br = —ikzA

Br = - 4 - * * .

(5.36)

co/kz We shall see in section 5.2.2 that both the prediction for the phase velocity of the dynamo wave and for the phase between Br and B^ have stringent consequences for the solar dynamo mechanism. Let us finally discuss the symmetry properties of the eigenmodes of (5.22) and (5.23). With the usual symmetry assumptions with respect to the equator z = 0 or colatitude 6 = 90°, a and VQ odd, vr and Q even, the eigenmode parities are either A even, B^ odd,

(dipole)

or A odd, Bcj) even,

(quadrupole).

As will be seen later, the dipole mode is the relevant one for most astrophysical dynamos. Since the pioneering work of Steenbeck et a\. (1966), the kinematic mean-field theory has become the standard tool for modeling dynamo action in various astrophysical objects. Choosing the appropriate geometry and profiles a(x), Q(x) based on the properties of the underlying convective turbulence and consistent with observations, one considers the dynamo mode with the lowest dynamo number D = Q C Q . This procedure gives

5.2 Mean-field MHD dynamo theory

161

an approximate picture of the spatial structure and the time variability of the magnetic field of the object. Though the conventional assumption of isotropy of the turbulence, which leads to the simplification a^ = ad^, is hardly ever really justified, anisotropy does not seem to change the dynamo behavior qualitatively. 5.2 Mean-field MHD dynamo theory If dynamo action occurs kinematic theory tells us that the magnetic field with a certain spatial structure corresponding to a particular dynamo mode grows exponentially in time. To describe the saturation of this process the reaction on the fluid through the Lorentz force must be included, such that the fluid motion adjusts self-consistently. In this section we will first discuss some phenomenological models which have been introduced in the framework of mean-field theory. We then give an overview of the present understanding of the solar dynamo, the prototypical application of mean-field theory, pointing out the difficulties of the theory which have emerged in the light of the recent helioseismological observations. 5.2.1 Nonlinear quenching processes and magnetic buoyancy Nonlinear effects enter the dynamo process in two different ways: through a finite-amplitude modification of the fluctuating quantities invalidating the quasi-linear approximation, which is called micro-feedback, and through the Lorentz force acting on the mean flow, called macro-feedback. The most important process in the first category is the nonlinear reduction of a, called a-quenching, which acts directly on the development of the mean magnetic field Bo in the induction equation (5.21). The physical origin of this process is the Alfven-wave effect, which is considered in more detail in section 5.4. Basically it means that, in the presence of a mean field, small-scale velocity and magnetic field fluctuations tend to be coupled, forming Alfven waves, (3.22), B ~ +v in the present normalization. This correlation, which competes with the local (in /c-space) decorrelating interactions of the small-scale fluctuations, is the stronger the larger the mean field. Since the electromotive force e = (v x B) vanishes for pure Alfven waves, e and in particular a, should decrease as the mean field grows due to the a-term, which ultimately leads to a saturated dynamo state. A rigorous theory of a-quenching has not yet been developed. Instead the process is modeled in a simple phenomenological way. Since the Alfven-wave effect depends only on the magnitude of the mean field, a

162

5 Dynamo theory

must be a decreasing function of Bo. The following ansatz is used in many applications (e.g., Riidiger, 1973; Ivanova & Ruzmaikin, 1977) a = ao/(l + CMBO),

(5.37)

where CM is a constant and small n, n = 0(1), implies a smooth reduction and large n an abrupt quenching.* To discuss the reaction on the mean flow, we consider the momentum equation, again in the incompressible limit, p( p is possible. If this process is dynamically accessible, only convective transport of B out of the production region can limit the field intensity, either by magnetic buoyancy or by poleward convection (as a strained rubber band slides on a smooth sphere). 5.2.2 The solar dynamo Let us briefly summarize the main observations concerning the magnetic field of the Sun (for more details see, e.g., Priest, 1984): (a) The most obvious manifestation of this field are sunspots, which usually occur in pairs in a latitude belt of ~ +30° about the equator. The two partners of a pair have opposite magnetic polarity; they appear to be connected by field lines, which emerge through the solar surface in one spot and submerge again in the partner spot. The orientation of the connecting field is essentially in azimuthal direction 0, with one orientation in the northern hemisphere and the reversed one in the southern (Hale's polarity law). (b) Apart from the azimuthal field associated with sunspots, there is a weaker poloidal component which is distributed on a grainy network all over the solar surface. Its average radial component also reverses sign from northern to southern hemisphere, hence has an approximate dipole character. (c) The most remarkable feature of the solar magnetic field is its regular time variation. Both azimuthal and poloidal field reverse sign every 11 years. This gives rise to the sunspot cycle, starting from a minimum number of spots (close to zero), proceeding to sunspot maximum and back to minimum. The full period including both polarities is 22 years. (d) The B(j) and Br components are in anti-phase, i.e., B^By < 0 for most of the time. While the field inside the sunspot latitude belt migrates toward the equator (which gives rise to the famous butterfly diagram, fig. 5.8), the more diffusely distributed field outside this belt migrates toward the poles. (e) Because of the very intermittent spatial distribution, the field intensity varies strongly. In sunspots the field reaches several kilogauss, which is about the limit given by cross-field equilibrium B2/8n ~ p

165

5.2 Mean-field MHD dynamo theory

V

40 30 20 10 0

-10 -20 -30 -40

i 111111 1875

1111111111111111111111

u

1880 1885 1890 1895 1900 1905 1910 f(year)

Fig. 5.8. Butterfly diagram (from Maunder, 1913). It shows the incidence of sunspots as a function of colatitude and time. A vertical bar is plotted for each degree and each 27-day period if, during this period, at least one sunspot is observed in this angular range. Thefiguredemonstrates the migration of sunspots toward the equator during the 11-year cycle. in the photosphere. Also in the thin flux tubes of the network called magnetic knots, where most of the remaining field is concentrated, the field seems to reach this magnitude. Since sunspots and network knots cover less than one percent of the solar surface, the spatially averaged field is much weaker, of the order of a few gauss. (f) The Sun possesses an extended convection zone. While the hot central part of the solar interior is in static equilibrium dominated by radiative transport, in the cooler outer shell opacity is increased by ionization and recombination processes, which gives rise to steeper temperature gradients and the onset of thermal convection (see section 5.3.1). The lithium depletion observed on the Sun suggests that the convection zone reaches down to the lithium-burning region, i.e., to about 0.7 of the solar radius. The visible manifestation of the turbulent convection are the granular structure of the photosphere and the associated motions. (g) The Sun is rotating with a mean period of 27 days; the rotation is not rigid, the surface value of the angular velocity being about 10 percent larger at the equator than at the pole. The variation of Q inside the convection zone, in particular the radial dependence, has recently been determined from helioseismology. Contrary to the

166

5 Dynamo theory original theoretical picture, Q-contour surfaces are not cylindrical but more disc-like.

(h) The rotation of the Sun is relatively slow and its magnetic field weak. Observations of the magnetic intensity and its time variation in other more rapidly rotating stars shows a close relation between the rotation rate and the mean magnetic field, (B) oc Q, which corroborates the dynamo interpretation. The time variation of the field essentially rules out a primordial origin (though the decay time of such a field in the non-turbulent solar interior would be of the same order as the age of the solar system and it is legitimate to ask what happened to the primordial field the Sun should have captured during its formation). The most plausible cause of the observed field is a dynamo process due to turbulent fluid motions in the convection zone. From the characteristics of sunspots and the poloidal field this should be an aQ-dynamo: differential rotation transforms a poloidal into a toroidal field (fig. 5.5), of which the sunspot field is the direct manifestation, while the helicity of the convective turbulence regenerates the poloidal field. A theory of the solar dynamo should account for the solar cycle time (at least give the correct order of magnitude), the equatorward sunspot drift and the phase relation between B^ and Br. Because of the extremely high Rayleigh number (for a definition of the Rayleigh number and other quantities see section 5.3.1), direct numerical simulations of the entire convection zone are out of the question and one has to resort to a mean-field approach characterized primarily by the dynamo coefficients a, rjr (for recent reviews of mean-field solar dynamo theory see Radler, 1990; Stix, 1991; Gilman, 1992; Brandenburg, 1994). Typical values of the velocity fluctuation amplitude and the correlation length in the photosphere are v ~ 1 km/s and / ~ 103 km. In the convection zone / increases with depth, while v decreases such that vl is roughly constant, hence by (5.25) riT~\vl~'Sx

102 km 2 /s.

(5.44)

By its very nature a is a more complicated effect, which is usually estimated by oc ~ ea£ll, (5.24). Here ea is a numerical factor, which will be discussed further below. Since / must be taken for the region of maximum dynamo action which, as we will see, is close to the bottom of (or even below) the convection layer, we assume ! ~ 3 x 10 4 km and, with fi~3x 10~6s~1, one obtains OZ~ 0.1 km/s. (5.45)

5.2 Mean-field MHD dynamo theory

167

Since the original paper by Parker (1955a), the fundamental concept in solar dynamo theory is a dynamo wave of stationary amplitude with appropriate frequency and propagation direction. The condition Re{y} = 0 in (5.34) and (5.35) gives the frequency (for negligible poloidal flow \p) co = T1'2 = rjTk2,

(5.46)

where k = 2n/L and L is the global scale of the dynamo process, which is of the order of the solar radius JR. Hence, using (5.44), the cycle time would be x ~ R2/2%r\T ^ 3 x 108 s ~ 10 yr (5.47) which, compared with the observed value T = 22 yr, is of the correct order of magnitude. The imaginary part of (5.34) or (5.35) gives a relation between a and the radial scale-length IQ of the differential rotation 3 r lnQ = Z^1, co = Y1'2 = ((xkRQ/ln)1/2, (5.33). Writing kR ~ 2n and a = eJQ, we find the relation

300

t1^)

(5.48)

co \eJJ Hence for large helicity ea ~ 1 the differential rotation needed for dynamo action is very low, IQ ~ 3 x 109 km corresponding to Rdr In Q ~ 10~~3 —10~4, which is not plausible. To avoid this conclusion the effective a should be strongly reduced compared with the estimate (5.24), i.e., ea ~ 10~2 or a ~ 1 m/s. The more stringent condition is posed by the observed equatorward migration, which corresponds to F < 0. (This can be seen directly from (5.35): since Re{y} = 0 requires the upper sign, a positive (poleward) meridional flow would reduce the frequency, hence the inherent migration of the dynamo wave is equatorward.) This implies a3 r Q < 0. As the arguments for a > 0 are quite strong (see section 5.1.4), this seems to require drQ < 0, i.e., Q increasing with depth (sub-rotation). (The alternative possibility a < 0, 3rQ > 0 (super-rotation) is ruled out also because in this case BrB^ > 0 from (5.36) contrary to the observation.) Observational results from helioseismology (fig. 5.9) show sub-rotation only below the poles, while below the equator there is a (weak) super-rotation. These properties of the differential rotation have also be obtained theoretically by solving the mean-field equation for Q (see, e.g., Kuker et a/., 1993). There seems now to be rather general agreement that the only way out of this dilemma is to locate the main dynamo action at the bottom of the convection zone or even in the convectively stable overshoot region below the convection zone. Fully three-dimensional numerical simulations of magnetoconvection, which will be discussed in section 5.3, reveal that in fact the magnetic field, mainly toroidal, is strongest in this region.

168

5 Dynamo theory

Fig. 5.9. Contours of the solar differential rotation Q (from Libbrecht, 1988). The field intensity depends on the mechanism by which the field is transported across the convection zone to the surface of the Sun. The process must be sufficiently rapid so that the Coriolis force does not tilt the field substantially out of the azimuthal direction. Basically two different mechanisms are conceivable. The field could be swept along with rapid updrafts. In this case the field intensity at the bottom would have to be a few kilogauss. If this transport process is not effective, magnetic buoyancy has to be the dominant mechanism, in which case the field has to be higher by at least one order of magnitude to provide the necessary buoyancy force. The question about the dominant dynamo process on the Sun is far from being decided. It might be that the mean-field theory itself is not the appropriate framework and that a more direct approach is required. 5.3 MHD theory of thermal convection 5.3.1 Thermal convection in a rotating sphere The engine for dynamo action in stars and planets is (very probably) convection of an electrically conducting fluid under the influence of the Coriolis force. Convection is usually driven thermally, but even if it has a different physical origin, as in the Earth's interior, the theoretical treatment is formally very similar. Let us first discuss the properties of thermal convection in a rotating system ignoring the magnetic field. Here a particular goal is to understand the generation of differential rotation which, as we have seen before, is of crucial importance for the dynamo process. We write the basic equations in a coordinate frame rotating with angular velocity ft: dtp + V • pv = 0, (5.49)

5.3 MHD theory of thermal convection dt\ + v • Vv = - - Vp + g + 2v x a + iV(O x r) 2 - V • 7i, P dtT + v • VT = -(y - 1)T V • v + K V 2 T + -(y-D-^. (5.53) dz p dz Relation (5.53) can easily be understood. Consider a fluid element of volume dV which, originally located at height z, is slightly displaced upward to z + dz. It is decompressed to the ambient pressure at the new height p(z + dz) = p(z) + p'dz, pf = dp/dz < 0. Since heat conduction is neglected, the change of state dp, dp is adiabatic S

P = ±5l = ^dz.

(5.54)

p yp y p If the density of the displaced fluid element is smaller than the surrounding density p + dp < p{z + dz) = p + p'dz, (5.55)

there is a buoyancy force pushing the fluid element still further up, i.e., the fluid is convectively unstable. Insertion of (5.54) into (5.55) and use of the relation p'/p = (pip) + Tf/T gives the instability condition, (5.53). Since fluid velocities are usually much smaller than the sound speed, the fluid motions can be taken as incompressible. It is therefore useful to simplify (5.49)-(5.52) by eliminating the sound wave. In the Boussinesq approximation one assumes, in addition, that the density perturbation is small compared to the mean density and that the wavelength is short compared to the density scale-height kLp > 1, L" 1 = |Vlnpo|, such that po can be considered constant. We introduce the smallness parameter e, writing pv2 = ep, i.e., v = el/2cs, where cs is the sound speed. The quantities p, T vary on the convective time-scale, hence dt ~ e ^2, and we write po + ep, T = Tb + eT. By contrast, p relaxes on the sound time-scale and hence is in instantaneous equilibrium, such that p = (dpo/dT)pT

= -ap0T,

(5.56)

170

5 Dynamo theory

where in an ideal gas the thermal expansion coefficient is a = T~l. To zeroth order, (5.50) gives the barometric equation (the centrifugal force, which is small in most cases of interest, can formally be incorporated in the total pressure) Vpo = gpo, (5.57) relating the global profiles of po and po- In order e, (5.49)-(5.52) reduce to

dt\ + v • Vv =

V • v = 0,

(5.58)

Vp - (xgT + 2\ x H + vV2v,

(5.59)

Po dtT + v • VT = KV2T + 0,

(5.60)

where we have used the equilibrium relation (5.57) in (5.59) and omitted the tildes. Equations (5.58)—(5.60) are called the Boussinesq equations. A somewhat more general approach is provided by the anelastic approximation (see, e.g., Gilman & Glatzmaier, 1981), which still eliminates the sound wave but relaxes the scale requirement by allowing convective wavelengths to be of the order of the equilibrium scales. The essential change consists of replacing the incompressibility condition, (5.58), by Po *V • (pov) = V • v + v • Vlnpo = 0,

(5.61)

which is similar to the quasi-incompressibility condition, (3.48), of the perpendicular velocity in a strong curved magnetic field. The physics contained in the Boussinesq equations becomes clearer when we write the equations in dimensionless form by introducing the normalizations x/d -* x,

tx/d2 -* t,

T/AT

-+ T,

(5.62)

where d is the extent of the convection zone, i.e., the thickness of the convectively unstable shell, and AT is the change of the temperature across the shell:* dty + v • Vv = —Vp - Pr Ra Te g + Pr Ta1/2y xen + Pr V2v, V 2 T,

(5.63) (5.64)

neglecting the heat source. The dimensionless parameters in (5.63) are the Rayleigh number

Ra=C^H, KV

* For a thick shell, the radial variation of g oc r has to be included.

(5.65)

5.3 MHD theory of thermal convection

171

the Prandtl number Pr = - ,

(5.66)

K

and the Taylor number

-(¥)

2

The Rayleigh number is the ratio of the convective and the diffusive heat flux, while the Taylor number is the square of the ratio of the viscous diffusion time and the rotation period. (Instead of the Taylor number, one often uses the Ekman number E = Ta~1/2) Equation (5.63) tells us that for high Taylor number the convective flux oc v decreases with increasing Q and hence the critical Rayleigh number for onset of convection must also increase. Balancing gravity with the Coriolis force in (5.63) and inserting v into (5.64) we obtain the scaling of the critical Rayleigh number Rac~ Ta1/2kd, (5.68) where k is the wavenumber of a convective mode. In a rotating slab with il || g, where k is the horizontal wavenumber, the critical wavenumber is (Chandrasekhar, 1961) kcd oc Tal/\ (5.69) i.e., the wavelength of the unstable modes decreases with increasing Q, as the Coriolis force tends to disrupt larger eddies. Insertion into (5.68) gives the relation Rac ~ Ta2/\ (5.70) Convection is a very efficient process of heat transport. Increasing Ra above Rac only increases the vigor of the convective motion, while the temperature profile remains close to the adiabate. For mildly supercritical conditions (Ra — Rac)/Rac 1, where no steady or at least approximately regular convection pattern is to be expected. In fact, numerical simulations in the regime Ra/Rac ~ 102 and at relatively high Taylor numbers Ta < 109 (Sun & Schubert, 1995) exhibit a strongly time-dependent velocity field with a few major but rather irregular vertical columns concentrated near the inner shell boundary and embedded in a bath of small-scale turbulent motions. At fixed Ta the characteristics of the mean azimuthal flow depends strongly on Ra, as illustrated in fig. 5.11, which gives contours of the azimuthal velocity v^. While for Ra/Rac = 25 the pattern is essentially that found in the mildly supercritical regime with super-rotation at the equator and sub-rotation in the polar region, the Ra/Rac = 50 case shows a more complicated pattern with alternating regions of super- and sub-rotation on the surface, similar to the mildly supercritical but high-Prandtl number regime (Zhang, 1992). Because of this rather sensitive dependence on Ta, Ra/Rac and Pr, predictions for the differential rotation in the Sun or the planets are not yet possible, since Ra and Ta values achievable in present-day direct simulations are far below those prevailing in these astrophysical objects. In the present simulations no tendency toward a disc-like differential rotation pattern, as expected for the Sun (fig. 5.9), is recognizable. 5.3.2 The self-consistent MHD dynamo The problem of dynamo action by convection in a rotating sphere has attracted much interest, mostly in the geophysical community, because of its significance for the geodynamo. Before discussing the latter more specifically in section 5.3.3, we want to understand the basic mechanisms of magnetic field amplification by thermal convection and the reaction on the

174

5 Dynamo theory

(a)

(b)

Fig. 5.11. Contours of longitudinally averaged azimuthal velocity v^; solid contours represent super-rotation, dashed contours sub-rotation, for Ta = 109 and Pr = 1. (a) Ra/Rac = 25, (b) Ra/Rac = 50 (from Sun & Schubert, 1995). fluid motion. For large Taylor number, where inertia and viscous effects can be neglected, the equation of motion reduces to the magnetostrophic approximation 2p0v x O = Vp - - j x B, (5.74) c generalizing the geostrophic balance (5.71). From this equation, Taylor (1963) has derived an important constraint on the Lorentz force satisfied in a magnetostrophic state. The conducting fluid is assumed to be confined in a spherical vessel. We introduce cylindrical coordinates r, (j>, z with the z-axis being the axis of rotation, and take the integral of the (^-component of (5.74) over a cylindrical surface r = ro up the the points z = +zo, where the cylinder pierces the spherical boundary. It can easily be shown that the left side vanishes, 2po /

Jr=ro

(v x

= 0,

hence we find the condition = 0. r=r0

(5.75)

5.3 MHD theory of thermal convection

175

(The buoyancy force, omitted in (5.74), does not change this condition, since it has no ^-component.) A state satisfying (5.75) is called a Taylor state in dynamo theory.* The physical meaning of this constraint can be understood by returning to the full equation of motion. Although for small velocities the inertia effect is in general much weaker than the Coriolis force, there are special motions v^r), where (v x fl)^ vanishes. If (5.75) is not satisfied, large azimuthal accelerations occur. This will, in turn, lead to strong magnetic stresses through a radial field Br, until the system settles in a state where (5.75) holds approximately with some torsional oscillations about this state. Hence it can be expected that a fully developed dynamo is close to a Taylor state. To study the actual dynamics, even the simplest case of steady convection can only be studied numerically by three-dimensional direct simulation. The first global computations of this kind were performed by Gilman & Miller (1981), Gilman (1983) and Glatzmaier (1985a,b). Detailed analysis of the dynamo action near critical conditions is due to Zhang & Busse (1988, 1989, 1990; see also Wicht & Busse, 1997). These authors determine the critical Rayleigh number for onset of dynamo action, which depends on the Taylor number and the magnetic Prandtl number PrM, or the magnetic Reynolds number PrM = v/rj,

ReM = Lv/rj.

(5.76)

It is pointed out that the Lorentz force of a finite dynamo field counteracts the effect of the Coriolis force and hence tends to soften the constraint imposed by the Taylor-Proudman theorem, (5.72). The dynamo field may thus increase the convection velocity and even allow subcritical dynamo action, i.e., a finite magnetic field below the critical Rayleigh number for linear dynamo action. However, because of the symmetry constraints imposed in the computations by Zhang & Busse it is difficult to draw conclusions for more strongly supercritical dynamo behavior. Only recently have fully three-dimensional numerical simulations relevant to this regime been performed by St. Pierre (1993), Kageyama et al (1995), Kagayama & Sato (1997), Glatzmaier & Roberts (1995, 1996, 1997), Kuang & Bloxham (1997). While a discussion of the results by the latter two groups is postponed to section 5.3.3, here we want to follow more closely the work by Kageyama et al., which appears to illustrate the basic dynamo physics particularly clearly. (The fact that these authors solve the fully compressible MHD equations instead of making the usual * This is a generalization of the static force-free field (2.40), to which a plasma tends to relax under the constraint of constant magnetic helicity, see section 2.4. Such a state V x B = aB with a = const is also called a Taylor state.

176

5 Dynamo theory

Boussinesq or anelastic approximations is insignificant, since the resulting Mach numbers remain small and hence compressibility effects are weak. The same goes for the unusual boundary condition on the magnetic field, the vanishing of the tangential component, applied in these studies.) The Rayleigh number is only mildly supercritical, Ra = 2 x 104, for the Taylor number chosen, Ta ~ 6 x 106, giving rise to a coherent convection pattern, which consists of a certain number of pairs of Taylor-Proudman columns and convective heat transport smaller than conductive transport. The Prandtl number is Pr = 1. Dynamo action depends on the magnetic Prandtl number, which varies in the interval 7 < PrM < 280. Increasing PrM above the threshold value PTMC — 13, a magnetic seed field is amplified or, more specifically, a magnetic eigenmode grows exponentially with the growth rate y oc \n(PrM/PrMc\

(5.77)

i.e., the 77-dependence is rather weak. The magnetic energy EM grows up to saturation and subsequently exhibits irregular oscillations about a mean value EM. For the largest value of PYM used the saturation magnetic energy is much larger than the kinetic energy, EM/EK ~ 30. While the kinetic energy remains almost constant being rather insensitive to the dynamo field, the convection pattern is strongly affected, as soon as EM exceeds EK. (Though a reaction by the Lorentz force could be expected, if EM and EK are comparable, the rather sharp transition at EM/EK precisely unity seems to be fortuitous in view of the complex structures of the magnetic and velocity fields.) The transition is illustrated in fig. 5.12, which shows the time evolution of the magnetic energy for different magnetic Prandtl numbers, and in fig. 5.13, displaying the azimuthal drift of the (pairs of) convection columns. As long as EM < EK, the reaction is too weak to be visible, the drifts are steady and retrograde (westward) (note also the short period at the bottom of the figures, t < 50 (in units of the sound-crossing time), corresponding to the linear phase of the convective instability, where the drift is prograde, see the discussion in section 5.3.1). If EM exceeds £ x , the drift reverses direction and the columns behave more irregularly, (d), becoming seemingly turbulent for EM > EK, (e). Case (c) is particularly interesting. After reaching EK briefly at t ~ 1100, EM falls back temporarily allowing a regular convection state to re-establish itself, if only with a higher mode number, until, after a further transitional phase, the system settles finally in the nonlinear MHD regime with EM > EK. Let us discuss the dynamo mechanism illustrated in fig. 5.14. The toroidal field is pulled into the gap between cyclonic and anticyclonic columns, where it is convected up close to the latter and down close to the former, thus generating a poloidal field component (a). In the

5.3 MHD theory of thermal convection

_ """'

177

20 14

7

Kinetic energy (basis case)

-17

1000

2000

Time

Fig. 5.12. Time evolution of the magnetic energy for eight different values of the magnetic Prandtl number PrM = v/rj (the viscosity remains fixed, v = 2.8 x 10~3). The kinetic energy is nearly identical in all cases with that of the nonmagnetic case, the horizontal line (from Kageyama et a/., 1995).

180°

360° 0°

Longitude

Fig. 5.13. Azimuthal hodographs of the convection columns in the equatorial plane for (a) PrM = 20, (b) 23, (c) 28, (d) 35, (e) 280 (from Kageyama et al, 1995).

178

5 Dynamo theory

(a)

cyclone

anti-cyclone

v\

A

(c)

Fig. 5.14. Schematic drawing of the generation of the poloidal field (from Kageyama & Sato, 1997). southern hemisphere the toroidal field and the flow directions are generating a poloidal field of the same direction as in the hemisphere (b). Reconnection finally gives an overall poloidal which may subsequently be deformed again into the toroidal

reversed, northern field (c), direction

5.3 MHD theory of thermal convection

Poyorcgrl J U U

1 fl UL

-4

-3

Normal

n r~| r~

i

in n

III I I i

-2

1 |

-1

179

r

J i

6

f(10 years) Fig. 5.15. Record of geomagnetic reversals over the last 4.5 x 106 yr (from Cox, 1969). by differential rotation. We see that this coherent process functions as an aQ-dynamo, the process (a) representing the "cyclonic event" (Parker, 1955a), which produces the a-effect in mean-field theory (cf. fig. 5.7).

5.3.3 The geodynamo The geomagnetic field is certainly the best-known example of a magnetic field in an astrophysical body. In fact, it is more than just an object of scientific curiosity, because of its importance of protecting the Earth's surface from the bombardment by the solar wind and cosmic radiation. Of the inner planets, only the Earth carries a sizable magnetic field with a mean surface value of about 0.5 gauss. The field has a rather simple predominantly dipole structure with the magnetic moment oriented roughly along the Earth's rotation axis. Though we probably know more about the interior of stars than about that of our own planet, there is now a fairly detailed and reliable picture of the Earth's internal structure resulting both from seismology and materials science. The main constituents are: the iron core, consisting of a solid inner core with a radius of 1250 km and a liquid outer core extending up to a radius of Re ~ 3500 km; and the solid mineral mantle covered by a thin crust. Only the core is electrically conducting with a magnetic diffusivity X ~ 2m 2 /s, which would give rise to a free decay time of any current flowing in the core x ~ R*/k ~ 105 years. Since paleomagnetic studies clearly indicate that the Earth carries a magnetic field since its early history, this field cannot be of primordial origin. Moreover the magnetization of minerals, in particular on the sea floor, shows that the geomagnetic field frequently reversed sign. The abrupt reversals occur, apparently randomly in time, after a period of almost steady field of average duration 2 x 105 yr (see, e.g., Bullard, 1968; Cox, 1969), fig. 5.15. A qualitatively similar behavior may result from disc dynamo models, as shown in fig. 5.16 obtained from Rikitake's system of two coupled disc dynamos, see fig. 5.2(b).

180

5 Dynamo theory

Fig. 5.16. Time development of X(T) CC I\ in Rikitake's disc dynamo system (see fig. 5.2(b)), showing a nonperiodic sequence of reversals (from Cook & Roberts, 1970).

One concludes that there must be a dynamo process responsible for the geomagnetic field and its observed dynamics, which can only be located in the liquid outer core. The chemical composition of the outer-core material is that of an alloy consisting mainly of iron but mixed with some lighter elements such as Si, S and O. The engine driving the convection has its origin in the slow cooling of the Earth, which leads to (i) a freezing of the iron component onto the inner core with the corresponding release of latent heat at the bottom of the convection zone, and (ii) a melting of iron out of the mantle, hence a cooling due to negative latent heat at the top of the convection zone.* In addition to this thermal buoyancy source (Verhoogen, 1961) there is an even stronger gravity source caused by the buoyancy force on the light constituents of the alloy set free by the solidification of the iron (Braginskii, 1963). The fact that the outer-core density is roughly uniform and the sound speed is by many orders higher than the convection velocity justifies the Boussinesq approximation in a numerical treatment of the geodynamo, in contrast to the highly stratified compressible solar convection zone (section 5.2.2). However, direct simulations using the true values of the dissipation coefficients are still far out of reach, because of the extreme values of the Taylor number Ta ~ 1028 and the magnetic Prandtl number PVM ~ 10~5 (the Prandtl number itself is Pr ~ 1), though these values should not be taken too literally, since the molecular kinematic viscosity is only very poorly known (the usual assumption is v ~ 10~~5 m 2 /s). In order to perform numerical simulations with the presently achievable spatial resolution, one uses an eddy viscosity and eddy thermal diffusivity resulting from smallscale / convective motions dv, v ~ re ~ (dv)2/l ~ 103 m 2 /s, which brings the Taylor number down to Ta ~ 1012. * One thus finds the unfamiliar situation that the system is freezing at the hotter bottom and melting at the cooler top. The reason is that the melting temperature of the core material increases with pressure faster than the temperature increase along the adiabate, i.e., the temperature profile across the convection zone.

5.3 MHD theory of thermal convection

181

Using these assumptions, several numerical studies of the geodynamo have been performed by Glatzmaier & Roberts (1995, 1996, 1997). In these simulations the system is strongly supercritical Ra/Rac ~ 103, and the resulting convection is highly irregular. Field penetration into the inner core (assumed with the same magnetic diffusivity as the outer core) and the reaction on the dynamics in the outer core are included. The authors find generation of strong magnetic fields with energy EM ~ 103EK greatly exceeding the convective kinetic energy and a dominant dipole structure outside the convection zone with surface intensity of the correct order of magnitude. The inner core rotates faster (in an eastward direction in the rotating frame) than the mantle, in agreement with seismic observations, accelerated by the magnetic coupling to the high-latitude poleward convection along the bottom of the outer core (Glatzmaier & Roberts, 1996a). Concerning the question of field reversal, however, no coherent picture has appeared as yet. While in their first paper Glatzmaier & Roberts found that the poloidal fields in the inner and the outer core are essentially opposite, their interaction giving rise to an inherent tendency toward reversal, which actually occurred once during the computed time of 4 x 104 yr, in their later studies they deduce that the field direction was very stable (and parallel to the core field). A weak point in the numerical work of Glatzmaier & Roberts is the use of ad hoc hyperdiffusivities (viscous, magnetic, thermal) v(/) oc v/3, / = order of the spherical harmonic, which corresponds to adding a fifthorder derivative with respect to the colatitude 0.* The effect of such a choice is difficult to assess, in particular because of the strongly anisotropic character of this damping. At the least it leads to reduced effective values of the Taylor number. Kuang & Bloxham (1997) have pointed out the effect of the velocity boundary conditions. Since the viscous sublayers are far too thin in the Earth's core to be resolvable by present computational means, these authors assume stress-free boundary conditions, the continuity of the tangential components dn\t = 0, instead of the physical no-slip conditions v = 0, thus effectively reducing the layer width to zero. They argue that this choice is more adequate than unrealistically broad viscous layers due to the no-slip conditions assumed in the Glatzmaier-Roberts computations, which overemphasize the effect of the viscosity. Kuang & Bloxham find in their simulations a stable dynamo behavior with EM > EK. The configuration is close to a Taylor state, hence the effects of viscosity and inertia are indeed weak. The field distribution is markedly different from that * Hyperdiffusivities are often used in numerical turbulence studies (see section 5.4) in order to remove dissipation from the medium scales shifting it to the smallest possible scales.

182

5 Dynamo theory

found by Glatzmaier & Roberts, which they can essentially reproduce by applying no-slip boundary conditions. It thus appears that the difference in these models is primarily caused by the choice of boundary conditions and less by the different values of the nondimensional parameters such as the Taylor number. 5.4 MHD turbulence In many convectively unstable systems Rayleigh numbers are far above critical, resulting in fully developed magnefohydrodynamic turbulence, since, because of the dynamo effect, magnetic fields are usually excited to finite intensity, reacting on the fluid motion. MHD turbulence occurs almost inevitably in plasmas in motion and is therefore a widespread phenomenon in astrophysics, for instance in the solar system, where we encounter turbulent magnetic fields in the convection zone and in the solar wind, and in accretion discs, where MHD turbulence is responsible for the angular momentum transport. It is interesting to note that the presence of a magnetic field may even enhance the tendency toward turbulence, for instance in accretion discs, where, in the absence of magnetic fields, no turbulence would be generated (Balbus & Hawley, 1991). In this section I restrict consideration mainly to some recent developments, in particular in 3D turbulence, referring for a broader introduction to chapter 7 of my previous book (Biskamp, 1993a). 5.4.1 hom*ogeneous MHD turbulence The notion of hom*ogeneous turbulence is a theoretical abstraction or idealization. It refers to the quasi-hom*ogeneous properties of the turbulence in a small open subsystem L of the global inhom*ogeneous system away from boundary layers. The effect of turbulent structures of size exceeding L is treated as a hom*ogeneous background in the subsystem. The convenient representation of the fields in hom*ogeneous turbulence is the decomposition into Fourier modes

v(x) = Y, vk e*"x, B(x) = ^2 B k elk x,

(5.78)

where the Fourier components have the physical meaning of field intensities at the scale / = k~l. If the turbulence can also be assumed isotropic, the theoretical description and the interpretation of experimental and simulation results is further simplified considerably. However, whereas in hydrodynamic turbulence the assumption of isotropy is usually appropriate in a coordinate system moving with the mean flow, in an electrically

5.4 MHD turbulence

183

conducting fluid a mean magnetic field Bo has a profound effect on the turbulent fluctuations. While the dynamics perpendicular to Bo may develop small-scale structures, spatial variations along Bo are smooth because of the stiffness of the magnetic field. Hence in the presence of a strong Bo, viz., in a low-/? plasma, turbulence can be isotropic only in the perpendicular plane, i.e., the turbulence is quasi-two-dimensional, whereas fully isotropic turbulence can only be generated in a high-/? plasma with no mean field. The ideal invariants are characteristic quantities of a turbulence field. Here the quadratic invariants are most important, since these are particularly robust or "rugged" surviving finite truncations of the Fourier series, (5.78). In incompressible MHD there are three such invariants, which have been introduced in section 3.1.3: two mixed kinetic-magnetic quantities, the energy* i2

,

R

i2\

lV^/

k

_+i2

, i,,-|2\ _

_

FK

F

k

and the cross-helicity K = J2 Vk • B_ k = \ ^ k

(|z+| 2 - |ZJT|2) ;

(5.80)

k

and a purely magnetic quantity, the magnetic helicity B_ k .

(5.81)

Here z- = v + B are the Elsasser fields, (3.11), the natural variables in incompressible MHD. In 2D, where the MHD equations are identical with the simple reduced equations (3.65) and (3.66), we also have three quadratic invariants, again energy and cross-helicity, but the magnetic helicity is now replaced by the mean-square magnetic potential

^

= ElVk|2-

(5.82)

k

The spectral densities Ek, Kk, Hk or H% satisfy detailed balance relations, i.e., are conserved for each triad of modes k, p, q with k + p + q = 0. Hence, when injected in a certain /c-band, these quantities propagate or "cascade" in /c-space, either to large k, which is called direct or normal cascade, or to small k, called inverse cascade. The latter gives rise to formation of large-scale structures and is an example of the general phenomenon * Following conventions in turbulence theory we call the energy E in this section instead of W.

184

5 Dynamo theory

Table 5.1. Ideal invariants and cascade directions. HK is the kinetic helicity and Q the enstrophy.

Theory MHD

Navier-Stokes

Ek Kk Hk

Ef

3D direct direct inverse direct direct

Ek Kk

HZ E?

2D direct direct inverse inverse direct

of self-organization. The cascade directions can be determined from the absolute equilibrium distributions, i.e., the spectral equilibria of Ek etc. of a truncated ideal system of Fourier modes (see, e.g., Fyfe & Montgomery, 1976). While Ek, Kk have direct cascades, the magnetic quantities H^ H% are inversely cascading. Hence cascade properties are identical in 2D and 3D MHD, which reflects the basic similarity of 2D and 3D MHD turbulence, in contrast to hydrodynamic (Navier-Stokes) turbulence, where the inverse energy cascade in 2D results in the well-known formation of large eddies, a behavior which differs fundamentally from 3D turbulence, where no inverse cascade occurs. The cascade properties of MHD and Navier-Stokes theory are given in Table 5.1. It should, however, be noted that the formal identity of the cascades in 2D and 3D MHD suggests a somewhat oversimplified picture. Since H is an independent quantity, while HV is strongly coupled to £, Ef = -k2H%, MHD turbulence is expected to be more complex in 3D than in 2D, for instance the dynamo effect occurs only in 3D.

5.4.2 Selective decay and energy decay laws When dissipation effects are switched on, the ideal invariants in MHD turbulence decay, but not all at the same rate. While the energy decay rate is fast, independent of the values of resistivity and viscosity, the decay of K and H or H^ is much slower. This behavior, called selective decay (see, e.g., Ting et al., 1986), has a strong influence on the global turbulence dynamics. On a time-scale shorter than the decay time of K, for instance, K can be considered constant and the state to which the turbulence decays is determined by the variational principle

d[E-aK]

=dh

-a

B =0, /VXV-BI

(5.83)

5.4 MHD turbulence

185

where a is a Lagrange multiplier. Independent variation with respect to v and to B gives the relations v = aB and B = av, hence v = +B.

(5.84)

Since in such a completely aligned or Alfvenic state the dynamics is turned off, as directly seen in the z- formulation, (3.12), this is also the final state apart from the very slow laminar diffusion. Strongly aligned turbulence is observed in the solar wind, and numerical simulations, both in 2D and 3D, demonstrate the general tendency toward an aligned state. Note that alignment is a local property not connected with large-scale structures, which is consistent with the direct cascade of the cross-helicity. The degree of alignment is measured by the velocity-magnetic field correlation

p = K/E.

For weak initial correlation the tendency toward alignment is, however, superseded by a competing process of self-organization resulting from the effective constancy of H, or Hw in 2D. As shown in section 2.4, this leads to the variational principle d[E — txH] = 0, the solution of which is a static force-free field, j = aB. In 2D, the corresponding variational principle S[E — oc2Hw] = 0 has a similar linear solution j = a2xp. While smooth laminar fields decay exponentially, Ek ~ e~2/i/c f, the decay of turbulent fields shows a power-law behavior E ~ t~n, independent of the Reynolds number. The exponent n is a characteristic parameter, which reflects the similarity properties of the turbulence. In hydrodynamic turbulence the exponent is difficult to predict, and no universal decay law seems to exist, as will be briefly discussed below. In MHD turbulence, however, the decay process is govered by self-organization due to selective decay, giving rise to definite values of n. The basic theory was developed by Batchelor (1969) for the enstrophy decay in 2D hydrodynamic turbulence and applied later to MHD turbulence by Hatori (1984). Let us first discuss the 2D case, where Hw = const, while the energy E = EK +EM decreases. The macroscopic or integral scale fo of the turbulence can be defined by H v = EMll

(5.85)

In addition, numerical simulations show that in 2D MHD turbulence the ratio F = EK /EM remains constant,* i.e., the turbulence decays in a selfsimilar way. The integral scale should be insensitive to the details of the * There has been some discussion on this point in the literature. Certain simulation runs show a decrease of F, see, e.g., Kinney et a/., (1995), though the decay is much slower than that of the energy. Experience, however, suggests caution when drawing conclusions from individual runs, since one finds a significant scatter with some cases exhibiting quite anomalous features. Over the last decade we have performed many simulations of decaying 2D MHD turbulence with widely different parameters and initial states. The general impression is that F does remain nearly constant (an example is shown in fig. 5.17), quite differently from the 3D case.

186

5 Dynamo theory

turbulence depending only on the macroscopic quantities characteristic of the turbulent state, the energy E and the energy dissipation rate e = —dE/df (when restricting consideration to weak initial alignment, the cross-helicity does not enter). Hence, by simple dimensional arguments, lo satifies the scaling relation h ~ E3/2/e.

(5.86)

Combining (5.85) and (5.86) yields j ^

~~HV=

const,

(5.87)

and thus, making use of F = const, e=- ^

- £2,

(5.88)

which has the solution E ~ (t - to)"1 -> r 1

for

t>t0.

(5.89)

Figure 5.17 presents results from a computer run of freely decaying 2D MHD turbulence with a resolution of 20482 modes, choosing an initial scale of the turbulence k^1 = lo much smaller than the system size, which allows the inverse cascade of Hw to evolve freely. Figure 5.17(a) shows e/E2, the ratio of the l.h.s. and the r.h.s. of (5.88), which is seen to be nearly constant, thereby confirming the decay law, (5.88). Also shown is the energy ratio F, fig. 5.17(b), which is indeed practically constant, and the ratio of viscous to resistive dissipation A = e^/e11 = jdf w2d2x/rj J j2d2x, fig. 5.17(c). (Variation of the magnetic Prandtl number PrM = fi/r] suggests that A approaches unity for PrM —> 0.) The constancy of A over the entire time-span reflects the structure of the dissipative eddies, viz., current and vorticity micro-sheets, which is consistent with the observation that, at the dissipative scales, kinetic and magnetic energy spectra are nearly equal, ^k

— ^k

-

While, in 2D, finite EM entails finite ffv, in 3D the magnetic helicity may take any value in the interval 0 < \H\ < H max , where Hm8LX is the maximum magnetic helicity for a given energy spectrum Ejf. Selective decay is only effective for H =fc 0. Since MHD turbulence is usually generated in rotating systems, where Coriolis and buoyancy forces naturally lead to twisted field lines, finite H should be typical. In this case we can define the integral scale lo by EMl0 = H = const. (5.90) * This is the conventional terminology. Strictly speaking, the energy decay rate is

-E~ldE/dt.

187

5.4 MHD turbulence =

20

40

60

80

20

100

t

40

EK/EM

60 t

80

100

**AAtA* 2 (in units of the large-eddy turnover time), when the turbulence has become fully developed, and their small scatter shows that the relation (5.91) is well satisfied. From the same runs we also see the rapid decay of F, fig. 5.18(b), in particular we find T(t) oc E(i). Inserting this behavior in (5.91) gives the simple relation for the energy decay in 3D turbulence, valid asymptotically when F < 1,

Tt

(5.92)

5.4 MHD turbulence

189

whence the asymptotic decay laws E ~ EM ~ r 1 / 2 ,

EK ~ r 1 .

(5.93)

Both decay laws are consistent with the numerical observations. The fact that the kinetic energy decays faster than the magnetic energy is consistent with the selective decay to a static force-free state. The decay of nonhelical MHD turbulence H ~ 0 is not governed by a selective decay process, which makes the situation more complicated, similar to that in (3D) hydrodynamic turbulence. There the classical approach starts from the Loitsianskii integral ££ = /0°° dll4(v(x + l)v(x))9 which is assumed to be constant (see, e.g., Monin & Yaglom, 1975, Vol. II). Writing S£ = EIQ and assuming again IQ ~ £ 3 / / 2 /e, one immediately obtains -dE/dt ~ E11/m with the asymptotic solution E ~ t~10/7. Since S£ can also be written in the form S£ = / dk k~5Ek, where Ek is the energy spectrum, S£ is determined by the spectrum at the largest scales, which in general depends on the way of turbulence generation and hence cannot be universal. Moreover, the invariance of the Loitsianskii integral has been questioned (see, for instance, Frisch, 1995). For MHD turbulence a similar quantity has been considered, J$?MHD = f dl I4(z±(x + l)z±(x)) (Galtier et a/., 1997). However, here, too, the invariance is questionable. Both closure theory (Lesieur & Schertzer, 1978) and direct numerical simulations (see, e.g., Hossain et a/., 1995) indicate E ~ t~l. It is interesting to note that the transition from E ~ r 4 to E ~ t~1//2 is found to occur at relatively small values of the magnetic helicity, H/Hm2iX ~ 0.2-0.3, hence the latter decay law should be more typical. 5.43 Spatial scaling properties Theoretical concepts Turbulence is characterized by a broad range of spatial scales exhibiting certain scaling properties. The scaling range, or inertial range, is defined by k > I > Id,

(5.94)

or in /c-space ko I > kn and kn> I > Id- Here we only discuss recent results concerning the direct cascade. The scaling properties are conveniently described by the statistics of the increments of the turbulent fields, dv\ = v(x + /) — v(x) etc., which for hom*ogeneous turbulence is independent of x. (The vector character of v, x, 1 is not important at this point.) For stationary conditions the energy decay rate e equals the energy injection rate em, which is also equal to the energy transfer rate et across the inertial range and to the energy dissipation rate e are constant. The scaling properties in 3D MHD turbulence have recently been studied by numerical simulations of decaying turbulence with resolutions of up to 5123 modes (Miiller & Biskamp, 2000). Figure 5.20 shows a scatter plot of the normalized energy spectrum compensated by k5^ obtained from an extended period of fully developed turbulence."'' The * The bottleneck effect designates the hump in the energy spectrum observed at the transition from the inertial to the dissipation range, where the energy flux from the former encounters some kind of a soft barrier, a "bottleneck", which the energy must cross before being dissipated. Mathematically the effect is caused by an overshoot in the Fourier transform resulting from the rather abrupt dissipative fall-off of the fluctuation amplitudes. While in 3D hydrodynamic turbulence this hump has a finite spectral extent, which is independent of the Reynolds number (see, e.g., Lohse & Miiller-Groeling, 1995), in 2D MHD it is found to affect an increasingly larger part of the inertial range. t In decaying turbulence one cannot simply average over the spectrum at different times because of the secular change of the system. Hence the spectrum must be normalized to eliminate the time variation of the global quantities, Ek(t) —> E(k), where k is normalized to the dissipation scale ld9 k = kid. In the Kolmogorov case ld = lK = (/j3/6)1/4> the Kolmogorov length, and normalization of the energy spectrum (5.98) results from simple dimensional analysis: E(k) = Ek/{eii5)i/4 = CF{k)k~5/\ (5.104)

194

5 Dynamo theory

1.00 -

0.10 r-

0.01 0.01

1.00

Fig. 5.20. Scatter plot of the normalized energy spectrum compensated by fc5/3 from a simulation run with 5123 modes, H/HmSiX = 0.6, K = 0, using normal diffusion. The dashed line indicates the IK spectrum ocfc~3/2,the dotted line the spectrum Ck~5/3 with C = 1.1 (after Muller & Biskamp, 2000). spectrum exhibits a clear scaling behavior over almost one decade, which is close to k~5/3, definitely steeper than the IK spectrum k~3/2 indicated by the dashed line. The dotted line corresponds to Ck~5^ with C = 1.7. In contrast to the energy spectrum in fig. 5.20, the structure functions from this simulation run do not exhibit a clear scaling range. (It is generally observed, both in experiments and in numerical simulations, that the structure functions have a shorter scaling range than the energy spectrum, in particular for higher orders, p > 2.) However, using hyperdiffusion, i.e., replacing the operator rjV2 by — ^(V 2 ) 2 and similarly for the viscosity, concentrates dissipation more strongly at the smallest scales and results in a longer inertial range. In fig. 5.21 the normalized time-averaged structure functions S+(T) = (\dzf\P)/Ep/2J= l/ld, are plotted for p = 3 and 4. We where F is a dimensionless function, F(0) = 1, which describes the spectrum in the dissipation range k ^ 1, and C is the Kolmogorov constant. Though the energy spectrum in 3D MHD seems to follow a /c~5/3 law, the standard normalization, (5.104), must be generalized to account for the dependence on H. A suitable form is (see Muller & Biskamp, 2000) E{k) = Ek/ \{efi5)l/4(v5A/eH)3/s]

= CF(k)k-5/\

(5.105)

where the dissipation scale is slightly changed compared with the Kolmogorov case,

5.4 MHD turbulence

195

10 3

:i

i

i

i

| i i i i

i

i

i

X^

10°: - (b)

\ Slope • 1 1 1 1 ' i'

\

~z-

" \ i

~~-

1111111

i

I\I

111-

100 1

10

j

1O- U

-

o.6

-

/

-

: 3

y

_

~^y

10 -2

r: / /I\\

1 1

10 2

1 1

1

1

1

1

1

1

10 3

Fig. 5.21. Log-log plots of the normalized structure functions S+(l) for (a) p = 3, (b) p = 4, from a simulation run with hyperdiffusion but otherwise the same parameters as in fig. 5.20. The insets show the logarithmic derivative, the horizontal dashed lines indicating the most probable values of the scaling exponents (from Miiller & Biskamp, 2000).

find in particular that £3 is close to the Kolmogorov value £3 = 1, a result which has recently also been derived analytically for certain triple products of bz\ in MHD turbulence (Politano & Pouquet, 1998a,b). Assuming £3 = 1 to hold exactly, allows us to obtain rather accurate values for the remaining scaling exponents by using the property of extended self-similarity (Benzi et al, 1993). Here the structure functions Sp are plotted not as functions of / but of S3, / = /(S3), which yields a significantly extended scaling range. Figure 5.22, summarizing the results, gives Cp for MHD turbulence in 3D (diamonds) and 2D (triangles), compared with the She-Leveque model £^L for hydrodynamic turbulence (continuous line) and the £pK, (5.103) (dotted line). (Note that the 2D values are only relative, Cp/C3> since here (3 is smaller than unity, as mentioned above.) The figure shows that the scaling exponents for 3D MHD

196

5 Dynamo theory

fa.

0.5 0.0

Fig. 5.22. Scaling exponents C+ for 3D MHD turbulence (diamonds) and relative exponents C^/C^ for 2D MHD turbulence (triangles). The continuous line is the She-Leveque model £^L, the dashed line the modified She-Leveque model ^ which fits the 3D MHD results very well, and the dotted line is the IK model (from Muller & Biskamp, 2000). turbulence lie between those for 2D MHD and those for hydrodynamic turbulence, hence MHD turbulence is less intermittent in 3D than in 2D, but more intermittent than hydrodynamic turbulence. We also see that the simulation results differ basically from the IK prediction £*K. Let us now see, how the numerically observed values for 3D MHD turbulence fit into the framework of the generalized She-Leveque formalism (5.102). The numerical results indicate that the fundamental scaling is Kolmogorov-like, i.e., g = 3. Visualizations of the small-scale turbulent structures clearly demonstrate the sheet-like character of the dissipative eddies, i.e., Co = 1. We thus obtain the intermittency model for 3D MHD turbulence (5.106) which in fact agrees well with the numerically observed scaling exponents, as seen in fig. 5.22. These findings change the conventional interpretation of the turbulence dynamics in MHD. Although the Alfven-wave effect is present, as clearly seen from the fact that the kinetic and magnetic energy spectra come close at small scales, the interaction of Alfven waves propagating along the local field does not seem to govern the dynamics, which is dominated instead by cross-field eddy-like motions, much as in hydrodynamic turbulence. Thus the Alfven-wave interaction, which is basically a weak turbulence process,

5.4 MHD turbulence

197

is only a subdominant effect, which taken separately would lead to a steeper energy spectrum Ek ~ k~2 (Galtier et a/., 1999). The IroshnikovKraichnan mechanism giving rise to the fc~3//2 spectrum does not seem to be active in a physical system. It is interesting to note the difference between 2D and fully 3D MHD turbulence. Though the asymptotic behavior in 2D is difficult to assess, as mentioned above, at least in the Reynolds number range considered in the 3D studies the 2D results seem to be consistent with the IK phenomenology, showing in particular a clear fc~3/2 energy spectrum (see, e.g., Biskamp & Welter, 1989). The difference between the 2D and the 3D behavior is not implausible. There are two competing dynamical effects in MHD turbulence - eddy motions and Alfven waves. Since in 2D hydrodynamic turbulence the eddy dynamics is weak, resulting in the steep k~3 energy spectrum, Alfven waves are likely to play an important role. In 3D, by contrast, eddy dynamics is much more vigorous and may thus dominate over the Alfven-wave effect. 5.4.4 hom*ogenous turbulent dynamo The statistical properties just described refer to turbulent systems in which the magnetic energy is of the same order as or larger than the kinetic energy, assuming, in particular, that the magnetic field is smooth, but its derivative, the current density, is distributed sparsely or intermittently. An example in nature is the turbulence in the solar wind, where the magnetic field has simply expanded from the space-filling quasi-static field in the corona. Conditions are, however, quite different for systems in which the turbulence is primarily kinetic, stirred by mechanical or thermal forces. Here the magnetic field is convected with the fluid and possibly amplified by dynamo action. Such a field is not space-filling but tends to be concentrated in thin flux tubes, similar to vortex tubes in hydrodynamic turbulence. This behavior can be understood from the structure of the equations. If the magnetic field is weak, the Lorentz force is negligible, such that vorticity and magnetic field follow the equations 5fB + v V B dtW + v V w

= =

B-Vv + ^V2B, w-Vv + /iV2w,

(5.107) (5.108)

which are formally identical for magnetic Prandtl number PrM = n/r\ = 1. Hence a similar structure of the solutions can be expected. For the inertial range, the analogy between w and B implies the magnetic energy spectrum Ef oc k2Ef ~ kl/\

(5.109)

assuming a Kolmogorov kinetic energy spectrum. We thus see that increases with /c, up to the dissipative cutoff, which implies that the mag-

198

5 Dynamo theory

netic field is concentrated in small-scale structures. The argument, (5.109), leading the the spectrum E^f ~ k1^ should, however, be considered with some caution. Since the magnetic energy is not an ideal invariant, the spectrum cannot be inferred from a cascade argument but is only based on the analogy with the vorticity. Nor does the analogy give information about the magnitude of EM. Will a weak seed field decay or grow in time? Starting from a smooth field distribution, the magnetic energy will always increase initially due to field line stretching, twisting and folding, as idealized in the VainshteinZeldovich process, where the magnetic field is rapidly concentrated in rope-like structures of decreasing diameter d, until the thinning is limited by resistive diffusion. Whether, subsequently, the field energy continues to grow or decays depends on the balance between convective driving and diffusion, which is measured by the magnetic Reynolds number ReM or, for given kinetic Reynolds number, on the magnetic Prandtl number, ReM — PrMRe. Only if ReM > RCM,C (or P^M > P^M,C) does true dynamo action occur, where ReM,c ~ 20-40 typically. Saturation is caused by the Lorentz force, but no clear answer has yet been given about the final magnetic energy level and the spectral distribution. Three-dimensional numerical simulations (Kida et a/., 1991) indicate that for the most interesting case of PYM > 1 the magnetic spectrum is essentially flat at small k where Ef > E$*9 while for large k where £ f < Ejj*9 the magnetic field is passively convected by the largescale motions, hence E^f ~ k~l up to the resistive cutoff. In this regime the energy dissipation is mainly ohmic, in spite of the fact that r\ < /i.

6 Noncollisional reconnection processes

We have seen in the previous chapters that, in the framework of resistive MHD, reconnection rates are in general low (~ rj1/2) becoming independent of r\ only when the system is fully turbulent. The resistive MHD approximation is, however, only valid for relatively dense and cool plasmas. In hot plasmas, which are of special interest in laboratory fusion devices and in space, resistivity is no longer the most important nonideal effect in Ohm's law. Instead R in (2.1) is dominated by noncollisional terms, which are associated with different intrinsic microscopic lengthscales: the ion and electron inertial lengths C/CDVJ, j = i,e; the ion and electron Larmor radii pj = vtj/ilji and the Debye length Xv = vte/cope. Here the plasma frequencies coPj, the gyrofrequencies |Q7, and the thermal velocities vtj are defined by ojpj = J4nne2/nij,

Qj = eB/mjC,

vtj = JkBTj/nij.

(6.1)

In a low-/? plasma we have typically the ordering c/copi > c/cOpe •> pi »

pe •> XD.

(6.2)

The particle density n is, properly speaking, the electron density ne or the ion density nu respectively. However, deviations from charge neutrality ne ^ Hi occur only on spatial scales of order Xv. Since most plasmas of interest in the present context are sufficiently dense, such that XD is small and co~el short compared to all other spatial and time-scales, we have (rii — ne)/ne < 1, so that we can set ne = rii = n. This approximation is called quasi-neutrality, since the net electric charge density is small but not identically zero. There is still an electrostatic field, which is, however, no longer calculated from Poisson's equation but from the quasi-neutrality condition V • j = 0. Only in chapter 7 will we discuss certain electrostatic microinstabilities, which involve fluctuations on the Debye scale. Effects 199

200

6 Noncollisional reconnection processes

associated with the intrinsic scales in (6.2) become more important than resistivity, when the resistive scale-length 8n ~ (T/T^) 1 / 2 , the width of a resistive current sheet, becomes smaller than one of these scale-lengths. For instance, electron inertia dominates over resistive diffusion in a current sheet, if /

'

A l / 2

which is independent of the collision time. Though ^ is much larger than ji±, this is compensated by the difference in the derivatives. To estimate the ratio of the ion viscous effects we use V|| ~ L7 1 , V_L ~ Lj1 fi± for QT > 1.* Neglecting the smaller mixed parallel-perpendicular terms nXZ9 the stress tensor contribution in the ion momentum equation (6.7) assumes the form: —(V • TT)_L = nmijiii±V2\i± — nm^ob x V2v;. (6.18) In the electron momentum equation (6.5) the stress tensor is usually neglected, since it is proportional to the electron mass. Though the parallel viscosity could still be formally larger than the friction term T||, the latter is more important since it describes the momentum exchange between electrons and ions, i.e., resistivity, instead of a mere relaxation of the parallel electron momentum gradient. However the perpendicular electron * It is worth mentioning that the expression p?Qj(^: p2e€le) ^ cT/eB has a certain reputation or notoriety as the Bohm diffusion coefficient (originally including a factor 1/16) in the theory of plasma transport.

204

6 Noncollisional reconnection processes

viscosity fie± is often kept, since it may be effectively enhanced compared with the negligibly small collisional value p2e ve in the presence of braided magnetic field lines (Furth et a/., 1973; Kaw et a/., 1979). The electron momentum equation (6.5) is usually written in the form of Ohm's law. With the approximations on ne and T and substituting ne\e = ne\t —j, the generalized Ohm's law reads E = —v,- x B + — j x B c nee

Vpe ne

-^(dtme + V • meye) + m + ^ ? V i v (6.19) ne e Compared with the simple form of Ohm's law in resistive MHD, E = —(l/c)v x B + jyj, with v = v,-, several additional terms appear in twofluid theory: the Hall term ( x j x B ) ; the electron pressure gradient; the electron inertia; and electron viscosity, which are nondissipative except the last one. Electron inertia and viscosity are only important in cases where the electron flux me is large, which occurs on scales small compared to the smallest ion scales, hence we can substitute \e ~ —]/ne in these terms. To illucidate the order of magnitude of the non-MHD effects let us write (6.19) in non-dimensional form using the normalisations x/L —• X,V/VA -> v,B/Bo —• B,n/no -> n,cE/vABo -> E,pe/(Bl/A%) -> peJ/(cBo/4nL) -+ j : E = - v ( x B + -i(j x B - VPe) + -£ J + ,,j -

2

mSJ

\,

(6.20)

where di = -^j,

de = -?-,

(6.21)

and r\2 is also called hyperresistivity. While the dissipative terms are known to cause reconnection, this is, in general, not obvious for the nondissipative effects. For instance the Hall term itself conserves field line topology, see section 2.2, the field being simply convected with the electron flow ye, but it may drastically enhance the efficiency of a weak dissipative reconnection process as will be discussed in detail in the following section. The electron pressure contribution in Ohm's law is usually small. Inserting E in Faraday's law gives a term x V n x Vpe, which is only important in a dense plasma, while it vanishes for a polytropic gas law p = p(n). (Here the underlying assumption is that collisions, though rare, are still frequent enough to make the pressure tensor isotropic.) The electron pressure

6.1 Two-fluid theory

205

gradient will, however, play an crucial role through the diamagnetic drift in a low-/? plasma. Substituting E in the ion momentum equation, neglecting electron inertia, and using the approximation (6.18) and ne(\i — \e) = j , gives nmi(dt\i + v; • Vv/) = - V f e + pe) + - j x B + nmKmVivi -

M

b x V2v,).

(6.23)

Since magnetic reconnection processes are rather insensitive to the details of the evolution of the temperatures T,-, Te, it is sufficient to treat the transport terms on the r.h.s. of (6.7) and (6.8) in a crude way. We neglect all heat fluxes except for the parallel electron one qe\\, which we take to be infinite such that V\\Te = 0. Also, the terms involving the stress tensors are neglected. Qe consists mainly of the Ohmic heating, while Qt contains the viscous dissipation. Summarizing our two-fluid model written in non-dimensional form with the notation v = v,-,/i = /i_u,\e = v — di)/n, we have: dty + v • Vv = —(V(pt + pe) - j x B) + fN2±\ - juob x v]y9

(6.24) (6.25) (6.26)

(dtPi + v VPi) +

- v = -MV±v) 2 .

(6.27)

In addition, the continuity equation for n and Faraday's law have to be satisfied. This set of equations conserves energy in the form (neglecting dissipation and boundary effects): W = { I [ nvl + Bl + -Z-j1 + 3(Pi + pe) dV = const. -./(•

(6.28)

The energy expression differs from the (compressible) MHD expression, (3.32), only by the electron inertia term. What happens to the other ideal MHD invariants in the two-fluid framework, most notably the magnetic helicity H = J A • BdVl We find up to resistive and boundary effects

f ~*e/B-i**K-OK*

(629)

Hence also in two-fluid theory magnetic helicity is conserved except for the usually very weak electron inertia effect, since the contribution from

206

6 Noncollisional reconnection processes

the electron pressure term in (6.25) vanishes on the assumption V\\Te = 0. On very small scales pe, which are dominated by electron inertia, we introduce, in section 6.2.1, a generalized helicity expression, which is conserved in the framework of electron magnetohydrodynamics. 6.2 High-/? whistler-mediated reconnection We first consider noncollisional reconnection in high-/? plasmas, /? ^ 1, which is typical for plasma parameters in the solar wind and the geomagnetic tail. In this case the two-fluid equations can be simplified by assuming the ion motion to be incompressible. (A rigorous justification of this approximation requires j8 > 1, such that the sound wave is much faster than, and hence decouples from, the Alfven wave governing the reconnection time-scale. However, compressibility effects are found to be weak already at moderately high fi •> 1). If V • v,- = 0, one can choose uniform density, whence the quasi-neutrality condition V • j = 0 requires that also the electrons are incompressible V • \e = 0. We can now eliminate the pressure by taking the curl of (6.24), which leaves us with two equations for the vorticity w = V x v and the magnetic field B, generalizing the incompressible MHD equations (3.9), (3.10): 3tw + v • Vw - w • Vv = B • Vj - j • VB - V x (V • TT,-), dtB + v VB - B • Vv = -dt(B

(6.30)

Vj - j VB)

(6.31) where the last term in (6.30) comprises the viscous effects. While the w equation remains essentially unchanged, there are two new non-ideal terms in the induction equation (6.31), originating from the Hall effect and electron inertia, which are associated with the scales d\ and de, respectively. To study the linear mode characteristics, we linearize these equations about a hom*ogeneous static equilibrium neglecting dissipation and making the usual Fourier ansatz for the perturbations B = Bi exp{—icot + z'k • x} etc., +k\\B Bi = 0 , (6.32) eo(l +k2d2e)Bi +/c|Bvi = idtk^Bk x B1?

(6.33)

k • vi = k • Bi = 0, which gives the dispersion relation (writing v\ instead of B2 for clarity)

[co 2 (l + d2ek2) - k2v2A]2 - co2d2k2k2v2A

= 0.

(6.34)

6.2 High-fi whistler-mediated reconnection

207

While for long wavelength kdi < 1 one recovers the shear Alfven wave co2 = kh>\, in the short-wavelength regime kdi > 1 the mode is dispersive djk2kh2A

,

called whistler.* In this regime the frequency no longer depends on the ion mass, since vAc/coPi = QiC2/cojt = Qec2/coje = cB/Anne.

The dispersion relation (6.35) has two different regimes depending on the magnitude of dek. Written in dimensional form, we have 2 C co ~ Q e —^-/CII/C, (Ope

CO ' CO —— < k < ——, c

c

(6.36)

the whistler proper, and

co ~ o A

k>^ ,

(6.37)

k c the electron cyclotron wave. The coupling of the two components of Bi (and of vi) by the cross-product on the r.h.s. of (6.33) makes the whistler circularly polarized, ^T- = +1,

(6.38)

right-hand for wave propagation in the direction of Bo, left-hand for propagation in the opposite direction. The sense of rotation is that of electron gyration. 6.2.1 The EMHD approximation While on large scales / > c/copt ions and electrons move essentially together, \t ^ \e = O(VA), V,- — \e = ]/ne = O(vAc/coPil), at small scales / < c/copt the ions are too heavy to follow the electrons, hence v,- < v g ~ —\/ne. Neglecting the ion velocity in this regime, the system is determined by the induction equation (6.31) depending only on B and j = V x B. With * The name results from the properties of such disturbances in the reception of radio waves. Whistlers are audio-frequency electromagnetic bursts originating from lightning flashes, which propagate along the ionospheric magnetic field. Since the group velocity increases with increasing frequency dco/dk oc ^Jco, (6.36), the received signal exhibits a descending audio tone.

208

6 Noncollisional reconnection processes

the usual normalization to the Alfven time %A — L/VA the electron inertia term reads x |

= d2e V x (da - dt] • Vj) = -d] [dtV2B - dtV x (j x V 2 B)]. (6.39)

If, instead of TA, we choose as unit time the characteristic time of the whistler %w = IAI&U which can also be written as TW = L2(opi/(cvA) = {Oedlr1

(6.40)

(note that xw is independent of both m,- and me\ (6.31) assumes the form dt{B - d2eV2B) + V x [j x (B - d2V2B)] = r\V2B - */2(V2)2B.

(6.41)

This equation is called electron magnetohydrodynamics (EMHD), since the electrons play a role similar to that played by the ions in ordinary MHD.* EMHD is simpler than MHD, consisting of only one vector equation, but it contains an intrinsic scale-length de. There are two ideal invariants in EMHD, the energy W = i f(B2 + d2ej2)d3x

(6.42)

and a generalized helicity

He = y ( A - d2j) • (B - d2eV2B)d3x.

(6.43)

6.2.2 Properties of the reconnection region As in the resistive case, we consider reconnection processes in the twodimensional approximation dz = 0, where magnetic field and current density can be written in the form B = ez x V\p + ez(Bz0 + b),

(6.44)

j = Vb x ez + e z V V

(6.45)

splitting the axial field into a mean field Bzo and a fluctuating part b. Reconnection occurs in a region of small B_L, which for stationary * EMHD has previously been discussed mostly in the context of fast plasma processes such as plasma switches and Z-pinches, see the reviews by Kingsep et al., (1990) and Gordeev et al, (1994). Since in such processes the density is strongly inhom*ogeneous, the constantdensity approximation implicit in (6.41) assuming \e = —i/en0 is not justified, as the term V(l/n) x (j x B) can play an important role. In the context of magnetic reconnection, however, electron inertia is usually more important than density inhom*ogeneity, hence (6.41) is appropriate, see, e.g., Bulanov et al, (1992).

6.2 High-fi whistler-mediated reconnection

209

conditions has the form of an X-point configuration. In the region \x\ < d\ around the X-point, where the ion velocity can be neglected, the system is described by the EMHD equation (6.41), which in 2D reads dt(xp - d2ej) + ye - V(xp - d2ej) = ij - rj2V2j, dt(b - d2ewe) + \e • V(b - d2ewe) + B ; V/ = f]we - >?2V2w,,

(6.46) (6.47)

\e = - j ± = ez x Vb, ; = jz = V2t/), we = W2b, hence the axial field fluctuation b plays the role of an electron stream function, which we therefore denote by (j)e in the following, b = (f)e. Equation (6.46) expresses the conservation of F = xp — d2ej, the electron canonical momentum in the direction of the ignorable coordinate, cf. (7.7). We also assume stationarity, which is valid as long as the flow remains stable. Instability and resulting turbulence will be treated in section 6.2.4. As in section 3.3.3, the coordinate system is chosen such that +x defines the inflow and ±y the outflow directions. For \x\ > de (6.46) and (6.47) reduce to E + \e • Vxp = 0, (6.48) B • W2xp = 0,

(6.49)

where E = dt\p is the reconnection rate. The equations have the similarity solution (6.50)

-fin

x + ay

(6.51) 2 x - ay which is identical with the corresponding MHD solution, (3.142) and (3.143), illustrated in fig. 3.15. Finite resistivity or hyperresistivity is only needed to smooth the solution on the separatrix x = +y. The scale parameter a allows a finite uniform current density jo = 1 — a2, so that the separatrix branches may intersect at any angle. We now show (Biskamp et a/., 1997) that the similarity solution valid for \x\ > de can be matched to the solution in the electron inertia-dominated layer, or simply the inertial layer, |x| < de. In the limit of small resistivity, the current sheet which forms in this region exhibits a complicated multi-scale structure, in particular the current density develops a cusplike singularity, resulting from continuous acceleration at the X-point, the stagnation point of the perpendicular electron flow. However, the singularity does not give a finite contribution to the integrated current (see also the discussion in section 6.3.3). Because of the cusp-like current profile the hyperresistivity r\i associated with electron viscosity can be a more efficient damping and dissipation mechanism than resistivity even for very

210

6 Noncollisional reconnection processes

Fig. 6.1. Integration domain in (6.53). small values of r\i. In the following order-of-magnitude discussion these substructures are irrelevant and will be neglected. The inertial layer consists of a current sheet of width b and length A. From the convection term in (6.46), it is seen that inertia becomes important if dxxp = By~ d2edxj ~ d2eBy/S\ hence the layer width is 5 ~ de.

(6.52)

Integrating (6.47), d]ye • Vwe = B • V/,

(6.53)

over a quadrant of the current layer and using Gauss' theorem to transform the surface integrals into line integrals over the boundary we have venwe dl = meu3A ~ d5j\ whereas the energy flux out of the layer is mainly kinetic m v3S ~ d ^> vB2S ~ d5^

using again that in our units me = d2, Bx ~ A and By ~ jde. We also see that the energy fluxes into and out of the layer are of the same order, which differs from the behavior in resistive MHD, (3.156), where the main part of the energy flux into the layer is dissipated. Relation (6.59) implies that A -> 0 for de -> 0. Since By - dl/3 -> 0, there is no flux pile-up in front of the layer. The reconnection rate E is independent of the small parameter de. In the region \x\ < d\9 £ is a free parameter of the electron flow, determined only by the free energy of the global configuration. This behavior differs markedly from the scaling of resistive current-sheet reconnection. It is true that outside the current layer a solution equivalent to (6.50) and (6.51) is permitted also in resistive MHD, but it cannot be matched to the diffusion region, which has macroscopic length A ~ L. The basic difference is that in resistive MHD the outflow (= downstream) velocity equals the upstream Alfven velocity, which remains finite, such that the reconnection rate is slow, u ~ VAS/A < VA, depending strongly on the small parameter r\9E ~ ?/1//2.

213

6.2 High-ft whistler-mediated reconnection

t=u

Fig. 6.3. Coalescence of two flux bundles in EMHD: contours of xp and (f)e at two different times, t\, t2\ £2 > h.

6.2.3 Coalescence of EMHD flux bundles The EMHD scaling predictions have been confirmed by a series of numerical simulations (Biskamp et ah, 1995). We consider the coalescence of two flux bundles located on the diagonal in a square box of edge size L = 2n, xp(x, y9t = O)=^2

exp{-r, 4 /4},

(6.62)

7=1,2

where rj = (x-Xj)2+(y-yj)\ *i = }>i = (7r/2)+0.6,x2 = yi = (3TT/2)-0.6. Figure 6.3 gives contour plots of xp and <j)e at two times during the reconnection process. The conspicuous feature is that the flux surfaces seem to be pulled into the reconnection region instead of being pushed against it as in the MHD case. This property is due to the flow pattern with streamlines converging (i.e., velocity increasing) toward the X-point, which agrees with the similarity solution (6.51). In resistive reconnection

6 Noncollisional reconnection processes

214

OJLF

0.7/.

0.3/. 0.3L

0.7L

0.3L 0.3L

0.7/.

Fig. 6.4. Reconnection region of EMHD flux bundle coalescence: (a) xp; (b) cj)e for de = 0.06; and (c) for de = 0.03. xp is essentially identical in both cases (from Biskamp et a/, 1995).

the flow is rather uniform, even diverging in front of the macroscopic current sheet, see fig. 3.16. The structure of the region around the X-point and its dependence on de are illustrated in fig. 6.4. The inertial layer appears as a small layer of high electron velocity (note that the strong flows along the separatrix are a consequence of the similarity solution (6.51)), but does not show up in the flux function xp, which exhibits a structureless X-point behavior. The layer shrinks both in length and width with decreasing de, which is found to be in quantitative agreement with the scaling laws derived above. A further blow-up of the reconnection region is shown in fig. 6.2, which demonstrates an important feature of noncollisional electron dynamics. Equation (6.46) implies that F = xp — d2ej is convected with the electron fluid. Since symmetry requires a stagnation point of the flow at the Xpoint, F would simply pile up in front of the X-point and would remain

6.2 High-fi whistler-mediated reconnection

215

topologically invariant in the absence of dissipation. This pile-up of F is clearly seen in fig. 6.2(b) making F nearly constant along the layer. The distortion of the contours of F in the outflow region, however, indicates that reconnection is taking place due to the small dissipation included in the computation. Since the current density j assumes a quasi-singular cusplike behavior at the X-point (see section 6.2.4), which produces a similar behavior of F, F may change topology even in the limit of vanishingly small dissipation. 6.2.4 Electron Kelvin-Helmholtz instability of the current layer In the derivation of the scaling laws in section 6.2.2 we have ignored the substructure of the current density in the inertial layer de. Instead of the bell-shaped cross-layer current profile typical of resistive MHD (see section 3.2.1), noncollisional current layers, because of free electron acceleration at the X-point, develop a cusp-like profile. Hence transverse gradient scales become very short, ls < de, such that we > j/de, the vorticity terms dominate in (6.47), d2eve • Vwe > B • Vy, and the equation reduces to the 2D Euler equation for the electron flow, dtwe + \e - Vwe = 0. As shown in section 4.6.2, the strong flow along the inertial layer may be Kelvin-Helmholtz unstable, if the parallel magnetic field is sufficiently weak. While a resistive current sheet is Kelvin-Helmholtz stable since the outflow velocity is sub-Alfvenic, v < VA, the magnetic field at an EMHD current layer is weak as discussed above, hence shear flow instability is expected to occur for sufficiently narrow structure ls, i.e., for sufficiently small dissipation coefficients. Figure 6.5 illustrates the onset of the electron Kelvin-Helmholtz instability (EKHI) and the generation of EMHD turbulence. Since the largest flow velocity occurs at the layer edges, turbulence appears mainly in these regions. (For still smaller dissipation coefficients than used in the simulation of fig. 6.5, which further reduce the sublayer width ls and hence increase the velocity shear, the entire layer breaks up leading to large regions of fully developed EMHD turbulence, which is found to exhibit many interesting properties, e.g., a Kolmogorov k~5^ energy spectrum, see Biskamp et a/., 1996, 1999). To this point, our consideration of the EMHD dynamics has been confined to 2D geometry. While this approximation is justified for a stable laminar flow configuration, it is no longer valid when the system is susceptible to the EKHI, where not only the azimuthal flow v but also the axial current density jz in the central part of the layer will

216

6 Noncollisional reconnection processes

(a)

(b)

(c)

Fig. 6.5. Simulation of electron Kelvin-Helmholtz instability giving rise EMHD to turbulence in the outflow regions of coalescing flux bundles, (a) xp, (b) (j>69 (c) h (d) we. drive the instability. Since the extent of the layer in the axial direction is macroscopic, much larger than in the azimuthal plane, there is no finitelength effect limiting the excitation of turbulence due to convection out of the unstable region, which implies that in 3D the layer should become turbulent more readily. The first 3D EMHD simulations of the break-up of a current layer were performed by Drake et al, (1994). Here, we discuss some more recent numerical studies, where the current layer is generated self-consistently during the coalescence of two flux bundles (Drake et al, 1997). This behavior is illustrated in fig. 6.6, showing a numerical simulation with the same initial conditions as in the 2D simulations of figs. 6.3, 6.4 but for somewhat larger de than used in the 2D runs. Shown are greyscale plots of

6.2 High-fi whistler-mediated reconnection

(a)

217

(b)

Fig. 6.6. 3D simulation of flux bundle coalescence illustrating the development with time of 3D EMHD turbulence for de = 0.2. Grayscale plots of axial current density jz in the azimuthal {x,y) plane and the (x,z) plane along the diagonal x = y across the current sheet, (a) t = 0.8, (b) t = 1.75 (from Drake et al, 1997). the axial current density in the x, y and the x, z plane. At time t = 0.8 (in units of the whistler time, (6.40)) the coalescence process is still laminar and practically two-dimensional, characterized by a smooth narrow current sheet between the two flux bundles. Subsequently turbulence sets in, generated mainly by EKHI of the axial current, which leads to a turbulent broadening of the entire layer as shown at t = 1.75. A 2D control run with the same parameters and initial conditions remains laminar - excitation of 2D turbulence as in fig. 6.5 would require a smaller dissipation coefficient - which demonstrates the importance of fully 3D geometry. The 3D simulation just discussed has been performed assuming that the mean axial field Bzo vanishes. Since the EKHI is stabilized by a parallel magnetic field, an axial field of finite strength Bzo > 1 strongly reduces the effect of the instability, as shown in section 4.6.2. It might be expected that the turbulent break-up of the current layer would result in a strongly enhanced reconnection rate. This is, however, not the case. The reason is that, even in a laminar state, the electron current layer does not control the reconnection speed, as we have shown in the preceding section. Electron inertia is needed only to break up the magnetic topology, but does not lead to magnetic flux pile-up which could slow down the reconnection process. The main effect of the turbulence is to enhance the energy dissipation rate. In the absence of EKHI the magnetic energy is mainly converted into kinetic energy of the electrons in the form of high-speed flows, while the dissipation rate is small, decreasing

218

6 Noncollisional reconnection processes

with decreasing diffusion coefficients. However, when the current layer becomes turbulent, a finite fraction of the magnetic energy is dissipated, this fraction being independent of the diffusion coefficients. This is a general property of EMHD turbulence (Biskamp et al, 1996), which is similar to hydrodynamic and MHD turbulence. In a real physical system, the dissipated energy would appear as electron thermal energy. 6.2.5 Ion-controlled reconnection dynamics Whereas in the region / < d,-, which is governed by the properties of the whistler, the ion motion can be neglected, i.e., the ions can be assumed infinitely heavy, the global dynamics of course depends on the finite ion inertia. It determines the reconnection rate £, to which the EMHD configuration around the X-point adjusts automatically. The first simulations of whistler-mediated reconnection in the coalescence of two flux bundles were performed by Mandt et al. (1994), where the ions were treated as macro-particles. Here, we discuss results from two-fluid studies (Biskamp et a/., 1995), which allow us to use higher spatial resolution. Figure 6.7 illustrates a typical state during the coalescence process, where both de(= 0.015) and d\(= 0.1) (corresponding to mass ratio mi/me = 44) are small compared with the diameter of a flux bundle L ~ 2, in contrast to the EMHD coalescence studies in section 6.2.3, where L < d\ is implied. The difference of the flow patterns in the vicinity of the X-point in clearly recognizable, electron streamlines converging toward the reconnection point, where ion streamlines are already deflected into the downstream cone. The magnitude of the inflow velocities along the main diagonal in fig. 6.7 are plotted in fig. 6.8. While on the MHD scale / > d{ both ions and electrons move closely together, the ions are slowed down on a scale / ~ d,-, whereas the electrons are further accelerated reaching a maximum speed at a much smaller distance / ~ de. A series of computer runs has been performed with different values of di and de (Biskamp et a/., 1997). The main result is that the maximum reconnection rate reached during the coalescence is almost independent of du in other words, the reconnection speed is quasi-Alfvenic. Consistent with this fast process, the ion flow does not form an extended layer. In fact, the ion velocity changes abruptly in a shock-like way from the lower inflow to the higher outflow values, a behavior reminiscent of Petschek's shock configuration (see section 3.3.4). Recent simulations of the tearing mode (Shay et a/., 1998), which have been performed for substantially larger system size as used before, seem to confirm the Alfvenic scaling. We will come back to these studies in section 7.6.3. The results of the incompressible theory are corroborated by numerical solutions of the fully compressible equations (6.24)-(6.27) in the case of

6.2 High-P whistler-mediated reconnection

219

271

Fig. 6.7. Simulation of flux bundle coalescence with d\ = 0.1 and de = 0.015. Plotted are contours of xp, the electron stream function 4>e and the ion stream function , 1. This is also true for the case of a moderate axial field Bzo ~ Bj_, though this field destroys the up-down, right-left symmetry of the incompressible dynamics due to the term BZQV • \e in the Bz equation, which makes the analysis of the reconnection region more complicated. As we have already seen in the EMHD case, three-dimensional effects are important making the current sheet turbulent, both in the central layer and along the separatrix, see fig. 6.9. Simulations using a 3D compressible model have been performed by Rogers et a\. (2000). Two different instabilities can be identified driving the turbulence, the electron Kelvin-Helmholtz instability (EKHI) discussed in section 6.2.4, and the lower-hybrid drift instability (LHDI), which will be treated in more detail in section 7.3. The EKHI, being essentially incompressible, generates primarily fluctuations of the current density yz, while the LHDI gives rise to perturbations of the plasma density. Both instabilities remain active in the presence of a moderate axial field. The turbulence does not further enhance the reconnection speed, which is already quasi-Alfvenic in the 2D laminar case, but increases both ion and electron dissipation. It should, however, be kept in mind that dissipation in these high-/? cases is, properly speaking, a kinetic process, since Larmor radii are large and particle orbits are

6.3 Low-fi noncollisional reconnection

221

-2

Fig. 6.9. Contours of current density jz in (a) a 2D simulation, (b) a fully 3D simulation (from Rogers et al, 2000). nonadiabatic in the reconnection region. As will be seen in section 7.6, finite particle excursions broaden the sheet structure, in particular for the electrons. 6.3 Low-/? noncollisional reconnection 6.3.1 The four-field and three-field models Most laboratory plasmas are embedded in a strong, nearly constant magnetic field Bo. As already discussed for the reduced-MHD approximation in section 3.1.4, in a low-/? plasma only the field B± generated by currents along Bo is an independent dynamic variable. To determine the small change of the main field 5B\\9 one uses the fact that for low jS the magnetosonic wave is the fastest mode, even faster than the whistler. Let us compare typical time-scales of both waves, %A = (vAk±)~1, (3.24), and for 1> (6-37): xw

since in the vicinity of a reconnection region k\\/k±_ will be very small, e.g.,fe||/fcj_~ c/copeR ~ 10~4 in a typical tokamak plasma. Hence SB\\ is

222

6 Noncollisional reconnection processes

determined by the instantaneous perpendicular pressure balance Vj_(p + B 2 /8TT) = 0, hence SB^ = -4ndp/B. With B o = Bez and B ± = ez x Vxp and Ez = dtxp the z-component of Ohm's law (6.19) becomes

'i+ni,

(6.63)

at

where e z x V4>e = ye± = vu. -]±/ne

= viJL

- e z x V(p,- + pe)

fin

or, in nondimensional form with the units introduced in (6.20) and assuming isothermal behavior,

Insertion into (6.63) gives, with \f± = ez x Vcpu dtxp - B • V = -di(Tt + T,)V|, Inn + d2edj/dt + qj,

(6.77)

dtw + V_L • Vw = V • (v*; • VV(/>) + B • V/' + )U±V2w,

(6.78)

.

(6.79)

6.3.2 Linear stability theory The linear stability theory of reconnecting modes, in particular of the internal kink mode, has received much attention in the literature. A general classification of the different collisionality regimes was developed by Drake & Lee (1977). For the simplified treatment presented here it suffices to distinguish between the collisional regime if the resistive diffusion layer exceeds the electron inertial length, dn = ^rj/co > de, and the collisionless regime in the opposite case 8n < de. Moreover, one has to distinguish between very low plasma pressure [i < 2mg/m,-, where ps < de, and relatively high pressure /? > 2me/m,-, where ps > de. Here ps = cs/Q>i and cs = \J{Ti + Te)/nii is the isothermal sound speed. Though a proper treatment for large ion gyroradius seems to require a kinetic approach,

6.3 Low-P noncollisional reconnection

225

which we postpone to chapter 7, the fluid theory turns out to be at least semi-quantitatively correct. We start by deriving the dispersion relation for linear modes in a hom*ogeneous plasma embedded in a uniform field B = Bez. Linearizing the four-field equations (6.63), (6.64), (6.75), (6.76) neglecting dissipation and making the Fourier ansatz \p = \p\ exp{—icot + k • x} etc., r Te)m = 0, = 0, con\ —^||^||i —d//c||fcj_tpi = 0, cov\\\ — (Tj + Te)fe||ni = 0,

(6.80) (6.81) (6.82) (6.83)

gives the dispersion relation

{co2 - k2c]) [co\\ + k{d2e) - k2v2A] - co2 k\v\ k2±P2 = 0,

(6.84)

where we have reintroduced the dimensional notation to elucidate the physics, B -+ vA, Tt + Te - • c2, d?(T,- + Te) -> p2. For long wavelength ^_LPS 1 and for ps > de the mode co2 ^ k?\^p2k^o2A is called the kinetic Alfven wave. It is generated by parallel electron compressibility, the r.h.s. in (6.79), which will be shown to act as a powerful reconnection mechanism. In the next step we include a density gradient, n'o = —no/Ln* in the x-direction and hence also pressure gradients p^ = Tjnf0. Since we are primarily interested in small scales kLn > 1, the change of the equilibrium density over a wavelength is negligible, i.e., when not differentiated, no can be taken constant. The density gradient introduces the diamagnetic drift frequencies co*7-, which lead to important modifications of the dispersion relation. Neglecting the sound-wave branch as indicated above, we restrict consideration to the three-field model equations (6.77)-(6.79), which are * Since in a plasma column the radial density gradient is usually negative, it is convenient to define Ln with the minus sign.

226

6 Noncollisional reconnection processes

linearized about an equilibrium state with vanishing electrostatic field \E = 0 and hence with the equilibrium mass flow m^v*/, (co - (o*e)\pi + fe||B0i - k\\(Ti + Te)dtm = -k\{ir\ + cod2e)xpu

(6.86)

oj(t)i+k\\B\pi = 0,

(6.87)

(co — co*i)n\ + kyfiQ^i = k\\k]_dixpi9

(6.88)

where in our isothermal model T

A

e = -ky - ^ n

ri0,

(6.89)

or in dimensional form eBLn'

o)*e = k

y

^-.

(6.90)

Thus by including diamagnetic effects and resistivity the dispersion relation (6.85) is generalized

1

"" CO(CO — CO*j)

1+ \

felpj

(O — (J0*e

)+ )

CO — CO*e

fci

(6.108)

i.e., the behavior of the linear instability continues into the nonlinear regime, which is faster than the algebraic growth, (4.101), w/ oc t2, of the resistive kink mode. In the preceding discussion quasi-stationary reconnection is implied, in particular E + v • VF = 0 from (6.104). However, this equation cannot be valid at the X-point, where the second term vanishes; an argument which is sometimes used against electron inertia as a true reconnection mechanism. In fact, (6.104) indicates, that the dynamics at the X-point is governed by electron acceleration along the z-direction. Reconnection occurs, but not in a quasi-stationary way. The process has been analyzed by Ottaviani & Porcelli (1993). The authors consider the double tearing mode in a plane configuration, specifically the periodic equilibrium xpo = cosx, similar to (4.74), but the results apply, mutatis mutandis, also to the cylindrical kink mode, the dynamics being very similar, as discussed in sections 4.2 and 4.3. The starting point is the assumption that the outflow is confined to a layer of constant width de and macroscopic length A ~ k~l > de, where k is the wavenumber of the mode, hence one can make the simple ansatz for the stream function 0 for x < 0. The following integral is independent of x:

f vj

=

[[X M =

g Jxo J-oo g Jxo g(*') which equals the absolute value of the uniform shift in the region outside the layer, as seen by choosing xo,x in this region, where \g\ = 1. The nonlinear regime is characterized by £Q> de. Different from the flow layer de, the current layer width d(t) shrinks with time which is defined by •* dx>

L

8

, ——- = — de In — = Co, de g(x') de

„

whence .

(6.U2)

This implies that the initially shallow F(XQ) profile inside de is squeezed into a narrow sublayer of width d, forming much stronger gradients. But this only affects the current density, whereas the flux function xp remains smooth, varying only on the scale de, as can be seen by integrating the equation d2eV2\p — \ p ~ d2ed2x\p —xp = —F using Green's function — (l/2d e )exp{—|x — x'\/de}, 1 c°° tp(x,y,t)~— / ^x-x'^F(x\y,t)dx\ 2ae J-oD

(6.113)

which shows that xp is smoothed out over the flow layer width de. Hence in the resonant layer |x| < de one has F(xo) * F0 - \x2F^ = F0- \{x - ^fFl

(6.114)

F(x0) ^ ip(xo) ~xp0- i(x - £o)2Vo>

( 6 - 115 )

since in the nonlinear regime ^ ^ ^o > de and j(xo) is smooth, such that

fd2e < Vff.

Along the layer, F ~ const (see also fig. 6.2(b)), hence the variation of the current density along the layer is

dj = dxp/d2e~-\{^/de)2xpl

(6.116)

To derive an equation for £o we make use of the vorticity equation, which we integrate over a quadrant of the layer, as indicated infig.6.1: j

fwd2x=

fB-Vjd2x.

(6.117)

232

6 Noncollisional reconnection processes

Here we are neglecting the vorticity advection term, which is small compared with dtw in this strongly nonstationary phase. The l.h.s. of (6.117) becomes approximately

using (6.109) and (6.110) and fcA ~ 1. Similarly, the r.h.s. of (6.117) becomes {

B- Vjftc = jli/H, = \ ^f- = i«J?g,

(6.119)

inserting (6.114) and (6.116). Combining the two expressions gives the desired equation for the nonlinear evolution of £o d%/df

= c^

(6.120)

where £o = €o/de, T = tyo, yo = deky)Q (cf. (6.95)), and c is a numerical factor of order unity. Equation (6.120) has the similarity solution

indicating explosive (= faster than exponential) growth with a singularity at the finite time to = O(JQ1). The solution (6.121) can, however, only be valid in the very early nonlinear stage, where £o ~ de, since further increase of £o would lead to extremely narrow widths 5 given by (6.112), which should easily be smoothed out by residual resistivity or, more importantly, by electron viscosity. Hence the substructure in the current layer will not be important in the main nonlinear phase, much in the same way it has been neglected in the EMHD analysis in section 6.2. Instead of the explosive behavior, the dynamics is governed by the quasi-stationary reconnection process leading to exponential island growth, (6.108). Figure 6.10 gives the current density from a numerical simulation of the electron inertia-driven kink mode with the parameters de = 0.01, rj = 10~6. The value of the resistivity is sufficiently small such that de^> 8n, i.e., the system is in the collisionless regime, where the growth rate is given by (6.95) and the sheet width is de. The current-density distribution resembles that in the resistive kink mode, see fig. 4.14(a). 6.4.2 Kinetic Alfven-wave-mediated reconnection In most plasmas of interest jS is sufficiently high, /? > me/mu such that the pressure term in (6.64) is more important than electron inertia, ps > de. In this case the structure of the reconnection region is quite different

6.4 Nonlinear noncollisional kink mode

233

Fig. 6.10. Grayscale plot of simulation of the current density in the nonlinear development of the electron inertia-driven kink mode showing an extended current sheet of width de. from the current sheet of macroscopic length characteristic of the low-/? electron inertia-dominated dynamics. Since the pressure term requires a perpendicular field component, B- Vn ~ Bxdxn, it is more effective, leading to faster reconnection, if Bx is large, i.e., if xp forms an X-point. To discuss the structure of the reconnection region it is useful to simplify (6.77) by neglecting the ion temperature and assuming that the equilibrium pressure gradient are small, such that diamagnetic effects can be neglected. Writing V\\j = B~XB • V/ makes clear that in this case n obeys the same equation as (dt/Bz)w. Since in the evolution of an instability the memory of the initial perturbation is rapidly lost, we can substitute n in Ohm's law, which leaves us with only two equations for xp and w, dtxp - B

= -p?B- Vw +

d2edj/dt,

(6.122)

dtw + v • Vw = B • V/, 2

(6.123) 2

since dfTe/B^ = p s in the units introduced in section 6.1. The p s effect gives rise to the kinetic Alfven wave, (6.85). (Electron inertia and residual

234

6 Noncollisional reconnection processes

dissipation are only important in the immediate vicinity of the X-point to allow a finite reconnection rate dtxp, since here the pressure term vanishes.) For stationary conditions and de < p s , (6.122) and (6.123) reduce to

E - B • V(0 - p?V2(/>) = 0,

(6.124)

v V V 2 ( / > - B - V ; = 0,

(6.125)

where E = dt\p = const. Compared with the resistive MHD case, (3.140) and (3.141), the p2 term implies a smoothing of the flow over a distance ps. For not too large E the plasma inertia term in (6.125) can be neglected to lowest approximation and the solution of B • V/ = 0 around the X-point is again xp = \{x2-a2y2l (6.126) while the stream function (f> satisfies the equation

x — ay

= f(x,y),

(6.127)

where f(x9y) solves the equation E — B • V/ = 0. Equation (6.127) has the formal solution 0(x) = / G(x - x')/(x')dV

(6-128)

with Green's function (6 129)

'

A numerical solution of (6.127) is plotted in fig. 6.11. For finite ps the singularity of 4> on the separatrix x = +ay is eliminated, and there is only a weak logarithmic singularity of w = V2, which is smoothed by residual resistivity or electron viscosity. Hence the effect of plasma inertia, the first term in (6.125), remains small, which allows the Xpoint configuration, (6.126), to persist also for finite reconnection rate, in contrast to the resistive case, where the inertia term enforces the formation of a macroscopic current sheet. It is worth mentioning that for ps = 0 the electron inertia term in (6.122) modifies only the current profile, but does not smooth the singular flow pattern. Hence in this case the plasma inertia term in (6.125) is not negligible and an X-point configuration is not possible. Instead a macroscopic current sheet is formed, see the preceding section, much in the same way as in the resistive case. Numerical simulation studies of pressure-mediated reconnection have been performed for the coalescence of two flux bundles (Kleva et al, 1995),

6.4 Nonlinear noncollisional kink mode

235

(b)

Fig. 6.11. Solution of (6.127): (a) flow pattern <j>; (b) vorticity V2. the cylindrical kink mode (Aydemir, 1992), and the double tearing mode, which differ only by the geometry of the system, the reconnection dynamics being very similar. Simulations confirm that the magnetic configuration exhibits an X-point with a downstream cone of finite angle. In the limit d2e < p], vorticity and current density form shock-like structures along the magnetic separatrix just as in the solution of (6.124), see fig. 6.1 l(b). Figure 6.12 illustrates a simulation run with de = 0.01, ps = 0.02 and dt = 0.2, i.e., nii/me = 400, which demonstrates the change of the currentdensity distribution compared with the corresponding electron inertiadriven case with ps = 0 shown in fig. 6.10. Reconnection now occurs in a micro-current sheet. Within the range of ps/L values amenable in the

236

6 Noncollisional reconnection processes

Fig. 6.12. Grayscale plot of current density from a kink mode simulation. simulations the reconnection rate E = dt\p is independent of ps during the early reconnection period, decreasing only weakly with ps during the fastest dynamic phase. It is found that for sufficiently small values of the resistivity the structure of the reconnection region and the reconnection rate are independent of r\. When de is increased, the length of the small central current sheet increases, but the same basic configuration persists even for dezips. Only if de > ps does electron inertia dominate and a macroscopic current sheet is formed. The fast reconnection dynamics in the pressure-dominated regime is reflected in the nonlinear evolution of the kink mode. Whereas in the electron inertia regime the kink mode effectively continues to grow exponentially in the nonlinear phase, the pressure-dominated mode exhibits an explosive time behavior, as was first observed numerically by Aydemir (1992). This behavior can be interpreted by a simple model. Since reconnection does not depend on the smallness parameters of the system p s , de, rj, i.e., it is essentially Alfvenic, the angle a of the outflow cone is

6.4 Nonlinear noncollisional kink mode

237

determined purely geometrically, ria~f,

(6.130)

where £ is the helical shift of the of the main magnetic axis, and the island width is w/ ~ 2£, see section 4.3.2. In addition, we use the continuity equation (6.131) As in the resistive M H D case the outflow velocity is determined by pressure balance and hence equals the upstream Alfven speed vo = vA~\^\?;/2,

(6.132)

while the inflow velocity is given by the geometry of the system no = £/2. Combining these relations we find

fx^2,

(6.133)

with £ = £/ri and T = f|y>ol = t/zA- Hence the helical shift grows explosively I = —!—.

(6.134)

TO"?

Note that we did not make use of Ohm's law. The only physical assumption regards the pressure balance in the layer, i.e., vo = vA, which is justified for processes slow compared to the compressional Alfven time-scale. The finite-time singularity in (6.134) should not be taken literally, since this self-similar solution is only valid in the early nonlinear phase of the kink mode evolution, much as the self-similar solutions (6.108) and (4.101). At large amplitude £ ~ r\ the change of the radius of the plasma core and of the helical flux distribution must be considered, which leads to saturation and final decay of the growth rate. 6.4.3 Influence of diamagnetic effects In the previous subsection the equilibrium pressure gradient has been neglected, thus omitting diamagnetic rotation. In fact, the electron pressure term in the ^-equation contains two distinctly different effects, which can clearly be seen when linearizing the term S(B • Vpe) = Bo • VdPe + SBrp'a.

(6.135)

The first term on the r.h.s., which gives rise to the kinetic Alfven wave, is responsible for the nonlinear acceleration of the growth, (6.134), discussed

238

6 Noncollisional reconnection processes

above. The second term introduces the diamagnetic effect, which implies that a perturbation propagates in the poloidal direction with the drift velocity v*e = cTe/(eBzLn) = pscs/Ln, which constitutes an efficient damping process by phase mixing xp and 0, as discussed in section 6.3.2, stabilizing the linear kink mode if co* > yo. For co* ~ yo the kink mode is still unstable, but is found to saturate at a finite island width. Let us therefore study in more detail the nonlinear evolution in this regime. Because of the MQ term in the continuity equation (6.79) we can no longer replace n by (di/Bz)w, but have to regard n and w as independent dynamic variables, i.e., go back to the full three-field equations (6.77)-(6.79), which we write in a slightly different form, neglecting ion temperature and introducing electron viscosity instead of resistivity: dtF + v • VF = -atfpB • Vn - /ueV2j,

(6.136)

dtw + v • Vw = B • V/ + /iV2w,

(6.137)

2

dtn + v • Vn = a,-B • V/ + DV n,

(6.138)

2

where F = \p - d ej, a,- = diB0/Bz = c/copiR, since q = riBz/RB0 = 1, j}p = 4nnoTe/BQ, R = major torus radius, Bo = poloidal field at r = r\. We discuss a numerical simulation study with the following parameters: de = 10" 2 , a*: = 5 x 10" 2 , iie = 10" 9 , fi = D = 10" 6 ; fip and Ln are varied, where Ln is the gradient scale of a parabolic density profile «o(r) = 1 — ir/Ln)2. These parameter values are chosen more for numerical convenience than for a realistic tokamak simulation. We only want to demonstrate the qualitative behavior. If Ln is sufficiently large, Ln > 1, i.e., for weak diamagnetic effects, the pressure term is found to be destabilizing nonlinearly. Figure 6.13(a) gives the growth rate defined by the kinetic energy, y(t) = (dEK /dt)/2EK, for Ln = 10 and several values of jBp. While according to (6.99) the linear growth rate yo changes only weakly with /} p , yo °c /V > the most conspicuous feature is the accelerated growth in the nonlinear phase, first noted by Aydemir (1992), which corresponds to the explosive solution (6.134). We have therefore plotted y/yo as function of yo£, which shows that the nonlinear enhancement of the growth rate increases with $v. If, on the other hand, Ln is small, Ln < 1, where the diamagnetic effects dominate, increasing fip reduces the nonlinear growth rate, fig. 6.13(b), which for sufficiently high f}p results in finite-amplitude saturation, f$p >, 2 in this case. Saturation occurs when the helical shift £ of the plasma core no longer points toward the X-point but is rotated out of phase, such that reconnection is no longer driven. From a series of simulation runs one obtains a semi-quantitative condition for saturation: co* > ymax, (6.139)

239

6.4 Nonlinear noncollisional kink mode

(a)

10

Fig. 6.13. Normalized growth rate of the nonlinear kink mode as function of normalized time: (a) Ln = 10; (b) Ln < 1 (from Biskamp & Sato, 1997). where co* is the linear mode frequency and ymax the maximum (nonlinear) growth rate. For co* < ymax full reconnection occurs, even if co* exceeds the linear growth rate. This condition generalizes the linear stability criterion, which relates co* to yo6.4.4 Criterion for fast reconnection In this subsection we want to generalize the discussion following (6.124) and (6.125), formulating a criterion to decide whether reconnection is fast, i.e., basically Alfvenic, or depends more or less strongly on the respective

240

6 Noncollisional reconnection processes

small parameter. For stationary conditions the process is governed by the equations E-BV$ = R, (6.140) v W 2 0 - B - V / = O,

(6.141)

where R comprises the nonideal effects in Ohm's law. We neglect viscosity in the momentum equation, since it is not important for the argument. In the ideal limit R = 0, the solution of (6.140) in an X-point configuration, (6.126), E x + ay (6.142) ~1 x — ay is singular on the separatrix x = +y, hence does not solve the momentum equation, where the inertia term would diverge, while the magnetic term vanishes. What are the conditions to be imposed on R, which is a functional of xp and (/>, to allow a solution of both (6.140) and (6.141) with a finite-angle X-point, which is the prerequisite for a fast process? We suggest the following criterion: Fast reconnection occurs, if R regularizes the flow (j) in Ohm's law for an X-point configuration. In this case the inertia term in (6.141) is small for sufficiently small (but finite) velocity, such that only weak currents (essentially at the separatrix) are required for B • Vj to balance the inertia term, and the X-point character of the magnetic field is preserved. If R does not regularize the flow, the configuration is completely changed. In general a macroscopic current sheet is formed, in accord with Syrovatskii's theory (see section 3.3.1), and the reconnection rate depends sensitively on R. The criterion can be satisfied if R contains higher-order derivatives of (j) as, for instance, in the p2 term in (6.122), which smoothes the basic singularity, (6.128). Only in the immediate vicinity of the X-point is a micro-current sheet set up, providing the actual reconnection mechanism due to residual resistivity or electron inertia. The length of this sheet remains small, shrinking with r\ or de. By contrast, the electron inertia term R = d2ev • V/ alone does not affect the singularity of 0, which arises from the inversion of the operator B • V, nor does a purely resistive R = rjj. In both cases extended current sheets are formed as discussed in detail in sections 6.4.1 and 3.3.3, respectively. It is interesting to consider the relative importance of the p2 and the d2e terms in (6.122). The reconnection process is fast, dominated by the p2 term even for ps < de. Only if ps < de does the d2e term control the dynamics with a macroscopic current sheet. In the high-/? case, where the ions are immobile at scales < du the ion or plasma stream function cp in Ohm's law is replaced by (pe in this region. Here electron inertia (section 6.2.2) and non-gyrotropic terms in

6.5 Sawtooth oscillation in tokamak plasmas

241

the electron pressure tensor (section 7.6.2) have a regularizing effect on (j)e similar to that exerted by the p] term on / in the low-/? case, which allows an X-point configuration to be established. The X-point solution, (6.126), (6.142), or the regularized flow, (6.128), are symmetric with respect to inflow and outflow, the angle of the downstream cone being a free parameter. This is true for sufficiently low velocity, where the inertia term in the momentum equation is unimportant. At higher velocities the term introduces an inflow-outflow asymmetry. Inertia limits the outflow velocity to the Alfven speed, as shown in Petschek's model, the basic configuration in the presence of a fast reconnection mechanism, which now couples the downstream angle to the inflow speed. 6.5 Sawtooth oscillation in tokamak plasmas The sawtooth oscillation is the classical paradigm of a reconnection process in a tokamak plasma. In my previous book (Biskamp, 1993a) I have given an extended review of the historical developments during the 20 years following the first observations by von Goeler et al (1974) and the first theoretical model by Kadomtsev (1975). Here, I will therefore concentrate mainly on the results obtained during the 1990s, which have contributed essential new pieces to the puzzle posed by this intriguing phenomenon. 6.5.1 Basic experimental observations The sawtooth oscillation is a relaxation process in a tokamak plasma corresponding to a periodic flattening of the temperature and the density profiles in the core of the discharge column. The conventional diagnostic tool for studying this phenomenon is an array of soft X-ray (SX) diodes measuring the emission along different chords across the plasma, as indicated schematically in fig. 6.14. The signal along a chord, passing through the central part of the plasma column is modulated by a periodic oscillation consisting of a slowly rising part, the rise phase, followed by a rapid drop, the sawtooth collapse or crash. The signal has the shape of a sawtooth, whence the name. For a sufficiently off-centered chord, the signal is inverted, with a sudden rise coinciding with the sawtooth collapse and a gradual subsequent decay. The cross-over from one type of signal to the other occurs at the radius where the safety factor q passes through unity, which is called the inversion radius. Sawtooth oscillations, or simply sawteeth, occur in practically all tokamak devices under various operational conditions. While the standard diagnostics is still via SX emission, which allows the highest time resolution but provides only an indirect measure of the electron temperature,

242

6 Noncollisional reconnection processes Movable soft X-ray detectors

Vacuum vessel

slot

(b)

Fig. 6.14. Observation of sawtooth oscillations by soft-X-ray diagnostics. Schematic drawings of (a) the experimental set-up, (b) the SX signals along chords (1) and (2).

Fig. 6.15. Schematic drawing of the m = n = 1 precursor oscillation superimposed on the m = n = 0 sawtooth signal, TO is the sawtooth period, TI the collapse time-scale. these observations are now supplemented by direct measurements of (i) the electron temperature using electron cyclotron emission (ECE), and (ii) the electron density using microwave interferometry. An indication of the basic mechanism responsible for the sawtooth relaxation is obtained from a precursor oscillation, which in many cases is observed to precede the collapse, fig. 6.15. While the main sawtooth signal is symmetric with an m = n = 0 mode signature, the precursor corresponds to a helical shift of the central part of the plasma column, i.e., an m = n = 1 mode, which gives rise to the observed oscillating signal owing to diamagnetic and plasma rotation. In contrast to the simple regular sawteeth characteristic of earlier small tokamak experiments with minor radius a ~ 10 cm, the saw-

6.5 Sawtooth oscillation in tokamak plasmas

243

(a) 8.4

f( S )

(b)

10.5

11.0

11.5

12.0 f(S)

Fig. 6.16. ECE measurements of the central electron temperature for different types of sawtooth oscillations obtained in the JET tokamak. (a) Normal sawteeth from an ohmic discharge (from Edwards et a/., 1986); (b) so-called giant sawteeth from an ion-cyclotron-heated plasma (after Campbell et al, 1987).

tooth phenomenon exhibits a more complex behavior in large-diameter tokamak plasmas with minor radius a ~ 102 cm. Depending on the operational mode - pure ohmic, neutral beam, or radio-frequency (mainly ion-cyclotron resonance) heating - the period, amplitude and shape of the SX or ECE signals may vary considerably, examples being given in fig. 6.16. Sawteeth often show a compound structure, where roughly half-way through the rise phase a minor or partial collapse occurs, as seen in fig. 6.16(b), which affects mainly the region close to the q = 1 radius, while the central energy loss is small. The sawtooth period is directly related to the confinement time TEI of the energy W\ inside the q = 1 radius defined by — =P,

(6.143)

where P is the heating power deposited inside this radius. The sawtooth collapse results in a certain fraction AW\ of the energy being ejected. This generates an imbalance between heating power and energy transport, from which the plasma recovers during the rise phase, dt

_

(6.144)

244

6 Noncollisional reconnection processes

Since W\ increases roughly linearly, at least for normal sawteeth as in fig. 6.16(a), the sawtooth period to is given approximately by T0 ~ T£i oc r\

(6.145)

assuming a diffusive energy loss during the quiescent period. In compound sawteeth the period is roughly doubled. Relation (6.145) indicates that the sawtooth period increases rapidly with increasing plasma radius. The collapse time-scale x\ is found to be essentially independent of the sawtooth period, increasing at most linearly with radius x\ oc r\, such that in the large tokamak devices the ratio of rise to collapse time is 103. 6.5.2 The safety factor profile Knowledge of the g-profile and its time variation during sawtooth activity is of fundamental importance for understanding and modeling the sawtooth phenomenon. Since the signature of the precursor as well as the fast dynamics during the collapse, which will be discussed below, clearly point to the kink mode as the basic agent responsible for the collapse, the central value qo = q(0) should be smaller than unity, giving rise to a resonant surface q(r\) = 1 at some finite radius r\. On the other hand, the resistive kink mode itself leads to complete reconnection of the internal helical flux, which relaxes the plasma into a state with q > 1 as outlined in section 4.3.3. Since in typical hot and large-diameter plasmas the resistive skin time is much longer than the sawtooth period, one could expect that after the first sawtooth the ^-profile remains essentially flat and close to unity inside the inversion radius. In earlier tokamak experiments the current profile and hence the qprofile were inferred from the condition of resistive equilibrium rjj = const using the measured electron temperature profile and a reasonable guess of the impurity ion distribution, which resulted in rather inaccurate qvalues in the crucial center part. During the 1990s, however, several diagnostic methods have been developed to measure the poloidal magnetic field directly using either the Faraday rotation effect (Soltwisch et a/., 1987) or the motional Stark effect (see, e.g., Levinton et al, 1993), which allow determination of the ^-profile with an error of less than 5%. Averaging over many equivalent sawteeth, making use of the nearperiodicity of the sawtooth oscillation, yields accurate information about the time behavior of qo, the central q value, see fig. 6.17. The surprising and now firmly established result is that qo remains distinctly below unity, typically jumping up from 0.7 to 0.8 during the collapse and slowly decaying back to 0.7 during the sawtooth rise phase. The ^-profile

6.5 Sawtooth oscillation in tokamak plasmas

245

0 0.04 0.08 0.12 Time from sawtooth crash (s)

Fig. 6.17. Average time behavior of the g-value on the magnetic axis during a sawtooth period (from Levinton et al, 1993).

1.5

/p=1.8MA Pb=10.5MW £ f =4.6T n =2.6x10 19 rrr 3

• Before sawtooth crash - statistical uncertainty

0.5 2.7

I

I

I

2.8

2.9

3.0

J_ 3.1

I

3.2

Major radius (m)

Fig. 6.18. The ^-profile measured before and after the sawtooth collapse for two separate sawteeth (after Levinton et a/., 1993).

before and after the collapse is shown in fig. 6.18, which differs fundamentally from the behavior of the resistive kink mode, see fig. 4.17(b). These experimental findings impose severe constraints on the theoretical interpretation of the sawtooth collapse dynamics, as we will see later.

246

6 Noncollisional reconnection processes 6.5.3 Stabilization and onset of sawtooth oscillations

One of the problems arising from the observed behavior of the central safety factor refers to the kink mode stability of the plasma column. While the ideal MHD mode requires finite plasma pressure, in particular must the poloidal /? exceed a threshold, fip > fic — 0.3 for stability (Bussac et a/., 1975), finite resistivity makes the mode unstable as soon as qo drops below unity. But the plasma is apparently stable during the the sawtooth rise phase, hence there must be some additional stabilizing mechanism, which is only overcome at the moment of the collapse. Moreover, there are plasma conditions which do not show sawtooth activity at all, though qo is measured to be clearly below unity. The simplest of such effects arise in two-fluid theory. As discussed in section 6.3.2, diamagnetic rotation has a distinctly stabilizing influence. The general dispersion relation of the kink mode has been derived by Zakharov & Rogers (1992). For marginal ideal stability l j j = 0 a simple threshold condition can be given (Levinton et a/., 1994). The kink mode is stable, if

Note that, contrary to the ideal kink mode, in two-fluid theory the pressure gradient is stabilizing, while the shear is destabilizing. For Tre = T( = 0 the criterion (6.146) corresponds to the stability condition co* > yo, where yo is the growth rate, (6.99), in the absence of diamagnetic rotation.* When applied to discharges in the TFTR tokamak, the criterion (6.146) is found to predict very well the presence or absence of sawteeth. Socalled supershot discharges are characterized by peaked pressure profiles and are usually sawtooth-free, while so-called L-mode discharges with flatter pressure profiles typically do show sawteeth. Figure 6.19 gives two representative cases, where either the l.h.s. of (6.146) exceeds r\q[ and no sawteeth are observed (a), or the l.h.s. is smaller than r\q[ and sawteeth are present (b). Figure 6.20 presents a scatter plot of data comprising many discharges with a considerable range of plasma current and heating power, which demonstrates the tight correlation of the stability criterion (6.146) with the appearance of sawteeth. The success of the criterion (6.146) is surprising and somewhat counterintuitive since, by general physical arguments, the pressure gradient constituting a free-energy reservoir should be destabilizing, whereas the shear obstructs plasma motions and should have a stabilizing influence. * If XH > 0, the linear kink mode is always unstable. However, analyzing the nonlinear behavior, it is found that the mode saturates at low amplitude if the diamagnetic frequency is somewhat above the XH = 0 threshold value (Rogers & Zakharov, 1995), such that the criterion remains, practically speaking, valid.

6.5 Sawtooth oscillation in tokamak

247

plasmas

1.0 0.8

"

No sawteeth :

..k = &V + °V>

(7-43)

the probability being proportional to the intensities ("occupation numbers" in quantum mechanical language) of the decay products. It is easy to see that the resonance condition, (7.43), can only be satisfied, if the waves have negative dispersion dco^/dk < 0 or if waves of different types interact, for instance two high-frequency Langmuir waves and one low-frequency ion-sound wave. The second process, nonlinear wave-particle interaction, results from the wave-particle resonances in the expressions in brackets in the last term in (7.42), co*k - co*kf = (k - k') • v, (7.44) the resonance of particle v with the beat wave k — k'. The process is the analog of stimulated Compton scattering, the inelastic scattering of a wave on a particle, hco^ — hco*k> = AE = Ap • v, where from momentum conservation the change of particle momentum is Ap = h(k — k'). In this way Langmuir waves, whose phase velocity cope/k is much higher than the thermal velocity, can interact nonlinearly with the bulk of the electron distribution. Since the last term in (7.42) is proportional to Ik, it gives rise to a nonlinear change of the growth rate, while resonant mode coupling acts as a sink (or possibly also a source) of modal energy.

270

7 Microscopic theory of magnetic reconnection

Weak turbulence theory is a consistent model for a weakly excited dispersive medium. It is, however, found that for many turbulent plasmas the nonlinear processes just discussed are rather insignificant compared with the quasi-linear effect. To go beyond the quasi-linear approximation, nonlinear corrections to the particle orbits are important. A mathematically consistent formalism is Kraichnan's direct interaction approximation (Kraichnan, 1959), which is, however, simply too unwieldy when applied to the Vlasov equation (Orszag & Kraichnan, 1967; Biskamp, 1968) to be of much benefit without severe ad hoc assumptions. It is therefore practicable to follow a route based more on physical reasoning. Dupree (1966) has introduced the concept of broadening of the wave-particle resonance (see also Dum & Dupree, 1970). The basic ingredient of the theory is a modification of the free-particle propagator, i.e., the ( k v — co^) resonance in the quasi-linear diffusion coefficient (7.33), by taking into account the effect of just this diffusion. In particular, one replaces the expression exp{—ik • (x — \t)} containing the free-particle orbit by the average (exp{—ik • x(t)}) over the true orbits, which is calculated in the diffusion approximation. In a one-dimensional statistically hom*ogeneous problem this is described by a diffusion in velocity space U-ikx(t)\ _

2

Q-ikx+ikvt-±k

Dvt3

(7 45)

where Dv is determined self-consistently by the nonlinear integral equation

Dv{v)

= 4 / ^ k2|^|2Re f mZj J 2%

Jo

generalizing the quasi-linear expression, (7.37). The magnitude of the resonance broadening can easily be estimated. For (Dv/k)1^ < Av, where At; is the range of unstable phase velocities, the effect is small, while in the opposite case (Dv/k)1^ > Av resonance broadening dominates over free-particle propagation and the time integral in (7.46) is approximately

Re f° f° dte^-^-^20^ Jo

~ \3{k2Dv)-ll\

(7.47)

Substitution in (7.46) yields Dv in terms of the fluctuation energy (£ 2 ) =

J(dk/2n)k2\(t>k\2,

Dv - {qj/mjf2 (E2)^/kl'\

k'2/3 = (1 /k2")

(7.48)

~ (qrf/mj)1/2,

(7.49)

and the diffusion velocity dv = (2V/c o ) 1/3 ~ (qj/mj)V2(E2)l'4/kj>/2

271

7.1 Vlasov theory and microinstabilities

2(^0)1/2

Fig. 7.2. Phase-space orbits of particles in the potential <po sin x. Areas of trapped particles are shaded. which is just the velocity range of a particle trapped by the potential