An extended MHD study of the 16 October 2015 MMS diffusion region crossing
J. M. TenBarge, J. Ng, J. Juno, L. Wang, A. H. Hakim, A. Bhattacharjee
mmanuscript submitted to
JGR-Space Physics
An extended MHD study of the 16 October 2015 MMSdiffusion region crossing
J. M. TenBarge , , J. Ng , , , J. Juno , L. Wang , , A. H. Hakim , , and A.Bhattacharjee , , Department of Astrophysical Sciences, Princeton University, Princeton, NJ, USA. Princeton Center for Heliophysics, Princeton Plasma Physics Laboratory, Princeton, NJ, USA. Department of Astronomy, University of Maryland, College Park, MD, USA. IREAP, University of Maryland, College Park, MD, USA. Princeton Plasma Physics Laboratory, Princeton, NJ, USA.
Key Points: • An extended MHD model evolving the full electron pressure tensor and captur-ing aspects of Landau damping is applied to the MMS Burch event • Many aspects of the fluid model agree well with in situ observations, includingstrong parallel electron heating near the current layer • The event study demonstrates the efficacy of the ten moment model for globalmagnetosphere modeling efforts
Corresponding author: J. M. TenBarge, [email protected] –1– a r X i v : . [ phy s i c s . s p ace - ph ] M a r anuscript submitted to JGR-Space Physics
Abstract
The Magnetospheric Multiscale (MMS) mission has given us unprecedented accessto high cadence particle and field data of magnetic reconnection at Earth’s magne-topause. MMS first passed very near an X-line on 16 October 2015, the Burch event,and has since observed multiple X-line crossings. Subsequent 3D particle-in-cell (PIC)modeling efforts of and comparison with the Burch event have revealed a host ofnovel physical insights concerning magnetic reconnection, turbulence induced particlemixing, and secondary instabilities. In this study, we employ the
Gkeyll simulationframework to study the Burch event with different classes of extended, multi-fluidmagnetohydrodynamics (MHD), including models that incorporate important kineticeffects, such as the electron pressure tensor, with physics-based closure relations de-signed to capture linear Landau damping. Such fluid modeling approaches are ableto capture different levels of kinetic physics in global simulations and are generallyless costly than fully kinetic PIC. We focus on the additional physics one can capturewith increasing levels of fluid closure refinement via comparison with MMS data andexisting PIC simulations.
The coupled sun-Earth system is a critically important problem for protectingmankind’s technological infrastructure from damaging space weather events, but plan-etary magnetospheres represent a multi-scale, multi-process physical system that isweakly collisional and coupled to the highly dynamic sun, placing great demand onmodeling efforts. Modern in situ spacecraft missions like the Magnetospheric Multi-scale (MMS) mission (Burch, Moore, et al., 2016) are providing field and particle dataat unprecedented temporal and spatial resolution, revolutionizing our understandingof the fundamental plasma processes occurring in Earth’s magnetosphere, such asmagnetic reconnection, shocks, and turbulence. Magnetic reconnection of the shockedplasma at Earth’s magnetosheath transports energy from the solar wind plasma intoEarth’s magnetosphere (Paschmann et al., 1979) and is a core process for understand-ing the sun-Earth coupling.Despite the impressive quality of in situ measurements, numerical modelling re-mains our best resource for understanding the global behavior of planetary magne-tospheres. Today, the three main numerical approaches are: 1) fully kinetic modelssuch as particle-in-cell (PIC) and continuum Vlasov; 2) hybrid models wherein onespecies, typically ions, is treated kinetically and the other as a fluid; and 3) fluidmodels, such as magnetohydrodynamics (MHD), Hall MHD, and various extensionsof MHD. Each of these modeling approaches has inherent advantages. For example,fully kinetic models are capable of capturing all of the fundamental plasma physicsprocesses occurring in a magnetosphere; however, such complete modeling comes atthe expense of computational complexity, making global simulations of Earth’s mag-netosphere too costly to perform. Hybrid approaches are attractive because of theirreduced computational cost relative to fully kinetic models, but the electron are typi-cally treated as a massless, isothermal fluid, which is problematic for sub-ion physicssuch as reconnection, which will be diffused by explicit or numerical resistivity at thesimulation grid scale. Additionally, one may be given to assume that since the ions arekinetic, they are treated fully; however, the kinetic ions are coupled to fluid electronsthrough Maxwell’s equations, which can have consequences at large-scales, such asfast magnetosonic modes being undamped (Parashar et al., 2014; Groˇselj et al., 2017).Fluid simulations are the most commonly employed approach to model global magne-tospheres due to their significantly reduced computational demand, but they assumethat the plasma is highly collisional and thus near thermodynamic equilibrium; neitheris a good assumption in the weakly collisional plasmas that permeate the solar windand magnetospheres. Fluid models can be improved by retaining progressively higher –2–anuscript submitted to
JGR-Space Physics fluid moments, with each additional moment retained by the closure correspondingto the fluid model more accurately representing the fully kinetic system. The closureemployed by the fluid model can also be tuned to better approximate kinetic effects.Magnetic reconnection plays a fundamental role in coupling the sun-Earth sys-tem; therefore, any fluid model employed in global magnetosphere modeling shouldminimally capture the physics of reconnection accurately. Reconnection is ultimatelycontrolled by electron demagnetization, which is governed by a combination of electroninertia and pressure agyrotropy (Vasyliunas, 1975). Based upon this condition, theminimum requirement for a fluid model for global modeling is that the electron closureretain finite mass and evolve the full pressure tensor, implying at least ten momentsbe retained for each species, n , u , and P ij , i.e., 1 + 3 + 6 moments. Since the plasmasof interest are weakly collisional, the closure for a ten moment model should incorpo-rate collisionless effects, such as Landau damping. The Hammett-Perkins (Hammett& Perkins, 1990) approach to closure is one such attempt to incorporate aspects oflinear Landau damping in a fluid model. The Hammett-Perkins closure is based upona linearization of the Vlasov equation, followed by integration over velocity space toobtain an expression for the unclosed fluid moment.The Gkeyll simulation framework includes both a two-fluid, five-moment modelthat achieves closure via truncation of moments beyond the isotropic pressure, p ,and a ten moment model with an approximation of the Hammett-Perkins closure forthe heat flux. The ten moment model has been well compared to PIC in the anti-parallel reconnection case (Wang et al., 2015) and island coalescence (Ng et al., 2015). Gkeyll has also been successfully employed to model the magnetosphere of the Jovianmoon Ganymede (Wang et al., 2018). To further evaluate the efficacy of the tenmoment model, we consider here an asymmetric reconnection event based upon the 16October 2015 MMS magnetosheath crossing (Burch, Torbert, et al., 2016). The largeasymmetry of this event provides a stringent test for the ten moment model in the twoand three dimensions that can be validated against in situ observations and existingPIC simulations of the event.The paper is organized as follows. In section 2, we introduce the
Gkeyll frame-work and discuss the closures for the two fluid five ( n , u , and p ) and ten moment( n , u , and P ij ) models that are part of Gkeyll . Section 3 presents the initial con-ditions derived from the Burch event and the closure parameters chosen for the tenmoment model. Simulation results for the ten and five moment models are providedin section 4. Section 5 places the results of the paper in the context of contemporarysimulation methods. Finally, in section 6, we present a summary of this paper.
The simulation framework we will use for this study is
Gkeyll , which includessolvers for continuum Vlasov-Maxwell (Juno et al., 2018), gyrokinetics (Shi et al.,2015), and five and ten moment two fluid equations (Hakim et al., 2006; Hakim, 2008;Wang et al., 2015). In this work, we focus on the five- and ten-moment two-fluidmodels.The fluid models evolve the electromagnetic fields through Maxwell’s equations ∇ × E = − ∂ B ∂t , (1) ∇ × B = µ J + 1 µ (cid:15) ∂ E ∂t , (2)where E and B are the electric and magnetic field, µ and (cid:15) are the permeability andpermittivity of free space, and J = (cid:80) s n s q s u s , with q s , n s , and u s representing the –3–anuscript submitted to JGR-Space Physics species charge, density, and velocity. In the following, we drop the species subscriptsfor notational elegance.The evolution equations for the moments are obtained via the standard procedureof taking moments of the Vlasov equation: ∂n∂t + ∂nu i ∂x i = 0 , (3) m ∂nu i ∂t + ∂ P ij ∂x j = nq ( E i + (cid:15) ijk u j B k ) , (4) ∂ P ij ∂t + ∂ Q ijk ∂x k = nqu [ i E j ] + qm (cid:15) [ ikl P kj ] B l , (5)where (cid:15) ijk is is the completely anti-symmetric Levi-Cevita tensor, which is definedto be ± P ij = m (cid:90) v i v j f d v = P ij + mnu i u j , (6)is the pressure tensor and Q ijk = m (cid:90) v i v j v k f d v = Q ijk + u [ i P jk ] − mnu i u j u k (7)is the heat-flux tensor, where P ij and Q ijk are the rest-frame pressure and heat-fluxtensors. At this level, the fluid equations are exact but an endless hierarchy, with the n th moment requiring knowledge of the n + 1st moment. One approach to achieve closure is a simple truncation by neglecting the heat-fluxand considering only isotropic pressures. We begin by taking the trace of equation (5)to obtain the exact energy equation ∂ E ∂t + 12 ∂ Q iik ∂x k = nq u · E , (8)where E = 12 P ii = 32 p + 12 mn u , (9) p = P ii /
3, 12 Q iik = q k + u k ( p + E ) + u i π ik , (10) q k = Q iik / π ij = P ij − pδ ij is the viscous stress tensor, and δ ij is the Kronecker delta function. Closure is obtained by setting q k = 0 and π ij = 0,reducing equation (8) to ∂ E ∂t + ∂∂x k [ u k ( p + E )] = nq u · E . (11)Thus, equations (3)-(4) and (11) constitute a closed system of five-moment equations, n , u , and p , equations for each species. Note that this model reduces formally to HallMHD in the limit m e → (cid:15) → –4–anuscript submitted to JGR-Space Physics
To retain additional kinetic physics relative to the five-moment model, we desirethe full pressure tensor, thus constituting a ten moment model. Therefore, we mustobtain an expression for ∂Q ijk /∂x k appearing implicitly in equation (5) to close thesystem. Additionally, since space plasmas tend to be weakly collisional, we seek acollisionless closure and thus choose to use a Hammett-Perkins (Hammett & Perkins,1990) inspired closure. The closure is based upon a linearization of the Vlasov equation,followed by integration over velocity space to obtain an expression for the unclosedfluid moment, which in this case is the heat flux. This approach leads to expressionsinvolving plasma response functions in Fourier space, which are approximated by an n -pole Pad´e approximant. Following Hammett and Perkins (1990), we choose to employa 3-pole approximation to the response function; however, since the approximationis in Fourier space, the closure is non-local along magnetic field lines. To reducethe computation demand of integrating along field lines and Fourier transforming, wefurther use a local approximation of the heat-flux divergence ∂Q ijk ∂x k = v t | k | ( P ij − pδ ij ) , (12)where v t = (cid:112) T /m is the thermal speed, and k is a typical wavenumber, which definesa scale below which collisionless damping is thought to occur. Equation (12) togetherwith equations (3)-(5) represent the closed ten moment model for n , u , and P ij . Notethat this closure closely resembles a collisional closure that scatters the system towardsisotropy with collision frequency ν = v t | k | . k = 0 would allow the pressure to evolvefreely, while the limit k → ∞ isotropizes the pressure tensor at all scales, approachingthe five-moment limit. The initial conditions for the study contained herein are motivated by the firstobserved Magnetospheric Multiscale (MMS) spacecraft crossing of a magnetic field X-line at Earth’s magnetosheath (Burch, Torbert, et al., 2016), hereafter referred to asthe Burch event. The initial magnetic field and temperature profiles have the followingform, f ( Q s , Q m ) = 0 . Q s + Q m + ( Q m − Q s ) tanh ( y/w )] , (13)where subscripts s and m correspond to magnetosheath and magnetosphere quantities, w = d is is the initial current sheet width, d is = c/ω pis is the asymptotic ion inertiallength on the sheath side, and ω pis is the sheath ion plasma frequency. For notationalsimplicity, we will use d i ≡ d is and Ω ci = Ω cis going forward, where Ω cis is the sheathgyrofrequency. The initial magnetic field profile is given by B x = f ( B sx , B mx ), where B mx /B sx = − .
7, with a uniform guide field B z /B sx = 0 . T em /T es = 12 . T im /T is = 5 .
63. The density profile is chosen to satisfy force balance across thesheet, which leads to n m /n s = 16 .
7. An initial vector potential perturbation of theform A z = ψ cos (2 πx/L x ) [1 − cos (4 πy/L y )] (14)is added to initiate reconnection, where ψ = 0 .
1. The simulation dimensions are( L x , L y ) = (40 . , . d i in 2D and ( L x , L y , L z ) = (40 . , . , . d i in 3D,with fully periodic boundary conditions in all cases. The number of grid points are( N x , N y ) = (2048 , N x , N y , N z ) = (1024 , , x = 0 . d es and∆ x = 0 . d es in 2 and 3D. A reduced mass ratio of m i /m e = 100 and c/v A ⊥ s = 500are used in all cases, where v A ⊥ s = B sx / √ µ m i n s is the asymptotic sheath Alfv´enspeed based upon the reconnecting field. To break the inherent symmetry of the initialconditions in 3D, noise with root-mean-square amplitude B noise = 0 . B sx is added –5–anuscript submitted to JGR-Space Physics
Figure 1. (Left) The reconnection rate for runs 2D10 (black) and 3D10 (red dashed). (Right)The relative change in energy for runs 2D10 (solid), 3D10 (dashed), and 3D5 (dash-dotted). E isthe total initial energy, E = E tot ( t = 0), and ∆ E ∗ = E ∗ ( t ) − E ∗ ( t = 0). to the first 20 Fourier modes in the x , y , and z directions of B x and B y . We will referto the 2D ten moment and 3D five and ten moment Gkeyll simulations as 2D10, 3D5,and 3D10, respectively.Prior
Gkeyll reconnection studies (Ng et al., 2015; Wang et al., 2015) havedemonstrated that k s d s ∼ s refers to species. However, a value of 1 /
10 providesslightly better performance for turbulence simulations without significantly degradingthe reconnection performance (Juno et al., 2018). Since we expect these event sim-ulations to be turbulent, we choose closure parameters for both ten moment modelsimulations to be k s d s = 1 /
10. Formally, this choice for k will permit pressureanisotropy and agyrotropy to develop to a scale l s = 2 π/k s (cid:39) d s before the tenmoment closure scatters the pressure toward isotropy. We first present the results of the Burch event study using the ten-moment modelin two and three dimensions. In figure 1 is presented the reconnection rates (leftpanel) and relative changes of energy components (right) for runs 2D10 (solid) and3D10 (dashed). The reconnection rate in the 3D case was computed by averaging B x and B y over the z -direction and integrating to construct an average vector potential, (cid:104) A z (cid:105) . In both 2 and 3D, the rate reaches a plateau at the canonically observed valueof cE/v A ⊥ s B xs ∼ .
1. The relative change of the energies is essentially identical in 2Dand 3D for the ten-moment model. Thus, in terms of energy exchange and reconnectionrates, the 2 and 3D ten-moment simulations are fundamentally similar.The evolution of the out-of-plane current density, J z , in the xy -plane for runs2D10 and 3D10 are presented in figures 2 and 3. For run 3D10, the current density istaken at z = 0. The early-time evolution for t Ω ci = 10 to 30 in 5 / Ω ci increments ispresented in the left panel of each figure. Little qualitative or quantitative differenceexists during this phase, the most significant difference being a moderately increasedcurrent density in 2D due to the existence of sharper magnetic field gradients. Theright panels of the figures presents the current density at times t Ω ci = 35 to 55.Beginning at t Ω ci = 35, the 2 and 3D evolution noticeably diverges. The source ofthe divergence becomes clear by examining the full 3D evolution of the current in run –6–anuscript submitted to JGR-Space Physics
Figure 2.
Slices of the out-of-plane current density, J z , for run 2D10 at times t Ω ci = 10 to 30(left) and 35 to 55 (right), with each frame separated by 5 / Ω ci . Figure 3.
Slices of the out-of-plane current density, J z , taken at z = 0 for run 3D10 at times t Ω ci = 10 to 30 (left) and 35 to 55 (right), with each frame separated by 5 / Ω ci . t Ω ci = 30 the current layer is predominantly laminar,while at t Ω ci = 50, the entire layer has become turbulent. In figure 5 is shown thecurrent in the yz -plane at the X-line, x = 30 d i , at the same times as figure 3. It isclear that in three dimensions turbulence plays an active role, broadening the currentlayers, thereby further reducing the local current density. The entire current layer isalso kink unstable with wavelength equal to the longest wavelength accessible in thedomain.The turbulence is generated by a field aligned shear flow instability with wave-length kd e (cid:39) t Ω ci = 30 to 50. The instability drives a low amplitude density per-turbation along the current layer that spreads into the high density, magnetosheath,side. By using the large temperature gradient across the layer as a tracer in the rightpanel of figure 6, it becomes apparent that the instability drives Kelvin-Helmholtzlike vorticies along the current layer that lead to enhanced mixing relative to the 2Dsimulation. This instability is active along the entire separatrix layer; however, it ismost significant in the vicinity of the X-line.The electron shear flow instability present in run 3D10 has been observed inPIC simulations of the Burch event (Le et al., 2017); however, the shear instability is –7–anuscript submitted to JGR-Space Physics (a) (b)(c) (d) x/d i
The full 3D evolution of the current density J z in run 3D10. The upper two panels,(a) and (b), correspond to t Ω ci = 30, and the lower two correspond t Ω ci = 50. preceded by the growth of the lower hybrid drift instability (LHDI) (Krall & Liewer,1971; Daughton, 2003), which feeds off of the large density gradient present in thisevent, but the LHDI is absent from this ten-moment simulation. Prior particle-in-cellsimulations of the Burch event (Price et al., 2016; Le et al., 2017) find that the LHDIplays a fundamental role in generating turbulence that both supports the reconnectionelectric field and enhances mixing and electron heating. Linear analysis of the ten-moment model performed in Ng, Hakim, Juno, and Bhattacharjee (2019) shows thatthe LHDI is indeed captured by the model, but the closure parameter for both species, k s , can significantly modify the growth rate of the instability. To recover kineticgrowth rates, the ion closure parameter must be approximately equal to the wavelengthof the fastest growing LHDI mode. Unfortunately, the fastest growing LHDI modeshave kρ e ∼ k i d i (cid:46)
1, andthe value chosen for this study. Therefore, the ten-moment closure for the parametersemployed in run 3D10 damps the LHDI, leading to a growth rate that is reduced byapproximately an order of magnitude relative to kinetic theory.The electron temperature anisotropy, T (cid:107) /T ⊥ , for runs 2D10 (left) and 3D10(right) are presented in figure 7 for times t Ω ci = 30 to 50. Significant electron temper-ature anisotropy develops in both two and three dimensions in the ten-moment model;however, the presence of shear driven turbulence and mixing in three dimensions en-hances the temperature anisotropy along the separatrices and significantly broadensthe region with T ⊥ > T (cid:107) in the vicinity of the X-line. Examining the structure ofthe temperature anisotropy in the yz -plane in the left panel of figure 8 for run 3D10 –8–anuscript submitted to JGR-Space Physics
Figure 5.
Slices of the current density J z taken at at the X-line, x = 30 d i , for run 3D10 attimes t Ω ci = 10 to 30 (left) and 35 to 55 (right), with each frame separated by 5 / Ω ci . Figure 6.
Slices of the electron density (left) and electron temperature (right) taken at theX-line, x = 30 d i , for run 3D10 at times t Ω ci = 30 to 50, with each frame separated by 5 / Ω ci . reveals that the most significant excess parallel heating occurs in regions where thecurrent layer has oscillated toward the magnetosphere, creating local depressions ofthe magnetic field amplitude.The largest value of the temperature anisotropy occurs in run 3D10 and is T (cid:107) /T ⊥ (cid:39) .
5, which is in close agreement with the value observed by MMS [see fig-ure 3I of Burch, Torbert, et al. (2016)]. However, the Burch event observation of T (cid:107) /T ⊥ > T (cid:107) /T ⊥ > T (cid:107) /T ⊥ < T (cid:107) (a & c) and T ⊥ (b & d) separately for runs 2D10 (a & b) and3D10 (c & d) at t Ω ci = 40 and z = 0. Overall, the electron heating is isotropic,with regions of anisotropic heating being due to isolated and small variations in thetemperature components. Three dimensional PIC simulations of the Burch event findthat the entire magnetosphere separatrix layer is heated in only the parallel directionwith T (cid:107) /T ⊥ ∼ − –9–anuscript submitted to JGR-Space Physics
Figure 7.
Slices of the electron temperature anisotropy, T (cid:107) /T ⊥ for run 2D10 (left) and 3D10(right), taken at z = 0, times t Ω ci = 30 to 50, with each frame separated by 5 / Ω ci . Figure 8. (Left) Slices of the electron temperature anisotropy, T (cid:107) /T ⊥ , for run 3D10 taken atthe X-line, y = 30 d i at times t Ω ci = 30 to 50, with each frame separated by 5 / Ω ci . (Right) Slicesof the the parallel (a & c) and perpendicular (b & d) electron temperatures for runs 2D10 (a &b) and 3D10 (c & d) at time t Ω ci = 40 and z = 0 for run 3D10. Since the ten-moment model evolves the full pressure tensor, we present in figure 9the electron pressure agyrotropy √ Q for runs 2D10 (left) and 3D10 (right), where Q = P + P + P P ⊥ + 2 P ⊥ P (cid:107) (15)is the measure of gyrotropy defined in Swisdak (2016). By definition, 0 ≤ Q ≤ Q = 0(1) corresponds to gyrotropy (maximal agyrotropy). The peak valuesof √ Q (cid:39) . .
55 occur in the vicinity of the X-point and X-Line in 2 and 3D,respectively, with the reduced peak value in run 3D10 due to turbulence broadeningthe current layers. The gyrotropy measure also highlights the separatricies, where √ Q (cid:39) .
05. The structure and magnitude of √ Q are generally consistent with resultsobserved in PIC simulations and in situ spacecraft observations of reconnection, e.g.,Swisdak (2016); Zenitani and Nagai (2016); Bourdin (2017); Genestreti et al. (2018).To further emphasize the importance of the pressure agyrotropy at the X-line,we next examine the components of the generalized Ohm’s law for the ten-momentmodel, which can be derived from the electron momentum equation, equation (4). At –10–anuscript submitted to JGR-Space Physics
Figure 9.
Slices of the electron pressure agyrotropy, √ Q for run 2D10 (left) and 3D10 (right),taken at z = 0, times t Ω ci = 30 to 50, with each frame separated by 5 / Ω ci . Figure 10.
The 1D components of the generalized Ohm’s law across the current layer for runs3D10 at time t Ω ci = 40 (left) and 3D5 at t Ω ci = 20 (right), both slices taken at z = 0. Thevertical dashed line in each figure corresponds to location of the magnetic field reversal. the X-line, only the guide field, z , aligned component of the Ohm’s law contributes E z = − ( u i × B ) z + 1 n e | e | ( j × B ) z − n e | e | (cid:18) ∂P exz ∂x + ∂P eyz ∂y + ∂P ezz ∂z (cid:19) − n e | e | (cid:18) ∂ρ e u ez ∂t + ∇ · ( ρ e u ez u e ) (cid:19) , (16)where e is the charge of an electron and ρ e = m e n e is the electron mass density. Theone dimensional components of the generalized Ohm’s law across the current layerfor run 3D10 at time t Ω ci = 40 and z = 0 are plotted in the left panel of figure 10,where the vertical dashed line indicates the location of the magnetic field reversal. Theonly component of the Ohm’s law that contributes significantly at the reversal is theoff-diagonal pressure tensor (dot dashed blue). Away from the X-line, the Hall term, j × B , (dashed green) is the primary contributor. Moving to the five-moment model, we begin by examining the relative change inenergy for run 3D5 presented in the right panel of figure 1 as dash-dotted lines. Theevolution of run 3D5 is notably different from the ten-moment simulations, with themagnetic field decaying rapidly, and its free energy going to ion and electron thermal –11–anuscript submitted to
JGR-Space Physics
Figure 11.
Slices of the out-of-plane current density, J z , taken at z = 0 for run 3D5 at times t Ω ci = 5 to 20 (left) and 25 to 40 (right), with each frame separated by 5 / Ω ci . energy. Also notable is the poor conservation of total energy (magenta) in this runthat onsets shortly following the initial decay of magnetic energy. Poor conservationof energy indicates that a significant amount of energy is reaching the grid scale, evenat this early stage of the simulation. The reconnection rate for run 3D5 is not plottedin this figure for reasons that will become apparent below.In figure 11 is plotted the evolution of the out-of-plane current density taken at z = 0 for run 3D5 at times t Ω ci = 5 to 20 (left) and 25 to 40 (right), with each frameseparated by 5 / Ω ci . The full three dimensional structure of the layer at times t Ω ci = 10and 30 are presented in figure 12. Clearly, the behavior of the five-moment model forthese highly asymmetric parameters differs significantly from that of the ten-momentmodel. The current layer is permeated by turbulence, which leads to substantialdiffusion and broadening across the entire layer, essentially destroying the currentlayer. By t Ω ci = 20, the entire domain has become turbulent. The strong turbulencepresent in this simulation is responsible for the transport of energy to the grid scalethat results in poor energy conservation. Additionally, the turbulence precludes onefrom reliably calculating the reconnection rate using conventional methods, since themagnetic X and O points are ambiguous.The origin of the turbulence in this simulation is the LHDI, in stark contrastto the ten-moment results. Evidence for this can be seen in figure 13, where theevolution of the electron density (left) and E z (right) are plotted in the yz -plane attimes t Ω ci = 0 to 20. Already, at t Ω ci = 5, evidence of the LHDI has become apparent,with density and electric field fluctuations extending from the the edge of the centralcurrent layer into the low density, magnetosphere, side with wavelength kρ e ∼
1. By t Ω ci = 10, the initially electrostatic modes have fully penetrated the current layer andfill the entire magnetosphere side.These results concerning the LHDI are also supported by linear stability analysisof the five-moment model presented in Ng et al. (2019). Ng et al. (2019) demonstratesthat the LHDI grows by a factor of two to three faster in the five-moment model thanin kinetic theory, thus leading to the explosive growth of the instability and subsequentturbulence seen in run 3D5.The generalized Ohm’s for the five-moment model differs from that for ten-moment model given in equation (16) only in that the the off-diagonal pressure tensorcomponents do not contribute. The one dimensional components of the generalizedOhm’s law across the current layer for run 3D5 at time t Ω ci = 20 and z = 0 are plotted –12–anuscript submitted to JGR-Space Physics (a) (b)(c) (d) x/d i
The full 3D evolution of the current density J z in run 3D5. The upper two panels,(a) and (b), correspond to t Ω ci = 10, and the lower two correspond t Ω ci = 30. in the right panel of figure 10, where the vertical dashed line indicates the locationof the magnetic field reversal. The primary contributor to the reconnection electricfield at the X-line is the Hall term (dashed green). Away from the X-line, the electronpressure gradient (lavender) is the primary contributor. The ten-moment model presented in section 2.2 represents the most sophisticatedfluid closure available that is numerically tractable for global magnetosphere model-ing. The closure employed by the ten-moment model captures many aspects of kinetictheory that are well beyond the capabilities of the most commonly applied model toglobal magnetospheres, MHD. MHD assumes the plasma is highly collisional, suchthat it is in local thermodynamic equilibrium; therefore, the only mechanism avail-able to break field lines and allow magnetic reconnection to proceed is resistivity, beit grid-scale diffusion or explicitly imposed. In either case, the resistivity employedexceeds the real value in weakly collisional planetary magnetospheres by orders ofmagnitude, thus altering significantly the location, evolution, and transport inducedby reconnection. Indeed, MHD models predict either extended Sweet-Parker layersor copious plasmoid production, with the evolution being determined entirely by theLundquist number rather than kinetic processes (Daughton & Roytershteyn, 2012).Additionally, although the ten-moment model has difficulty capturing the LHDI, theLHDI is completely ordered out of MHD. –13–anuscript submitted to
JGR-Space Physics
Figure 13.
Slices of the electron density (left) and E z (right) taken at the X-line, x = 30 d i ,for run 3D5 at times t Ω ci = 0 to 20. Fluid models more sophisticated than MHD exist, such as Hall MHD or two-fluidmodels with isothermal closures, but they are rarely employed in the context of globalmodelling due to computational constraints. However, both Hall MHD and traditionaltwo-fluid models also neglect the electron kinetic effects that permit reconnection toproceed naturally, relying on resistivity to break field lines. As demonstrated by thefive-moment model in section 4.2, retaining finite electron inertia effects without fullyevolving the pressure tensor can lead to anomalously large growth rates for instabilitiesthat are common to global magnetospheres, such as the LHDI.As computational resources are expanding, hybrid approaches are gaining in-creased focus, because they treat the ions kinetically; however, hybrid models suffermany of the same issues outlined above. Since the electrons are treated as a fluid, theyalso require a closure model, which is that of a massless, isothermal fluid, althoughfinite electron mass is occasionally retained. Therefore, reconnection still proceeds viaresistivity, despite having fully kinetic ions. The electron closure also feeds back onthe ions through Maxwell’s equations, leading to consequences at large-scales, suchas undamped fast magnetosonic modes (Groˇselj et al., 2017) and electron equation ofstate dependent energy transport (Parashar et al., 2014). Further, electron instabil-ities such as the LHDI remain problematic in hybrid models, because they occur onelectron kinetic scales and often involve electron phase space resonances.Despite the overall success of the application of the ten-moment model to theBurch event presented in Section 4.1, it does display some weaknesses relative tofully kinetic, particle-in-cell simulations. Specifically, the ten-moment closure fails toproperly capture the LHDI and the important role it plays in heating and transportingplasma across the separatrix. The source of the disagreement follows directly from thelocal closure chosen in section 2.2. Ng et al. (2019) demonstrates that a non-local heatflux closure is able to accurately recreate the LHDI, where the local heat flux closurepresented in equation (12) is replaced withˆ q ijk = − i v t | k | χk [ i ˆ T jk ] , (17)where Q ijk = n ˜ q ijk , ˆ ∗ represents the Fourier transform, and χ is an adjustable con-stant. The | k | factor appearing in the denominator makes this closure non-local inreal space, drastically increasing the computational expense and making its utilityfor global modelling limited. A gradient-based closure to the ten-moment model hasbeen introduced by Allmann-Rahn, Trost, and Grauer (2018), where equation (12) is –14–anuscript submitted to JGR-Space Physics replaced with ∂Q ijk ∂x k = − v t | k | ∇ ( P ij − pδ ij ) , (18)and k is an adjustable constant. Allmann-Rahn et al. (2018) find that the gradientclosure better agrees with Vlasov and PIC simulations than the local closure presentedin equation (12), with minimal additional computational expense for typical choicesof k . Therefore, the gradient closure appears to be a good alternative to the localclosure for global modelling, but further study to determine how well the gradientclosure captures the LHDI is necessary. We have performed an event study of the Burch, Torbert, et al. (2016) 16 Oc-tober 2015 MMS diffusion region crossing using the two-fluid, multi-moment modelscontained in the
Gkeyll framework. The ten-moment model evolves the full pressuretensor and is closed with a local Hammett-Perkins model for the heat flux divergence.The closure is designed to capture linear Landau damping that sets in at a fixed,isotropic, and species dependent wavenumber, k s . The five-moment model achievesclosure by truncating the moment hierarchy at the isotropic pressure for each species,setting all higher moments to zero. The ten-moment model was applied to the Burchevent in two and three dimensions, while the five-moment model was applied only inthree dimensions.The ten-moment model was found to accurately capture aspects of the Burchevent that are beyond the capabilities of conventional fluid approaches, such as MHD,Hall MHD, or two-fluid models with isothermal closures. Inclusion of the full electronpressure tensor, electron inertia, and Hall terms permits the electrons to demagnetizenaturally at the X-line. The pressure agyrotropy structure and amplitude developedby the ten-moment model is consistent with existing PIC simulations of reconnection.The ten-moment model also accurately reproduces the same magnitude of electrontemperature anisotropy as observed in situ with MMS. However, the spatial locationof the temperature anisotropy is on the magnetosheath side of the X-line; whereas, insitu observations (Burch, Torbert, et al., 2016) and PIC simulations (Le et al., 2017)find it to be on the low density, magnetosphere side. The discrepancy is due to lowerhybrid drift instability (LHDI) induced heating and transport in the PIC simulation.The chosen closure parameter for the ion ten-moment model, k i d i = 1 /
10, is toofar in scale from the fastest growing LHDI mode and is thus anomalously damped.Choosing a larger value for the closure parameter, k i d e ∼
1, can reproduce the growthof the LHDI accurately but would prevent reconnection from proceeding correctly,since the electron pressure would be scattered toward isotropy at a much faster rate,precluding electron agyrotropy from forming. Electron shear flow and kink instabilitiesdo develop in the layer in three dimensional ten-moment
Gkeyll simulation, which arealso observed in PIC simulations (Le et al., 2017). The shear flow instability generatessignificant turbulence and increases the transport across the current layer relative tothe two dimensional
Gkeyll simulations of the Burch event.In contrast to the ten-moment model, the LHDI grows more rapidly in the five-moment model than predicted by kinetic theory. The explosive growth of the instabilityand concomitant turbulence rapidly destroy the current layer and drives energy to thegrid scale, where it is diffused away and lost from the simulation. Due to the explosivegrowth of the LHDI, the five-moment model was found to poorly model the Burchevent.The ten-moment model in
Gkeyll is a simple but effective closure for capturingreconnection dynamics, which are fundamental for understanding the coupled sun-Earth system and other planetary magnetospheres, both within and without the solar –15–anuscript submitted to
JGR-Space Physics system. The model well captures many of the kinetic aspects of the Burch event,demonstrating the efficacy of the ten-moment model for performing global magneto-sphere simulations. However, the local, isotropic closure does anomalously damp smallscale modes that play an important role in mixing and heating the plasma at the mag-netopause. Therefore, additional improvements in the closure are possible, such asthe gradient based closure introduced by Allmann-Rahn et al. (2018) and discussed insection 5. These improvements are necessary to attain maximal fidelity of global fluidmodels of planetary magnetospheres and will be explored in future works.
Acknowledgments
The authors are grateful for fruitful discussions with Marc Swisdak, Ari Le, MikeShay, and Paul Cassak. This research was supported by NSF grants AGS-1338944and AGS-1622306. This research also used resources of the National Energy ResearchScientific Computing Center (NERSC), a U.S. Department of Energy Office of ScienceUser Facility operated under Contract No. DE-AC02-05CH11231. Simulation dataused in this manuscript is available upon request.
References
Allmann-Rahn, F., Trost, T., & Grauer, R. (2018). Temperature gradient drivenheat flux closure in fluid simulations of collisionless reconnection.
J. PlasmaPhys. , (3).Bourdin, P.-A. (2017, September). Catalog of fine-structured electron veloc-ity distribution functions - Part 1: Antiparallel magnetic-field reconnection(Geospace Environmental Modeling case). Ann. Geophys. , , 1051-1067. doi:10.5194/angeo-35-1051-2017Burch, J. L., Moore, T. E., Torbert, R. B., & Giles, B. L. (2016, March). Magne-tospheric Multiscale Overview and Science Objectives. Space Sci. Rev. , , 5-21. doi: 10.1007/s11214-015-0164-9Burch, J. L., Torbert, R. B., Phan, T. D., Chen, L.-J., Moore, T. E., Ergun, R. E.,. . . Chandler, M. (2016, June). Electron-scale measurements of magneticreconnection in space. Sci. , , aaf2939. doi: 10.1126/science.aaf2939Daughton, W. (2003, August). Electromagnetic properties of the lower-hybrid driftinstability in a thin current sheet. Phys. Plasmas , , 3103-3119. doi: 10.1063/1.1594724Daughton, W., & Roytershteyn, V. (2012, November). Emerging Parameter SpaceMap of Magnetic Reconnection in Collisional and Kinetic Regimes. SpaceSci. Rev. , , 271-282. doi: 10.1007/s11214-011-9766-zGenestreti, K., Varsani, A., Burch, J., Cassak, P., Torbert, R., Nakamura, R., . . .others (2018). Mms observation of asymmetric reconnection supported by 3-delectron pressure divergence. J. Geophys. Res. , (3), 1806–1821.Groˇselj, D., Cerri, S. S., Ba˜n´on Navarro, A., Willmott, C., Told, D., Loureiro, N. F.,. . . Jenko, F. (2017, September). Fully Kinetic versus Reduced-kineticModeling of Collisionless Plasma Turbulence. Astrophys. J. , , 28. doi:10.3847/1538-4357/aa894dHakim, A. (2008). Extended mhd modelling with the ten-moment equations. J. Fu-sion Energy , (1-2), 36-43. Retrieved from http://dx.doi.org/10.1007/s10894-007-9116-z doi: 10.1007/s10894-007-9116-zHakim, A., Loverich, J., & Shumlak, U. (2006, November). A high resolution wavepropagation scheme for ideal Two-Fluid plasma equations. J. Comp. Phys. , , 418-442. doi: 10.1016/j.jcp.2006.03.036Hammett, G. W., & Perkins, F. W. (1990, June). Fluid moment models for Lan-dau damping with application to the ion-temperature-gradient instability. Phys. Rev. Lett. , , 3019-3022. doi: 10.1103/PhysRevLett.64.3019 –16–anuscript submitted to JGR-Space Physics
Juno, J., Hakim, A., TenBarge, J., Shi, E., & Dorland, W. (2018, January). Dis-continuous Galerkin algorithms for fully kinetic plasmas.
J. Comp. Phys. , ,110-147. doi: 10.1016/j.jcp.2017.10.009Krall, N. A., & Liewer, P. C. (1971, November). Low-Frequency Instabilities in Mag-netic Pulses. Phys. Rev. A , , 2094-2103. doi: 10.1103/PhysRevA.4.2094Le, A., Daughton, W., Chen, L.-J., & Egedal, J. (2017, March). Enhanced electronmixing and heating in 3-D asymmetric reconnection at the Earth’s magne-topause. Geophys. Res. Lett. , , 2096-2104. doi: 10.1002/2017GL072522Ng, J., Hakim, A., Juno, J., & Bhattacharjee, A. (2019). Drift instabilities in thincurrent sheets using a two fluid model with pressure tensor effects. J. Geo-phys. Res. . (Submitted)Ng, J., Huang, Y.-M., Hakim, A., Bhattacharjee, A., Stanier, A., Daughton, W., . . .Germaschewski, K. (2015, November). The island coalescence problem: scal-ing of reconnection in extended fluid models including higher-order moments.
Phys. Plasmas , , 112104.Parashar, T. N., Vasquez, B. J., & Markovskii, S. A. (2014, February). The roleof electron equation of state in heating partition of protons in a collisionlessplasma. Phys. Plasmas , (2), 022301. doi: 10.1063/1.4863422Paschmann, G., Papamastorakis, I., Sckopke, N., Haerendel, G., Sonnerup, B. U. O.,Bame, S. J., . . . Elphic, R. C. (1979, November). Plasma acceleration at theearth’s magnetopause - Evidence for reconnection. Nature , , 243-246. doi:10.1038/282243a0Price, L., Swisdak, M., Drake, J. F., Cassak, P. A., Dahlin, J. T., & Ergun, R. E.(2016, June). The effects of turbulence on three-dimensional magnetic re-connection at the magnetopause. Geophys. Res. Lett. , , 6020-6027. doi:10.1002/2016GL069578Romero, H., Ganguli, G., Lee, Y. C., & Palmadesso, P. J. (1992, July). Electron-ionhybrid instabilities driven by velocity shear in a magnetized plasma. Phys. Flu-ids B , , 1708-1723. doi: 10.1063/1.860028Shi, E. L., Hakim, A. H., & Hammett, G. W. (2015, February). A gyrokineticone-dimensional scrape-off layer model of an edge-localized mode heat pulse. Phys. Plasmas , (2), 022504. doi: 10.1063/1.4907160Srinivasan, B., & Shumlak, U. (2011, September). Analytical and computationalstudy of the ideal full two-fluid plasma model and asymptotic approxima-tions for Hall-magnetohydrodynamics. Phys. Plasmas , (9), 092113. doi:10.1063/1.3640811Swisdak, M. (2016, January). Quantifying gyrotropy in magnetic reconnection. Geo-phys. Res. Lett. , , 43-49. doi: 10.1002/2015GL066980Vasyliunas, V. M. (1975, February). Theoretical models of magnetic field line merg-ing. I. Reviews of Geophysics and Space Physics , , 303-336.Wang, L., Germaschewski, K., Hakim, A., Dong, C., Raeder, J., & Bhattacharjee,A. (2018, April). Electron Physics in 3-D Two-Fluid 10-Moment Model-ing of Ganymede’s Magnetosphere. J. Geophys. Res. , , 2815-2830. doi:10.1002/2017JA024761Wang, L., Hakim, A. H., Bhattacharjee, A., & Germaschewski, K. (2015, January).Comparison of multi-fluid moment models with particle-in-cell simulationsof collisionless magnetic reconnection. Phys. Plasmas , (1), 012108. doi:10.1063/1.4906063Zenitani, S., & Nagai, T. (2016, October). Particle dynamics in the electron currentlayer in collisionless magnetic reconnection. Phys. Plasmas , (10), 102102.doi: 10.1063/1.4963008(10), 102102.doi: 10.1063/1.4963008