This book provides a comprehensive introduction to the theory of magnetic field line reconnection, now a major subject in plasma physics. The book focuses on the various reconnection mechanisms dominating magnetic processes under the different plasma conditions encountered in astrophysical systems and in laboratory devices. The book consists of two major parts: the first deals with the classical resistive approach, while the second presents an overview of weakly collisional or collisionless plasmas. Applications primarily concern astrophysical phenomena and dynamo theory, with emphasis on the solar dynamo and the geodynamo, as well as on magnetospheric substorms, the most spectacular reconnection events in the magnetospheric plasma. The theoretical procedures and results also apply directly to reconnection processes in laboratory plasmas, in particular the sawtooth phenomenon in tokamaks. The book will be of value to graduate students and researchers interested in magnetic processes both in astrophysical and laboratory plasma physics. Dieter Biskamp received his Ph.D. from the University of Munich. Following a postdoctoral period at the MaxPlanckInstitute for Astrophysics he worked at the Space Research Institute in Frascati and became senior research scientist at the MaxPlanckInstitute for Plasma Physics in 1972. Since 1981 he has been head of the General Theory Group and since 1995 has been head of the Nonlinear Plasma Dynamics Group. In 1979 he was visiting professor at the University of Texas and in 1995 COE visiting professor at the National Institute for Fusion Science in Nagoya. His scientific activities cover many areas of plasma physics, in particular magnetohydrodynamics and reconnection theory. He is the author of the book Nonlinear Magnetohydrodynamics.
CAMBRIDGE MONOGRAPHS ON PLASMA PHYSICS General Editors: M. G. Haines, K. I. Hopcraft, I. H. Hutchinson, C. M. Surko and K. Schindler 1. D. Biskamp Nonlinear Magnetohydrodynamics 2. H. R. Griem Principles of Plasma Spectroscopy 3. D. Biskamp Magnetic Reconnection in Plasmas
Magnetic Reconnection in Plasmas Dieter Biskamp MaxPlanckInstitute for Plasma Physics, Garching
CAMBRIDGE UNIVERSITY PRESS Cambridge, New York, Melbourne, Madrid, Cape Town, Singapore, Sao Paulo Cambridge University Press The Edinburgh Building, Cambridge CB2 2RU, UK Published in the United States of America by Cambridge University Press, New York www. Cambridge. org Information on this title: www.cambridge.org/9780521582889 © Cambridge University Press 2000 This publication is in copyright. Subject to statutory exception and to the provisions of relevant collective licensing agreements, no reproduction of any part may take place without the written permission of Cambridge University Press. First published 2000 This digitally printed first paperback version 2005 A catalogue recordfor this publication is available from the British Library Library of Congress Cataloguing in Publication data Biskamp, D. Magnetic reconnection in plasmas / Dieter Biskamp. p. cm. Includes bibliographical references and index. ISBN 0 521 58288 1 (hb) 1. Magnetic reconnection. 2. Plasma (Ionized gases) I. Title. QC718.5.M3B53 2000 538\6dc21 99087680 CIP ISBN13 9780521582889 hardback ISBN10 0521582881 hardback ISBN13 9780521020367 paperback ISBN10 0521020360 paperback
Basic kinematic concepts Magnetic flux and magnetic helicity Conservation of magnetic topology Conditions for magnetic reconnection Conservation of magnetic helicity in reconnection
Current sheets Resistive MHD 3.1.1 The MHD equations 3.1.2 Equilibrium and linear modes 3.1.3 Conservation laws 3.1.4 Reduced MHD models Resistive current sheets 3.2.1 Currentsheet formation 3.2.2 Basic properties of resistive current sheets 3.2.3 Refined theory of the currentsheet structure Driven currentsheet reconnection 3.3.1 Syrovatskii's theory of current sheets 3.3.2 Dynamic structure of Y points 3.3.3 Scaling laws of stationary currentsheet reconnection 3.3.4 Petschek's slowshock model
Resistive instabilities The resistive tearing mode 4.1.1 Linear tearing instability
Contents 4.1.2 Smallamplitude nonlinear behavior 4.1.3 Saturation of the tearing mode 4.1.4 Effect of dynamic resistivity 4.1.5 Neoclassical tearing mode The double tearing mode The resistive kink mode 4.3.1 Resistive kink instability 4.3.2 Smallamplitude nonlinear evolution 4.3.3 Final state of the resistive kink mode Coalescence instability Pressuredriven instabilities 4.5.1 Interchange and ballooning modes 4.5.2 Nonlinear evolution of pressuredriven modes 4.5.3 Finite pressure effects on tearing modes Shear flow instability 4.6.1 Shear flow instability in neutral fluids 4.6.2 Effect of magnetic field on the KelvinHelmholtz instability 4.6.3 Instability of a magnetized jet Instability of a resistive current sheet 4.7.1 Threshold condition for tearing instability 4.7.2 Plasmoids
Dynamo theory Kinematic dynamo theory 5.1.1 Nonexistence theorems and special dynamo solutions 5.1.2 The rope dynamo 5.1.3 Meanfield electrodynamics 5.1.4 a2 and aQdynamos 5.1.5 Free dynamo modes Meanfield MHD dynamo theory 5.2.1 Nonlinear quenching processes and magnetic buoyancy 5.2.2 The solar dynamo MHD theory of thermal convection 5.3.1 Thermal convection in a rotating sphere 5.3.2 The selfconsistent MHD dynamo 5.3.3 The geodynamo MHD turbulence 5.4.1 hom*ogeneous MHD turbulence 5.4.2 Selective decay and energy decay laws 5.4.3 Spatial scaling properties 5.4.4 hom*ogenous turbulent dynamo
Noncollisional reconnection processes Twofluid theory High/? whistlermediated reconnection 6.2.1 The EMHD approximation 6.2.2 Properties of the reconnection region 6.2.3 Coalescence of EMHD flux bundles 6.2.4 Electron KelvinHelmholtz instability of the current layer 6.2.5 Ioncontrolled reconnection dynamics Low/? noncollisional reconnection 6.3.1 The fourfield and threefield models 6.3.2 Linear stability theory Nonlinear noncollisional kink mode 6.4.1 Electron inertiadominated reconnection 6.4.2 Kinetic Alfvenwavemediated reconnection 6.4.3 Influence of diamagnetic effects 6.4.4 Criterion for fast reconnection Sawtooth oscillation in tokamak plasmas 6.5.1 Basic experimental observations 6.5.2 The safety factor profile 6.5.3 Stabilization and onset of sawtooth oscillations 6.5.4 Observations of the collapse dynamics 6.5.5 Theoretical interpretation of the sawtooth collapse Laboratory reconnection experiments 6.6.1 The UCLA magnetic reconnection experiment 6.6.2 The PPPL magnetic reconnection experiment 6.6.3 The Tokyo University reconnection experiment
Microscopic theory of magnetic reconnection Vlasov theory and microinstabilities 7.1.1 Linear Vlasov theory 7.1.2 Quasilinear theory 7.1.3 Mode coupling, resonance broadening, and particle trapping 7.1.4 Turbulent transport coefficients Ionsound instability 7.2.1 Linear stability characteristics 7.2.2 Nonlinear saturation and longtime behavior Lowerhybriddrift instability (LHDI) Whistler anisotropy instability The collisionless tearing mode 7.5.1 Linear stability theory 7.5.2 Nonlinear saturation
Contents 7.5.3 The ion tearing mode Particle simulation of coUisionless reconnection 7.6.1 Particle simulation methods 7.6.2 CoUisionless electron dynamics 7.6.3 CoUisionless ion dynamics 7.6.4 GEM Magnetic Reconnection Challenge 7.6.5 Threedimensional simulations
Magnetospheric substorms The structure of the magnetosphere 8.1.1 The solar wind 8.1.2 Magnetopause and bow shock 8.1.3 The internal structure of the magnetosphere Magnetospheric convection Magnetopause reconnection Magnetospheric substorms 8.4.1 Substormrelated observations 8.4.2 MHD modeling of substorms 8.4.3 CoUisionless reconnection in the magnetotail 8.4.4 Nonreconnective substorm model 8.4.5 Particle acceleration by magnetotail reconnection
This book is in a sense a sequel to my previous book Nonlinear Magnetohydrodynamics, which contained a chapter on magnetic reconnection. Judging from many discussions it appeared that it was this chapter that was particularly appreciated. The plan to write a full monograph on this topic actually took a concrete shape during a stay at the National Institute for Fusion Science at Nagoya, where I found the time to work out the basic conception of the book. It became clear that resistive theory, to which most of the previous work was restricted, including that chapter of my previous book, covers only a particular aspect of this multifaceted subject and not even the most interesting one, in view of the various applications, both in fusion plasma devices and in astrophysical plasmas, where collisionless effects tend to dominate over resistivity. While resistive reconnection theory had reached a certain level of maturity and completion about a decade ago (few theories are really complete before becoming obsolete), the understanding of collisionless reconnection processes has shown a rapid development during the past five years or so. The book therefore consists of two main parts, chapters 35 deal with resistive theory, while chapters 68 give an overview of the present understanding of collisionless reconnection processes. I mainly emphasize the reconnection mechanisms, which operate under the different plasma conditions, to explain the apparent paradox that formally very weak effects in Ohm's law account for the rapid dynamic timescales suggested by the observations. Applications concern primarily astrophysical phenomena. Chapter 5 introduces dynamo theory, considering in some detail the generation of the solar and the geomagnetic field, while chapter 8 deals with magnetospheric substorms, the most important reconnection process in the Earth's magnetosphere. Both chapters are rather autonomous and can be read independently of the remainder of the book. Concerning laboratory plas
mas I resume the discussion of the sawtooth phenomenon considered in some length in Nonlinear Magnetohydrodynamics, giving an update of the experimental and theoretical situation. I also discuss briefly several laboratory experiments designed specifically to study magnetic reconnection physics. It is a pleasure to express my gratitude to the many colleagues with whom I enjoyed fruitful and illuminating discussions on the topics of this book, in particular Jim Drake, who taught me the importance of coUisionless reconnection, and Wolfgang Baumjohann, Michael Hesse and Manfred Scholer for introducing me to the realm of magnetospheric physics. I also acknowledge the financial support by the COE programme of Monbusho and the kind hospitality of the National Institute for Fusion Science with special thanks to Tetsuya Sato. Finally I would like to thank Brian Watts for his painstaking copyediting of the manuscript.
Since the early 1950s, when magnetohydrodynamics  MHD in short became an established theory and along with it the concept of a "frozenin" magnetic field within an electrically conducting fluid, the problem of how magnetic field energy could be released in such a fluid has been generally acknowledged. In the early days the major impetus came from solar physics. Estimates readily showed that the energies associated with eruptive processes, notably flares, can only be stored in the coronal magnetic field, all other energy sources being by far too weak. On the other hand the high temperature in the corona, which makes the coronal plasma a particularly good electrical conductor, appeared to preclude any fast magnetic change involving diffusion. For a coronal electron temperature Te ~ 106 K the magnetic diffusivity is rj ~ 10 4 cm 2 /s, hence field diffusion in a region of diameter L ~ 104 km as typically involved in a flare would require a timescale zn = L2/r] ~ 1014 s, whereas the observed flash phase of a flare takes less than ~ 103 s. It had, however, soon been realized that the discrepancy is not quite as bad as this. Contrary to magnetic diffusion in a solid conductor, a fluid is stirred into motion by the change of the magnetic field. As it carries along the frozenin field, it may generate steep field gradients typically located in sheetlike structures, and hence lead to much shorter diffusion times. By the early 1960s the concept of resistive current sheets, or SweetParker sheets named after their major protagonists, had become standard knowledge. The transverse scale of such a sheet is 5n ~ ( T ^ ) 1 / 2 , where the Alfven time XA is a typical MHD timescale. Though magnetic diffusion across the sheet is very fast, 82/rj ~ XA ~ 1 s, the overall dynamics is now limited by the rate of convective field transport toward the sheet. A simple backoftheenvelope analysis shows that the time for magnetic energy release in the presence of a SweetParker sheet is TSP ~ (^A^rj)1^2 ~ 10 7 s, which is much shorter than the global dif1
fusion time of 10 14 s, but still many orders of magnitude longer than observed. A breakthrough occurred in 1964 at a symposium on the physics of solar flares, when Petschek presented a mechanism which appeared to allow very fast energy release, depending only logarithmically on rj. The model was soon generally accepted, forming the basis for the theory during the following two decades, during which more and more refined configurations were proposed. However Petschek's model is, properly speaking, not a genuine theory of fast reconnection, but the MHD configuration set up in the presence of an efficient reconnection mechanism. Petschek and his successors assumed, explicitly or tacitly, that resistivity provides such a mechanism. It was only during the 1980s it became gradually clear that Petschek's configuration is not valid in resistive MHD where, instead, extended current sheets are formed. It thus appeared that theory was thrown back to the prePetschek period of relatively slow SweetParker sheet reconnection. The situation seemed to converge with the school of thought that had developed independently in the eastern hemisphere starting from Syrovatskii's elegant theory of current sheets. A possible way to reach faster reconnection speeds was to give up the usual assumption of quasistationarity by allowing MHD turbulence to develop. In many applications, however, such macroscopic turbulence is not observed. A different approach to the problem was followed in magnetospheric physics. Except for the ionosphere, the magnetospheric plasma is so dilute that Coulomb collisions are, practically speaking, absent and hence classical resistivity vanishes. Can magnetic reconnection occur in a collisionless plasma? The problem had already been realized by Dungey in the early 1960s in his theory of magnetospheric convection. The usual approach was to assume the action of smallscale turbulence excited by some microinstability, where the scattering of electrons by charge fluctuations has a similar effect as Coulomb collisions. Various candidates for such effective or anomalous resistivity were discussed, though none was found to be fully satisfactory to account for reconnection in the magnetotail. The question of fast quasicollisionless reconnection became also relevant for hot fusion plasmas, in particular in tokamaks, where much theoretical effort was made to explain the characteristics of the sawtooth phenomenon, an internal relaxation oscillation. It was realized that in Ohm's law nonideal terms other than resistivity play the dominant role. Combining the results from both magnetospheric and tokamak plasma theory achieved in recent years we now seem to be close to a consistent picture of collisionless reconnection. Such processes seem to be significantly more efficient than resistive diffusion and allow fast quasiAlfvenic reconnection velocities.
Fig. 1.1. Elementary picture of field line reconnection: relaxation of magnetic tension and plasma acceleration by a local change of field line connectivity. The reader may have noticed that the term magnetic reconnection has slipped into our vocabulary without introduction. There is hardly a term in plasma physics which exhibits more scents, facets and also ambiguities, and which at times seems to be used even with a touch of magic. The basic picture is that of two field lines being frozen in and carried along with the fluid, until they come close to each other at some point where, due to weak nonideal effects in Ohm's law, they are cut and reconnected in a different way. (The alternative term "field line merging" with essentially the same meaning, which was popular in the 1970s, is gradually falling from use.) Simple as it may look this picture already contains important characteristic features of the reconnection process. Consider the field line configuration shown in fig. 1.1 where the upper and lower endpoints move around such as, for instance, the footpoints of coronal fields anchored in the photosphere. Physically the close encounter of the field lines implies that magnetic field gradients become locally strong, thus enhancing the formally weak nonideal process in Ohm's law. Hence reconnection is a localized process, in contrast to the slow global field diffusion and dissipation by the omnipresent, but weak, resistivity. The local enhancement of the nonideal effects determines the timescale of the reconnection process, while the overall magnetic energy dissipation remains small owing to the localization of the process. The crucial point to note is that after local rearrangement of the magnetic connectivity, the global magnetic configuration may relax by reducing magnetic tension and accelerating the fluid in the way indicated in fig. 1.1, often called "the magnetic slingshot effect". Hence reconnection allows rapid conversion of magnetic into kinetic energy, which is the essence of impulsive events
in magnetized plasmas. In fact the observation of fast collimated plasma flow is regarded as the main indication of the presence of magnetic reconnection, for instance during the emergence of new flux in the solar corona or the reconnection at the Earth's magnetopause. There are two different aspects when considering such reconnection events. The first one concerns the global properties of the system. The magnetic configuration determines the magnitude of the event, i.e., the energy, which can be released. This is related to the complexity of the magnetic field  the twist and shear of the field lines  which by reconnection relaxes into a structurally simpler state. It also determines the location of the reconnection process, which is not simple to predict in a complicated threedimensional configuration. The second aspect refers to the local properties in the reconnection region  the physics of the actual reconnection mechanism  to which the major part of the book is devoted. Here one considers: (a) Timescale T. It depends of course on the free magnetic energy, being the shorter the larger the latter, simply because higher velocities are generated. But x is also expected to depend more or less strongly on the nonideal process R in Ohm's law. More specifically T ~ L/v, where L is the spatial extent of the configuration which is affected by the relaxation and v some average velocity generated during the process. According to the rule given above, we have v oc VA, where vA = SB/^J4np is the Alfven velocity corresponding to the change 5B of the magnetic field. Hence we can write T ~ [L/vA]f(R) = TAf(R).
The function f(R) represents the influence of the particular reconnection mechanism, which dominates for the given plasma parameters, for instance / = ( T ^ / T ^ ) 1 / 2 for the SweetParker process. Determination or at least estimation of f(R) is a major topic of this book. (b) Energy partition. The magnetic energy released by reconnection in an eruptive event is transformed into (i) bulk plasma motion often generating a strong shock, the blast wave, (ii) electron and ion heating and (iii) acceleration of a certain number of ions and electrons to high superthermal energies. Energy partition depends on the reconnection process. Until recently this point has received little attention in reconnection theory. A fundamental question concerns the efficiency of particle acceleration, in particular whether the observed high particle energies are primarily generated at the reconnection site or by acceleration at the shock.
(c) Threshold conditions. In general, a certain amount of free energy has to be accumulated before rapid relaxation sets in. A particular facet is to explain the trigger mechanism responsible for the sudden onset of energy release often observed. The book consists of three parts: an introductory one, chapter 2, where the more formal, conceptual aspects of magnetic reconnection are discussed; a second part, chapters 35, which deals with resistive theory; and a third part, chapters 68, which treats collisionless reconnection. While the first two chapters in each of the two main parts are devoted to the basic theory, the last chapters introduce a particular application  dynamo theory and magnetospheric substorms, respectively. Chapter 2 gives an introduction to the more formal aspects of reconnection theory. These are basic kinematic concepts referring to the topological properties of the magnetic field, where the dynamics of the embedding conducting fluid need not to be specified. Only Ohm's law and Faraday's law enter this discussion, and give rise to two conservation laws: the conservation of magnetic flux, from which derives in particular the physical meaning of field lines as infinitesimally thin flux tubes; and the conservation of magnetic helicity, which is a measure of the structural complexity of the configuration, the interconnection of the individual field lines. We then define magnetic topology and derive a criterion for topology conservation, which is related to but more general than magnetic flux conservation. A question which has evoked considerable discussion concerns the definition of reconnection; in particular, what are the criteria to decide whether, in an evolving magnetic configuration, reconnection actually takes place. While this point can easily be decided for a symmetric configuration where a separatrix is defined, general threedimensional systems require individual consideration which refers to the local field line connectivity rather than their global topology. Finally, the problem of helicity conservation in reconnection is considered. We define the concepts of twist, knottedness and linkage and show that, while the helicities of the constituent flux tubes lose their meaning, the total helicity is in fact conserved, which is the basis of Taylor's theory of relaxation to a linear forcefree state. Chapters 3 and 4 deal with resistive reconnection theory: chapter 3 is focused on the characteristic stationary configuration, the resistive current sheet; chapter 4 considers the dynamics arising in resistive MHD, instabilities and their nonlinear evolution. We start, in chapter 3, with a cursory derivation and discussion of the MHD equations and their conservation properties and introduce a reduced set of MHD equations which, by eliminating the fast compressional effects, is particularly suited for treating the typically slower reconnection processes. Resistive recon
nection occurs in current sheets. We analyze the selfsimilar process of sheet formation and the structure of the stationary resistive sheets. The properties of currentsheet reconnection are conveniently studied in the form of driven reconnection, i.e., in an open system of a configuration with reversing magnetic field and boundary conditions enforcing stationary inflow and outflow of plasma and magnetic flux. A simple and elegant analytical approach has been developed by Syrovatskii, where current sheets appear as branch cuts of a complex function. This theory is confirmed by numerical solution of the resistive MHD equations, which yields the currentsheet scaling properties in terms of the resistivity and the (enforced) reconnection rate. The simulations demonstrate that, in the framework of resistive MHD, Petschek's slowshock solution is not valid. Petschek's configuration, which we briefly discuss, is only set up in the presence of a reconnection mechanism more efficient than resistivity. In chapter 4 we discuss in some detail the various types of instabilities which occur in a resistive MHD system driven by gradients of current density, pressure or velocity. For each mode we first outline the linear stability characteristics and subsequently study its observable effects, i.e., nonlinear evolution and saturation. The prototype of a currentdriven resistive instability is the tearing mode, which leads to a configuration containing magnetic islands. Though the growth rate is rather small, the perturbation may reach a substantial amplitude or island width. Analytical expressions are derived for the nonlinear growth and for the saturation level. Nonlinear evolution is actually a slow magnetic diffusion process, where no current sheet is formed. Hence not only the value of the resistivity at the reconnection point, but also the resistivity profile in the island affect the nonlinear properties. The neoclassical tearing mode, which may be destabilized in a toroidal plasma at low collisionality, is due to a similar nonlocal effect. There are other paradigms of currentdriven instabilities, which give rise to faster dynamics and thus to current sheets. In the double tearing mode, magnetic islands grow at neighboring resonant surfaces mutually enforcing the reconnection process. In a plasma column we find the resistive kink mode, a quasirigid helical shift of the central part of the column, which is energetically favored. Here a single magnetic island boosts itself due to the geometry of the configuration. Finally, two magnetic islands or flux bundles having parallel currents attract each other and coalesce, their magnetic fields reconnecting across a current sheet in between. Pressuredriven modes arise due to field line curvature, which is equivalent to a gravitational force and may thus drive a RayleighTaylortype instability. Here finite resistivity can destabilize modes by weakening the effect of magnetic shear. Nonlinearly, however, the mode is selffocusing, thus essentially avoiding the stabilizing shear effect. Thus reconnection
is not required, which makes the nonlinear evolution of this socalled ballooning instability particularly violent. A third source of free energy is a sheared flow driving the KelvinHelmholtz instability. Intense shear flows occur along current sheets, but the KelvinHelmholtz instability is, in general, stabilized by the parallel magnetic field adjacent to the sheet. The current sheet may, however, break up due to the tearing mode, if the sheet is sufficiently thin. Chapter 5 gives an introduction to dynamo theory. Devoting a full chapter to this topic may at first sight appear somewhat surprising, since the dynamo effect  the amplification and sustainment of magnetic fields by motions in a conducting fluid  is usually not directly associated with a reconnection process. In particular the basic model of field amplification by stretching, twisting and folding of flux tubes does not involve reconnection. However, an important aspect of the dynamo problem is to account for largescale fields, which have to be generated by reconnecting the naturally amplified smallscale flux tubes. We start with kinematic dynamo theory, where the nonlinear reaction of the field through the Lorentz force is neglected. Here the convenient tool is the meanfield approximation, which allows a linear analysis in terms of dynamo modes for a phenomenological description of the magnetic fields in astrophysical objects such as stars and galaxies. We briefly review the theory of the solar dynamo, which has again become an open subject since recent observations, especially from helioseismology, have shattered what was believed to be a rather solid building. The fundamental process of planetary and stellar dynamo action is convection in a rotating fluid shell. Numerical computations of model systems have demonstrated the basic mechanism. Though direct dynamical simulations of the Earth's liquid core are out of reach, simulations using phenomenological eddy diffusion coefficients have already revealed important properties of the geodynamo. The convection driving dynamo action is highly turbulent in most systems of interest. Chapter 5 therefore presents also a brief overview of MHD turbulence and turbulent magnetoconvection. I concentrate on recent developments obtained mainly from highresolution threedimensional numerical simulations, while referring for a broader introduction to chapter 7 in my book Nonlinear Magnetohydrodynamics. Conventional theory has considered resistivity as the natural, if not the only possible, mechanism for magnetic reconnection. As a result, in weakly collisional plasmas such as in the solar corona or in a tokamak discharge the discrepancy between observed and predicted timescales is particularly large caused by the inherent inefficiency of the SweetParker process. In reality, however, in such systems noncollisional effects in Ohm's law are more important than resistivity, which no longer controls the reconnection
dynamics. The simplest framework to deal with these effects is a twofluid model, generalizing the onefluid MHD theory. Since in this approach the presence of some collisions is still required to justify the use of isotropic pressures and to smooth weak flow singularities, we call reconnection in this regime, which is treated in chapter 6, noncollisional instead of collisionless. It is useful to distinguish between high and low/? systems. In the high/? case, ions and electrons decouple at scales below the ion inertia scale c/copi, where the dynamics is described by the approximation of electron magnetohydrodynamics (EMHD). It is derived analytically and verified by numerical simulations, that the reconnection rate does not depend on the electron scalelength c/cope. Also, the dependence on the ion scalelength c/copt is found to be weak, no macroscopic current sheet being formed. In the low/? case, the strong guide field tightly couples electron and ion motions, apart from diamagnetic drifts, which give rise to the smallness parameter p s , the effective ion Larmor radius. As in the high/? case, there is no macroscopic current sheet and reconnection is quasiAlfvenic. The sawtooth oscillation in a tokamak plasma is a classical reconnection paradigm. I give a brief review of the present understanding of this phenomenon in chapter 6, again referring for a more detailed introduction to chapter 8 of my previous book. For truly collisionless plasmas, such as in the solar wind, a kinetic description is required, which is presented in chapter 7. Since the usual approach to reconnection phenomena occurs from the fluid side, I presume that many readers are less familiar with the kinetic plasma aspects and will appreciate a somewhat broader introduction to microscopic theory. We therefore start the chapter with an outline of linear Vlasov theory and the different nonlinear effects, quasilinear relaxation, mode coupling and particle trapping. We then discuss several microinstabilities and their nonlinear properties, which may give rise to anomalous transport effects in Ohm's law: the ionsound instability; the lowerhybriddrift instability; and the anisotropydriven whistler instability. The collisionless tearing mode in a neutral sheet has played a special role in collisionless reconnection theory, in particular in the presence of a weak normal field component. After a longstanding debate it is now generally believed that the mode is stable, so that the breakup of such a collisionless configuration requires a finite perturbation generating an Xpoint. Since collisionless reconnection is studied by particle simulation, we briefly describe the different numerical approaches. Simulation results indicate that, different from the resistive regime and similar to the behavior in the high/? noncollisional approximation, reconnection rates are independent of the electron physics and do not seem to depend on ion inertia, either, i.e., reconnection speeds scale with the Alfven velocity.
Probably the most important application of collisionless reconnection is the magnetospheric substorm, which is considered in chapter 8. We first briefly describe the structure of the magnetosphere, which results from the interaction of the solar wind with the geomagnetic field. Dungey's model of magnetospheric convection gives a qualitative account of the reconnection processes at the front and in the tail and their dependence on the orientation of the interplanetary magnetic field. Observations indicate that tail reconnection occurs mainly in the form of a relaxation oscillation, a sequence of substorms. Combining observational and numerical results, the most probable substorm model is a fast reconnection of the tail magnetic field generating a plasmoid which is ejected tailward, though a nonreconnective process involving a ballooningtype instability cannot be excluded. Finally, the mechanism of particle acceleration in the magnetotail by the reconnection electric field is outlined. Concerning notation and coordinate systems, I have tried to compromise between the demands of uniformity throughout the book and the use of notation which has become standard in a particular field. To a large extent the individual chapters are autonomous, but a certain amount of crossreferencing between chapters is certainly useful. The equations are primarily written in Gaussian units which, in my view, are more natural in collisionless theory than are SI units. But, as soon as the basic equations are specified, I introduce convenient normalizations depending on the specific problem to write the equations in dimensionless form.
In this chapter we introduce some formal properties of the magnetic field embedded in an electrically conducting fluid, which play a fundamental role in reconnection theory. Conservation of magnetic flux leads to the concept of magnetic field lines as infinitesimally thin flux tubes, which are frozen in to the fluid. A slightly different concept is that of magnetic topology or connectivity. Since this is a property of the magnetic field alone, topology conservation is valid under more general conditions than flux conservation. We then try to give a comprehensive definition of reconnection. While, in two dimensions, reconnection can be defined globally in terms of plasma flows across a separatrix surface, such surfaces, separating regions of different global fieldline topology, no longer exist in a general threedimensional configuration. Here, reconnection can only be defined by a local change of fieldline connectivity, which is due to a nonvanishing parallel electric field. Finally we consider the behavior of magnetic helicity under the influence of reconnection. Though the helicities of the individual flux tubes involved in the reconnection process lose their meaning, the total helicity of the system is conserved. 2.1 Magnetic flux and magnetic helicity The coupling between an electrically conducting fluid and the magnetic field immersed in the fluid is described by the generalized Ohm's law E + V x B = R, (2.1) c where v(x, t) is the fluid velocity and R comprises the different socalled nonideal effects of the plasma, dissipative or nondissipative, which are discussed in more detail in section 6.1. Here, we only note that in most cases of interest R is formally very small. The left side of (2.1) is the 10
electric field E' in the rest frame of the fluid element. Neglecting R simply means that E' = 0, which is called the ideal Ohm's law E+
x B = 0.
c (We shall see, in section 6.1, that Ohm's law is just the equation of motion of the electron fluid.) Insertion into Faraday's law ^B =  c V x E
which has the form of a continuity equation for the vector field B. In integral form (2.4) states the conservation of the magnetic flux
through an arbitrary surface F(t) bounded by a closed curve l(t) moving with the fluid. Taking the surface integral of (2.4) one obtains
Fu(x,to,to) = x.
The flow is topology conserving if the points that form a field line at time to remain on the same line for any t > to. Consider two points x,y on a field line at time to, x = F#(x,O,£o)? y = F 5 (x,s,£o), which are transported * To simplify notations, B has the dimension of a velocity, i.e., should be considered as the Alfven velocity v^ = B/^/4np, p— mass density.
.
Fig. 2.1. Field line conservation of a vector field B mapped by the flow function FM (from Hornig & Schindler, 1996). by u to x', y' at time t, hence x' = FM(x, t, to), y7 = FM(y, r, to) or, inserting
yf = Fu[FB(x,s,to),t,to].
The condition for topology conservation is that x', y are again connected by a field line, i.e., by the mapping y! = FB(x\s\t) or, inserting x', yf
FB [FM(x, t, ^o)? s\ t\.
(Note that sf may differ from s, since parallel flows are admitted, hence sf(s, t) is some arbitrary function.) From (2.19) and (2.20) we obtain the commutation relation (seefig.2.1) =
FB[Fu(x,t,to\s\t].
Equation (2.21) can be transformed into a more familiar differential form. Differentiation with respect to both s and t tFu
taken at s = 0, t = to, and use of the definitions of F B , F M , (2.17) and (2.18), gives x,t,to),tW}
\t=t0.
Carrying out the differentiations, using dsii(F) = 3SF Vu = B Vu, yields B Vu = dtB + u VB + B dtds^.
Hence we obtain the condition dtB + u VB  B Vu = Bp.
Fig. 2.3. Plasmoid formation in a weakly threedimensional configuration. Unfortunately this concept holds only in the strictly 2D situation, since the separatrix is structurally unstable to slight deviations from 2D geometry. Consider the generalization of the process shown in fig. 2.2 to a plasmoid of finite axial extent in the presence of a weak axial field, fig. 2.3. Obviously reconnection takes place in a way similar to that in the 2D case, but we can no longer distinguish different global field line topologies, all field lines coming from and ending up on the left part of the configuration. There is no separatrix surface anymore. An equivalent of the 2D separatrix for the 3D case can be defined in the presence of points where B = 0, called magnetic nulls. These isolated points are structurally stable, a slight change of the configuration just shifts them by a short amount. Magnetic nulls have therefore attracted some attention (f*ckao et a/, 1975; Green, 1988; Lau & Finn, 1990). Let us briefly discuss the topological structure of the field in the presence of nulls. In the vicinity of a null, which we assume to be located at x = 0, B has the form
Because of V • B = 0 the trace of the (real) matrix ay vanishes. The matrix has three eigenvalues, which we can use to classify the different types of nulls. First, consider the case where all eigenvalues are real. Because of the trace condition either two are positive and one negative, which we call ^4type, or two are negative and one positive, which we call Btype. The eigenvectors of the two equalsign eigenvalues span a plane locally, which can be continued into a global surface 2 by following the field lines in this plane. Following the field lines along the third eigenvector defines a curve y in space. The field lines near an ^4type null are shown in fig. 2.4, the surface HA dividing space into two regions not connected by field lines, all field lines in one region converging into a bundle around y&. A Btype null has the same topology with the field line directions reversed. A and I?type nulls are generalizations of a twodimensional Xpoint, to
Fig. 2.4. Field lines near an ,4type magnetic null (after Lau & Finn, 1990).
Fig. 2.5. Topology of a magnetic configuration including one v4type and one Btype null (after Lau & Finn, 1990). which they degenerate if one of the two equalsign eigenvalues vanishes. These nulls are therefore expected to be important in the context of reconnection. Nulls with one real and two conjugate complex eigenvalues are called Stype. They are generalizations of a twodimensional 0point, to which they degenerate if the real eigenvalue vanishes (in this case the remaining two are purely imaginary). A magnetic configuration with two nulls, one .4type and one 5type, is of special interest. The relative positions of the surfaces ^A^B and the curves yA,7B are illustrated in fig. 2.5. Since in £# all field lines are directed away from B and are collected in the bundle around JA, the line an JA bounds YJB> d conversely all field lines in Z^ must originate in the bundle y#, i.e., JB bounds Z^. Hence Z^,Sj5 are semiinfinite halfsheets intersecting in a line connecting A and B, called a nullnull line. A topologically equivalent configuration arises when superimposing a constant field on a dipole field, for instance the Earth's magnetic field together with the interplanetary
field. Since the surfaces E^, Z^ are not pierced by field lines they serve as separators, the threedimensional generalization of the twodimensional separatrix, dividing space into regions of topologically different field lines. In this case, Vasyliunas' definition of reconnection is still applicable. Studies of the threedimensional dynamics of interacting twisted flux tubes performed by Lau & Finn (1996) confirm this concept of reconnection across a separator. In the case of two tubes with antiparallel field and opposite twist, i.e., parallel current, the evolution proceeds in two phases. In the first phase, the (transient) presence of two magnetic nulls leads to formation of a closed Xline (nullnull line) encircling the flux tubes. In the second phase, field lines from both tubes are reconnected, forming loops thereafter. However, few magnetic configurations in laboratory and space do contain magnetic nulls. Moreover, numerical simulations of the dynamics of general 3D fields with nulls (Politano et a/., 1995) show that magnetic nulls do not seem to play a special role. The typical loci of reconnection characterized by sheetlike structures of quasisingular current density, are found to be uncorrelated to the presence of nulls, the structures being locally similar to that in 2D discussed in detail in chapter 3. Hence reconnection occurs typically at finite B. One can therefore in general not refer to the global field line behavior to define reconnection as done in Vasyliunas' criterion. Instead, one has to consider directly the local breakdown of flux conservation resulting in a change of the field line connection, as suggested by Axford (1984), i.e., a violation of the topology conservation discussed in section 2.2. Reconnection occurs if there are points which are originally located on the same field line but end up on different field lines. This concept has been elaborated by Schindler et al (1988) and Hesse & Schindler (1988). It clearly applies also to the 3D plasmoid formation illustrated in fig. 2.3. (Since the effect of reconnection should be global, we exclude the case of a purely local change of the magnetic configuration by requiring that the points considered remain outside the small diffusion region, where R is important). A formal criterion for reconnection is obtained from (2.30): B x (V x S) ^ 0.
Here S = R\\ B/B is the parallel component of the nonideal part in Ohm's law, since the perpendicular component can always be incorporated in the general velocity u in (2.29) by writing R i = u ' x B , u = v  u / . Hence only R\\ = E\\ can cause reconnection. But not all R\\ actually do so. Clearly an electrostatic field does not, since for S = V(/> the left side of (2.35) vanishes. Even in the standard resistive case R = f/j (see section 3.1), topology is conserved for special fields, for instance for a linear forcefree field j = /iB, \i = const, and uniform resistivity. Practically speaking these
exceptions are irrelevant because of the selfconsistent dynamics of v in the diffusion region, as is discussed in subsequent chapters. One should also note that although R± does not directly lead to reconnection, it may strongly increase the reconnection efficiency, for instance the Hall term in high/? noncollisional reconnection considered in section 6.2. 2.4 Conservation of magnetic helicity in reconnection As shown in section 2.1 magnetic helicity is conserved for idealflows,i.e., processes where the ideal Ohm's law (2.2) holds. For more general motions, H is not strictly conserved, not even in the general topologyconserving case (2.25). Straightforward evaluation gives, assuming for simplicity that the system is bounded by a perfectly conducting vessel for which the tangential electric field component vanishes, Et = 0, f
which vanishes only in the ideal case a = Vv. However the helicity change is in general very slow, of the order of the global diffusion rate. Using (2.11), again for Et = 0 on the boundary, m
we find that dH/dt is small, since E\\ = R\\ is usually small and may at most be finite in the small regions where reconnection takes place. By contrast the total energy, which is also ideally conserved, changes more rapidly. Neglecting viscosity and pressure (low/? approximation) we have
(where p = the mass density), which may be larger than O(R), since j may become large locally. (A more detailed discussion of the energy balance is given in section 3.1.) Using the disparity of relaxation rates, Taylor (1974, 1986) has conjectured (see also Berger, 1984) that a magnetic configuration relaxes through a dynamic, generally turbulent phase to the lowest energy state consistent with the conservation of the initial value of the magnetic helicity. The final state is obtained from the variational principle S[WocH]=S
where a is a Lagrange multiplier. Variation with respect to A gives the equation V x B = 87iaB, (2.40) while variation with respect to v gives v = 0. Hence the final state is a static constanta forcefree magnetic field, where a is determined by the value of the helicity. In general, the relaxation process requires reconnection of field lines resulting in a simpler, less complex magnetic field structure. Since magnetic helicity is a measure of the linkedness and knottedness of the magnetic field, breaking such linkages seems to involve also a change of helicity. We will show or at least make plausible that while the individual helicities of the constituent flux ropes become meaningless, since the latter are merged and reorganized and hence lose their identity, the sum of their helicities, nevertheless, is conserved. To analyze the topological structure of a magnetic configuration (to be definite we restrict consideration to closed configurations with a perfectly conducting boundary surface) it is useful to divide the volume into small subregions bounded by flux surfaces, the simplest case being that of a toroidal configuration with nested flux surfaces or magnetic surfaces, as they are usually called. Such a surface is generated by following a field line around the torus. In general, field lines are infinitely long such that a single line will span the entire surface, which is called an irrational surface. Intermixed are special surfaces, called rational surfaces, where field lines are closed after a certain number of turns. Such surfaces are simply defined as limiting cases of the neighboring irrational ones. The subregions we consider are thin toroidal shells between two magnetic surfaces. The helicity of such a shell is AH
Here § ds and / d\ indicate the line integrals the short way and the long way around the torus, d¥s and d¥i are the surface elements with normals along these directions, — A\pp and A\pt are the poloidal and toroidal fluxes in the shell. xpt = § Ap  ds is the toroidal flux in the torus volume bounded by the shell and \pp = § At • d\ is the poloidal flux threading the hole of the shell surface. The geometry is illustrated in fig. 2.6. Note that increasing the radius of a toroidal flux surface xpt increases, while xpp decreases, whence the minus sign in the first term of (2.41). There is a simple interpretation
W= const.
Fig. 2.6. Toroidal geometry defining the surface elements d¥s and of the two terms in AH. The first term accounts for the linkage between the flux in the shell volume AV with the fields inside the shell surface, the second one for the linkage with the fields outside. The helicity of the torus can be written in the following way r
with T = —dxpp/dxpt and 0 = total toroidal flux, which is obtained by partial integration of the second term in (2.41) ft f° dxpt / xppdxpt= VP^—d\pp = Jo J(j>p dxpp "
f° f^ xptdxpp= T\ptdxpt. J^p * Jo
Here the boundary contributions vanish, since at the lower boundary xpt vanishes, at the upper xpp, the field outside the torus being assumed zero. The quantity T is the (in general fractional) number of times a field line winds around the torus in the poloidal direction for one turn in the toroidal direction. In the theory of toroidal MHD equilibria T is usually called the rotational transform denoted by i (pronounced iota), while in tokamak theory the inverse T~l is called the safety factor q. For a uniformly twisted torus with twist T = const we have the helicity H = Tcj)2
Fig. 2.7. Measuring helicity. A crossover in (a) gives a unit of positive helicity, in (b) a unit of negative helicity.
If the torus crosssection is sufficiently small, we are dealing with a thin closed flux tube or a flux rope. In this case the variation of T in the tube becomes unimportant and we can use expression (2.44) with T representing the average twist. Let us now discuss the connection between the internal fieldline twist T and the external distortion or kink K of the flux tube. A plane closed flux tube with twist T = +1 may be kinked into a figure8 configuration with no internal twist. To see this, unkink an untwisted figure8 tube (as for instance shown in fig. 2.7) with a field line drawn on its upper side into a plane torus to find that this field line now winds around the torus just once corresponding to T = +1 depending on the orientation of the field with respect to the kink. Hence we see that the helicity of a flux tube consists of an internal part Hj and an external one HK, H =
The decomposition is, however, not topologically invariant since, as we have just seen, one can be converted into the other. Let us assume that an arbitrarily kinked and knotted flux tube has zero twist, HT = 0, if one can draw a field line on top of the tube* along its entire length. There is a simple way of measuring the helicity of an untwisted tube. We have seen that untwisting a tube by kinking it into a figure8 shape * Strictly speaking on the plane projection of the tube. To simplify the discussion, we avoid the complicated discussion of nonplanar flux tubes, see, e.g., Moffatt & Ricca (1992).
Fig. 2.8. Helicity of a trefoil knot (after Berger & Field, 1984).
Fig. 2.9. Two linked flux tubes. introduces a crossover. There are two types of such crossovers as indicated infig.2.7, which we call positive,fig.2.7(a), and negative,fig.2.7(b). Then the helicity of such a tube is H = HK=$2(N+N),
where N+ and N are the numbers of positive and negative crossovers. Consider for instance in untwisted trefoil knot, fig. 2.8(a), which according to relation (2.46) has helicity H = 30 2 . This can be proved by reconnecting the knot in the way indicated infig.2.8(b,c). The final state consists of three kinked tubes and two unkinked ones, all nontwisted. While the former ones contribute 4>2 each, the helicity contributions of the latter are, of course, zero. However, how do we know that the reconnection does not change the helicity? We postpone the answer to this question for a while and introduce first the concept of linkedness. Consider two flux tubes with fluxes 0i, (/>2, interlinked as shown infig.2.9, untwisted for simplicity. The total helicity can easily be computed as
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 SweetParker 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 socalled Ypoints. We then discuss the scaling laws of stationary currentsheet 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 twofluid 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, pr = ~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 timescales are usually long compared to the dynamics described by the MHD equation (3.1), such that the system changes quasistatically 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: VPB
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. Onedimensional 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 fieldline 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 reversedfield 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 kv = 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\ + <W4o\ ((02/k)2 > max{i;,4,cs}. Here, the lefthand 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 (lefthand 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 quasistatic change of equilibrium. For parallel propagation (righthand 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 superAlvenic plasma stream such as the Earth's bow shock. In reconnection theory, slowmode 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 ~ pvVv, we find
Vt,~4f
c2s L where L is a typical gradient scalelength. 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 quasiincompressible, 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#~)'dFri
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 frozenin 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 crosshelicity K = f \BdV. IV(t) Jv(t)
(3.33)
The crosshelicity 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: •+VVV>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 quasineutrality 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 quasiincompressibility. 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
nc^.
(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=
foA\
fl = meA\
(3.86)
corresponding to exponential flattening of the Xpoint 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 onedimensional, 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 sheetlike structure only the crosssheet 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 finitetime 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 xderivative 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., threedimensional Euler flows, FTS seems to be generic (except for highsymmetry 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 sheetlike eventually, even if initially the magnetic field is weak (Grauer & Marliani, 2000). The absence of a finitetime singularity seems to weaken the relevance of the similarity solution (3.86) since, for long times, spatial regions far away from the Xpoint, 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 nonuniversal effects remain small. Though the similarity solutions given above provide strong evidence for rapid currentsheet formation, their validity is restricted to the center region of the sheet, the immediate vicinity of the Xpoint, 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 OrszagTang vortex (Orszag & Tang, 1979). Four current sheets develop at the Xpoints 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 upperright 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 
(395)
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 = CQ2V/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) 7oo 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 quasistationary resistive current sheet, called a SweetParker 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 quasionedimensional 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 threedimensional 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 selfsimilar
53
3.2 Resistive current sheets
Fig. 3.5. Schematic drawing of a SweetParker 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 xaxis. 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 midplane 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 onedimensionality 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 lengthscale L. In this case Ohm's law,
E  —xB = r\B\ JLJ
B = \p\x\
3 Current sheets
58
Fig. 3.6. Theoretical currentsheet profiles: (3.125) (full), seen2* (dashed), Gaussian (dotdashed).
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 selfsimilar 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 selfsimilar behavior in the ideal phase of currentsheet 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 currentsheet reconnection
59
3.3 Driven currentsheet 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 welldefined 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 Lsystem to the Asystem occurs through the boundary conditions imposed on the Lsystem, which change on the global timescale ~ 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 twodimensional 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 frozenin 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 frozenin 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 currentsheet 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 CauchyRiemann 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 secondorder or Xtype neutral points. (Note that there is no nonsingular Otype 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 Xpoints, as illustrated in fig. 3.7. The change of topology occurs at these secondary Xpoints, 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 Xpoint. (a) Initial nonsingular configuration; (b) effect of an induced singular current line in the original Xpoint, leading to a fictitious Opoint and two adjacent Xpoints. 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 Xpoints arising by adding the line current to the original Xpoint, (3.133). As indicated in fig. 3.7(b), the cut passes through this point connecting the new Xpoints, which drift further apart as I(t) increases. Addition of the current sheet splits the original Xpoint into two Ypoints. 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 yaxis 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) =laz
(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 currentsheet
(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 currentsheet 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 yaxis, while the current sheet continues along the yaxis, 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 frozenin 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(09y) =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 Xpoint 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 Xline 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 Xline and hence to a current sheet, an example being given in fig. 3.9. Let us briefly discuss the case of a fully threedimensional magnetic field. Imagine a configuration characterized by two magnetic nulls, one xp =
3.3 Driven currentsheet
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).
^4type and one 5type (see fig. 2.5). If there is a nonvanishing electric field E along the nullnull line AB, a current sheet is expected to form. It appears all we have to do is to split the nullnull 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 threedimensional 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 ydirection, (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 halfribbons, 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 ,4type nulls and B\ B" being Btype 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 Ypoints We have seen in the preceding section that a singular Xpoint, 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 superAlfvenic 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 subAlfvenic speed. In fact, simulations show that the flow in the downstream region is clearly subAlfvenic 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 yaxis 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 currentsheet 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 quasisingular 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 jplots) 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 currentsheet 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 flowblocking shocklike structure. In the limit *;—•() a hierarchy of higherorder current sheets is generated with a selfsimilar scaling behavior, as drawn schematically in fig. 3.13.
3.3 Driven currentsheet reconnection
69
Fig. 3.13. Schematic drawing of the selfsimilar hierarchy of current sheets in the diffusion layer: thick line = positive j ; thin line = negative j ; open circle = quasisingularity at sheet endpoint. These properties of the diffusion layer are consistent with Syrovatskii's currentsheet model, the multicurrentsheet edge behavior representing the dynamically resolved singularity predicted in Syrovatskii's quasistatic 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 nonstationary 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 currentsheet 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 nonseparable, i.e., fully twodimensional.
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 updown and rightleft 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 xaxis and the yaxis, these conditions follow from the imposed symmetry, \p and j being symmetric, (j) and w antisymmetric, at the upper (x = Lx) and the lefthand (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 currentsheet 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 hyperbolicelliptic type. Moreover, the equations are nonlinear and twodimensional, 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 currentsheet 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 pileup) 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 pileup, 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 currentsheet 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 quasistatically because of changes in the boundary conditions and magnetic diffusion. By solving these equations numerically it is found that reconnection always occurs at a finiteangle Xpoint and that the reconnection rate seems to remain finite for Y\ —> 0. No indication of currentsheet 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 currentsheet 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 constantt/; 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 semiquantitative. 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 constantt/; 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 constanty; property breaks down for k ~ 3/a2, where A' ~ 3~{. A more general analysis without the constanttp 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 fieldline 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 (427)
= 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)[lnq(r)/m]J
= 0
The shortwavelength 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  ^ xvo = 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 Xpoint 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 Xpoint. 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 currentdensity 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 selfstabilization process of an unstable tearing mode when growing to finite amplitude. The process may, however, destabilize other modes owing to the currentprofile 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 highrc tearing modes with n up to 20 (Hallatschek et al, 1998). By currentprofile tailoring around the major resonant surfaces completely tearing modestable 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 timescale 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 Xpoint and the Opoint, where V\p = 0, so that the advection term vanishes: —(xpxxp0) = rjxjx  rjojo
(4.51)
96
4 Resistive instabilities
Since \pxxp0 = 2xpu jx  jo  2xpi A'/wj and r\xr\oobtains approximately w/ ~ 77 A' + jrwi.
2xpidr\/dxp, one (4.52)
The drj/dxpterm is related to the aterm 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 Xpoint 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 tokamaklike 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 0point, 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 Xpoint 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 Xpoint. 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 77profiles across the Opoint (full line) and the Xpoint (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 lowfield 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 bootstrapcurrent 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 phasespace 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 fieldline 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 crossfield 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 phasespace 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 constanttp 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 longwavelength 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 firstorder 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{xxsr
[
}
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 Ffunction, 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 (469) Hence in the MHDstable 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)
xx 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 nonreconnected 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 SweetParker 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.\ = lKx s ^ lvolft/2.
(4.98)
(Note that in our normalization B* = VA 1, Lp = po/pfo = pressure scaleheight, one has Bi • Vjo < Bo • V/i, hence the kinkmode 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_,fcrefer to the direction of the magnetic field. For K pointing along the pressure gradient, Kp^ > 0, convective instability arises, if the linebending 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, fieldline 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 socalled 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 inverseaspectratio e = r/R < 1 and circular crosssection 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 crosssection. In a reversedfield pinch, however, where the toroidal field is much weaker and hence the fieldline 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 Pressuredriven 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
CorrLs.
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 largeaspectratio tokamak of circular crosssection, Mercier's criterion becomes (Shafranov & Yurchenko, 1968) K/r^lgV^^r
(4134)
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 linebending effect. The dispersion relation in the local approximation, (4.118), now reads
(4135)
^Tinhpo 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, ^
(4136)
such that the ballooning equation (4.131) becomes > = 0.
(4.137)
4.5 Pressuredriven 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 scalelengths 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 (CorreaRestrepo, 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 nondimensional 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 (CorreaRestrepo, 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 Pressuredriven 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 longwavelength tearing modes, smallscale pressuredriven 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 RayleighTaylor 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 pressuredriven modes The nonlinear dynamics of pressuredriven modes is quite different from that of currentdriven 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,* pressuredriven instabilities tend to accelerate nonlinearly by a selffocusing process, which avoids reconnection. This behavior has been analyzed by Cowley & Artun (1997) for the simple model of a plane stratified plasma with a RayleighTaylorunstable density profile g • Vp < 0. Gravity is chosen in the negative xdirection g = — gex. The system is stabilized by a horizontal magnetic field B = Bo(x)ez, which is linetied 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 linetying 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 resistivitygoverned 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 fingerlike 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 twodimensional 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 quasilinear effect, which is stabilizing, while the last term acts as a nonlinear drive leading to explosive growth. Ignoring the quasilinear term, the xdependence becomes parametric. Integration over y and further simplification leads us to the following basic equation (4.147) where the overline indicates the yaverage. 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 finitetime singularity to, Ay ~ J(to — t)/tfQ, as can easily be seen. Hence an initially smooth function £(y9t = 0) develops a fingerlike 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 Pressuredriven 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 crosssection 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 fieldline 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 fingerlike structures just discussed. The dynamics found in the RayleighTaylorunstable 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 threedimensional 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 threedimensional quasiequilibrium. Where the helical shift is toward the torus outside, a steep pressure gradient is generated, on which highn ballooning modes become unstable, exhibiting explosive nonlinear growth with a fingerlike 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 selffocusing 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 smallscale 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/? circularcrosssection, largeaspectratio tokamak, where the dispersion relation reads (Glasser et al, 1975)
r()v 5/4 /
A' = 2TC—^—^
1DR^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 smallscale 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 shearflow 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 pinchtearing 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 (dashdotted) (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 grayscale 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 KelvinHelmholtz mode is stable. Here, the pinchtearing mode remains unstable leading to magnetic islands which, due to their advective motion, assume a droplike 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 KelvinHelmholtz 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. Grayscale plots of the vorticity illustrating the nonlinear evolution of the Bickley jet. (a) Initial state; (b) nonlinear state of the nonmagnetized 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 nonuniformity 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 SweetParker 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 ydependence 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,
(4183)
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 onehalf, F/yo > 0.5. Hence the sheet is tearingmode 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 tearingunstable SweetParker sheet. The figure shows tpcontours 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 zdirection instead of the ydirection. 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 SweetParker 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 SweetParker sheet the more typical case is that of one or a few isolated islands, which are convected along the sheet, growing in a selfsimilar 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 xpprofile 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 smallscale 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 steadystate 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 wellknown Lorenz model, a truncation of the equations of twodimensional 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 antidynamo theorems. A major step from the proofofprinciple models to a theoretical tool for practical applications has been the introduction of meanfield electrodynamics (MFE). Here the fluid motions are assumed smallscale 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 largescale 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 aeffect may give rise to field amplification. MFE is the framework for most theoretical studies of stellar or galactic dynamos. While the basic mechanisms of smallscale field amplification, namely flux tube stretching and twisting, do not require reconnection, the generation of largescale observable fields from smallscale 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 bestknown 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 smallscale fluid motion will be modified by the magnetic force even at still low mean magnetic energy, which reduces very efficiently their contribution to the aeffect, such that progressively larger turbulent scales become involved. This process requires a selfconsistent MHD description. Section 5.2 introduces selfconsistency in the framework of meanfield 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. HighRayleigh 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 twodimensional 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 steadystate 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) \xBd\=
(5.7)
150
5 Dynamo theory
Since
J>E'dl= hxEdF
=  fdtBd¥ = 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 antidynamo 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 steadystate helical dynamo solution has been given by Lortz (1968). Such helical fields have a certain relevance for the dynamo effect in the reversedfield 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 socalled 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 forcefree magnetic field j = XR). A special case of a twodimensional ABC dynamo was first analyzed by Roberts (1972). Helicity is the essential requirement for dynamo action, as is obvious from the aeffect 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 proofofprinciple 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 crosssection F carrying the flux (/>. Stretching the tube to twice its original length reduces the crosssection by onehalf 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 crosssection 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 stretchtwistfold process does not involve reconnection, it occurs on the fastest possible timescale  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 twodimensional geometry, by contrast, only stretching can occur. Whereas a flux tube of the original crosssection 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 antiparallel, hence only the magnetic energy not the flux has increased. By repeating the process, the crosssections of the unidirectionalfield filaments shrink, until resistive diffusion leads to field annihilation, prohibiting further growth of the energy, which ultimately decays, consistent with the antidynamo theorem. The rope dynamo is the prototype of a socalled 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 largescale 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 Meanfield 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 selfconsistency. In most systems of practical relevance, both in astrophysical and laboratory plasmas, the fluid motions responsible for dynamo action are caused by some instability buoyancydriven in stars and planets or currentdriven in the reversedfield 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 largescale magnetic fields, these motions generate. This suggests a twoscale statistical approach called meanfield 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 VainshteinZeldovich 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 smallscale 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 quasilinear 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 timeintegral is along the Lagrangian orbit x' = x(tf) and the weak spacetime 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 firstorder 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 aeffect 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 aterm in the I n equation, which now gives the steadystate 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 aQdynamo, which carries a toroidal field much larger than the poloidal field. An aQdynamo is therefore also
5.1 Kinematic dynamo theory
159
characterized as "strongfield dynamo", compared with an a 2 dynamo, where both fields have similar magnitude, which is called "weakfield 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 afitype. 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{ikx} 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 aQdynamo, where the aterm 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 socalled dynamo waves. For sufficiently small poloidal velocity \p, the wave propagates in the +kdirection, 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—F1//2. For \p = 0 the dynamo wave (upper sign in (5.35)) now travels in the —kdirection. 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 meanfield 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 Meanfield 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 Meanfield 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 selfconsistently. In this section we will first discuss some phenomenological models which have been introduced in the framework of meanfield theory. We then give an overview of the present understanding of the solar dynamo, the prototypical application of meanfield 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 finiteamplitude modification of the fluctuating quantities invalidating the quasilinear approximation, which is called microfeedback, and through the Lorentz force acting on the mean flow, called macrofeedback. The most important process in the first category is the nonlinear reduction of a, called aquenching, 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 Alfvenwave effect, which is considered in more detail in section 5.4. Basically it means that, in the presence of a mean field, smallscale 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 /cspace) decorrelating interactions of the smallscale 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 aterm, which ultimately leads to a saturated dynamo state. A rigorous theory of aquenching has not yet been developed. Instead the process is modeled in a simple phenomenological way. Since the Alfvenwave 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 antiphase, 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 crossfield equilibrium B2/8n ~ p
165
5.2 Meanfield 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 27day 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 11year 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 lithiumburning 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, Qcontour surfaces are not cylindrical but more disclike.
(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 nonturbulent 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 aQdynamo: 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 meanfield approach characterized primarily by the dynamo coefficients a, rjr (for recent reviews of meanfield 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 Meanfield 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 scalelength 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 (subrotation). (The alternative possibility a < 0, 3rQ > 0 (superrotation) 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 subrotation only below the poles, while below the equator there is a (weak) superrotation. These properties of the differential rotation have also be obtained theoretically by solving the meanfield 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 threedimensional 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 meanfield 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 + (yD^. (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 scaleheight 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 timescale, hence dt ~ e ^2, and we write po + ep, T = Tb + eT. By contrast, p relaxes on the sound timescale 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 quasiincompressibility 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 timedependent velocity field with a few major but rather irregular vertical columns concentrated near the inner shell boundary and embedded in a bath of smallscale 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 superrotation at the equator and subrotation in the polar region, the Ra/Rac = 50 case shows a more complicated pattern with alternating regions of super and subrotation on the surface, similar to the mildly supercritical but highPrandtl 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 presentday direct simulations are far below those prevailing in these astrophysical objects. In the present simulations no tendency toward a disclike differential rotation pattern, as expected for the Sun (fig. 5.9), is recognizable. 5.3.2 The selfconsistent 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 superrotation, dashed contours subrotation, 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 zaxis 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 threedimensional 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 TaylorProudman 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 threedimensional 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 forcefree 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 TaylorProudman 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 77dependence 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 soundcrossing 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 reestablish 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
anticyclone
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 aQdynamo, the process (a) representing the "cyclonic event" (Parker, 1955a), which produces the aeffect in meanfield theory (cf. fig. 5.7).
5.3.3 The geodynamo The geomagnetic field is certainly the bestknown 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 outercore 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 outercore 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 highlatitude 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 stressfree boundary conditions, the continuity of the tangential components dn\t = 0, instead of the physical noslip 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 noslip conditions assumed in the GlatzmaierRoberts 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 noslip 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 quasihom*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 smallscale 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 quasitwodimensional, 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 kineticmagnetic quantities, the energy* i2
,
R
i2\
lV^/
k
_+i2
, i,,2\ _
_
FK
F
k
and the crosshelicity K = J2 Vk • B_ k = \ ^ k
(z+ 2  ZJT2) ;
(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 crosshelicity, but the magnetic helicity is now replaced by the meansquare magnetic potential
^
= ElVk2
(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 /cband, these quantities propagate or "cascade" in /cspace, 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 largescale 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
NavierStokes
Ek Kk Hk
Ef
3D direct direct inverse direct direct
Ek Kk
HZ E?
2D direct direct inverse inverse direct
of selforganization. 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 (NavierStokes) turbulence, where the inverse energy cascade in 2D results in the wellknown formation of large eddies, a behavior which differs fundamentally from 3D turbulence, where no inverse cascade occurs. The cascade properties of MHD and NavierStokes 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 timescale 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[EaK]
=dh
a
B =0, /VXVBI
(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 largescale structures, which is consistent with the direct cascade of the crosshelicity. The degree of alignment is measured by the velocitymagnetic field correlation
p = K/E.
For weak initial correlation the tendency toward alignment is, however, superseded by a competing process of selforganization 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 forcefree 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 powerlaw 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 selforganization 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 crosshelicity 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 timespan reflects the structure of the dissipative eddies, viz., current and vorticity microsheets, 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 largeeddy 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 forcefree 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.20.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 /cspace 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 falloff 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 & MiillerGroeling, 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 timeaveraged 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)k5/\
(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. Loglog 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 selfsimilarity (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 SheLeveque 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 SheLeveque model £^L, the dashed line the modified SheLeveque 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 SheLeveque formalism (5.102). The numerical results indicate that the fundamental scaling is Kolmogorovlike, i.e., g = 3. Visualizations of the smallscale turbulent structures clearly demonstrate the sheetlike 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 Alfvenwave 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 crossfield eddylike motions, much as in hydrodynamic turbulence. Thus the Alfvenwave 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 Alfvenwave 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 spacefilling quasistatic 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 spacefilling 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
= =
BVv + ^V2B, wVv + /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 smallscale 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 ropelike 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 ~ 2040 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. Threedimensional 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 timescales, we have (rii — ne)/ne < 1, so that we can set ne = rii = n. This approximation is called quasineutrality, 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 quasineutrality 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 scalelength 8n ~ (T/T^) 1 / 2 , the width of a resistive current sheet, becomes smaller than one of these scalelengths. 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 parallelperpendicular 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 nonMHD effects let us write (6.19) in nondimensional 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 Twofluid 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 twofluid model written in nondimensional 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 + Zj1 + 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 twofluid framework, most notably the magnetic helicity H = J A • BdVl We find up to resistive and boundary effects
f ~*e/Bi**KOK*
(629)
Hence also in twofluid 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/? whistlermediated 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 twofluid 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 timescale. 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 quasineutrality 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 nonideal 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 +/cBvi = 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 Highfi whistlermediated reconnection
207
While for long wavelength kdi < 1 one recovers the shear Alfven wave co2 = kh>\, in the shortwavelength 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 crossproduct on the r.h.s. of (6.33) makes the whistler circularly polarized, ^T = +1,
(6.38)
righthand for wave propagation in the direction of Bo, lefthand 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 audiofrequency 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 scalelength 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 Zpinches, 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 Highfi whistlermediated reconnection
209
conditions has the form of an Xpoint configuration. In the region \x\ < d\ around the Xpoint, 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 inertiadominated 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 multiscale structure, in particular the current density develops a cusplike singularity, resulting from continuous acceleration at the Xpoint, 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 cusplike 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 orderofmagnitude 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 pileup 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 currentsheet 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 Highft whistlermediated 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 = (xXj)2+(yyj)\ *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 Xpoint, 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 Xpoint 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 Xpoint 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 blowup 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 Xpoint and would remain
6.2 Highfi whistlermediated reconnection
215
topologically invariant in the absence of dissipation. This pileup 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 quasisingular cusplike behavior at the Xpoint (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 KelvinHelmholtz 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 bellshaped crosslayer current profile typical of resistive MHD (see section 3.2.1), noncollisional current layers, because of free electron acceleration at the Xpoint, develop a cusplike 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 KelvinHelmholtz unstable, if the parallel magnetic field is sufficiently weak. While a resistive current sheet is KelvinHelmholtz stable since the outflow velocity is subAlfvenic, 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 KelvinHelmholtz 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 KelvinHelmholtz 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 breakup 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 selfconsistently 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 Highfi whistlermediated 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 twodimensional, 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 breakup 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 pileup 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 highspeed 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 Ioncontrolled 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 Xpoint adjusts automatically. The first simulations of whistlermediated reconnection in the coalescence of two flux bundles were performed by Mandt et al. (1994), where the ions were treated as macroparticles. Here, we discuss results from twofluid 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 Xpoint 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 quasiAlfvenic. Consistent with this fast process, the ion flow does not form an extended layer. In fact, the ion velocity changes abruptly in a shocklike 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 HighP whistlermediated 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 updown, rightleft 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, threedimensional 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 KelvinHelmholtz instability (EKHI) discussed in section 6.2.4, and the lowerhybrid 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 quasiAlfvenic 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 Lowfi 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 fourfield and threefield models Most laboratory plasmas are embedded in a strong, nearly constant magnetic field Bo. As already discussed for the reducedMHD 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 timescales of both waves, %A = (vAk±)~1, (3.24), and for 1> (637): 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 zcomponent 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 LowP noncollisional reconnection
225
which we postpone to chapter 7, the fluid theory turns out to be at least semiquantitatively 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 fourfield 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//cfcj_tpi = 0, cov\\\ — (Tj + Te)feni = 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 xdirection 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 soundwave branch as indicated above, we restrict consideration to the threefield 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 + feB0i  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 quasistationary reconnection is implied, in particular E + v • VF = 0 from (6.104). However, this equation cannot be valid at the Xpoint, 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 Xpoint is governed by electron acceleration along the zdirection. Reconnection occurs, but not in a quasistationary 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 Joo 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)~— / ^xx'^F(x\y,t)dx\ 2ae JoD
(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=
fBVjd2x.
(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 quasistationary reconnection process leading to exponential island growth, (6.108). Figure 6.10 gives the current density from a numerical simulation of the electron inertiadriven 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 currentdensity distribution resembles that in the resistive kink mode, see fig. 4.14(a). 6.4.2 Kinetic Alfvenwavemediated 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 inertiadriven kink mode showing an extended current sheet of width de. from the current sheet of macroscopic length characteristic of the low/? electron inertiadominated 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 Xpoint. 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 Xpoint 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 Xpoint is again xp = \{x2a2y2l (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
(6128)
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 Xpoint 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 pressuremediated 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 Xpoint with a downstream cone of finite angle. In the limit d2e < p], vorticity and current density form shocklike 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 microcurrent 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 pressuredominated 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 pressuredominated 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 = fy>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 timescale. The finitetime singularity in (6.134) should not be taken literally, since this selfsimilar solution is only valid in the early nonlinear phase of the kink mode evolution, much as the selfsimilar 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 threefield 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 finiteamplitude saturation, f$p >, 2 in this case. Saturation occurs when the helical shift £ of the plasma core no longer points toward the Xpoint but is rotated out of phase, such that reconnection is no longer driven. From a series of simulation runs one obtains a semiquantitative 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 EBV$ = 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 Xpoint 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 finiteangle Xpoint, 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 Xpoint 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 Xpoint 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 higherorder 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 Xpoint is a microcurrent 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 nongyrotropic 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 Xpoint configuration to be established. The Xpoint 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 inflowoutflow 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 Xray (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 offcentered chord, the signal is inverted, with a sudden rise coinciding with the sawtooth collapse and a gradual subsequent decay. The crossover 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 Xray detectors
Vacuum vessel
slot
(b)
Fig. 6.14. Observation of sawtooth oscillations by softXray diagnostics. Schematic drawings of (a) the experimental setup, (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 timescale. 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) socalled giant sawteeth from an ioncyclotronheated plasma (after Campbell et al, 1987).
tooth phenomenon exhibits a more complex behavior in largediameter tokamak plasmas with minor radius a ~ 102 cm. Depending on the operational mode  pure ohmic, neutral beam, or radiofrequency (mainly ioncyclotron 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 halfway 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 timescale 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 gprofile 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 largediameter 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 gvalue 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 twofluid 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 twofluid 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 sawtoothfree, while socalled Lmode 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 freeenergy 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>
(743)
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 highfrequency Langmuir waves and one lowfrequency ionsound wave. The second process, nonlinear waveparticle interaction, results from the waveparticle 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 quasilinear effect. To go beyond the quasilinear 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 waveparticle resonance (see also Dum & Dupree, 1970). The basic ingredient of the theory is a modification of the freeparticle propagator, i.e., the ( k v — co^) resonance in the quasilinear 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 freeparticle orbit by the average (exp{—ik • x(t)}) over the true orbits, which is calculated in the diffusion approximation. In a onedimensional statistically hom*ogeneous problem this is described by a diffusion in velocity space Uikx(t)\ _
2
Qikx+ikvt±k
Dvt3
(7 45)
where Dv is determined selfconsistently by the nonlinear integral equation
Dv{v)
= 4 / ^ k2^2Re f mZj J 2%
Jo
generalizing the quasilinear 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 freeparticle 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. Phasespace 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