Investigation of Kinetic Physics of Magnetic Reconnection under a Finite Guide Field by Xiang Lu A dissertation submitted to the Graduate Faculty of Auburn University in partial ful?llment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama December 12, 2011 Keywords: Plasmas, Magnetic reconnection, Particle simulation Copyright 2011 by Xiang Lu Approved by Yu Lin, Chair, Professor of Physics Stephen Knowlton, Professor of Physics Edward Thomas, Jr., Professor of Physics Abstract Magnetic reconnection, the process by which the magnetic energy is converted into heat and ?ow energy, is believed to be a very important process in space and lab- oratory plasmas. Among various reconnection models, the Petschek model suggests a promising explanation for the fast time scale phenomena caused by the reconnection. For the antiparallel (zero guide ?eld) cases, the Petschek model contains a small di?u- sion region and four standing slow shocks that bound the out?ow region. In contrast to the pure antiparallel case, theoretical and simulation studies show that in general cases with a ?nite ?eld, the Petschek model is modi?ed with complicated wave struc- tures. Such models, however, only correspond to the two-dimensional (2D) magnetic reconnection with a long, extended X-line. Two remaining issues in the Petschek model are still poorly understood. One is the onset mechanism for the reconnection, and the other is the wave properties in the out?ow region of general three-dimensional (3D) reconnection. Tearing mode instabilities are thought to play a fundamental role to trigger the reconnection. In this thesis, we extend the recently-developed gyroki- netic electron and fully kinetic ion (GeFi) simulation model to nonuniform plasmas, and use it to study tearing mode instabilities under a ?nite guide ?eld. By remov- ing the rapid gyro-motions of electrons in the calculation, the main advantage of the current GeFi model is that the realistic mass ratio between electrons and ions is al- lowed to be employed in the simulation code. Through an eigenmode analysis, the improved GeFi model is benchmarked against the linear theory of the tearing instabil- ity. Furthermore, our simulation results show that the ion kinetics play an important role in the linear growth rate of the tearing mode instability as well as in the non- linear saturation of tearing modes. For the cases with multiple tearing modes, the ii interactions between them lead to the coalescence of neighboring magnetic islands, and form larger saturated magnetic islands, which has been suggested as a necessary condition to trigger a fast reconnection. In order to explore the properties of low frequency waves generated by the reconnection under a ?nite guide ?eld, a 3D hy- brid simulation is carried out for the large-scale structure of the reconnection layer. For the case with an in?nitely long X-line, quasi-steady rotational discontinuities are formed behind a leading plasma bugle that propagates away from the reconnection site. Field-aligned structures are also observed in the transition region between the steady layer and the leading bulge. For the cases with a limited extent of X-line with length < 30di, the perturbations caused by the reconnection propagate along the magnetic ?eld line, and the wavefront of propagation has the properties of kinetic Alfv?en wave. There are no evident discontinuities that form to bound the out?ow region. Due to the propagation of the waves, a layer of accelerated plasmas is found beyond the extent of the X-line. The simulation indicates that the wave structures in the reconnection are greatly modi?ed from that in the 2D Petschek reconnection model when the X-line has a ?nite length. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 An Overview of Magnetic Reconnection . . . . . . . . . . . . . . . . . . . 1 1.1 Steady State Reconnection Models . . . . . . . . . . . . . . . . . . . 3 1.1.1 Sweet-Parker Model . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.2 Petschek Model . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Tearing Mode Instability . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.1 Resistive Tearing Mode . . . . . . . . . . . . . . . . . . . . . . 7 1.2.2 Collisionless Tearing Modes . . . . . . . . . . . . . . . . . . . 10 1.3 The Structure of Reconnection Layer . . . . . . . . . . . . . . . . . . 13 1.3.1 Summary of MHD Waves and Discontinuities . . . . . . . . . 13 1.3.2 Previous Research on the Structure of Reconnection Layer . . 16 1.4 Numerical Models of Magnetic Reconnection . . . . . . . . . . . . . . 18 1.5 Motivations and Objectives of the Thesis . . . . . . . . . . . . . . . . 21 2 Improved GeFi Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.1 Formulations of the GeFi Model . . . . . . . . . . . . . . . . . . . . . 24 2.1.1 Kinetic Equations of Particles . . . . . . . . . . . . . . . . . . 24 2.1.2 Maxwell Equations of GeFi Model . . . . . . . . . . . . . . . . 27 2.2 Simulation Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.1 The ?f Method . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.2 Advance Particles . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.3 Operators < .. > and ?.. . . . . . . . . . . . . . . . . . . . . . . 35 iv 2.2.4 Field Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.3 Benchmark for Uniform Plasmas . . . . . . . . . . . . . . . . . . . . . 37 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3 Simulations of Tearing Mode Instabilities under a Finite Guide Field . . 40 3.1 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Eigenmode Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 Benchmark for a Simple Current sheet . . . . . . . . . . . . . . . . . 45 3.4 Linear Simulation of the Tearing Mode in the Harris current Sheet . . 47 3.5 Nonlinear Simulation of Tearing Modes in the Harris current Sheet . 59 3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4 Generation of Low Frequency Electromagnetic Waves by Magnetic Recon- nection under a Finite Guide Field . . . . . . . . . . . . . . . . . . . . . 67 4.1 3D Hybrid Simulation Model . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Properties of Kinetic Alfv?en Waves . . . . . . . . . . . . . . . . . . . 70 4.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.3.1 Case with an In?nite X-line: 2?0 = ? . . . . . . . . . . . . . 71 4.3.2 Case with a Finite X-line: 2?0 = 8di . . . . . . . . . . . . . . . 77 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 v List of Figures 1.1 A schematic representation of reconnection is shown. Magnetic ?eld lines become reconnected in the di?usion region and form a X-point. Two islands are generated and bounded by the separatrix. Picture is cited from http://www.aldebaran.cz/astrofyzika/plazma/reconnection.html. . 2 1.2 Sweet-Parker model of reconnection is shown with a steady di?usion region of length 2L and width 2l. vi is the in?ow speed, and vo the out?ow speed. Picture is cited from http://www.scholarpedia.org/article/reconnection. 4 1.3 Petschek model of reconnection is shown with a small di?usion region. The out?ow region is bounded by two pair slow shocks. Picture is cited from http://www.scholarpedia.org/article/reconnection. . . . . . . . . . 5 1.4 Magnetic ?eld-lines in the vicinity of a magnetic island generated by tear- ing modes is shown. A magnetic island is bounded by separatrix. Picture is cited from http://farside.ph.utexas.edu/teaching/plasma/lectures1. . . 8 1.5 A simple slab current sheet is shown with antiparallel magnetic compo- nents Bx and a strong guide ?eld By. Picture is cited from Drake and Lee [1977a]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.6 Directions ofkandBnear the center of the current sheet is shown. Picture is cited from Drake and Lee [1977a]. . . . . . . . . . . . . . . . . . . . . 12 1.7 Phase velocity and group velocity of MHD waves. Picture is cited from http://solar.physics.montana.edu/magara/Research/Topics/MHDwaves. 14 vi 2.1 Coordinates transformation from particle phase space to guiding-center phase space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2 ?b??2?b for Bx0By0 = 1 (black) and Bx0By0 = 0.5 (red). The horizontal line is the normalized z-axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3 The four points gyro-averaging scheme is shown where b is the local mag- netic direction and r1 is a point on the gyro-circle. . . . . . . . . . . . . 36 2.4 The top plot : a comparison between the dispersion relations of ?Bz for the fast magnetosonic/whistler branch obtained from the GeFi simulation (open dots), and the corresponding analytical linear dispersion relation based on the fully kinetic mode (solid lines). The bottom plot : a compar- ison between the dispersion relations of ?By for the shear Alfv?en/kinetic Alfv?en mode branch obtained from the GeFi simulation and the analytical theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.1 The magnetic ?eld (left, and the red solid line on the right) and particle density (black solid line on the right) of the Harris current sheet with a tanh and sech pro?le along Z direction, respectively. The magnetic ?eld consists of two components, Bx0 and By0. Particle density is indicated by the black line, where the dotted line means the background particle density. 41 3.2 Contour of Ay and eigenfunction of ?Ay with Lz = 10.0?i, Lx = 10.0?i, L = 0.5?i, By0/Bx0 = 10, Ti = Te, and C1 = ?0.14vthe. . . . . . . . . . . 45 3.3 Relationship between the linear growth rate and Bx0/By0 with Lz = 10.0?i, Lx = 10.0?i, L = 0.5?i, Ti = Te, and C1 = ?0.14vthe . . . . . . . . 46 3.4 The linear growth rate as a function of kxL (top) and L/?s (bottom) with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 (top), and kxL = 0.41 (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 vii 3.5 The contours of perturbations magnetic ?eld and the corresponding eigen- functions with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 kxL = 0.5. . . 49 3.6 The e?ect of inhomogeneity with cases Bx0/By0 = 0.05 (left column) and Bx0/By0 = 1 (right column), for Ti/Te = 0.01 and L = 0.5?s. The top row : eigenfunctions ?Ay (solid lines from GeFi simulation, dashed lines from eigenmode calculation; the middle row : the contour of ?Ay from simulations; the bottom row : the contour of ?Ay with ?H?B turned o?. . 50 3.7 The growth rate as a function of Ti/Te for three cases L = ?s,1.5?s,3.0?s with Bx0/By0 = 0.05, kxL = 0.5. . . . . . . . . . . . . . . . . . . . . . . 51 3.8 The growth rate as a function of kxL for three di?erent guiding ?elds Bx0/By0 = 5.0, 3.0, 1.0 with Ti/Te = 1.0, L = 0.5?s. . . . . . . . . . . . . 53 3.9 The growth rate as a function of L with Ti/Te = 1.0, Bx0/By0 = 5.0, kxL = 0.41. ? : cases with ions; ? : cases without ions; solid lines: the growth rate ratio between cases with ions and cases without ions. . . . . 55 3.10 The growth rate as a function of Ti/Te with L = 0.5?s, Bx0/By0 = 3.0, kxL = 0.41. ? : with electrostatic e?ect; ? : without electrostatic e?ect; the solid line : the growth rate ratio between them. . . . . . . . . . . . . 56 3.11 The contours of typical perturbation structure with L = 0.5?s, Bx0/By0 = 3.0, kxL = 0.41, Ti/Te = 2.0. . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.12 The comparison among tearing mode growth rate as a function of Bx0/By0 obtained from GeFi simulation (open circles), the eigenmode calculation (solid line) and Drake and Lee analytical asymptotic matching Eq.3.23 (dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 viii 3.13 The time evolution of the magnitude of eigenfunction A? (left) and the contour plot of Ay at time t = 20??1i (right) with L = 0.5?s, Bx0/By0 = 10.0, Ti/Te = 1.0, and kxL = 0.41. . . . . . . . . . . . . . . . . . . . . . 61 3.14 The time evolutions of the magnitude of eigenfunction ne (top left), ni (top right), A? (bottom left), and the contour plot of Ay at time t = 20??1i (bottom right) including ions. . . . . . . . . . . . . . . . . . . . . . . . . 62 3.15 Contours of vector potential Ay at di?erent times for the case with four unstable modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.16 The time evolution of the magnitude A? for the four modes. . . . . . . . 65 4.1 The contour plots of Bx, By, Bz, Vx, Vy, Vz, T, N, and con?guration of magnetic ?eld lines. The solid line corresponds to the position at z = 40di, along which there appears a quasi-steady structure. The solid line at z = 40di indicates the boundary of the leading bugle. The solid line z = 40di shows the transit boundary between the plasma bugle and the steady state layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2 The spatial pro?les of Bx, By, Bz, N, Vx, Vy, Vz, T?, T?, J?, and B at z = 40di and hodograms of the tangential magnetic ?eld (Bz-By). . . . . 73 4.3 Top ?gures: the contours of variations of Bx and N along the y direction in the x-z plane at y = 0; bottom ?gures: the contours of variations of Bx and N along the y direction in the y-z plane at x = 4di. The solid line along z = 52di in the top left ?gure indicates the boundary of the leading bugle. The solid lines in the top ?gures indicate the y-z plane at x = 4di; the solid lines in the bottom ?gures show the local initial magnetic ?eld direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 ix 4.4 The contours of variations along the y direction of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ey, J?, E? in the y-z plane at x = 3di when t = 150??1i . . . . 76 4.5 The con?guration of magnetic ?eld lines near the current sheet at t = 50??1i (the top left) and t = 50??1i (the top right), and the contour plots of Bx (the bottom left) and Vx (the bottom right) in the xz plane of y = 25di at t = 50?i. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.6 The relationship between time t and the distance L from the center of the di?usion region (0,0,0) to the position of wave front. . . . . . . . . . . . 79 4.7 The contour plots of Bx and density N in the x-z plane of y = 16di at t = 80?i. The red dot stands for the position (4, 16, -20). . . . . . . . . 80 4.8 The time evolution of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ez, E?, B, Ex?b, B?, and the hodogram of Bx ?Bx?b at the position (4, 16, -20) as indicated by the red dot in Fig.4.7. . . . . . . . . . . . . . . . . . . . . . 81 4.9 The contour plots of Vz at x = 0 (left) and x = 2di (right) at t = 80??1i . 82 4.10 The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane of y = 0 when t = 150??1i . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.11 The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane of y = 20di when t = 150??1i . . . . . . . . . . . . . . . . . . . . . . . . 85 x List of Tables 1.1 Structure of reconnection layer under various conditions [Lin and Lee, 1994a]. SS: slow shock, SE: slow expansion wave, RD: rotational discon- tinuity, CD: contact discontinuity, IS: intermediate shock, TIDS: time- dependent intermediate shock. . . . . . . . . . . . . . . . . . . . . . . . . 17 xi Chapter 1 An Overview of Magnetic Reconnection Magnetic reconnection is a process that changes the topology of magnetic ?eld lines. It is believed to take place between two plasma regions with antiparallel mag- netic ?eld components. During magnetic reconnection, sheared magnetic ?eld lines will be e?ectively annihilated through magnetic di?usion, and then magnetic energy is transferred into plasma heat and kinetic energy. Magnetic reconnection provides the free energy for phenomena such as solar ?ares, magnetosphere substorms, the tokamak sawtooth collapse, and the ion heating in reversed-?eld pinch [Yamada et al., 2010]. Fig.1.1 is a schematic representation of magnetic reconnection. antipar- allel magnetic ?eld lines break at the di?usion region and form an X-point at the center. Two magnetic islands are generated and bounded by the separatrix. The in?ow plasmas are accelerated and ?ow out from the di?usion region. In a plasma characterized by the ideal MHD Ohm?s law E+V?B = 0, where V is the velocity, E and B are the vector electric and magnetic ?elds, magnetic ?eld lines follow the plasma motion as if they were ?frozen in? the ?ow. This property conserves the ?eld line topology, and thus reconnection can never happen. Varieties of mechanisms can lead to a violation the ideal Ohm?s law. Resistivity is the most commonly cited dissipation that causes reconnection. Models based on resistivity include the Sweet-Parker model [Parker, 1957; Sweet, 1958], in which the extent of the dissipation layer in the out?ow direction is much bigger than the width of the layer, and the Petschek model [Petschek, 1964], in which the dissipation layer is highly localized in both the in?ow and out?ow direction. Theoretical calculations show the magnetic reconnection rate, a direct measure of the speed of the mixing process, of 1 Figure 1.1: A schematic representation of reconnection is shown. Magnetic ?eld lines become reconnected in the di?usion region and form a X-point. Two islands are generated and bounded by the separatrix. Picture is cited from http://www.aldebaran.cz/astrofyzika/plazma/reconnection.html. the Sweet-Parker model is too slow to explain the energy release processes that are thought to be caused by the magnetic reconnection. The Petschek model contains a set of four MHD discontinuities which in the compressible ?ow may be identi?ed as slow mode shocks across which the plasma is accelerated into the outgoing ?ow regions. Instead of dissipating magnetic energy, the di?usion region in the Petschek model mainly plays a role as the perturbation source of slow modes. The Petschek- type reconnection is the ?rst fast reconnection model, in which the reconnection rate is fast enough to explain observed dynamic times of phenomena caused by the reconnection. To validate the Petschek model, two questions should be answered. One is what is the physics in the dissipation region, i.e., the triggering mechanism of reconnection; the other is what is the structure of out?ow region referred as the reconnection layer. The tearing mode instability is one candidate to trigger the magnetic reconnection in 2 both the resistive MHD regime [Furth et al., 1963, 1973] and the collisionless kinetic regime [Drake & Lee, 1977a, 1977b]. On the other hand, theoretical and numerical studies [Lin and Lee, 1994a, 1994b] showed that the structure of the reconnection layer is not as simple as predicted by Petschek but varies with di?erent parameters. This thesis is composed of two parts, which addresses the kinetic physics of magnetic reconnection on two di?erent spatial and temporal scales, from the di?usion region to the large-scale reconnection layer. First, we are going to improve the exiting gyrokinetic electron and fully kinetic ion simulation (GeFi) model [Lin et al., 2005a; Wang et al., 2010], and use it to study the tearing mode instability. Second, a three- dimensional (3D) hybrid model [Lin et al., 2005b; Lin et al., 2010] will be used to investigate the electromagnetic waves generated by the reconnection. In the following, we ?rst give a brief review of the steady state reconnection models, the tearing mode instability, and the structure of the reconnection layer, is given. In the end, the results of this thesis are summarized. 1.1 Steady State Reconnection Models 1.1.1 Sweet-Parker Model The Sweet-Parker model is sketched in Fig.1.2. Sweet [1957] and Parker [1958] determined the speed vi with which ?eld lines are carried into a steady di?usion region of length 2L and width 2l. First of all, for a steady state, the convective in?ow of the magnetic ?eld is balanced by the outward di?usion, so that vi = ?l, (1.1) where ? is the magnetic di?usivity. 3 Figure 1.2: Sweet-Parker model of reconnection is shown with a steady di?usion region of length 2L and width 2l. vi is the in?ow speed, and vo the out?ow speed. Picture is cited from http://www.scholarpedia.org/article/reconnection. Second, the rate 4?Lvi at which mass is entering the sheet must equal the rate 4?lvo at which it is leaving, so that, if the density is uniform, Lvi = lvo, (1.2) where vo is the out?ow speed. The width l may now be eliminated from Eqs.(1.1) and (1.2) to give the in?ow speed as v2i = ?vo/L. If the plasma is accelerated along the sheet by a Lorentz ?j ? ?B force, the out?ow speed vo is the Alfv?en speed at the in?ow, namely, vo = vAi = Bi/?4?? and the reconnection rate is vi = vAiR mi 1/2 (1.3) where the magnetic Reynolds number Rmi = LvAi/?, or in dimensionless form Mi = R?1/2mi , with Mi = vi/vAi being the in?ow Alfv?en Mach number. In the Sweet-Parker mechanism, we identify the sheet length L with the global length and Rmi therefore with the global magnetic Reynolds number. Since in practice 4 Rmi ? 1, the reconnection rate is very small: for instance, in the solar corona where Rmi lies between 106 and 1012, the ?elds reconnect at between 10?3 and 10?6 of the Alfv?en speed. Note the phenomena caused by the reconnection such as solar ?ares usually have the Alfv?en time scale ?A = 1/vA where vA is the local Alfv?en speed. Important extra e?ects are the compressibility, which slows down the out?ow when ?o > ?i, and a pressure gradient along the di?usion region, which can also slow the out?ow when the out?ow pressure po is large enough. 1.1.2 Petschek Model Figure 1.3: Petschek model of reconnection is shown with a small di?usion re- gion. The out?ow region is bounded by two pair slow shocks. Picture is cited from http://www.scholarpedia.org/article/reconnection. The con?guration of the Petschek model is sketched in Fig.1.3. Two waves (slow mode shocks) stand in the ?ow on either side of the out?ow region, where the direction of B reverses, marking the boundaries of the plasma out?ow regions. A small di?usion region still exists, but has much smaller length compared with that in the Sweet-Parker model. The upstream plasma ?ows into the shock where the 5 magnetic energy is converted into the energy of out?owing plasma jets; the outward wave propagation velocity in the ?z direction balances the plasma in?ow velocity in a steady state. The reconnection rate should be equal to the shock propagation velocity which is independent of the plasma resistivity and can be much larger than theSweet-Parkermodel. Atthesametime, themagneticdi?usionisstillthedominant mechanism of energy conversion in the di?usion region near the center z = 0 of the reconnection site. The match of reconnection rate between the di?usion region and the shock region leads to a very narrow di?usion length. In the switch-o? limit of the slow shock, the in?ow speed can be estimated as vzi = vphi = |Bzi|?4?? i , (1.4) where vphi is the phase speed of slow shocks, Bzi is the z component of magnetic ?elds in the in?ow region and is approximated as the normal ?eld of the shocks. The subscripts ?i? and ?o? stand for the in?ow (upstream) and the out?ow (downstream). Note that the switch-o? limit corresponds to the strongest slow shock, across which the tangential magnetic ?elds change to 0. To simplify the calculation, we assume that the plasma is incompressible, i.e., ? = constant. The MHD mass momentum of conservation laws give, vzix = vxo?(x), (1.5) d dx(?v 2 xo?(x)) = BziBi 4? , (1.6) where Bi is the magnitude of magnetic ?eld in the in?ow region. Insert Eq (1.5) into Eq (1.6), and do a little algebra, d dx( x2 ?(x)) = 1 M2Ai Bzi Bi = Bzi Bi , (1.7) 6 where MAi the upstream Mach number and MAi = 1 for the switch-o? limit of slow shocks. After integrating Eq.1.7, we have x ?(x) = Bzi Bi . (1.8) Since the downstream velocity can be expressed as vxo = vAi, (1.9) we have vzi ? BziBi vAi. It has indicated BziBi ? 0.1 by simulations. Therefore the reconnection rate is vziv Ai ? 0.1, much faster than in the Sweet-Parker model. There has been some controversy whether the Petschek model describes the cor- rect physics. It gives a description of reconnection that is only very weakly dependent on the properties of the reconnection region and thus on the resistivity. Numerical simulations [Biskamp, 1986], however, showed that with a spatially constant resistiv- ity the Petschek-type reconnection does not occur. In order to have the Petschek-type reconnection, the resistivity should be enhanced in the di?usion region [Lottermoser et al., 1997]. For collisionless plasmas as in most of the space environments, the de- termination of the anomalous ? in the di?usion region is a hot topic in reconnection research. Current driven instabilities such as the modi?ed two-stream instability [Wu et al., 1983] the low-hybrid drift instability [Daughton, 2003; Ricci et al., 2005] may play a key role in the formation of anomalous resistivity. 1.2 Tearing Mode Instability 1.2.1 Resistive Tearing Mode In the resistive (i.e., collisional) plasma, a current sheet tends to di?use at a slow rate with a time-scale of ?d = l2/?, where 2l is the width of the current sheet 7 and ? = ???1 is the magnetic di?usivity. During the process of magnetic di?usion, magnetic energy is converted into heat at the same slow rate by the Ohm?s law. As pointed out in the ?rst section, however, the magnitude of ?d is often far too large to explain the time scale of dynamical physical processes. Furth [1963] ?rst found that the di?usion can drive an instability with a time scale much faster than the di?usion time, referred as the tearing mode instability. Considering the perturbation as the one in Fig.1.4, the tension force of reconnected ?eld lines tends to pull new loops of ?eld up and down away from the X-point, while the magnetic pressure gradient tends to push plasma in from the sides towards the X-point. The ?eld lines at the sides are curved and so possess a restoring magnetic tension force. Figure 1.4: Magnetic ?eld-lines in the vicinity of a magnetic island generated by tearing modes is shown. A magnetic island is bounded by separatrix. Picture is cited from http://farside.ph.utexas.edu/teaching/plasma/lectures1. 8 Detailed theoretical studies showed [Furth, 1963; Wesson, 1966] that the tearing mode occurs with a wavelength greater than the width of the sheet (kl < 1). The smallest allowable wavelength grows in a time ?3/5d ?2/5A , whereas the longest wave- length has the fastest growth rate ?, ? = 1?? d?A . (1.10) The presence of a guiding ?eld does not a?ect the existence of tearing modes. At any particular location, the maximum growth of tearing modes occurs where the wave vector k is perpendicular to the equilibrium ?eld B0, that is, k?B0 = 0. (1.11) With the addition of a uniform guiding ?eld, this condition can be satis?ed at the center of the current sheet. One may have noticed that the typical time scale of tearing modes ?t = ??d?A is in general much larger than the di?usion time ?d, but much smaller than the Alfv?en time ?A, or the typical fast reconnection time scale, i.e., ?A << ?t << ?d. So the tearing mode is still not fast enough to explain the phenomena that are thought to be caused by the reconnection. The geospace environment magnetic (GEM) reconnection challenge was to aim to discover the physics of fast reconnection [Birn et al., 2001]. They found that nonlinear evolution of a large magnetic island in a current sheet can lead to a fast reconnection with the rate of the order ??1A . The Hall terms in the generalized or kinetic Ohm?s law were believed to be important. For high ? plasmas it was clari?ed that the whistler wave plays a key role in the fast reconnection physics [Pritchett, 2001; Birn et al., 2001]. At lower ? plasmas in the presence of a guide ?eld, it was 9 shown that the kinetic Alfv?en wave drives the fast reconnection [Ricci, 2008; Rogers, 2001]. However, the GEM reconnection challenge avoided the issue on the triggering of the fast reconnection, or more precisely on how a magnetic island forms initially with the size comparable to the system?s spatial scale. One may argue that the role of tearing modes may present a mechanism that form a large magnetic island. Un- fortunately, the tearing mode often saturates at a small level, and the width of the magnetic island is not large enough to trigger the fast reconnection [Rutherford, 1973; Daughton and Karimabadi, 2005]. However, it is also found that a chain of magnetic islands created by multiple tearing modes can lead to an instability, namely the coa- lescence [Finn and Kaw, 1977; Pritchett, 1992]. During the coalescence, neighboring islands attract each other due to the parallel currents inside islands and then form larger magnetic islands [Karimabadi et al., 2005a; Pritchett, 2005]. 1.2.2 Collisionless Tearing Modes Plasmas in the cosmic and tokamak environments are generally collisionless. So the resistive MHD approach to the tearing mode may be invalid. Drake and Lee [1977a] applied a kinetic theory to the tearing mode instability for a current sheet with a strong guide magnetic ?eld (perpendicular to the antiparallel components of the ?eld). They found a perturbed current produced by the induced electric ?eld around the central layer can drive the tearing mode instability, and thereby cause the magnetic reconnection. To estimate the growth rate of the collisionless tearing mode, the simple slab current con?guration as in Fig.1.5 was studied. First, consider a perturbation with the wave number kx, and the growth rate ?. The perturbed magnetic ?eld can be 10 Figure 1.5: A simple slab current sheet is shown with antiparallel magnetic compo- nents Bx and a strong guide ?eld By. Picture is cited from Drake and Lee [1977a]. expressed with the vector potential Ay as Bz = ikxAy, (1.12) Bx = ?Ay?z . (1.13) Then the induced electric ?eld can be written as Ey = ??Ay. (1.14) It has a component parallel to the local magnetic ?eld E|| = Eycos(?), where sin(?) = B0x/B0 = k||/kx(see Fig.1.6). The parallel electric ?eld will generate a large parallel current at z = 0 where k?B0 = 0. (1.15) 11 Figure 1.6: Directions of k and B near the center of the current sheet is shown. Picture is cited from Drake and Lee [1977a]. Note that the induced current is restricted inside a narrow region of |z| < ?. After integrating the z component of the Ampere?s equation, we have [Bx]??? = 4?c ?jy, (1.16) which can be rewritten as [?Ay?z ]??? = ?4?c ?jy. (1.17) So ?Ay?z is not continuous across the induced current layer. Outside of the current layer, E|| = 0, and the ideal MHD equations was used. A simple expression is given as ?? = 1A y [?Ay?x ]???. (1.18) Then combining the above two equations, one has ??Ay = ?4?c ?jy. (1.19) 12 The induced current is ?y = ?n0e 2Ay me . (1.20) According to the cherenkov resonances condition, the growth rate should be ? = k|| ?vthe = kxvthe?l s (1.21) where ls is the magnetic ?eld shear length, and vthe is the electron thermal velocity. A little algebra gives the growth rate as ? = ? ?k xvthe lsd2e , (1.22) where d?1e = c/?pe is the electron skin depth. 1.3 The Structure of Reconnection Layer 1.3.1 Summary of MHD Waves and Discontinuities Startting from the ideal MHD equations, one can easily deduce the dispersion relationship of the linear MHD waves, which are summarized below [Landau and Lifshitz., 1960]. (a) Fast mode: CF = 12 ? [(C2s +C2A) + ? (C2s +C2A)?4C2sC2I]; (b) Intermediate mode (Alfv?en wave): CI = CAcos(?); (c) Slow mode: CS = 12 ? [(C2s +C2A)? ? (C2s +C2A)?4C2sC2I]; where C denotes the wave phase speed, and the subscript F denotes the fast mode, I the intermediate mode, and S the slow mode. Cs = ??P ? is the sound speed and ? is the speci?c heat ratio. CA = B?4?? is the Alfv?en speed and ? is the wave propagation 13 angle with respect to the background magnetic ?eld. Fig.1.7 shows the dispersion- relation of MHD waves. Phase velocity is the propagation speed of a single-mode wave while group velocity is the envelope speed of the waves, i.e., the propagation speed of a group of single-mode waves. From Fig.1.7, one can tell that the fast mode is isotropic, propagating around all angles with respect to the background magnetic ?eld, while the propagations of the intermediate mode and the slow mode concentrate in the direction parallel to the background magnetic ?eld. Figure 1.7: Phase velocity and group velocity of MHD waves. Picture is cited from http://solar.physics.montana.edu/magara/Research/Topics/MHDwaves. An MHD discontinuity can be considered as a thin transition region between two uniform stationary plasma regions across which the magnetic ?eld, plasma den- sity, pressure and ?ow velocity may have a signi?cant jump. The jump conditions are obtained from the conservation laws of the MHD formulations referred to as the Rankine-Hugoniot (RH) jump conditions. It is well known that four types of MHD discontinuities exist: contact discontinuity, tangential discontinuity, rotational dis- continuity and MHD shocks. The jump conditions of di?erent discontinuities are summarized below [Hans and Poedts, 2004]. 14 (1) Contact discontinuity: vn = 0, Bn ?= 0, [Bt] = 0, [vt] = 0, [?] = 0, [P] = 0. A contact discontinuity separates two regions with di?erent plasma densities but the pressure, magnetic ?eld, and plasma velocity are continuous. The component of magnetic ?eld normal to the discontinuity front, Bn, is not equal to zero. The contact discontinuity is a non-propagating structure. (2) Tangential discontinuity: vn = 0, Bn = 0, [Bt] ?= 0, [vt] ?= 0, [?] ?= 0, [P + B28? ] = 0. A tangential discontinuity links two regions with di?erent tangential magnetic ?elds. The component of magnetic ?elds normal to the discontinuity front Bn is equal to zero. No ?ow normal to the discontinuity is allowed. Across the tangential discontinuity the total pressure is balanced , while the plasma densities and tangential plasma ?ow velocities on the two sides of the discontinuity are di?erent. (3) Rotational discontinuity: vn = Bn?4??, [vt] = [Bt]?4??, [?] = 0, [P] = 0. A rotational discontinuity links two regions with di?erent tangential magnetic ?eld and a ?nite normal magnetic ?eld component, and the normal ?ow is allowed. The rotational discontinuity is associated with the propagating nonlinear Alfv?en mode structure through which the normal component of the plasma ?ow velocity relative to the discontinuity is constant and equal to the normal Alfv?en speed, which is vAn = Bn? 4??. Across a rotational discontinuity, the plasma density, the pressure, the normal magnetic ?eld and the normal plasma ?ow velocity remain unchanged. MHD shocks can be further divided into three types: fast shocks, intermediate shocks, and slow shocks. They are formed by the steepening of MHD fast waves, intermediate waves, and slow waves, respectively. The jump conditions of MHD shocks can be described as the following [Priest and Forbes, 2006] : 15 (1) For a fast shock vn1 > CF1 and CI2 < Vn2 < CF2, the tangential magnetic ?eld does not change direction, and [?] > 0,[P] > 0,[Bt] > 0,[vn] < 0. Across the fast shock, the plasma density, the pressure, the magnitude of the tangential magnetic ?eld increase. The normal component of plasma ?ow velocity decreases. (2) For an intermediate shock vn1 > CI1 and vn2 < CI2, the tangential magnetic ?eld changes the direction by 1800, and [?] > 0,[P] > 0,[Bt] ?= 0,[vn] < 0. Across an intermediate shock, the plasm density and the pressure increase. The normal component of the plasma ?ow velocity decreases. The magnitude of the tangential magnetic ?elds is changed. (3) For a slow shock with Cs < vn1 < CI1 and vn2 < Cs2, the tangential magnetic ?eld does not change the direction, and [?] > 0,[P] > 0,[Bt] < 0,[vn] < 0. Across a slow shock, the plasma density and pressure increase. The normal component of the plasma ?ow velocity and the magnitude of the tangential magnetic ?eld decrease. MHD theory, however, does not give the internal structure of discontinuities. With the particle kinetic e?ects in a collisionless plasma, the discontinuities discussed above have a ?nite width on the order of several ion radii [Swift and Lee, 1983; Richter and Sholer, 1989, 1991]. Moreover, MHD theory fails to explain the dissipation mech- anism across MHD shocks in collisionless plasmas. Numerous numerical experiments have con?rmed that the dissipation of MHD shocks in collisionless plasmas arise from the self consistent wave-particle interactions inside of shocks [Sato, 1979; Ugai, 1984; Lee and Lee, 1991]. 1.3.2 Previous Research on the Structure of Reconnection Layer Based on the ideal MHD formulation, the structure of the reconnection layer was studied by solving the one-dimensional (1D) Riemann problem [Heyn et al., 1988; Biernat et al., 1989; Lin and Lee, 1994a]. In this initial value problem, a nonzero normal component of the magnetic ?eld is superposed on a 1D current sheet, 16 By = 0 By ?= 0 Symmetric SS + SS? RD + SS + SS? + RD? Asymmetric RD + SE + CD + SS? RD + SE + CD + SS? + RD? Table 1.1: Structure of reconnection layer under various conditions [Lin and Lee, 1994a]. SS: slow shock, SE: slow expansion wave, RD: rotational discontinuity, CD: contact discontinuity, IS: intermediate shock, TIDS: time-dependent intermediate shock. which breaks the equilibrium of the current sheet. The initial sheet then evolves to a layer with several MHD discontinuities, shocks. The temporal development of this layer is viewed as an approximation of the reconnection layer at increasing distances from the X-line. These studies indicate that ?ve discontinuities and expansion waves may develop in a reconnection layer, including rotational discontinuities, intermediate shocks, slow shocks, slow expansion waves, and a contact discontinuity. The rotational discontinuities and intermediate shocks change the direction of the magnetic ?eld, the slow shocks and slow expansion waves change the magnitude of magnetic ?eld and plasma density, and the central contact discontinuity matches di?erent plasma densities on the two sides. In addition, two fast expansion waves are also present in the solution of the Riemann problem, but they quickly propagate away from the reconnection layer. The results of the Riemann problem are summarized in Table.1.1. The MHD simulations have been carried out for the resistive structure of recon- nection layers. Using a 1D resistive MHD code, Lin et al. [1992], Lin and Lee [1994a], Ugai and Shimizu [1995] have simulated the Reimann problem for the evolution of an initial current system. It is found that the rotational discontinuities are replaced by time-dependent intermediate shocks. To deal with the multi-dimensional recon- nection con?guration, two-dimensional (2D) MHD simulations have been carried out to study the steady-state reconnection. The MHD discontinuities predicted from the 1D simulations are found to develop in the 2D reconnection layer [Shi et al., 1990; Sholer, 1989; Yan, 1992]. Slow shocks in the Petschek-type symmetric reconnection 17 layer and intermediate shocks in the reconnection layer with asymmetric magnetic ?eld and plasma density have been obtained in the cases with By = 0. In the cases with By ?= 0, rotational discontinuities are replaced by time-dependent intermediate shocks that have been found. MHD models and simulations, however, do not account for the particle kinetic e?ects in collisionless plasmas. In an e?ort to include the ion kinetic e?ects, both 1D hybrid simulations and 2D hybrid simulations, in which ions are treated as discrete particles and electrons are treated as a massless ?uid, of the reconnection layers have been performed. 1D hybrid simulations [Lin and Lee, 1995] indicate that the contact discontinuity disappears because of the mixing of ions along the ?eld lines; the slow shock and slow expansion waves are modi?ed. The kinetic structure of the reconnection layer have also been investigated by 2D large-scale hybrid simulations [Lin and Swift, 1996; Lottermoser et al., 1998; Kruassvarban and Omidi, 1995]. These hybrid simulations have shown that the structures of MHD discontinuities and shocks are modi?ed in collisionless plasmas.. Note that previous simulations of the reconnection layer are limited to the 1D or 2D geometry. In these cases, the normal direction to the fronts of low-frequency electromagnetic waves are mainly along the direction of the normal magnetic ?eld, whichisalmostperpendiculartothecurrentsheet. Intherealspaceplasmas, however, magnetic reconnection is frequently observed to have a 3D structure, for which the length of the X-line has a short or ?nite extent. To further understand the physics of the reconnection, we will address the 3D e?ects of the reconnection layer in this thesis. 1.4 Numerical Models of Magnetic Reconnection Numerical simulations have proven to be a powerful tool to understand the physics of magnetic reconnection. Simulation codes can be rougthly divided into 18 four categories: (1) MHD codes, (2) fully-kinetic particle codes, (3) hybrid codes that treat ions as particles and electrons as a massless ?uid, and (4) gyrokinetic codes that treat charged particles as gyrokinetic particles. In the MHD codes, MHD equations are solved numerically with ?nite di?erence or ?nite elements. Since MHD simulations do not contain the individual particle dynamics, they only describe the large-scale dynamics of plasmas. However, the MHD model may still give a correct description of physics on the large-scale bulk dynamics of collisionless plasmas. There are two reasons for that. First, in the presence of a magnetic ?eld particles undergo gyro-motions around the local magnetic ?eld, and thus the motions across ?eld lines are prevented. Second, along the magnetic ?eld direction numerous short-scale wave-particle interactions tend to impede the motion of charged particles. MHD simulations of the magnetic reconnection are mainly focused on the global structure of reconnection. It usually has a much larger spatial scale than di or ?i, and much larger time scale than 1/?pi or 1/?i, where di is the ion inertial length, ?i is the ion gyroradius, ?pi is the ion plasma frequency, and ?i is the ion gyrofrequency. To incorporate the local dynamics due to the separation between the electron and ion masses, the so-called Hall MHD codes and two ?uid MHD codes have been developed. Nevertheless, kinetic e?ects are again missing in these models. Full-particle simulations, in which electrons as well as ions are treated as fully kinetic particles, have been utilized for decades to investigate the kinetic behavior of reconnection [Terasawa, 1981; Leboeuf et al., 1982; Hoshino, 1987; Ding and Lee, 1990; Mandt et al., 1994; Cai and Lee, 1997; Hesse and Winske, 1998; Shay and Drake, 1998; Shay et al., 2001, 2007; Pritchett, 2001; Pritchet and Coroniti, 2004; Swisdak et al., 2005; Daughton et al., 2006; Karimabadi et al., 2007]. In these codes, the particle-in-cell (PIC) technique is often used in which the simulated ?macroparticles? represent many plasma particles. Nevertheless due to large-scale separations both in 19 time (from electron plasma oscillation to reconnection time) and in space (from the Debye length to the system size), various approximations or compromises have to be used. Most of full-particle simulations have employed unrealistically low ion-to- electron mass ratio mi/me, and also limited the domain to a few or ten?s of ?i, and the simulation time much less than the global Alfv?en time scale. Hybrid simulation is another kinetic approach often used in reconnection studies, in which ions are treated as particles, but electrons are treated as a massless resistive ?uid. In general, hybrid codes do not resolve the small spatial and time scales asso- ciated with the electron mass, and have been used to simulate large-scale structures of reconnection associated with the ion dynamics [Lin and Swift, 1996; Nakamura and Scholer, 2000; Petschek, 1964; Lin and Lee, 1994a; Lin and Wang, 2006]. Note that the triggering mechanism of collisionless reconnection, due to the electron kinetic e?ects, is not included in hybrid models. An alternative particle model, called gyrokinetic model, has been broadly used in the fusion research. For the low frequency physics in strong magnetized plasmas such as a tokamak plasma, the gyromotion can be averaged out. The nonlinear gyrokinetic Vlasov-Maxwell equations have been derived by Fireman and Chen [1982] and ?rst applied in gyrokinetic simulations by Lee [1983]. Since the fast gyromotion of charged particles have been removed, larger spatial and time steps are allowed in gyrokinetic simulations than in traditional particle simulations. In the presence of a guide ?eld, particles are strongly magnetized during the magnetic reconnection. Therefore the gyrokinetic treatment can also be applied to the reconnection research. Tearing mode instabilities have been simulated using the gyrokinetic model. Wan et al. [2007] treated ions as gyrokinetic particles and electrons as drift kinetic particles, while both ions and electrons are gyrokinetic particles in Ricci et al.?s work [2004]. But the gyrokinetic and drift kinetic treatment of ions are not suitable to the reconnection physics since the whistler and lower-hybrid regime is beyond the model capability. 20 Recently an innovative fully nonlinear new particle simulation model has been de- veloped, in which the electron dynamics is handled by the gyrokinetic (GK) equations and the ions are treated as fully kinetic (FK) particles. In the gyrokinetic electron and fully kinetic ion (GeFi) simulation, the rapid electron cyclotron motion is re- moved, while ?nite electron Larmor radii, wave-particle interactions, and o?-diagonal components of the electron pressure tensor are retained, as well as the fully kinetic ion physics. This treatment results in a larger time step and allows the realistic mass ratio mi/me in the simulation of magnetic reconnection. The computation e?ciency can be signi?cantly improved as compared with that of the full-particle codes. The model is particularly suitable for plasma dynamics with wave frequencies ? < elec- tron gyrofrequency ?e, and for problems in which the wave modes range from Alfv?en waves to lower-hybrid/whistler waves. The GeFi scheme has been benchmarked in the uniform plasma [Lin et al., 2005] and important new results have been obtained for current sheet instabilities in the presence of a ?nite guide ?eld [Wang et al., 2008; Yoon et al., 2008]. The original GeFi model, however, does not include the e?ects of nonuniformity of the background magnetic ?eld and plasmas. In this thesis, we are going to improve the existing GeFi codes for a current sheet, and study the collisionless tearing mode instability with the realistic mass ratio me/mi = 1/1836. 1.5 Motivations and Objectives of the Thesis As described above, full particle simulations includes the kinetic electron dy- namics but cannot properly handle the realistic mass ratio and the large spatial and temporal scales associated with the reconnection layers, while hybrid simulations cannot model the triggering physics of reconnection in a self-consistent manner. Our approach: ? We utilize the GeFi simulation model, in which the electron dynamics is handled by the gyrokinetic equations and the ions are treated as fully kinetic particles. 21 ? In the GeFi model, the rapid electron cyclotron motion is removed, while ?nite electron Larmor radii, wave-particle interactions, and o?-diagonal components of the electron pressure tensor (necessary for the reconnection physics) are re- tained. ? This treatment results in a larger time step and allows us to treat the realistic mass ratio mi/me. ? The model is particularly suitable for plasma dynamics with wave frequencies < ?e, and wave numbers k? << k?. Analytical theories on collisionless reconnection are very limited to cases with either a zero guide ?eld or a very large guide ?eld, and neglect many kinetic e?ect. On the other hand, because of its low frequency properties, the exiting fully particle simulations on tearing modes often employ mi/me < 100. Our approach: ? Although the GeFi model is promising to handle the electron and ion scales with the realistic mass ratio in one code essential for the reconnection research, it needs to be improved to include the e?ects of the background inhomogeneity because the wavelength of the tearing mode is comparable to the current sheet width. We will modify the GeFi model. ? We derive the eigenmode equations of the tearing modes based on the drift kinetic approximation. By solving them, we benchmark the GeFi model in the nonuniform plasmas. ? Using the GeFi model, various kinetic e?ects on the linear growth of tearing modes are investigated. ? Nonlinear evolutions of both a single tearing mode and multiple tearing modes are studied by the GeFi simulations. 22 While the full development of the GeFi model for the nonlinear fast reconnection is still underway, we also address in this thesis the 3D physics of large-scale kinetic structure of reconnection with the hybrid model. In the Petschek model, the recon- nection layer is bounded by shocks and discontinuities. Previous simulations show the ion kinetic e?ect plays an important role on the structure of reconnection layer. But these simulations are limited to 1D or 2D con?gurations. The large-amplitude wave/discontinuity structures in the general cases with a ?nite X-line length have not been studied. Our approach: ? We extend the previous hybrid simulations to 3D and investigate the 3D physics of the reconnection layer. ? Using the 3D hybrid code, the ?nite length X-line e?ects on the generation of low frequency electromagnetic waves are studied. 23 Chapter 2 Improved GeFi Model Magnetic reconnection involves a wide range of spatial and temporal scales, and is usually triggered by micro instabilities. The reconnection changes the topology of magnetic ?elds on a global scale. Aiming at incorporating multi scale physics of reconnection in one code, the gyrokinetic electron and fully kinetic ion (GeFi) particle simulation scheme has been developed [Lin et al., 2005; Wang et al., 2008]. While the original scheme has been benchmarked for uniform plasmas, further improvement of GeFi is still required for the simulation of magnetic reconnection, in which the interested wavelengths are usually comparable to the scale length of the current sheet nonuniformity. In this chapter, a description of the improved GeFi model is given ?rst. Then I will discuss the algorithm of the GeFi numerical scheme. Finally, the improved model will be validated for uniform plasmas. 2.1 Formulations of the GeFi Model 2.1.1 Kinetic Equations of Particles In the GeFi model, ions are treated as fully kinetic particles, and electrons are modeled by gyrokinetic particles. The Vlasov equation is still valid for ions, while for electrons we transform the phase space variables to guiding-center variables. Di?erent from the original GeFi scheme in Lin et al. [2005], both ions and electrons are advanced by the scalar and vector potential ? and A. 1.Ions 24 In the limit of collisionless plasmas, the time evolution of ion distribution function is totally determined by the Vlasov equation: ?fi ?t + dxi dt ? ?fi ?xi + dpi dt ? ?fi ?pi = 0, (2.1) where pi is the ion canonical momentum. For a singe ion, its motion follows dxi dt = vi = (pi ?qA/c)/mi, (2.2) dpi dt = ?q ??(??xi ?A/c). (2.3) 2. Electrons The particle phase space variables (x, v, t) can be related to the guiding-center phase space variables (R, ?,?,?, t) through the following guiding-center transforma- tion R = x+ ??1v?e? (2.4) ? = mv? 2 2B (2.5) ? = 12mv2 +q? = 12m|p? ecA|2 +q?, (2.6) where ? is the particle gyrofrequency, ? is the magnetic moment, ? represents the Hamiltonian particle energy in the electromagnetic ?eld, ? is the gyro-angle. e? and e?, respectively, are the parallel and perpendicular units vectors relative to the equi- librium magnetic ?eld. The relationship between R and x is illustrated in Figure.2.1, where ?=v?/?=v?e?/? indicates the gyroradius, v?=v?(e1 sin?+e2 cos?). The units vectors e1 and e2 together with e? form an orthogonal system, in which e? = e1?e2. The set of guiding-center variables (R, ?,?,?, t) is equivalent to (R, ?,v?,?, 25 Figure 2.1: Coordinates transformation from particle phase space to guiding-center phase space. t) , so the Vlasov equation can be expressed as ?F ?t + dR dt ? ?F ?R + d? dt ? ?F ?? + dv? dt ? ?F ?v? + d? dt ? ?F ?? = 0. (2.7) In the case that the magnetic ?eld varies slowly in time, the magnetic moment ? is an adiabatic invariant, thus it can be treated as a constant in the guiding-center phase space. To get the electron guiding-center distribution function Fge, one could average Fge around the gyro-motion [Lee, 1983; Frieman and Chen, 1982]. Noting Fge is independent of the gyro-angle ?, the Vlasov equation of guiding-centers can be simpli?ed as ?Fge ?t + dRe dt ? ?Fge ?Re + dve? dt ? ?Fge ?ve? = 0. (2.8) We can rewrite the above equation with the electron generalized momentum pe? as ?Fge ?t + dRe dt ? ?Fge ?Re + dpe? dt ? ?Fge ?pe? = 0. (2.9) 26 To determine a single electron motion, one should take the gyro-averaging of the Lorenz force, and get dRe dt = ve?b ? + c qeB0b0 ?[qe < ?? ? > +??B0] (2.10) dpe? dt = ?b ? ?[qe < ??? > +??B0]. (2.11) where b? = b0 + (ve?/?e)b0 ? (b0 ? ?)b0, ?? = ?? ? ve ? ?A/c, where b0 is the local background magnetic direction, ?? and ?A are the perturbed scalar and vector potentials, respectively [Frieman and Chen, 1982; Lin et al., 2005]. The operator < .. > represents the gyro-averaging, which is carried out numerically on a discretized gyro-orbit in real space. 2.1.2 Maxwell Equations of GeFi Model Once the distribution function is given, physical quantities such as the charged particle density and current density are known by integrating the distribution function in the v space. Then combining with Maxwell?s equations, the plasma behavior is self-consistently determined. Integrating the ion distribution function fi over the velocity moments, the num- ber density ni and current density ji are ni = ? fid3pi, (2.12) ji = qi ? vifid3pi = (qi/mi) ? pifid3pi ?q2i niA/mic. (2.13) To proceed with the electron distribution, let us introduce the gyrokinetic order- ing for electrons ? ?e ? ?e L ? k??e ? ?B B ? O(?), (2.14) 27 k??e ? 1, (2.15) where ? is the mode frequency of interest, L is the macroscopic system length, and k? is the component of the wave vector in the parallel direction with respect to the equilibrium magnetic ?eld B0 [Rutherford and Frieman, 1968; Taylor and Hastie, 1968]. The total magnetic ?eld B consists of the equilibrium magnetic ?eld B0 and the perturbed magnetic ?eld ?B, B = B0+?B. ? is a small parameter. The interested perpendicular wavelength ranges from the scale of electron Larmor radius to the global scale length. In our GeFi model, the parallel direction is de?ned as b0 = B0/B0. The electron gyro-averaging charge density and current in the guiding-center phase space are < ne >= ? Fed3v, (2.16) < je? >= ? pe?Fedp, (2.17) To deduce Poisson?s equation of the GeFi model, we start from the general form ?2??? = ?4?(qini +qene), (2.18) where we have assumed ?2? >> ?2?. Note ne is the electron density rather than the guiding-center density, which can be calculated from the guiding-center distribution by ne = qem e ? d3v(? ?Fge/?w)[??? < ?? > +1c < v? ??A >]+ < Ne >, (2.19) where w = v2/2 [Lin et al., 2005]. By taking |?e??| < 1, Poisson?s equation becomes ?? ?[(1 + ?? 2 pe ??2e )????] + 4?qe?ne ?B ?B? = ?R? +?H?, (2.20) 28 where ?R? and ?H? are expressed as ?R? = ?4?[qi?ni +qe < ?Ne >] (2.21) ?H? = 4?qe[4??ne?B2c j??A??( ?ne?B2)?(?A? ?B)], (2.22) and ?B is the background magnetic ?eld average over the spatial nd temporal scales of wave perturbations. Compared with equation (6) in Lin et al. [2005], the appearance of ?H? is due to the usually higher order terms associated with the inhomogeneities in the background magnetic ?eld and density. To solve Poisson?s equation, ?B? must be known. In the low frequency limit ? << ?e, the electron force keeps balanced in the perpendicular direction. Using this approximation, ?B? can be obtained from the equation below ?neqe???? = ?SB +??[ 14?(1 +??e) ?B?B?] +?H?B +??P?, (2.23) where ?SB = ?? < ?Peg > +?[ji ?B/c] (2.24) ?H?B = 14??[(???e)?(?A? ?B)???4?jc ??A]??HB (2.25) ?HB = 14?{?[(B??)B]???B2/2} (2.26) ??P? = ?[??( P0cB?? e ????)] (2.27) ??e = ??e/2 (2.28) ?? = (1 +? ? e) ?B 4? ?B? + mc2 4?qe??(? ? e??). (2.29) Compared with equation (20) in Lin et al. [2005], the newly added term ?H?B is again due to the usually higher order contribution from the inhomogeneities in the background magnetic ?eld and density. 29 In the GeFi model, A0 is determined by the initial condition, and the perturba- tion ?A is solved at each simulation step. We decompose ?A into three components ?A = A??b + A2e2 + Azez, where e2 = ?b?ez and ?b = ?B/ ?B. Then they obey the following equations (?2 ? ? 2 pe c2 ? ?2pi c2 + ?b??2?b)A? + (?b??2e2)A2 = ?4?c (?ji?+ < ?je? >), (2.30) (?2 +e2 ??2e2)A2 + (e2 ??2?b)A? = ?4?c ?j2, (2.31) ?2Az = ?4?c ?jz, (2.32) where ?j2 = c4?(???B??b)?e2 (2.33) ?jz = c4?(???B??b)?ez. (2.34) Compared with the original model [Lin et al., 2005], terms such as ?b??2?bA? are added in the improved model to include the inhomogeneities of the background magnetic ?eld direction. Di?erent from the original scheme [Lin et al., 2005], the above improved model includes the e?ects of nonuniformity in the background magnetic ?eld direction b and density ne. To appreciate the e?ects of the terms associated with the variation of the background magnetic ?eld direction, in the following we compare the ?rst and the fourth term in the Eq.2.30, ?2A? and ?b??2?bA?. Considering a the Harris current sheet with a sheet normal in the z direction and the antiparallel ?eld component in the x direction, the magnetic ?eld pro?le with a spatial coordinate normalized to the 30 current sheet half-width can be expressed as Bx = Bx0 tanh2(z) (2.35) By = By0 (2.36) Bz = 0. (2.37) In most cases, we can estimate the scale length of A? as the half current width, and therefore the the ?rst term can be estimated as ?2A? ? A?. The pro?le of ?b??2?b is shown in Fig.2.2. The peak values ?b??2?b for the cases with Bx0By0 = 1 and Bx0By0 = 0.5 are located in the center of the current sheet and have values around 1.0 and 0.22 repectively. It is seen that By0 ? Bx0, the background nonuniformity cannot be neglected. Figure 2.2: ?b??2?b for Bx0By0 = 1 (black) and Bx0By0 = 0.5 (red). The horizontal line is the normalized z-axis. 31 2.2 Simulation Algorithms In the GeFi simulation program, full particle method and ?f method have been used. Full particle method is a standard approach for most particle in cell codes, while ?f method is a newly developed tech to reduce the computational variations [Lee, 1987; Parker and Lee, 1993; Chen and Parker, 2003]. 2.2.1 The ?f Method The ?f method is developed to get rid of the background noise caused by the ?nite number of particles. Since only the perturbed distribution function is evolved in the ?f scheme, thermal ?uctuations associated with the background particles are removed. To illustrate the idea of the ?f method, we ?rst linearize the ion and electron Vlasov equations. For ions, it is d?fi dt + qi mi(????+vi ?(???A)) ?fi0 ?t = 0, (2.38) and for electrons, it is d?Fge dt + [ve0??b ? +?ve?b? 0 + cb0 qeB0 ?(qe < ??? ? >)]? ?Fge0 ?R +[?b? ?(??B0) +b?0 ?(qe < ???? >)]? ?Fge0?p e? = 0, (2.39) where?b? = (?ve?/?e0)b0?(b0??)b0?(ve0???e/?2e0)b0?(b0??)b0,?ve? = ?qe?A?/(mec). The distribution function fe and Fge, the potentials A and ? have been divided into the equilibrium part with subscript 0 and the perturbed part with the notation ?, f = f0 +?f (2.40) ? = ?0 +?? (2.41) 32 A = A0 +?A. (2.42) The linearized Vlasov equation above can be simply written in a more general form as d?f dt = ? df0 dt . (2.43) Now de?ne the weight carried by each particle as w = ?ff . (2.44) Since dfdt = 0, we have dw dt = d?f dt f . (2.45) Using f = f0 +?f, one can get 1 f = 1 f0(1? df f ) = 1 f0(1?w). (2.46) Note that in the above approximation, only the ?rst order is kept. It is easy to show that dw dt = ?(1?w) 1 f0 df0 dt . (2.47) Since f0 is usually given in an analytical form, the time evolution of particle weights can be calculated by Eq.2.47. 2.2.2 Advance Particles To advance particles in the code, we need to numerically integrate the ordinary di?erence equations with the following form dZ dt = G[Z,t], (2.48) 33 where Z represents the motion quantities, such as R(t),p?(t) and w(t) of particles in the simulation. There are many methods available for the integration of an ODE, both explicit and implicit. Explicit schemes are straight forward, but should satisfy the Courant condition ?max?t < 1 in order to assure the numerical stability, where ?max refers to the highest frequency in the simulation [Hockney and Eastwood, 1981]. Implicit schemes are generally stable, but require large computation time to solve linear algebra equations. In the improved GeFi model, two implicit methods have been employed, the predictor-corrector method and the fourth order Runge-Kutta method [Burden and Faires, 2005]. In the predictor-corrector method, we need to know not only the value Z at time n but also the value Z at time n-1, then Z at time n+1 can be calculated in two steps: Predictor : Z(n+1) = Z(n?1) + 2G(n)[Z(n)]?t Corrector : Z(n+1) = Z(n) + 12[G(n+1)[Z(n+1)] +G(n)[Z(n)]]?t. (2.49) The fourth-order Runge-Kutta method needs the value Z at time n and involves four steps as Zn+1 = Zn + 16?t(k1 + 2k2 + 2k3 +k4) (2.50) k1 = G(tn,Zn) (2.51) k2 = G(tn + 12?t,Zn + 12?tk1) (2.52) k3 = G(tn + 12?t,Zn + 12?tk2) (2.53) k4 = G(tn +?t,Zn +?tk3). (2.54) Note the predictor-corrector method has the second order accuracy while the Runge- Kutta method above is of the forth order. 34 2.2.3 Operators < .. > and ?.. Two operators < .. > and ?.. have been seen in the gyrokinetic equations (Eq.2.20, Eq.2.23, Eqs.2.30-32). < .. > means the average around the particle gyro-motion, while ?.. stands for the background ?eld averaged over the spatial and temporal scales ofperturbations. Inthissection, Iamgoingtointroducehowtorealizetheseoperators numerically. We use the so-called four points gyro-averaging scheme to realize < .. > numeri- cally [Lee, 1983; Lee, 1987]. The basic idea is to ?nd four points uniformly distributed on the gyro-circle of a guiding-center and then take average of ?eld values at these four points. To ?nd the positions of these four points, one can follow the steps as indicated by Fig.2.3. We let e1 = b?ez where b is the local magnetic ?eld direction and ez is the units vector along the z direction. (If e1 = 0, one can just try ex or ey instead.) The ?rst point r1 is located along the direction e1 at a distance ?e from the guiding-center. Since the magnetic moment ? is conserved in the gyrokinetic model, the gyroradius can be easily calculated as ?e = ? 2?me B /qe. (2.55) In the opposite direction ofe1, another pointr2 can be located. Now makee3 = b?e1. Repeating the above procedure, one can ?nd all four points. Then the gyro-averaging ?eld quantities are calculated as < f >= 14(f(r1) +f(r2) +f(r3) +f(r4). (2.56) To realize ?.., we should store the data of ?elds before the current time step. However, how many previous steps should be kept and how to average them is still a 35 Figure 2.3: The four points gyro-averaging scheme is shown where b is the local magnetic direction and r1 is a point on the gyro-circle. ongoing issue. For the linear simulation or the cases with ?B << B, we can simply treat ?.. as the initial value of it. 2.2.4 Field Solver As is seen, Maxwell?s equations in the GeFi model are complicated partial dif- ference equations (PDE) and mutually dependent on each other. To solve them, iterations are done simultaneously at each time step to obtain ??, ?B?, and ?A. At each iteration, one has to solve the nonlinear Poison?s equation with the following form ?2f(x,y,z)?g(x,y,z)f(x,y,z) = C(x,y,z) (2.57) Two methods have been employed to handle the above task. The ?rst one is ?nite di?erences, which is straight forward, but requires much computational time solving 36 the linear algebra equations, especially for 2D and 3D cases. The second one is related to our particular boundary conditions. In our simulations, periodic boundaries are adopted in the x and y directions, while conductor boundaries are used in the z direction. If g(x,y,z) is only a function of z, i.e., g(x,y,z) = g(z), we are allowed to take the fourier transform on the x and y variables and have ?2f(kx,ky,z)?g(z)f(kx,ky,z) = C(kx,ky,z). (2.58) Such 1D PDE is solved quickly by ?nite di?erences. If g(x,y,z) is not as simple as above, we divide g(x,y,z) into two parts like g(x,y,z) = g0(z)+g1(x,y,z). Note g0(z) is determined by the initial condition that is independent on the y and z coordinates, whereas g1(x,y,z) is treated as a perturbation. Rewrite Eq.2.57 as ?2f(x,y,z)?g0(z)f(x,y,z) = C(x,y,z) +g1(x,y,z)f(x,y,z) (2.59) Eq.(2.59) is solved by iteration with the same method of Eq.(2.58). Since usually g1 << g0, it converges very fast. 2.3 Benchmark for Uniform Plasmas In this section, we are going to benchmark the modi?ed GeFi scheme for a one dimensional uniform system [Lin et al., 2011]. The wave vector k is assumed to be along x. The background magnetic ?eld is in the x-z plane and is allowed to point in various directions relative to k. The top plot of Fig.2.4 presents a comparison between the dispersion relations of ?Bz for the fast magnetosonic/whistler branch obtained from the particle-in-cell GeFi simulation, shown as open dots, and the corresponding analytical linear dispersion relation based on the fully kinetic mode, shown as solid lines. Cases with ?e = ?i = 0.04, mi/me = 1836, and k?/k? = 0.2,0.06, and 0 are plotted. The linear ?uctuations in the simulation are due to random noises, as 37 in usual particle simulations. In the case with k? = 0, the electromagnetic mode approaches the quasi electrostatic lower hybrid mode, and the frequency ?/?i = ?LH/?i = ??i?e/?i = ? mi/me = 42.8, where ?LH is the lower-hybrid frequency and ?i is the ion gyrofrequency of the background plasma. Figure 2.4: The top plot : a comparison between the dispersion relations of ?Bz for the fast magnetosonic/whistler branch obtained from the GeFi simulation (open dots), and the corresponding analytical linear dispersion relation based on the fully kinetic mode (solid lines). The bottom plot : a comparison between the dispersion relations of ?By for the shear Alfv?en/kinetic Alfv?en mode branch obtained from the GeFi simulation and the analytical theory. The bottom plot of Fig.2.4 shows a comparison between the dispersion relations of ?By for the shear Alfv?en/kinetic Alfv?en mode branch obtained from the GeFi simulation and the analytical theory for k?/k? = 0.06. The analytical solution of the MHD shear Alfv?en mode, ?/?i = k?vA/?i = (k?/k) ? 2/?i(k?i) = 0.42, is also shown as the dashed line. 38 It is seen that GeFi simulation results are in excellent agreement with the theo- retical analysis for k? << k?. The benchmark has been performed for cases with ?e and ?i ranging from O(10?2) to O(1). 2.4 Summary In summary, our GeFi particle simulation model, in which the electrons are treated as gyrokinetic particles and ions are treated as fully kinetic particles, has been improved and modi?ed to allow the the existence of modes with wavelengths on the same scale of the background nonuniformity. With fast electron gyro-motion and Langmuir oscillations removed from the dynamics, the GeFi model could be a powerful candidate to solve physics with realistic mass ratio mi/me is a global system. The linearized GeFi scheme is benchmarked for uniform plasmas, and the sim- ulation results using the ?f method show that the model can accurately resolve the physics ranging from Alfv?en waves to lower-hybrid/whister waves, for k? << k? and ? << ?e. 39 Chapter 3 Simulations of Tearing Mode Instabilities under a Finite Guide Field In this chapter, we use the improved GeFi model to simulate the tearing mode instability. The initial conditions are described in the ?rst section. The asymptotic approximation of the Drake-Lee theory [1977a] and the eigenmode theory are pre- sented next. Simulation results for a current sheet in a simple slab geometry are shown in the third section. The current sheet con?guration that we use here has been studied by various people, and can be used to benchmark the GeFi model for the re- connection physics. In the last two sections, linear and nonlinear simulation results of tearing mode instabilities for symmetric the Harris current sheets with various guide ?elds are presented and discussed, repectively. 3.1 Simulation Model the Harris current sheet is an equilibrium solution of the Maxwell-Vlasov system, which has been extensively used in the magnetic reconnection research [Harris, 1962]. The equilibrium magnetic ?eld consists of two components, the uniform guide ?eld By0 and the tanh pro?le of antiparallel component Bx, B = ?xBx0 tanh(z/L) + ?yBy0, (3.1) where L is the half-width of the Harris current sheet. The particle density is described by n = nHsech2(z/L) +nb0, (3.2) 40 Figure 3.1: The magnetic ?eld (left, and the red solid line on the right) and particle density (black solid line on the right) of the Harris current sheet with a tanh and sech pro?le along Z direction, respectively. The magnetic ?eld consists of two components, Bx0 and By0. Particle density is indicated by the black line, where the dotted line means the background particle density. where nb0 is the background density and nH is the peak density in the current sheet. The peak density is obtained from the total pressure balance nH(Ti +Te) = 18?B2x0, (3.3) Here we have assumed that the temperatures for ions and electrons are constant everywhere. The equilibrium velocity distributions for the background ions and elec- trons are Maxwellian. The equilibrium velocity distributions for ions and electrons of the Harris current sheet in the GeFi model are given by [Lin et al., 2005; Wang et al., 2008;] fHi = nh0(2?T i/mi)3/2 e?mi[v2x+(vy?vdi)2+v2z]/2Ti ?e?vdiqiAy(z)/Ti (3.4) fHe = nHsech 2(z/L) (2?Te/me)3/2 exp(? B2xmev2de B22Te )exp(1? mevde Te?e dBx dz ?) 41 ?exp{? 12T e [2?B +me(v? ? vdeBGB )2]}. (3.5) We focus on the instabilities in the x-z plane, where the tearing instability is included, while the current driven instabilities such as the kink instability, the lower hybrid instability are excluded. The normalization follows: the magnetic ?eld is in the units of B0 = ? B2x0 +B2y0, the spatial length is in the units of ?e, the time step is in the units of ??1e , particle density is in the units of n0, where n0 is calculated from ?0 = n0Te/(B20/8?) and ?0 is an input parameter. To benchmark our GeFi model, we also do the calculations for a di?erent current sheet geometry that has been used in the previous tearing mode studies [Drake and Lee, 1977; Katanuma and Kamimura, 1980; Wang et al., 2005;]. In this current sheet, the current pro?le along a strong guide ?eld By0 is given Jy0(z) = ?eneu0(z) = C1enee?z 2 L2 . (3.6) In addition to the guide ?eld By0, the current Jy0 generates an antiparallel ?eld Bx0(z) = 12eneC1?0L??Erf(zL), (3.7) where the Erf(z) is the error a function. In the following, I refer to this model as ?the Drake-Lee current sheet?. Di?erent from the Harris current sheet, the Drake-Lee current sheet is not an equilibrium con?guration, and ions are initialized without the drift motion. This current sheet is particularly suitable for the laboratory plasmas, in which there is usually a strong guide ?eld. 42 3.2 Eigenmode Analysis The energy exchange between tearing modes and charged particles takes place e?ectively in the region where the condition of the Cherenkov resonance is satis?ed: ? = k?vthe. (3.8) Since the tearing mode is a low frequency phenomenon, the condition Eq.(3.8) is mainly satis?ed in the region where k? ? 0. If the gyroradius of charged particle is much smaller than the current sheet width, the drift kinetic approximation is valid. My benchmark of the tearing instability against the analytical theory will be conducted for this simpli?ed situation [Hoshino, 1988; Wan et al., 2005]. The linearized drift kinetic equation is given as ??fj ?t +v? Bx0 B0 ??fj ?x +v? ??fj ?y = [?v? ?Bz B0 1 n(z) dn(z) dz ? dv? dt 2(vy ?Vdj) v2thj ]f0j, (3.9) where j stands for species of particles and uj means the drift speed of ?uid j [Hoshino, 1988; Wan et al., 2005]. Note that in the above equation the electrostatic term has been neglected, since its e?ect has proven to be small for the tearing mode. We can rewrite ?Bz and dv?dt in terms of ?Ay: ?Bz = ??Ay?x , (3.10) dv? dt = qj mj ??Ay ?t , (3.11) Writing ?Ay=?Ay(z)ei(kxx??t) and using ? ?x = ikx, ? ?y = 0, ? ?t = ?i?. (3.12) 43 Eq.(3.9) becomes ?fj = 1 ? ? v?Bx0kxB0 ?( v?kxB 0nj(z) dnj(z) dz + ? mj (vy ?Vdj) v2thj )foj?Ay. (3.13) ?Ay is induced by the perturbed current ?Jz through Ampere?s law, ??Ay = ??0?Jz = ??0 ? j=i,e qj ? dvyvy?fj (3.14) Combining the above two equations, we have a second order ordinary di?erential equation [ ddz?(z) ddz ?k2x?(z)? ? j=i,e Vj(z)]?Ay = 0, (3.15) where ?(z) = 1 + ? j=1,e ?2PjV 2dj c2?2j , (3.16) Vj(z) = ?2? 2 Pj c2 (? 2[1 +?z(??)] +? Vdj vthj(1 + kxBx0Vdj ?B0 )), (3.17) ? = ?B0k xBx0vthj , (3.18) ?? = ?(1? kxBx0Vdj?B 0 ), (3.19) and Z(??) is the plasma dispersion a function. Eq.(3.15) is solved for eigenfunctions by following the usual shooting method, with the boundary condition ?Ay = 0 at the z-boundaries on both ends. A factor of 2 di?erence is noticed between the expression V5j in Hoshino?s equation(20) and the corresponding expression Eq.(3.15) for Vj in our eigenmode analysis. Following the same procedure, Wan et al. [2005] has conducted the eigenmode analysis for the Drake-Lee current sheet. One should keep in mind that the drift kinetic approach is valid only when k??i << 1. 44 3.3 Benchmark for a Simple Current sheet Figure 3.2: Contour of Ay and eigenfunction of ?Ay with Lz = 10.0?i, Lx = 10.0?i, L = 0.5?i, By0/Bx0 = 10, Ti = Te, and C1 = ?0.14vthe. Fig.3.2 shows the contour plot of Ay and the corresponding eigenfunction of ?Ay obtained from the nonlinear simulation with Lz = 10.0?i, Lx = 10.0?i, L = 0.5?i, By0/Bx0 = 10, Ti = Te, and C1 = ?0.14vthe. A magnetic island is generated at the center with an almost symmetric structure, representing a signature of the tearing mode. ?Ay has a double-peak near the center, which is consistent with the tearing mode unstable condition ?? > 0. Recall the de?nition of ?? ?? = (??Ay(?)?z ? ??Ay(??)?z )/?Ay(0) (3.20) where ? is the width of the singular layer. The linear growth rate as a function of Bx0/By0 is shown in Fig.3.3 with the same parameter above except for a varying Bx0/By0. ?Fixed ions? means ions are only advanced by the initial ?elds, meaning the motion of ions is kept unperturbed. 45 Figure 3.3: Relationship between the linear growth rate and Bx0/By0 with Lz = 10.0?i, Lx = 10.0?i, L = 0.5?i, Ti = Te, and C1 = ?0.14vthe 46 Again the linear growth rate agrees the eigenmode theory well. Two other conclusions can be reached immediately: the growth rate is increasing with Bx0/By0, and the ions contribute little to the growing of tearing modes. According to the linear theory, the linear growth rate is given by ? = ? ?k xvthe lsd2e , (3.21) where ls is the magnetic ?eld shear length, vthe is the electron thermal velocity. Since ls = By0?Bx0/?z ? By0/Bx0, we have ? ? Bx0/By0. Since the current is carried totally by electrons for the Drake-Lee current, ions contribute little to the perturbed current based on ?Jy = ?j=i,e(nj0?vjy + vj0y?nj0), where the ratio of the ?rst term between ions and electrons is ni0?viy/ne0?vey = me/mi, and the second term of ions is 0. 3.4 Linear Simulation of the Tearing Mode in the Harris current Sheet the Harris current sheet is the ?rst equilibrium model of one-dimensional current sheet, which has been widely used in the investigation of reconnection. Note that since the Drake-Lee current sheet is not in equilibrium, it can only be applied to cases with By0/Bx0 >> 1. In this section, I show the linear simulation results of tearing modes in the Harris current sheet. In order to benchmark the GeFi simulation scheme against the existing eigen- mode analysis, GeFi simulations with ? = 0, k??i ? 0, ky = 0 and ?B? = 0 are carried out. These parameters are consistent with the assumptions of eigenmode analysis, where ions and electrons are both treated as drift kinetic particles requiring k??i << 1. Since the electrostatic, compressible modes are neglected, we have ? = 0, ?B? = 0. Fig.3.4 shows the linear growth rate as a function of kxL (top) and L/?s (bottom) with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 (top), and kxL = 0.41 (bottom), where ?s is de?ned as ?s = ? mi/me?e which is the ion gyroradius when Ti/Te = 1. Again the simulation results agree with the eigenmode analysis well. In 47 the top ?gure, we can see the growth rate is maximized around kxL = 0.4 and tends to be 0 for kxL > 1, both of which represent the characteristics of tearing modes. The bottom ?gure depicts that the growth rate is decreasing with the current width with a relationship ? ? L?2.4. Figure 3.4: The linear growth rate as a function of kxL (top) and L/?s (bottom) with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 (top), and kxL = 0.41 (bottom). The contours of the magnetic ?eld perturbations Bx, By, Bz from the eigenmode theory and the GeFi simulation are shown in the left Fig.3.5 with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75, kxL = 0.5, while the corresponding eigenfunctions are shown on the right. Little di?erence exists between the simulations and the eigenmode analysis. The Bx, By perturbations are symmetric, while Bz shifts in the opposite directions on each side of the current sheet. As discussed in the chapter2, the nonuniform terms associated with the magnetic ?eld direction b and the electron density ne have been added into the improved GeFi 48 Figure 3.5: The contours of perturbations magnetic ?eld and the corresponding eigen- functions with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 kxL = 0.5. 49 Figure 3.6: The e?ect of inhomogeneity with cases Bx0/By0 = 0.05 (left column) and Bx0/By0 = 1 (right column), for Ti/Te = 0.01 and L = 0.5?s. The top row : eigenfunctions ?Ay (solid lines from GeFi simulation, dashed lines from eigenmode calculation; the middle row : the contour of ?Ay from simulations; the bottom row : the contour of ?Ay with ?H?B turned o?. model. As is seen in Fig.2.2, such an e?ect becomes more profound as By0 gets smaller. The top row of Fig.3.6 shows the absolute values of eigenfunctions of ?Ay obtained from simulations with Bx0/By0 = 0.05 and Bx0/By0 = 1, for Ti/Te = 0.01 and L = 0.5?s. The solid line in each plot shows the simulation results based on the 50 improved GeFi scheme, while the dashed line corresponding to the run in which ?H?B is turned o?. The corresponding contours of the eigenfunctions of ?Ay obtained from the GeFi simulation are shown in the middle row of Fig.3.6, while the runs with ?H?B turned o? are plotted in the bottom row. Clearly, while the di?erence between the results with and without the ?H?B is negligible for Bx0/By0 = 0.05, the corresponding eigenmode structures of ?Ay are, however, signi?cantly di?erent in the case with Bx0/By0 = 1.0. With ?H?B terms included, the tearing mode shows ?ner structure in the current sheet. The growth rates of tearing modes for cases Bx0/By0 = 1.0 with and without the ?H?B term are, respectively, 0.098 and 0.128. Further simulation indicates that the ?H?B term in general cannot be ignored for linear tearing modes when By0 ? or < Bx0. Figure 3.7: The growth rate as a function of Ti/Te for three cases L = ?s,1.5?s,3.0?s with Bx0/By0 = 0.05, kxL = 0.5. 51 The above benchmark work has been limited to the extremely small tempera- ture ratio between ions and electrons Ti/Te = 0.01 in order to satisfy the drift kinetic ion approximation requirement k??i << 1. Since k? ? L?1 (L ? ?S for our appli- cations) and ?i = ?s ? Ti/Te, the parameters k??i are evaluated as ? 0.1 and ? 1, respectively, for cases with Ti/Te = 0.01 and for cases with the more realistic tem- perature ratio Ti/Te = 1.0. It is expected that as the ratio Ti/Te increases, the drift kinetic ion approximation is not valid any more, and therefore the agreement between the eigenmode theory and the GeFi simulation cannot hold. The growth rate as a function of Ti/Te is plotted in Fig.3.7, where three cases with L = ?s,1.5?s,3.0?s for Bx0/By0 = 0.05 and kxL = 0.5 are presented. The di?erence between the eigen- mode theory and the GeFi simulation becomes obvious after the value Ti/Te ex- ceeds 0.5,1.2,3.0, corresponding to k??i ? 0.7,0.73,0.67, repectively, for cases with L = ?s,1.5?s,3.0?s. Our simulation results here are qualitatively consistent with the previous study by Ricci et al. [2005], in which it is pointed out both the ?nite electron Larmor radius e?ect and the ion Larmor radius e?ect are potentially signi?cant in the calculation of the growth rate of tearing modes. To explore the tearing mode properties, the realistic temperature ratio Ti/Te ? 1 should be adopted in the simulations. Fig.3.8 displays the linear growth rate as a function of kxL with Ti/Te = 1.0, L = 0.5?s for three di?erent guide ?elds Bx0/By0 = 5.0, 3.0, 1.0. The linear growth rate is rapidly reduced as the guide ?eld deceases and kxL corresponding to the peak growth rate shifts to the short wavelength a bit. It seems not to make sense that the linear growth rate is bigger for larger guide ?elds. Nevertheless, we ?nd such di?erence comes from the way of our normalization. Based on the Drake-Lee kinetic theory, Karimabadi et al. [2005] deduced an asymptotic formulation of the growth rate of tearing modes in the Harris current sheet as ?/?i = 2?2mim e (1 + TiT e )1?k 2 xL 2 ?? Bx0B y0 (?eL)3. (3.22) 52 In the above equation the magnetic ?eld is normalized to the antiparallel magnetic ?eld Bx0, di?erent from our normalization to the asymptotic total magnetic ?eld B0. Transforming it to our units, we have ?/?i = 2?2mim e (1 + TiT e )1?k 2 xL 2 ?? (Bx0B y0 + By0B x0 )(?eL)3. (3.23) ? ? (Bx0By0 + By0Bx0), which is increasing with By0Bx0 when By0Bx0 > 1. Figure 3.8: The growth rate as a function of kxL for three di?erent guiding ?elds Bx0/By0 = 5.0, 3.0, 1.0 with Ti/Te = 1.0, L = 0.5?s. 53 The ion kinetic e?ect in the tearing mode is examined in Fig.3.9 with Ti/Te = 1.0, Bx0/By0 = 5.0, kxL = 0.41 for di?erent current widthes. ? represents cases of ?with ions?, while ? represents cases of ?without ions? that has the same meaning of ??xed ions? we discussed above, and the line stands for the ratio of growth rate between them. Two conclusions can be reached immediately. One is that the growth rate decreases rapidly as the current sheet gets wider. The other is that ions play a signi?cant role in the tearing mode, raising the growth rate almost by 50%. Daughton et al. [2005] has demonstrated ions can directly resonant with tearing modes and contribute to around 30% in?uence in the growth rate of tearing modes for thin current sheets. The compressibility of ions can also contribute to the tearing mode. Our simulations show that this term vdi?ni contributes to 10% of the total perturbed current ?j. Finite ion Larmor radius e?ect is another potentially important factor for cases with k??i ? 1. One may have noticed that in the previous section our simulations have shown that the ion e?ect is negligible in the Drake-Lee current sheet. Such inconsistency is due to that in the Drake-Lee current sheet the current is totally carried by electrons, while in the Harris current sheet the ratio between the ion current and the electron current is determined by ji/je = Ti/Te. Because of the low frequency property of tearing modes, the electrostatic e?ect was not considered in the kinetic theory by Drake and Lee [1977a]. Nevertheless, Hoshino [1988] pointed out that the electrostatic ?eld in?uences the induced current, and modi?es the linear growth rate. His analysis shows that in the presence of a guide ?eld, the electrostatic ?eld generated by the compressibility of electrons will increase the growth rate, but the ion e?ect is ignored due to its large Larmor radius. However, previous studies and our simulations have shown that the ?nite ion Larmor radius e?ect and the compressibility of ions are potentially important in the calculation of the growth rate, so it is necessary to investigate the electrostatic e?ect using the GeFi simulation. Fig.3.10 displays the growth rate as a function of Ti/Te with L = 0.5?s, 54 Figure 3.9: The growth rate as a function of L with Ti/Te = 1.0, Bx0/By0 = 5.0, kxL = 0.41. ? : cases with ions; ? : cases without ions; solid lines: the growth rate ratio between cases with ions and cases without ions. Bx0/By0 = 3.0, kxL = 0.41. ? represents cases of ?with electrostatic e?ect?, ? represents cases of ?without electrostatic e?ect?, and the line stands for the growth rate ratio between them. We see that the growth rate is increasing with Ti/Te and the in?uence of electrostatic e?ect depends on the value of Ti/Te. Electrostatic e?ect reduces the growth rate when Ti/Te < 4, but raises the growth rate when Ti/Te > 4. The ion Larmor radius becomes larger as Ti/Te increases, and then Hoshino?s assumption is satis?ed. In Fig.3.11, the contours of typical perturbation structures are plotted with L = 0.5?s, Bx0/By0 = 3.0, kxL = 0.41, Ti/Te = 2.0, where Bx, By, Bz, ni,, ne, ?, A?, Uix, Uiy, Uiz, Uex, Uey, Uez, and je? are repectively three components of the magnetic ?elds, the ion density, the electron density, the electrostatic potential, the parallel component of the vector potential, three components of the ion ?uid velocity, three components of the electron ?uid velocity, and the parallel component of the electron 55 Figure 3.10: The growth rate as a function of Ti/Te with L = 0.5?s, Bx0/By0 = 3.0, kxL = 0.41. ? : with electrostatic e?ect; ? : without electrostatic e?ect; the solid line : the growth rate ratio between them. current. An important quantity in the study of the tearing instability is the singular layer width?, which is de?ned as the region where electrons can free stream along the direction of the perturbed electric ?eld and Landau resonate with the wave. From the contour of Uey, we can see that the singular layer?extends to the current sheet width L with the scale of ?i, and therefore the Drake and Lee kinetic theory assumption ? << L is clearly violated for our parameters. Considering the current parameters ?e = 0.01 and ?i = 0.02, we have de/?e = 1/??e = 10, and di/?i = 1/??e = 7. Since the characteristic thickness of the current sheet is comparable to de and ?i, two factors can contribute to broaden the singular layer: one is the ?nite electron skin depth e?ect, and the other is the ?nite ion Larmor radius e?ect. The linear stage of collisionless tearing modes contains important information of the out-of-plane perturbation By. In the limit of a weak guide ?eld, the perturbation has a quadrupole structure, which has been identi?ed as a key signature of the Hall 56 Figure 3.11: The contours of typical perturbation structure with L = 0.5?s, Bx0/By0 = 3.0, kxL = 0.41, Ti/Te = 2.0. 57 e?ect in fully nonlinear fast reconnection [Sonnerup, 1974; Terasawa, 1983; Karam- abadi et al., 2004;]. Simulations by Daughton et al. [2005] show that the addition of a guide ?eld complicates the out-of-plane perturbation and compresses the spatial width from the ion inertial scale down to the electron kinetic scale. Our simulation shows that the perturbation of By in Fig.3.11 is not of a quadrupole structure but is located in the center of the current sheet with a certain asymmetry. Since di >> ?i ? L for our parameters, the Hall e?ect is excluded form the formation of the perturbation of By. Examinations on the contours of Uix and Uex indicate that the ion and electron motions follow di?erent spatial scales and thus cause the perturbation of By. Since electrons are strongly magnetized, the electron motion is mainly guided along the magnetic ?eld. Due to the di?erent directions of magnetic ?eld below and above the center of current sheet, Uex, the projection of the electron ?uid velocity on the x axis, shows a quadrupole structure and has a peak value around the O-type neutral point. Since the ion Larmor radius is comparable to the current sheet width, ions are able to move across the ?eld lines and be accelerated directly by magnetic tension forces. The structure of Uix is symmetric with respect to the X-type neutral point with a ?/2 phase di?erence compared with Uex. The separation between Uex and Uix forms an in-plane current, and consequently generates the out-of-plane perturbation By. Nonlinear simulations by Ricci et al. [2005] show the asymmetry of Uex results in the compression of electron ?owing along the parallel direction, and leads to a density asymmetry across the dissipation region. Our conclusions here are consistent with the discussion by Drake and Shay [edited by Birn and Priest, 2007]. The comparison among the growth rate of tearing modes as a function of Bx0/By0 obtained from the GeFi simulation (open circles), the eigenmode calculation (solid line) and the Drake and Lee analytical asymptotic matching [1977a] Eq.3.23 (dashed line) is plotted in Fig.3.12. Results based on the asymptotic matching Eq.3.23 devi- ate signi?cantly from the eigenmode calculation and the GeFi simulation, although 58 the general trend of ? versus By0 is consistent with them. This deviation can be understood, since in the cases shown in Fig.3.12 the thickness of current sheet L is comparable to de at the sheet center and thus the conditions assumed in the asymp- totic matching analysis break down. In addition, the ion e?ects are ignored in the matching theory. Figure 3.12: The comparison among tearing mode growth rate as a function of Bx0/By0 obtained from GeFi simulation (open circles), the eigenmode calculation (solid line) and Drake and Lee analytical asymptotic matching Eq.3.23 (dashed line). 3.5 Nonlinear Simulation of Tearing Modes in the Harris current Sheet In this section, we perform nonlinear simulations to study the saturation of tear- ing modes. It has been suggested that when the magnetic island reaches the spatial order of the singular layer width, the electron orbits are altered. This results in the saturation of tearing modes. Two cases, one without ions and the other with ions, are shown ?rst. In these two cases, only parallel component of the vector potential (i.e., A? is dropped due to weak magnetic compression) is kept in the simulation model, 59 which is consistent with the previous theoretical analysis. In the end, the coalescence of multiple tearing modes is studied. Fig.3.13 shows the time evolution of the magnitude of eigenfunction A? (left) and the contour plot of Ay at time t = 20.0 (right). The parameters here are L = 0.5?s, Bx0/By0 = 10.0, Ti/Te = 1.0, and kxL = 0.41. This implies that only the corresponding fastest growing mode is kept in the system. The tearing mode is in the linear growing stage before t = 9??1i with the growth rate ? = 0.83, which is almost the same as ? = 0.812 obtained from the linear simulation. After t = 9??1i , the tearing mode saturates and oscillates with a frequency ? = 2?/8 = 0.785?i. Wan et al. [2005] deduced a relationship between the growth rate ? and the oscillation frequency ? as ?/? = 1.22. Clearly our result ?/? = 0.94 is di?erent from their prediction, for their analysis still relies on de << L. At t = 20??1i , the width of the saturation island is w = 0.44?s, smaller than 2L = 1.0. Considering ? = 0.005, we have de/?e = 14.14 and de/?i = 0.33, so the saturation width of the magnetic island is on the order of de. Based on the Drake and Lee theory, the singular layer thickness in the Harris current sheet can be approximated as [Karamabadi et al., 2005] ? = 1?? 1?(kxL) 2 kxL ?2e L B0 Bx0(1 +Ti/Te). (3.24) For our parameters, the above equation gives ? = 1.1?e, which is much smaller than the saturation island width w obtained from the simulation. Clearly the ?nite electron skin length e?ect that is neglected in the Drake and Lee theory can broaden the singular layer and yield a larger saturation amplitude. The results of simulation with the same parameters above but including the ion motion are depicted in Fig.3.14. It shows the time evolutions of the magnitude of eigenfunctions ne (top left), ni (top right), A? (bottom left), and the contour plot of Ay at time t = 20??1i (bottom right). Two main di?erences can be observed 60 Figure 3.13: The time evolution of the magnitude of eigenfunction A? (left) and the contour plot of Ay at time t = 20??1i (right) with L = 0.5?s, Bx0/By0 = 10.0, Ti/Te = 1.0, and kxL = 0.41. by comparing with the case without ions: one is the faster oscillation frequency ? = 1.1?i; the other is the wider saturation island w = 1.05?i ? 2L. By including the ion in?uence that was neglected in the work of Drake and Lee, it has been suggested that much larger saturation amplitudes comparable to the ion gyroradius could be obtained [Galeev, 1988; Kuznetstova and Zelenyi, 1990]. The strong correlations between the time evolution of ?A? and ni in Fig.3.14 proves that ions play a signi?cant role in the saturation mechanism. The nonlinear evolution of multiple tearing modes is shown in Fig.3.15. Four unstable modes are kept in the domain, with the wave number kxL = 0.11 (mode1), kxL = 0.22 (mode2), kxL = 0.33 (mode3), and kxL = 0.44 (mode4) respectively. Note that the wavelength of mode1 is the same as the system length in the x direction. Since mode4 corresponds to the wavelength with roughly the largest growth rate, four magnetic islands are formed at the initial stage. After t = 8.72??1i , the islands start 61 Figure 3.14: The time evolutions of the magnitude of eigenfunction ne (top left), ni (top right), A? (bottom left), and the contour plot of Ay at time t = 20??1i (bottom right) including ions. 62 to coalesce due to the attraction between currents, and two wider islands become dominated at t = 13.08??1i . Two much smaller plasmoids become evident at t = 15.26??1i . The formation of such secondary island may be due to the thin current near the x-point [Daughton et al., 2005; Karamababidi et al., 2005]. The coalescence of two big islands continues to compress the smaller plasmoid and triggers a forced reconnection. At t = 17.44??1i , a very narrow reconnection layer can be seen between two large islands. At t = 19.62??1i , the system relaxes back to the state that two islands are dominated. Fig.3.16 shows the time evolution of the magnitude of eigenfunction A? for the four modes. Initially only the magnitude of mode4 corresponding to the fastest grow- ing tearing mode is increasing. Around t = 9??1i , the coalescence of the islands transfer the energy from mode4 to mode2. At t = 14??1i , all modes except for mode2 start to increase rapidly. This is due to the generation of small islands near the X-lines as shown in Fig.3.16 at t = 15.26??1i . The magnitude of all modes are deceasing dra- matically after t = 17??1i , because the energy is released by the forced reconnection as discussed above. 3.6 Summary In this chapter, the GeFi simulation model with the realistic electron/ion mass ratio is employed to study tearing mode instabilities under a ?nite guide ?eld. The results are summarized below: 1. The GeFi simulation model is benchmarked with tearing mode instabilities in a simple current sheet. Our results agree with the eigenmode analysis and the previous studies. 2. With the drift kinetic approximation, the eigenmode theory for the Harris cur- rent sheet is deduced. For cases with Ti << Te when the drift kinetic approximation is valid, the GeFi simulation results agree with the eigenmode analysis. 63 Figure 3.15: Contours of vector potential Ay at di?erent times for the case with four unstable modes. 64 Figure 3.16: The time evolution of the magnitude A? for the four modes. 65 3. The inhomogeneity terms a?ect the structure of tearing modes, and cannot be neglected when By0 ? Bx0. 4. For the temperature ratio Ti/Te ? 1, the simulation studies reveal the re- lationship between the growth rate of tearing modes and kxL, Bx0/By0, Ti/Te, and L/?s. 5. The electrostatic e?ect is negligible on the calculation of the growth rate of tearing modes. 6. In the general cases in which the ion gyroradius is comparable to the current sheet width, the ion e?ects are found to play a signi?cant role on the calculation of the growth rate of tearing modes. In such cases, the Drake and Lee theory [1977a] and the eignemode analysis based on the drift kinetic assumption are not invalid. 7. The nonlinear GeFi simulation shows that the ion e?ect can broaden the width of saturated magnetic islands from the electron dynamic spatial scale to the ion gyroradius. 8. The case with multiple tearing modes shows that the coalescence of magnetic islands can form a larger saturation island than the case of a single tearing mode. During the process of coalescence, small magnetic islands are formed around the X point. The small-scale magnetic reconnection can occur when the small magnetic islands are compressed by the large magnetic islands. This process results in fast energy release. 66 Chapter 4 Generation of Low Frequency Electromagnetic Waves by Magnetic Reconnection under a Finite Guide Field In the collisionless plasma, magnetic reconnection is usually triggered by mi- cro instabilities. The change of topology of the magnetic ?eld lines and and the perturbations in the di?usion region will lead to the generation of low frequency elec- tromagnetic waves. Previous 2D simulations of steady-state reconnection have found these waves to be present in the form of MHD discontinuities, as also observed in the magnetosphere as well as the interplanetary space [Sonnerup et al., 1981; Phan et al., 2006]. In this Chapter, we investigate the 3D structure of the reconnection layer. The parameters chosen will be suitable for the Earth?s magnetosphere. It is found that kinetic Alfv?en waves are present in the 3D reconnection. First the simulation model is introduced. Then a brief review of the kinetic Alfv?en wave is given. Our results show the kinetic e?ects play an important role in the 3D reconnection geometry. The simulation results will be presented in the third section. 4.1 3D Hybrid Simulation Model In the hybrid code, ions are treated as discrete particles and electrons are treated as a massless ?uid. The code was developed by Swift[1995], which has been used in a general curvilinear coordinate system to model the Earth?s magnetosphere. In this study we use a Cartesian coordinate system to model the generation of low frequency electromagnetic waves. 67 In the hybrid code, quasi charge neutrality is assumed. The equation for the ion motion is given by mi qi dvi dt = E+vi ?B??(Vi ?Ve) (4.1) where vi is the ion particle velocity, E is the electric ?eld in units of ion acceleration, B is the magnetic ?eld in units of the ion gyrofrequency, ? is the collision frequency which is used to model the resistivity at the X-line, and vi and Ve are the bulk of ?ow velocities of ions and electrons, repectively. The electron momentum equation is written in the form E = ?Ve ?B??(Ve ?Vi) +?Pe, (4.2) where Pe is the electron pressure. The electron ?ow speed is evaluated from Ampere?s law, Ve = Vi ? ??B?N , (4.3) where ? = 4?e2/mic2, e is the electron charge, mi is the ion mass, and N is the particle number density. The magnetic ?eld is advanced in time using Faraday?s law ?B ?t = ???E, (4.4) where the electric ?eld is calculated from Eq.4.2. In the simulation the ion particle velocity is updated at half time steps with a second-order accuracy. The magnetic ?eld and the particle positions are advanced at the integer time step. The time step is chosen to satisfy the Courant condition with respect to the whistler mode. In the simulation, the magnetic ?eld is advanced 10 time steps for every time step the ion velocities are advanced. Here, the magnetosphere plasmas are considered, in particular for the symmetric reconnection as usually occurring in the magnetostail. The simulation domain is a cube with x being normal to the current sheet and z being the antiparallel magnetic 68 ?eld direction. Di?erent with the previous chapters, we choose the z direction as the antiparallel direction, which is generally used in the magnetospause research. The initial pro?le of the z component magnetic ?eld is given by Bz(x) = Bz0(x)tanh(x/?) (4.5) where Bz0 is the z-component magnetic ?eld in the lobes and ? is the half-width of initial current sheet. The pro?le of the y-component magnetic ?eld, i.e., the guide ?eld in the lobes, is assumed to be By = By0. (4.6) The initial x-component magnetic ?eld is Bx = 0. The pro?le of ion thermal pressure is determined from the total pressure balance across the current sheet. A constant, isotropic temperature with Ti = Te is assumed, where Ti is the ion temperature and Te the electron temperature. The ion number density is thus given by N(x) = N0(1 + 1? 0 [1?B2z0/B20tanh2(x/?)?B2y0/B20]) (4.7) where B0 = (B2z0 +B2y0)12 is the total magnetic ?eld in each lobe. We study the evolution of a spontaneous reconnection [Scholer, 1989]. The dy- namics are not driven. A constant resistivity in time is imposed in the system to trigger the reconnection. The resistivity is modeled through a collisional term in the ion equation of motion and the electron momentum equation. The collision frequency is assumed to be ? = ?0exp(?(x2 +z2)/?20)exp(?y2/?20) (4.8) where ?0 ? 0.1?2?i and ?i is the ion gyrofrequency in the lobes. We de?ne ?0 as the half-length of X-line. 69 In our simulation, the time is in units of ??1i . For the results shown in this thesis the magnetic ?eld is expressed in units of B0, the ion number density in units of N0, and the temperature in units of the lobe temperature T0, where the subscript ?0? represents the quantities in the lobes. The velocity is normalized to the lobe Alfv?en speed VA0, and the spatial coordinate is normalized to di the lobe ion inertial length. The conductor boundary is used in the x direction, , while the periodic bound- aries are used in both the z and the y directions. To limit the boundary e?ect, the simulation stops when the waves propagate near the boundaries. The domain size X ?Y ?Z is 128di ?128di ?256di with the grid size ?x = 0.5di, ?y = 2.0di, and ?z = 1.0di. Almost 100 particles per cell are placed in the lobe region. 4.2 Properties of Kinetic Alfv en Waves The Alfv?en wave is one of the most important wave modes in magnetized plas- mas. It exists due to the balance between the magnetic ?eld tension and ion inertia. These waves were ?rst theoretically predicted by H.Alfv?en and are now called shear Alfv?en waves because the perturbations of the magnetic ?eld are perpendicular to the ambient ?eld and wave vector. When the perpendicular wavelengths are comparable to the ion gyroradius ?i, the properties of shear Alfv?en waves need to be modi?ed. In a plasma with intermediate beta values (me/mi < ?e << 1)), the modi?ed waves are usually refereed as kinetic Alfv?en waves. Starting from the two ?uid-MHD equations, one can deduce the dispersion rela- tionship of kinetic Alfv?en waves as [Stix, 1992] ?2 = k2?V 2A(1 + (1 + TiT e )k2??2i). (4.9) Obviously the kinetic Alfv?en wave is dispersive, while the shear Alfv?en wave is non- dispersive. 70 Another important di?erence between the shear Alfv?en wave and the kinetic Alfv?en wave is that the kinetic Alfv?en wave is compressible and generates an electric ?eld parallel to the ambient magnetic ?eld by E? = ? Teen e ??ne. (4.10) It is suggested that such electric ?eld plays an important role in the energy transport processes and particle heating/acceleration of space plasmas. 4.3 Simulation Results Two typical cases with By0 = 0.5Bz0 are presented in this section. One possesses an in?nitely long X-line, and the other possesses with a ?nite X-line. The cases with an in?nite X-line corresponds to ky ? 0. Since the propagation of waves is almost limited in the x-z plane, results similar to the 2D simulations are expected although some di?erences are found to exist. For the case with a ?nite X-line, the spatial e?ects along the y direction are found to lead to ky ?= 0, which results in the propagation of waves very di?erent from those in 2D simulations. 4.3.1 Case with an In nite X-line: 2?0 = ? Fig.4.1 shows the contour plots of the magnetic ?eld Bx, By, Bz, the ion ?uid velocity Vx, Vy, Vz, the temperature T, the particle number density N, and the con?guration of magnetic ?eld lines in the x-z plane located aty = 0 whent = 150??1i . A plasma bulge and a trailing quasi-steady reconnection layer are formed in the out?ow region of magnetic reconnection, as shown in the magnetic ?eld and plasma plots in Fig.4.1. The magnetic ?eld decreases and the plasma temperature increases from the lobe to the center of the reconnection layer. Since ky ? 0, the ?eld and 71 plasma properties in the x-z planes of other y positions resemble the results shown here. Figure 4.1: The contour plots of Bx, By, Bz, Vx, Vy, Vz, T, N, and con?guration of magnetic ?eld lines. The solid line corresponds to the position at z = 40di, along which there appears a quasi-steady structure. The solid line at z = 40di indicates the boundary of the leading bugle. The solid line z = 40di shows the transit boundary between the plasma bugle and the steady state layer. 72 Figure 4.2: The spatial pro?les of Bx, By, Bz, N, Vx, Vy, Vz, T?, T?, J?, and B at z = 40di and hodograms of the tangential magnetic ?eld (Bz-By). 73 Fig.4.2 depicts the spatial pro?les of the magnetic ?eld Bx, By, Bz, density N, the ion ?ow velocity Vx, Vy, Vz, the parallel temperature T?, the perpendicular temperature T?, the parallel current density J?, and the magnitude of magnetic ?eld B at z = 40di, and the hodograms of the tangential magnetic ?eld (Bz-By). Note that z = 40di is located in the steady reconnection layer as indicated in Fig.4.1. We can see that two discontinuities bound the steady reconnection layer with the wave front at x = ?8.5di and the width of ? 2.3di. Across the discontinuities, the plasma is accelerated so that Vz is changed from ? 0 in the upstream to ? 0.75VA0 in the downstream corresponding to the change of the Alfv?en speed ? 0.7VA0. The hodogram shows that the tangential ?eld rotates by an angle of around? ?/2 through each discontinuity. These pair of discontinuities are roughly rotational discontinuities. Notice that the magnetic ?eld strength slightly decreases and the ion density increases across the rotational discontinuity. The changes in the magnetic ?eld and ion density at the rotational discontinuity are due to the increase of the parallel temperature T? > T? because of the streaming ions along the magnetic ?eld lines. The slow shock, however, is too weak to be identi?ed in the simulation because of T? > T?. All results shown here are consistent with the previous 2D simulations. To study the 3D e?ect, we examine the variations along the y direction. The variation of a certain quantity f is calculated as f(x,y,z) ? ?f(x,z), where ?f(x,z) is the average along the y direction at ?xed positions (x,z). The top two ?gures in Fig.4.3 present the contour plots of variations of Bx and N along the y direction in the x-z plane at y = 0, and the bottom two ?gures present the contour plots of the corresponding variations in the y-z plane at x = 4di as indicated by the solid line in the top ?gures. At time t = 150??1i , one can see that the variations are signi?cant in the transition region between the leading bulge and the quasi-steady structure. In the y-z plane at x = 4di the structure of the variation of Bx stretches out along the initial local magnetic direction b. Since the x direction is perpendicular to the initial 74 magnetic ?eld, the variation of Bx represents the shear waves, while the variation of N represents the compressible waves. Figure 4.3: Top ?gures: the contours of variations of Bx and N along the y direction in the x-z plane at y = 0; bottom ?gures: the contours of variations of Bx and N along the y direction in the y-z plane at x = 4di. The solid line along z = 52di in the top left ?gure indicates the boundary of the leading bugle. The solid lines in the top ?gures indicate the y-z plane at x = 4di; the solid lines in the bottom ?gures show the local initial magnetic ?eld direction. Fig.4.4 shows the contours of variations along the y direction of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ey, J?, E? in the y-z plane at x = 3di when t = 150??1i . In these 75 structures, the shear components Bx and Vx, Ey, and the parallel current density J? are well correlated. In Fig.4.3, one may have noticed that the spatial scale of the variation of Bx is ? 3di or ? 10?i, considering to ?i = 0.1. Since kx?i ? 2?10?i?i ? 0.62, certain kinetic e?ect are expected, which is justi?ed by the strong correlation between the parallel electric ?eld E? and Bx. Figure 4.4: The contours of variations along the y direction of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ey, J?, E? in the y-z plane at x = 3di when t = 150??1i . 76 4.3.2 Case with a Finite X-line: 2?0 = 8di In reality, the length of X-line cannot be in?nite. Shay et al. [2010] have shown that 3D structures develop in a spontaneous reconnection, corresponding to spatially isolated sites of reconnection with scale lengths around 10di. In the following, we show a case in which the length of the X-line is set as 2?0 = 8di to study the ?nite length X-line e?ect. The top two ?gures in Fig.4.5 show the con?gurations of magnetic ?eld lines near the current sheet at t = 50??1i (left) and t = 80??1i (right). At t = 50??1i , some magnetic ?eld lines have been reconnected. Since the length of the X-line is ?nite, the process of reconnection is limited in a ?nite region. At t = 80??1i , the reconnected ?eld lines are dragged towards ?z directions by the accelerated plasmas and are twisted at the same time. It is due to that the large-amplitude perturbations generated by the reconnection and their consequent propagations change the local topology of magnetic ?eld lines. The contour plots of Bx (left) and Vx (right) in the x-z plane at y = 25di of t = 80??1i are presented in the bottom two plots. We can see that the shear perturbations of Bx and Vx are evident in the plane far from the reconnection region, and they are well correlated. The perturbation regions are associated with the reconnected ?eld lines, while the ?eld lines outside of the perturbations stay unreconnected. It indicates that the perturbations caused by the reconnection propagate into the sites far from the di?usion region along the magnetic ?eld lines. We also notice that the perturbations in the x-z plane of y = 25di are just located at one side of the current sheet. Crossing the current sheet along the red line indicated in the top left plot, we observe that the ?eld lines are connected on the side of current sheet and unconnected on the other side. Such asymmetry is simply due to the ?eld line geometry geometry of this ?component reconnection?, i.e., the reconnection with a ?nite guide ?eld. 77 Figure 4.5: The con?guration of magnetic ?eld lines near the current sheet at t = 50??1i (the top left) and t = 50??1i (the top right), and the contour plots of Bx (the bottom left) and Vx (the bottom right) in the xz plane of y = 25di at t = 50?i. 78 To investigate the properties of propagating waves, we now examine the propaga- tion speed of the wave front. Fig.4.6 shows the relationship between the time t and the distance L from the center of the di?usion region (0,0,0) to the position of the wave front. Obviously t and L are linearly correlated with a linear ?t of L = 0.89t?11.99. Therefore the propagation speed is 0.89VA0 that is close to the lobe Alfv?en speed VA0. Note that the propagation is mainly located outside the edge of the current sheet, where the Alfv?en speed is slight smaller than the lobe Alfv?en speed. Figure 4.6: The relationship between time t and the distance L from the center of the di?usion region (0,0,0) to the position of wave front. In Fig.4.7, the contour plots of Bx and N in the x-z plane of y = 16di at t = 80?i are presented. The contours of Bx exhibit clear structure of perturbations, while the density perturbations are also observed at the same position. Since the perturbations propagate from the di?usion region that has the characteristic scale length ? ?i, it makes sense that the spacial scale length of the perturbation Bx is ? 3di, corresponding to kx?i ? 0.62. We expect the waves to have certain kinetic e?ects. 79 Figure 4.7: The contour plots of Bx and density N in the x-z plane of y = 16di at t = 80?i. The red dot stands for the position (4, 16, -20). Fig.4.8 shows the time evolutions of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ez, E?, B, Ex?b, B?, and the hodogram of Bx ?Bx?b at the position (x,y,z) = (4,16,?20), as indicated by the red dot in Fig.4.7. Among them, B? is calculated by B ? b, where b is the initial magnetic ?eld direction, and the subscript x ? b stands for the direction perpendicular to x and b. We see that the transverse components Bx, Vx, and Ey, Ez are well correlated. The ratio Ex?b/Bx is ? 0.8VA0 close to the lobe Alfv?en speed. During the whole time, the magnitude of the magnetic ?eld B stays almost constant, as well as the parallel component of the magnetic ?eld B?. The wave appears an Alfv?en mode. Moreover, the spatial scale length of perturbations along the x direction is comparable to ?i, and E? is correlated with the perturbation of the local density that develops after the arrival of these perturbations. The slight increase of density is also due to the ion kinetic e?ect. Recall that due to the ion kinetic e?ect, the Alfv?en wave becomes compressible and generates an electric ?eld parallel to the local magnetic ?eld. Since b points into the paper in the hodogram of Bx?Bx?b, the waves rotate along the counterclockwise direction and have a lefthand polarization. 80 Figure 4.8: The time evolution of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ez, E?, B, Ex?b, B?, and the hodogram of Bx ?Bx?b at the position (4, 16, -20) as indicated by the red dot in Fig.4.7. 81 The contour plots of Vz at x = 0 (left) and x = 2di (right) are shown in Fig.4.9. In the center of the current sheet, the plasma is accelerated by the tension fore caused the reconnected ?eld lines, and the accelerated plasma is mainly towards ?z directions. Outside the current sheet, the Alfv?en wave propagates along the ?eld lines, and the plasma is accelerated through the wave during its propagation. Therefore the accelerated plasma is mainly along the ?eld lines as shown in Fig.4.8. Figure 4.9: The contour plots of Vz at x = 0 (left) and x = 2di (right) at t = 80??1i . Fig.4.10 shows the contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane of y = 0. Note the x-z plane of y = 0di is in the center of X-line. We can see that the resulting structures here are very di?erent from the 3D run with an in?nitely long X-line. No clear leading bulge exists during the reconnection, although an accelerated plasma layer is formed near the center of current sheet where the parallel heating is observed. Moreover, we do not ?nd the discontinuity fronts on the edges of the reconnection layer, as in 2D simulations and 3D simulations for cases with an in?nite X-line. In 3D simulation with the limited X-line, the component Bx is a shear perturbation to the local magnetic ?eld. In such case with ky ?= 0, such 82 3D perturbations propagate with a group velocity along the ?eld lines. Therefore, no steady discontinuities or shocks are formed. Figure 4.10: The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane of y = 0 when t = 150??1i . Fig.4.11 shows the contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane of y = 20di. The propagations of perturbations generated by the reconnection 83 accelerate the local plasma and form a structure similar to that in the x-z plane of y = 0di, although no reconnection actually occurs in this x-z plane. Across the current sheet, the layers are asymmetric due to the reason discussed already. Again, the parallel heating occurs in the layers. The enhancement of density is found in the layers because of the acceleration of the plasma. Cases with various lengths of X-line are also simulated. It is found that the critical X-line length for the existence of steady-discontinuities is ? 30di. 4.4 Summary We have used the 3D hybrid simulation model to study the generation of low frequency waves by the reconnection under a ?nite guide ?eld. The simulation results are summarized below. For the case with an in?nite X-line: 1. 3D simulation results are similar to those of 2D cases. Two quasi-steady rota- tional discontinuities are formed behind a plasma bulge. The plasmas are accelerated across the rotational discontinuities. 2. The contour plots show that there are perturbation structures in the transition region between the leading bulge and the steady discontinuities. The transverse perturbations lead to the presence of short-wavelength ?eld-aligned structures, in which E? ?= 0 due to the ion kinetic e?ects. For the cases with a ?nite X-line: 3. It is found that the perturbations caused by the reconnection propagate along the local magnetic ?eld lines. 4. The wavefront of perturbations propagates with the local Alfv?en speed. The shear components dominate the wavefront of perturbation which are seen as locations far from the reconnection region, and the ratio between the corresponding transverse electric ?eld and magnetic ?eld components is almost equal to the local Alfv?en speed. 84 Figure 4.11: The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane of y = 20di when t = 150??1i 85 Thesewavesareidenti?edaskineticAlfv?enwaves, inwhichthewavenumberkx?i ? 1. Our simulation indicates that the ?nite X-line length leads to a structure of the reconnection layer that is very di?erent from the cases with an in?nite X-line, in which steady discontinuity fronts are dominant. 5. Parallel electric ?eld E? is generated in the kinetic Alfv?en waves, and the structure of E? has a strong correlation with the structures of the shear components. 6. The critical length of X-line that leads to the Petschek-type steady shocks or discontinuities is found to be ? 30di. Our study indicates that a ?nite extent of X-line alters the wave structure of the reconnection layer. The generation of kinetic Alfv?en waves in the general cases of reconnection may provide an important mechanism for the particle heating and acceleration as these waves propagate along ?eld lines to the ionosphere. 86 Chapter 5 Summary Magnetic ?eld reconnection provides an e?cient mechanism for the transfer of magnetic energy to plasma heat and kinetic energy. The kinetic physics of the fast reconnection in collisionless plasmas, however, is still poorly understood due to the complex plasma dynamics that is largely caused by the disparate spatial and tem- poral scales between the electron and ions. The 2D Petschek reconnection model is the ?rst fast reconnection model, which provides a mechanism of reconnection that is fast enough to explain the phenomena caused by the reconnection. In order to fully understand the validity of the Petschek model, two basic questions should be answered. One is how the fast reconnection is triggered by micro instabilities, and the other is how the reconnection layer will be modi?ed or altered in the general 3D geometries in reality. To tackle the problem due to the vast di?erence between the electron and ion scales and investigate the physics of reconnection from the micro to macro scales, we use the newly developed GeFi particle scheme to simulate the reconnection in a 2D current sheet with a ?nite guide ?eld and a realistic mass ratio. In this thesis, the GeFi model of Lin et al. [2005] has been improved for the nonuniform current sheets. Tearing mode instabilities are believed to play a fundamental role to trigger the fast reconnection. As a ?rst step to apply the GeFi model to the reconnection research, we have developed a linear eigenmode theory of the tearing instability to benchmark the improved GeFi scheme. The GeFi results have also been compared with the asymptotic matching theory of Drake and Lee [1977a]. Furthermore, the kinetic e?ects that are missing in previous theoretical studies of the tearing mode 87 have been examined through the GeFi simulations. After the code benchmark, we have then carried out linear and nonlinear GeFi simulations of the tearing instability in magnetic reconnection. While the application of the GeFi model to magnetic reconnection has produced interesting initial results, the complete investigation of 3D nonlinear magnetic re- connection using the GeFi model is still beyond the capability of the current GeFi scheme. To explore the properties of low-frequency waves generated by magnetic re- connection, a hybrid code has been used to simulate the 3D structure of reconnection layer under a ?nite guide ?eld. The main results of this thesis are summarized: The improved GeFi model: 1. Our GeFi particle simulation model, in which the electrons are treated as gyrokinetic particles and ions are treated as fully kinetic particles, has been improved and modi?ed to allow the the existence of modes with wavelengths on the same scale of the background nonuniformity. 2. The linearized GeFi scheme is benchmarked for uniform plasmas, and the simulation results using the ?f method show that the model can accurately resolve the physics ranging from Alfv?en waves to lower-hybrid/whister waves, for k? << k? and ? << ?e. The tearing mode instability under a ?nite guide ?eld: 1. The GeFi simulation model is benchmarked with tearing mode instabilities in a simple current sheet and the Harris current sheet. Our results agree with the eigenmode analysis and the previous studies under simpli?ed conditions of these the- oretical models. The inhomogeneity terms a?ect the structure of tearing modes, and cannot be neglected when By ? Bx. 2. In the general cases in which the ion gyroradius is comparable to the current sheet width, the ion kinetic e?ects are found to play a signi?cant role on the growth 88 rate of tearing modes. In such cases, the Drake and Lee theory [1977a] and the eignemode analysis based on the drift kinetic assumption are invalid. 3. The nonlinear GeFi simulation shows that the ion e?ect can broaden the width of saturated magnetic islands from the electron dynamic spatial scale to the ion gyroradius. The case with multiple tearing modes shows that the coalescence of magnetic islands can form a larger saturation island than the case of a single tearing mode. The generation of low frequency waves by reconnection under a ?nite guide ?eld: 1. For the case with a in?nitely long X-line, 3D simulation results are similar to those of 2D cases. In addition to the steady rotational discontinuities in the out?ow region, wave perturbations are also found between between the leading plasma bulge and the steady discontinuities. The transverse perturbations lead to the presence of short-wavelength ?eld-aligned structures, in which E? ?= 0 due to the ion kinetic e?ects. 2. For the case with a ?nite extent of X-line, it is found that the perturbations caused by the reconnection propagate along the local magnetic ?eld lines. The wave- front of perturbations propagates with the local Alfv?en speed. The shear components dominate the wavefront of perturbation which are seen at locations far from the re- connection region. These waves are identi?ed as kinetic Alfv?en waves, in which the wave number k??i ? 1. Parallel electric ?eld E? is generated in the kinetic Alfv?en waves. Our simulation indicates that the ?nite X-line length leads to a structure of the reconnection layer that is very di?erent from the cases with an in?nite X-line as as in the Petschek model, in which steady discontinuity fronts are dominant. 89 Bibliography [1] Biernat, H.K., Heyn, M.F., Semenov, V.S., Farrugia, C.J., The structure of reconnection layers: Application to the Earths Magnetopause, J.Geophys.Res., 94, 287, 1989 [2] Birn, J., et al., Geospace Environmental Modeling(GEM) magnetic reconnection challenge, J.Geophys.Res., 106, 3715, 2001 [3] Birn, J., E. R. Priest, Reconnection of Magnetic Fields: Magnetohydrodynamics and Collisionless Theory and Observations, Cambridge, NewYork, 2007 [4] Biskamp, D., Magnetic reconnection via current sheets, Phys.Fluids, 29, 1520, 1986 [5] Burden, R.L., Faires, J.D., Numerical Analysis , Brooks/Cole, 2005 [6] Chen, Y., Parker, S.E., A ?f particle method for gyrokinetic simulations with kinetic electrons and electromagnetic perturbations, J.Comput.Phys, 189, 463, 2003 [7] Drake, J.F. Lee, Y.C., Kinetic theory of tearing mode instabilities, Phys.Fluids, 20, 1341, 1977a [8] Drake, J.F. Lee, Y.C., Nonlinear evolution of collisionless and semicllisional tear- ing modes, Phys.Rev.Lett. 39, 453, 1977b [9] Daughton, W., Electromagnetic properties of the lower-hybrid drift instability in a thin current sheet, Phys. Plasmas 10, 3103, 2003 [10] Daughton, W., Karimabadi, H., Kinetic theurgy of collisionless tearing at the magnetopause, J.Geophys.Res., 109, A0317, 2005 [11] Fireman, E.A., Chen., L., Nonlinear gyrokinetic equations for low-frequency elec- tromagnetic waves in general plasma equilibria, Phys.Fluids, 25, 502, 1982 [12] Furth, H.P., Killeen, J., Rosenbluth, M.N., Finite resistivity instabilities of a sheet pinch, Phys.Fluid, 6, 459, 1963 [13] Furth, H.P., Rutherford, P.H., Selberg,H., Tearing mode in the cylindrical toka- mak, Phys.Fluid, 16, 1054, 1973 90 [14] Heyn, M.F., Biernat, H.K., Semenov, V.S., The structure of magnetic reconnec- tion layer, J.Geophys.Res., 90, 1781, 1981 [15] Hans, J.P., Poedts, S. Principles of Magnetohydrodynamics: With Applications to Laboratory and Astrophysical Plasmas, Cambridge Press., 2004 [16] Harris, E.G., On a plasma sheath separating regions of one directional magnetic ?eld, Nuovo Cimento 23, 115, 1961 [17] Hockney, R.W., Eastwood, J.W., Computer Simulation Using Particles, Mcgraw- Hill, New York, 1981 [18] Hoshino, M., The electrostatic e?ect for the collisionless tearing mode, J.Geophys.Res., 92, 7368, 1988 [19] Landau, L.D., Lifshitz, E.M., Electrodynamics of Continuous Media, Pergamon Press., 1960 [20] Katanuma, I., Kamimura, T., Simulation studies of the collisionless tearing in- stabilities Phys.Fluid, 23, 2500, 1980 [21] Karimabadi, H., Daughton, W., Quest, B.K., Physics of saturation collisionless tearing mode as a a function of guide ?eld, J.Geophys.Res., 110, A010749, 2005 [22] Karimabadi, H., Huba, J., KraussVarban, D., On the generation and structure of the quadrupole magnetic ?eld in the reconnection process, Geophys.Rev.Lett. 31, L070806, 2004 [23] KrasussVarban, D., Omidi, N., large-scale hybrid simulations of the magnetotail during reconnection, Geophys.Res.Lett., 22, 3271, 1995 [24] Karimabadi, H., Daughton, W., Quest, K.B., Antiparallel versus component merging at the magnetopause: Current bifurcation and intermittent reconnec- tion, J.Geophys.Res., 10, A03213, 2005b [25] Karimabadi, H., Daughton, W., Quest, K.B., Physics of saturation of collisionless tearing mode as a a function of guide ?eld, J.Geophys.Res., 110, A03214, 2005c [26] KrasussVarban, D., Omidi, N., large-scale hybrid simulations of the magnetotail during reconnection, Geophys.Res.Lett., 22, 3271, 1995 [27] Lee, W.W., Gyrokinectic approach in particle simulation, Phys.Fluids, 26, 556, 1983 [28] Lee, W.W., Gyrokinectic particle simulation model, J.Comput.Phys., 72, 243, 1987 [29] Lin, Y., Lee, L.C., Structure of reconnection layers in the magnetosphere, Space.Sci.Rev., 65, 59, 1994a 91 [30] Lin, Y., Lee, L.C., Reconnection layer at the ?ank mangetopause in the presence of shear ?ow, Geophys.Res.Lett., 21, 59, 1994b [31] Lin, Y., Wang, X.Y., Chen, L., Lin, Z.H., A gyrokinetic electron and fully kinetic ion plasma simulation model, Plasma Phys.Control.Fusion, 47, 657, 2005 [32] Lin, Y., Wang, X.Y., Three-dimensional global hybrid simulation of dayside dynamics associated with the quasi-parallel bow shock, J.Geophys.Res., 110, A12216, 2005 [33] Lin, Y., Wang, X.Y., Three-dimensional global hybrid simulation of dayside dynamics associated with the quasi-parallel bow shock, J.Geophys.Res., 115, A04208, 2010 [34] Lin, Y., Lee, L.C., Chaos and ion heating in a slow shock, Geophys.Res.Lett, 18, 1615, 1991 [35] Lin, Y., Lee, L.C., Kennel, C.F., The role of intermediate shocks in magnetic reconnection, Geophys.Res.Lett, 19, 229, 1992 [36] Lin, Y., Lee, L.C., Simulation study of the Riemann problem associated with magnetotail reconnection, J.Geophys.Res., 100, 19,277, 1995 [37] Lin, Y., Swift, D.W., A two dimensional hybird simulation of the magnetotaill reconnection, J.Geophys.Res., 101, 19,859, 1996 [38] , Lin, Y., Wang, X.Y., Chen, L., Lu, X., and Kong, W., An improved gyrokinetic electron and fully kinetic ion particle simulation scheme: benchmark with a linear tearing mode, Plasma Phys. Control. Fusion., 53, 054013, 2011 [39] Lottermoser, R.F., Scholar, M., Undriven magnetic reconnection in magnetohy- drodynamics and Hall magnetohydrodynamics, J.Geophys.Res., 102, A3, 4875, 1997 [40] Lottermoser, R.F., Scholar, M., Matthews, A.P.,Ion kinetic e?ects in magnetic reconnection, J.Geophys.Res., 103, 4547, 1998 [41] Masaaki, Y., Kulsrud, R., Ji, H.T., Magnetic Reconnection, Rivews of Modern Physics, 82, 603, 2010. [42] Mandt, M.E., Denton, R.E., and Drake, J.F., Transition to whistler mediated magnetic reconnection, Geophys.Res.Lett, 21, 73, 1994 [43] Parker, E.N., Sweet?s mechanism for merging magnetic ?eld in conducting ?uids, J.Geophys.Res., 62, 509, 1957 [44] Parker, S.E., Lee, W.W., A fully nonlinear characteristic method for gyrokinetic simulation, Phys.Fluids B, 5, 77, 1993 92 [45] Petschek, H.E., Magnetic ?eld annihilation, in AAS-NASA Symposium on the Physics of Solar Flares, NASA Spec. Publ., SP-50, 425, 1964 [46] Priest, E., Forbes, T., Magnetic Reconnection: MHD Theory And Applications, Cambridge Press., 2006 [47] Prtchett, E., Collisionless magnetic reconnection in a three open system, J.Geophys.Res., 106, 3783, 2001 [48] Prtchett, E., Onset of saturation of guide-?eld magnetic reconnection, Phys.Plasmas, 12, 062301, 2005 [49] Richter, P., Sholer, M., on the stability of rotational discontinuities, J.Geophys.Res., 16, 1257, 1989 [50] Richter, P., Sholer, M., on the stability of rotational discontinuities with large phase angles, Adv.Space.Res., 11, 111, 1991 [51] Ricci, P., Brackbill, U. B., Doughton, W., and Lapenta, G., New role of the lower-hybrid drift instability in the magnetic reconnection, Phys. Plasmas 12, 5501, 2005 [52] Ricci, P., Brackbill, U. B., Doughton, W., and Lapenta, G., Collisionless mag- netic reconnection in the presence of a guide ?eld, Phys. Plasmas 11, 4102, 2008 [53] Rogers, B.N., Denton, R.E., Drake, J.F., Shay, M.A., The role of dispersive waves in collisionless magnetic reconnection, Phys.Rev.Lett., 97, 195004, 2001 [54] Rutherford, P.H., Frieman, E.A., Drift instabilities in general magnetic ?eld con?guration, Phys.Fluids, 11, 569, 1968 [55] Rutherford, P.H.,Nonlinear growth of tearing mode, Phys.Fluids, 16, 1903, 1973 [56] Sweet, P.A., The neutral point theory of solar ?ares, in Electromagnetic Phe- nomena in Cosmic Physics, 135, 1958. [57] Swift, D.W., Lee, L.C., Rotational discontinuities and the structure of the mag- netopause, J.Geophys.Res., 88, 111, 1983 [58] Swift, D.W., Use of hybrid code to model the Earth?s magnetosphere, Geo- phys.Res.Lett., 22, 311, 1995 [59] Sato, T.,Strong plasma acceleration by slow shocks resulting from magnetic re- connection, J.Geophys.Res., 84, 7177, 1979 [60] Shi, Y., Lee, L.C., Structure of the reconnection layer at the dayside magne- topause, Planet.Space.Sci., 38, 437, 1990 [61] Scholer, M., Undriven magnetic reconnection in an isolated current sheet, J.Geophys.Res., 94, 15099, 1989 93 [62] Sonnerup, B., Magnetic ?led reconnection, in Solar System Plasma Physics, 46, Elsevier, New York, 1974 [63] Taylor, J.B., Hastie, R.J., Stability of general plasma equilibria-I. Formal theory, Phys.Fluids, 10, 479, 1968 [64] Terasawa, T., Hall current e?ect on tearing mode stability, Geophys.Res.Lett., 10, 475, 1983 [65] Ugai, M., Self-consistent development of fast reconnection with anomalous plasma resistivity, Plasmas.Phys.Control.Fusion, 26, 1549, 1984 [66] Ugai, M., Shumizu, T., Computer studies on the fast reconnection mechanism in a sheared ?eld geometry, Phys.Plasmas, 1, 296, 1993 [67] Wang, X.Y., Lin, Y., Chen, L., Lin, Z.H., A particle simulation of current sheet instabilities under ?nite guide ?eld, Phys.Plasmas 15, 072103, 2008 [68] Wu, C.S., Zhou, Y.M., Tsai, S.T., and Gao, S.C., A kinetic cross ?eld streaming instability, Phys. Fluids 26, 1259,1983 [69] Wesson, J.A., Finite resisitivity instability of a sheet pinch, Nuclear Fusion, 18, 130,1966 [70] Wan, W., Chen, Y., and Parker, S.E., Gyrokinetic ?f simulation of the collision- less and semi-collisional tearing mode instability, Phys.Plasmas 12, 012311,2005 [71] Yan, M., Lee, L.C., Priest, E.R., Fast magnetic reconnection with small shock angles, J.Geophys.Res., 97, 8277, 1992 94