Collisionless Reconnection in the Large Guide Field Regime: Gyrokinetic Versus Particle-in-Cell Simulations
J. M. TenBarge, W. Daughton, H. Karimabadi, G. G. Howes, W. Dorland
CCollisionless Reconnection in the Large Guide Field Regime: GyrokineticVersus Particle-in-Cell Simulations
J. M. TenBarge, a) W. Daughton, H. Karimabadi,
3, 4
G. G. Howes, and W. Dorland IREAP, University of Maryland, College Park, MD 20742 Los Alamos National Laboratory, Los Alamos, NM 87545 SciberQuest, Del Mar, CA 92014 Department of Electrical and Computer Engineering, UCSD, La Jolla, CA 92093 Department of Physics and Astronomy, University of Iowa, Iowa City, IA 52242 (Dated: 16 October 2018)
Results of the first validation of large guide field, B g /δB (cid:29)
1, gyrokinetic simulations of magnetic recon-nection at a fusion and solar corona relevant β i = 0 .
01 and solar wind relevant β i = 1 are presented, where δB is the reconnecting field. Particle-in-cell (PIC) simulations scan a wide range of guide magnetic fieldstrength to test for convergence to the gyrokinetic limit. The gyrokinetic simulations display a high degreeof morphological symmetry, to which the PIC simulations converge when β i B g /δB (cid:38) B g /δB (cid:29) Introduction —Magnetic reconnection is a ubiquitousphenomenon that likely plays a significant role in thetransport of energy in the solar corona, solar wind,geospace environment, and magnetic confinement fusiondevices such as ITER . Many of these environmentsare in the strongly magnetized, weak shear limit, whereparticle-in-cell (PIC) simulations become challenging.The gyrokinetic (GK) system of equations has a longhistory in the magnetic confinement fusion community,e.g., , and has been recently applied to space and astro-physical plasmas, e.g., . Fundamentally, GKs assumeslinear frequencies are small compared to the cyclotronfrequency and that the plasma is strongly magnetized.Applying this ordering to the Vlasov-Maxwell system re-sults in the GK equations.Gyrokinetics has been used to model magnetic recon-nection directly, e.g., , and indirectly in simulationsof magnetic confinement plasmas; however, the limita-tions of GKs regarding reconnection have not been thor-oughly explored. Specifically, reconnection events ofteninvolve important physics that may violate the GK order-ings, such as large amplitude, high frequency fluctuationsand instabilities.In this letter, we perform the first detailed compar-ison between 2.5D magnetic reconnection as simulatedby a GK code, AstroGK, and fully electromagnetic PICcodes, NPIC and VPIC, to determine if the full kinetic,PIC, description converges to GKs in the limit of strongguide magnetic field. To address this question, we initial-ize as nearly as possible identical simulations in PIC and a) Electronic mail: [email protected]
AstroGK and examine both a fusion and solar corona rel-evant β i = 8 πn i T i /B g = 0 .
01 and solar wind relevant β i = 1. In both cases, the ratio of the guide to equilib-rium reconnection magnetic field is scanned in PIC andcomparisons are made between the simulation results.Reconnection rates, relative energy transport, and overallmagnitudes are found to be well predicted by GKs, evenin the weaker guide field cases examined. The morphol-ogy of the PIC reconnection results are found to convergeto GKs in the limit of strong guide field, but the conver-gence has a linear β i dependence. The correct predictionof overall magnitudes by GKs implies that the reconnec-tion results scale linearly with guide field strength. Thisis an important result since it implies that a single sim-ulation can represent a range of large guide fields whenproperly scaled. Gyrokinetics —The GK system of equations is basedupon a multi-scale expansion in which small-scale fluc-tuations vary on the linear ( ω ) time-scale while macro-scopic (background) quantities vary on the slow trans-port time-scale . The basic assumptions of this systemare: weak coupling, low frequency compared to the cy-clotron frequency, strong magnetization, and small fluc-tuations. We define the fundamental ordering parameterof GK to be (cid:15) = ω/ Ω ci (cid:28) δF/F ∼ δB/B g ∼ δU s /v ts ∼ (cid:15) ,where δ quantities represent the fluctuating portion, B g is the equilibrium, guide magnetic field, Ω cs = q s B g /m s c is the gyrofrequency, and v ts = (cid:112) T s /m s is the ther-mal speed. When applied to the Vlasov-Maxwell systemof equations, this ordering results in the GK equations.The ordering averages over the fast cyclotron motion ofthe particles and thus removes all cyclotron physics andthe fast magnetosonic wave, reduces the phase space di- a r X i v : . [ phy s i c s . p l a s m - ph ] F e b mensionality from six to five, and describes the evolutionof rings rather than particles. The ordering also enforcescharge neutrality and magnetic moment conservation inthe absence of collisions. Note that this ordering retainsfinite Larmor radius (FLR) effects, non-linear physics,and collisionless wave-particle damping. The details ofthis system of fully non-linear equations can be foundin . As noted above, background quantities such as F are assumed to vary on the slow, transport time-scale,which is one order higher in (cid:15) than typically retained inthe GK system. Therefore, this GK formulation repre-sents a δF system, meaning the background distributionof particles is not evolved. Code Descriptions —The Astrophysical Gyrokineticscode (AstroGK) has been chosen to represent the GKaspect of this study. AstroGK is the slab geometry sib-ling of the toroidal geometry code GS2 , which hasbeen used extensively by the fusion community. As-troGK is an Eulerian continuum code with triply periodicboundary conditions that solves the five-dimensional,electromagnetic gyroaveraged Vlasov-Maxwell equations.The background distribution functions for both speciesare Maxwellian. To obviate the numerical complexityof calculating asymptotically small values, all fluctuat-ing quantities are normalized by (cid:15) to make them O (1),e.g., δB AGK = δB/B g (cid:15) . Note that (cid:15) is neither a fixednor user chosen parameter. Due to the GK ordering andthe normalization employed in AstroGK, a value for (cid:15) must be chosen to make contact with reality and PIC. Anatural choice is (cid:15) = δB /B g ≡ / ˆ B g , where δB is theupstream magnitude of the reconnecting component ofthe magnetic field. For the purposes of this comparison,AstroGK is run in the fully collisionless limit.The first principles electromagnetic relativistic kineticPIC codes NPIC and VPIC have been chosen torepresent the full, unordered kinetic system. NPIC usesan algorithm described in detail by Forslund and is atwo spatial dimension Lagrangian particle code that in-tegrates the full Maxwell-Boltzmann system of equationsin 5D phase space and employs a semi-implicit scheme forimproved efficiency and better noise properties. VPIC isa full, 6D Lagrangian particle code that employs state-of-the art computational techniques to maximize efficiencyand scaling. No explicit collision operator is applied toeither PIC code for the runs contained herein; however,the particle noise and coarse graining inherent to PICsimulations effectively smooths velocity space in a man-ner similar to collisions. Simulation Details —For this study, both systems arerepresented by 2.5D simulations of an electron-ion plasmain the x - y plane, with a guide magnetic field in theout-of-plane direction, ˆ z . A reduced mass ratio of m i /m e = 25 is used, and all simulations initialize T i = T e = T and n i = n e = n . Both systems em-ploy a fully periodic domain in the x - y plane, requir-ing that two current sheets be present in all simula-tions. We here focus on one of the current sheets. For β i = 0 .
01 (1), both simulations use simulation domains L x = L y = 40 π (20 π ) ρ i and initial current sheet halfwidths w = 2 (1) ρ i , where ρ i = v ti / Ω ci is the ion gyro-radius. The parameters employed for the PIC runs areas follows for β i = 0 .
01 (1): ˆ B g = 5 , , ,
50 (1 , , n x = n y = n = 1024 , , , n ppc =1000 (2000 , , ω pe / Ω ce = 0 . n ppc is the number of super particles per cell and ω pe = (cid:112) πn e q e /m e . The AstroGK simulations employ n =1024 and 1280 with velocity space grids ( n λ , n E ) =(32 ,
16) and (16 ,
8) for β i = 0 .
01 and 1 respectively, where n λ and n E are pitch angle and energy grids. To keep β i fixed, the PIC simulations fix the magnitude of B g and reduce δB to achieve the effect of increasing guidefield strength. This approach necessarily implies that thesignal-to-noise ratio (SNR) becomes increasingly poor asˆ B g is increased. Note that β i < m i /m e for β i = 0 . .The initial magnetic geometry employed for this com-parison is a force-free equilibrium of the form B = δB tanh ( y/w )ˆ y + (cid:113) B g + δB cosh − ( y/w )ˆ z . (1)The ions are stationary and the electrons have a net driftvelocity δU e consistent with Ampere’s law. This systemhas recently been studied in detail by Liu et al. .A 1% GEM-like perturbation to B x is added tothe force-free equilibrium to initiate reconnection. SincePIC simulations inherently suffer from particle noise, thefirst twenty non-zero Fourier modes in ˆ x and ˆ y of theAstroGK runs have been populated with uniform ran-dom noise whose spectra and amplitude approximate theinitial noise present in the parallel vector potential ofthe PIC simulations—approximately 0 .
2% of the equilib-rium.
Reconnection Rate —We begin our analysis by exam-ining the reconnection rates of the runs, which are pre-sented in Figure 1. Time is normalized to the in-plane Alfv´en crossing time, τ A = tv Ay /w , where v Ay = δB / √ πn m i . The reconnection rate is fast and agreeswell for all runs, but the rate does show a β i dependence.As noted in , reconnection remains fast in fully kineticsystems, even in the regime of β i < m i /m e , where slowreconnection is observed in two fluid models. To confirmthat the observed reconnection rate is indeed fast for the β i = 0 .
01 case, a scaling study was performed using As-troGK wherein the box size was halved and doubled rel-ative to the simulation presented herein while all otherparameters were held fixed. The peak reconnection rateswere found to be 0 .
10 and 0 .
12 for the half and doublesize domains, which is comparable to the observed rate of0 .
088 for the nominal run. A more detailed study span-ning a wider range of simulation domains and plasmaparameters will be presented elsewhere. The sharp gra-dients and differences in reconnection rate apparent forthe β i = 1 case for τ A (cid:38)
80 are due to the current sheetsthinning and elongating to the point of being unstable to c δ E z / v A y δ B τ A β i . B g β i . B g β i . B g β i . B g β i . c δ E z / v A y δ B τ A β i B g β i B g β i B g β i FIG. 1: Reconnection rate versus time for the, a), β i = 0 .
01 and, b), β i = 1 simulations. Black representsAstroGK and red, green, blue, and cyan representincreasing values of ˆ B g for the PIC simulations.the formation of secondary islands, whose developmentis highly sensitive to noise in the simulation. Thus, mod-erate disagreement in the reconnection rate at these latetimes is not unexpected. Morphology —For the purpose of making detailed com-parisons of the structure and magnitude of relevant phys-ical quantities, we will focus on times τ A = 80 for β i = 0 .
01 and the times prior to secondary island for-mation for β i = 1 respectively—secondary islands format τ A (cid:39) , , , and 73 for ˆ B g = 1 , ,
10 and AstroGKrespectively. These times have been chosen because re-connection is well developed and the morphology is un-affected by secondary islands. Figures 2 and 3 presentthe out of plane current, J z , from PIC and AstroGK atthe indicated times. To decrease the significance of thenoise in the high ˆ B g , β i = 1 runs, the current has beenaveraged over ∆ τ A = 0 . B g = 5 and 10 respec-tively. These figures were produced by converting thePIC results to AstroGK units, which linearly normalizesthe PIC results by the strength of ˆ B g . Immediately clearfrom Figure 2 is that the morphology converges with in-creasing ˆ B g toward the highly symmetric AstroGK resultfor β i = 0 .
01. In contrast, the β i = 1 PIC simulationspresented in Figure 3 display a high degree of morpho-logical convergence even at ˆ B g = 1, although some asym-metry in the current is visible for ˆ B g = 1.The physical explanation for the morphological differ-ences follows from pressure balance across the currentsheet, δP tot + B g δB z / π + δB / π = const . δP tot hasbeen observed to maintain quadrupolar (odd) symmetryacross reconnection layers for nonzero B g , e.g., , whilethe final term on the LHS is manifestly even. Therefore,the balance of these two terms determines the symmetry −505 a) −505 b) −505 c) −505 d) −20 0 20−505 e) y/ ρ i x / ρ i −100102030405060 FIG. 2: Contours of J z normalized to AstroGK units at τ A = 80 for β i = 0 .
01. a)-d) are PIC simulations withˆ B g = 5 , , , and 50 respectively and e) is AstroGK.of δB z , and the density is anti-correlated to δB z alongthe separatrices, which determines the structure of thecurrent layer in the same region. In the limit B g (cid:29) δB ,it is reasonable to assume the ordering δP/P ∼ δB /B g .Therefore, δB oddz δB evenz ∼ β i B g δB . (2)Thus, the β i = 0 .
01 set of PIC simulations exhibit sig-nificant but decreasing asymmetry up to ˆ B g ∼
50, whereequation (2) is O (1), while the β i = 1 simulations ex-hibit symmetry for ˆ B g >
1, and the equations simulatedin AstroGK order equation (2) formally infinite.
Magnitude —As can be seen in Figures 2 and 3, themagnitude of the normalized J z agrees very well withthe results from AstroGK. This agreement in magni-tude extends across all quantities, including magneticfield strengths, outflow velocities, temperatures, andanisotropies. Table I summarizes this finding by examin-ing the ratio of the AstroGK normalized PIC quantitiesto the same AstroGK quantity. The root mean square(RMS) value of each quantity is found over the same re-gions presented in Figures 2 and 3. The temperatureanisotropy, T ⊥ /T (cid:107) , is computed with respect to the totalmagnetic field. To construct this quantity in GKs, thebackground temperature and guide field are assumed tobe 1 /(cid:15) . For low values of ˆ B g , this construction can leadto cases wherein T (cid:107) → β i . Bg
5. All ratios in Table I agree relatively well,aside from the higher ˆ B g runs, whose values are domi-nated by small-scale noise fluctuations. AstroGK doesconsistently underestimate some quantities such as den- −202 a) −202 b) −202 c) −10 −5 0 5 10−202 d) y/ ρ i x / ρ i −0.500.511.5 FIG. 3: Contours of J z normalized to AstroGK units attimes prior to secondary island formation and centeredon the x point for β i = 1. a)-c) are PIC simulationswith ˆ B g = 1 , , and 10 respectively and d) is AstroGK.TABLE I: Ratios of the RMS of each value from PIC tothe RMS AstroGK value over the regions presented inFigures 2 and 3. Run J z v iz v ey B x n e ( T ⊥ /T (cid:107) ) i ( T ⊥ /T (cid:107) ) e β i . B g . . .
90 1 . .
92 0 .
97 0 . β i . B g
10 1 . . .
94 0 .
93 0 .
98 1 . . β i . B g
20 1 . . . .
85 1 . . . β i . B g
50 1 . . . .
95 2 . . . β i B g . . . . . .
84 0 . β i B g . . . . . .
99 1 . β i B g
10 1 . .
98 1 . .
99 1 . . . sity for β i = 1 and v iz ; the reason for this disagreementwill be explored in more detail in a forthcoming paper.There are two significant implications stemming from theoverall good agreement: 1) The linear normalization ofthe PIC results by the strength of ˆ B g implies that quan-tities associated with reconnection scale linearly with ˆ B g in the ˆ B g (cid:29) B g and β i implies that AstroGK accurately predicts such phenom-ena as suprathermal outflow velocities and O (1) densityand magnetic fluctuations, predictions that violate theordering assumptions of GKs. Energy —For the PIC simulations, E = (cid:82) d r /V (cid:0) B / π + (cid:80) s m s n s δU s / P s / (cid:1) is awell conserved quantity, where all bulk velocity, U , is assumed to be fluctuating and the small contributionfrom the electric field has been neglected. However, theconservation of this expression for the energy assumesa collisional system, which the effective collisionalityinherent to PIC simulations provides. Since AstroGKis a continuum code and being run in the collisionlesslimit for this comparison, the above form of the energyis not a well conserved quantity. This is because energyin velocity space cascades toward sharper gradientsthrough linear and non-linear phase space mixing, whichimplies that energy is transferred to progressively highermoments as the system evolves unless collisions smoothvelocity space. The relevant conserved quantity in GKsis the generalized energy W = (cid:90) d r V (cid:34) δB π + (cid:88) s (cid:90) d v T s δf s F s (cid:35) = (cid:90) d r V (cid:34) δB π + (cid:88) s m s n s δU s + (cid:90) d v T s δ ˜ f s F s (cid:35) , (3)where δ ˜ f s is the portion of the perturbed distributionthat does not contribute to the bulk velocity. We willtreat the δ ˜ f s term as equivalent to the perturbed portionof the 3 P s / (cid:15) to compare theevolution of the energy in each system, we ex-amine only the fluctuating energy in PIC by re-moving B g and the initial background P s , δE = (cid:82) d r /V (cid:0) δB / π + (cid:80) s m s n s δU s / δP s / (cid:1) . Thisquantity is well conserved for most of the PIC simula-tions but is conserved only to the 10 −
20% level for allvalues of ˆ B g at β i = 1 and ˆ B g = 50 for β i = 0 .
01. Thepoor conservation is due to poor SNR in these runs andappears primarily in the electron bulk kinetic and ther-mal energies.The evolution of the fluctuating magnetic, bulk kinetic,and thermal energies from their respective initial energiesare plotted in Figure 4 for all of the simulations. All en-ergies are normalized to the total fluctuating energy atthe beginning of each simulation, δE . Aside from thehigh ˆ B g run having poor SNR, the evolution of the mag-netic, bulk kinetic, and ion thermal energies show excel-lent agreement across all simulations for which ˆ B g > β i = 1 is dominated by numerical heating of theentire simulation domain. Note that although the rela-tive energy change matches well across all runs, the freeenergy in each simulation scales with ˆ B − g . Therefore,the large ˆ B g cases have significantly less energy availableto accelerate and heat particles. Conclusions —In the large guide field limit, we haveshown that the fully kinetic particle-in-cell simulationsconverge to gyrokinetic results for 2.5D magnetic recon-nection. Morphological convergence was shown to require β i B g /δB (cid:38) B g /δB (cid:29)
1, implying much strongerguide fields are needed at low β i to achieve convergence,where δB is the reconnecting field. Reconnection rates, ∆ E B / δ E , β i = 0 . τ A ∆ E U / δ E , β i = 0 . τ A ∆ E P / δ E , β i = 0 . τ A ∆ E B / δ E , β i = 1 d) τ A ∆ E U / δ E , β i = 1e) τ A ∆ E P / δ E , β i = 1 f ) τ A FIG. 4: Plots of the change in the, a) and d), magnetic,b) and e), bulk kinetic and, c) and f), thermal energiesfrom their initial values for, a)-c), β i = 0 .
01 and, d)-f), β i = 1. Solid lines indicate ion quantities and dashedare electron. All energies are normalized to the totalinitial energy in the system. Color scheme is the sameas Figure 1.relative energy conversion from magnetic to bulk kineticand thermal energies, and appropriately normalized over-all magnitudes were found to match well for all values of β i and B g /δB >
1. The observed amplitude scalingimplies that gyrokinetics is capable of making accuratepredictions well outside its formal regime of applicabil-ity, and that magnetic reconnection in the large guidefield limit produces quantities that scale linearly withthe guide field strength, which implies that a single sim-ulation can be scaled to represent a range of guide fields.This study validates the use of gyrokinetics as alternativemeans of studying reconnection in the strongly magne-tized kinetic regime. Future work will examine these sim-ulations in greater detail, including their temporal evo-lution, comparison to linear theory, and development ofsecondary instabilities.Acknowledgments—The authors thank James Drakefor helpful discussions. This work was supported bythe US DOE grant DEFG0293ER54197, NSF CAREERAGS-1054061, NASA grant NNH11CC65C, and NASA’sHeliophysics Theory Program. This work used the Ex-treme Science and Engineering Discovery Environment(XSEDE), which is supported by NSF grant numbersPHY090084 and OCI-1053575, with runs performed onKraken at the National Institute for Computational Sci-ences, Stampede at the Texas Advanced Computing Cen-ter. Additional simulations were performed on Pleiadesprovided by NASA’s HEC Program and with Los AlamosInstitutional Computing resources. E. G. Zweibel and M. Yamada, Ann. Rev. Astron. Astrophys. , 291 (2009). M. Yamada, R. Kulsrud, and H. Ji, Reviews of Modern Physics , 603 (2010). H. Ji and W. Daughton, Phys. Plasmas , 111207 (2011),arXiv:1109.0756 [astro-ph.IM]. P. H. Rutherford and E. A. Frieman, Phys. Fluids , 569 (1968). J. B. Taylor and R. J. Hastie, Plasma Phys. , 479 (1968). T. M. Antonsen, Jr. and B. Lane, Phys. Fluids , 1205 (1980). E. A. Frieman and L. Chen, Phys. Fluids , 502 (1982). G. G. Howes, S. C. Cowley, W. Dorland, G. W. Hammett,E. Quataert, and A. A. Schekochihin, Astrophys. J. , 590(2006), astro-ph/0511812. G. G. Howes, J. M. TenBarge, W. Dorland, E. Quataert, A. A.Schekochihin, R. Numata, and T. Tatsuno, Phys. Rev. Lett. , 035004 (2011), arXiv:1104.0877 [astro-ph.SR]. J. M. TenBarge and G. G. Howes, Astrophys. J. Lett. , L27(2013), arXiv:1304.2958 [physics.plasm-ph]. B. N. Rogers, S. Kobayashi, P. Ricci, W. Dorland, J. Drake, andT. Tatsuno, Phys. Plasmas , 092110 (2007). A. Perona, L.-G. Eriksson, and D. Grasso, Phys. Plasmas ,042104 (2010). R. Numata, G. G. Howes, T. Tatsuno, M. Barnes, and W. Dor-land, J. Comp. Phys. , 9347 (2010). R. Numata, W. Dorland, G. G. Howes, N. F. Loureiro, B. N.Rogers, and T. Tatsuno, Phys. Plasmas , 112106 (2011),arXiv:1107.5842 [physics.plasm-ph]. M. J. Pueschel, F. Jenko, D. Told, and J. B¨uchner, Phys. Plas-mas , 112102 (2011). A. Zocco and A. A. Schekochihin, Phys. Plasmas , 102309(2011), arXiv:1104.4622 [physics.plasm-ph]. A. Perona, D. Borgogno, and D. Grasso, Journal of PhysicsConference Series , 012019 (2012). O. Zacharias, R. Kleiber, and R. Hatzky, Journal of PhysicsConference Series , 012026 (2012). N. F. Loureiro, A. A. Schekochihin, and A. Zocco,Phys. Rev. Lett. , 025002 (2013), arXiv:1301.0338[physics.plasm-ph]. M. Kotschenreuther, G. Rewoldt, and W. M. Tang,Comp. Phys. Comm. , 128 (1995). W. Dorland, F. Jenko, M. Kotschenreuther, and B. N. Rogers,Phys. Rev. Lett. , 5579 (2000). W. Daughton, Phys. Plasmas , 3103 (2003). K. J. Bowers, B. J. Albright, L. Yin, B. Bergen, and T. J. T.Kwan, Phys. Plasmas , 055703 (2008). K. J. Bowers, B. J. Albright, L. Yin, W. Daughton, V. Royter-shteyn, B. Bergen, and T. J. T. Kwan, Journal of Physics: Con-ference Series , 012055 (2009). D. W. Forslund, Space Sci. Rev. , 3 (1985). B. N. Rogers, R. E. Denton, J. F. Drake, and M. A. Shay,Phys. Rev. Lett. , 195004 (2001). Y.-H. Liu, W. Daughton, H. Karimabadi, H. Li, and V. Royter-shteyn, Phys. Rev. Lett. , 265004 (2013). Y.-H. Liu, W. Daughton, H. Karimabadi, H. Li, and S. P. Gary,Phys. Plasmas (2014), in press. J. Birn, J. F. Drake, M. A. Shay, B. N. Rogers, R. E. Denton,M. Hesse, M. Kuznetsova, Z. W. Ma, A. Bhattacharjee, A. Otto,and P. L. Pritchett, J. Geophys. Res. , 3715 (2001). R. G. Kleva, J. F. Drake, and F. L. Waelbroeck, Phys. Plasmas , 23 (1995). B. N. Rogers, R. E. Denton, and J. F. Drake, J. Geophys. Res. , 1111 (2003). A. A. Schekochihin, S. C. Cowley, W. Dorland, G. W. Hammett,G. G. Howes, E. Quataert, and T. Tatsuno, Astrophys. J. Supp.182