Simulations of anti-parallel reconnection using a nonlocal heat flux closure

Magnetic reconnection is a change in topology of the magnetic field lines in a plasma,¹ often followed by the conversion of stored magnetic energy to the kinetic energy of accelerated particles. It is believed to play an important role in many laboratory and astrophysical plasma processes, including sawtooth crashes in tokamaks, solar flares, magnetic substorms in the Earth's magnetosphere, and coronal mass ejections.^2–6 

For reconnection to take place, the motion of the plasma must decouple from the magnetic field lines. In collisionless environments such as the magnetosphere, this takes place in the electron diffusion region, and due to the kinetic scales involved, it cannot be described purely by resistive magnetohydrodynamic (MHD) models such as the Sweet–Parker model.^7,8 For large scale global numerical studies of the magnetosphere, which are crucial for accurate space weather prediction, the use of fluid models is common due to the computational cost involved in simulating the large domains involved. Thus, the development of extended fluid models incorporating collisionless effects is desirable for capturing the physics of reconnection in large systems.

The major issue with fluid models is the closure of the hierarchy of fluid equations obtained by taking moments of the kinetic equation. While increasing the number of moments allows more information about the distribution function to be captured,⁹ a choice must be made as to when to close the equations and how to treat the highest moment. Most closures of the moment equations involve an expansion in small Knudsen number in the collisional regime, and attempts have been made to use these models to simulate magnetic reconnection.^10–13 

However, in the collisionless case, kinetic effects such as phase mixing are present and cannot be captured by these fluid closures.^14,15 Instead, approximations can be made by introducing collisionless damping coefficients by expressing the highest moments in terms of lower moments. Unlike the collisional case, the closure coefficients can be complex and are nonlocal in space.^14,15 While these models have been used extensively in studies of fusion plasmas,^16,17 they have not been applied to studies of magnetic reconnection.

In this paper we describe the self-consistent implementation of a fluid closure using the ten-moment equation system without any adjustable parameters. We evolve the full pressure tensor for both electrons and ions, which is necessary to balance the reconnection electric field in collisionless reconnection^18,19 and close for the heat flux using a three-dimensional extension of the Hammett–Perkins closure.¹⁴ We perform simulations of the island coalescence problem, in which the existing Hall MHD fluid model does not reproduce the kinetic result due to the absence of important ion physics.^20,21 We demonstrate here with the collisionless fluid model that the wider ion diffusion region can be reproduced without the fine tuning of a free parameter as was previously required.²¹ While the strong system-size scaling of the reconnection rate in kinetic and hybrid (kinetic ion, fluid electron) simulations is not completely reproduced for larger scales, the new closure successfully captures the ion-scale structures such as agyrotropy which are not accounted for consistently by earlier fluid models.

The ten-moment equations are derived by taking moments of the kinetic equation

where $Pij$ and $Qijk$ are the second and third moments of the distribution function

and the square brackets denote a sum over permutations of the indices (e.g., $u[iEj]=uiEj+ujEi$⁠). The third moment tensor $Qijk$ can be written in terms of the heat flux tensor $qijk≡m∫(vi−ui)(vj−uj)(vk−uk)fd3v$

Earlier multi-fluid models of reconnection involving a pressure tensor have employed a relaxation to local isotropy in the form $∂iqijk=nvt|k0|(Tij−T0δij)$⁠. Here, $vt=2T/m$ is the thermal velocity of the associated species and k₀ is a free parameter which effectively allows deviations from isotropy at length scales less than $1/|k0|$⁠.^10,11,13,21 It has been shown that simulation results are sensitive to the choice of k₀ for each species^13,21 and such free parameters are present in earlier efforts as well.^10,11,13

In order to close the equations in the collisionless limit, we use a three-dimensional extension of the Hammett–Perkins closure, which can be expressed as follows for both electrons and ions¹⁴ 

where $q̂ijk$ in Fourier space is $q̃ijk$ and is calculated as

Here, $T̃jk$ is the Fourier transform of the deviation of the local temperature tensor from the mean. The $1/|k|$ scaling makes this a non-local closure when expressed in real space^15,16 and provides a 1 to 3 pole Padé approximation of various components of the dielectric tensor. The coefficient $χ=4/9π$ is the best fit value for the diagonal q_iii component and reduces to the closure in Refs. 14 and 15 in the 1-D limit. This is an unmagnetised closure (especially relevant to ions in this geometry in which there is no guide field) which approximates linear phase mixing, allowing the wavenumber-dependent damping of spurious short-wavelength oscillations which are present in higher moment fluid models,^15,22 and to our knowledge, is the first of its kind in closure studies of magnetic reconnection.

The moment equations coupled to Maxwell's equations are implemented in the finite-volume version of the Gkeyll code, which uses a high-resolution wave propagation method for the hyperbolic part of the equations and a point implicit method for the source terms.^23,24 The closure is evaluated using an interface to the parallel FFTW library.²⁵ 

We first illustrate the effects of the closure in Fig. 1(a), which shows the evolution of the electric field amplitude for a Langmuir wave with kλ_de = 0.35 in a 1-D two-fluid simulation, where λ_de is the electron Debye length. In the simulation, the measured frequency and damping rate, in terms of the electron plasma frequency ω_pe, are ω = 1.16ω_pe and γ = 0.044ω_pe, respectively, compared to the exact solution from the linear Vlasov dispersion relation ω = 1.22ω_pe, γ = 0.034ω_pe. Here we note that the measured values are consistent with those calculated using the Padé approximant of the plasma dispersion function, and improved agreement with the Vlasov solution requires evolving higher fluid moments.^14–16 In contrast, the simulation without any explicit dissipation is qualitatively different and shows no damping as expected.

For magnetic reconnection, while the Harris sheet geometry is a standard test problem,²⁶ and has been successfully simulated using this closure (not shown), we focus on the island coalescence geometry in this paper as it has been found to be sensitive to ion rather than electron dynamics,^20,21 for which this closure would be a better approximation. This is a Fadeev equilibrium with initial conditions²⁷ 

Here B₀ is the x-component of the magnetic field upstream of the layer, ϵ controls the island size, and λ is the half width of the current sheet. We use the same physical parameters as described in previous studies,^20,21,28 with ϵ = 0.4, which corresponds to an island half-width of approximately 1.2λ, and background density n_b = 0.2n₀. The system size is L_x × L_z= 4πλ × 2πλ, with periodic boundary conditions in the x direction. Conducting walls for fields and reflecting walls for particles are used in the z direction. We use mass ratio m_i/m_e= 25, electron thermal speed v_t,e/c = 0.35, and T_i= T_e = T, and the value of T is set by the upstream equilibrium condition β = 1. The ratio of electron plasma frequency to gyrofrequency is ω_pe/Ω_ce = 2. A 10% initial perturbation in the magnetic field is applied to initiate merging. The kinetic (particle-in-cell) results in this paper are described in Ref. 20 and additional simulations for the λ/d_i = 5, 10 cases were performed using the particle-in-cell code PSC (particle simulation code), which is described in Ref. 29.

In MHD, multi-fluid and kinetic models, the initial phase in which the islands approach each other is driven by the ideal MHD coalescence instability,³⁰ and differences between the models only appear after the central current sheet begins to form.^20,21 This can be seen in Fig. 1(b), which shows the normalised reconnection rate for simulations using the new closure and earlier kinetic, ten-moment, and Hall MHD simulations.^20,21 Here the reconnection rate $ER=(1/B′VA′)∂ψ/∂t$ is found in the same manner as in Refs. 20, 21, and 28, where $B′$ and $VA′$ are calculated using the maximum magnetic field between the centres of the two islands at t = 0. The flux within an island ψ is defined as the difference between A_y at the X- and O-points.

Because reconnection in this system is bursty, it is useful to consider the average as well as maximum reconnection rates.^20,21,28 In Figs. 1(c) and 1(d), the scaling of the maximum and average (over 1.5 global Alfvén times) reconnection rates with the system size is shown. With the nonlocal closure, there is a stronger negative scaling of the average reconnection rate $∝(λ/di)−0.45$ for system sizes up to λ/d_i= 25 compared to $(λ/di)−0.2$ for both Hall MHD and ten moment simulations with a local closure. However, this is still weaker than the results of Ref. 20 where the rate scales like $∝(λ/di)−0.8$ and $(λ/di)−0.65$ for kinetic and hybrid models, respectively. The discrepancies for the ten moment systems are partially explained by the enhancement of the reconnection rate due to the formation of secondary islands in smaller systems (λ/d_i = 15 for the nonlocal closure, λ/d_i = 25 for the local closure). In the presence of multiple secondary islands, the location of the major X-point is determined by the maximum A_y along the line x = 0. The flux is measured between this point and the O-point at the centre of either of the initial islands.

Figure 2 shows a comparison between the kinetic, ten-moment, and Hall MHD simulations with λ = 5d_i at t = t_A, which is close to the time of maximum reconnection rate. Due to the coalescence geometry [Eq. (6)], the inflow is in the x (horizontal) direction while the outflow is in the z (vertical) direction, which is opposite to the usual Harris sheet. As can be seen in the bottom panel, the width of the current sheet in the ten-moment simulation is comparable to that in the kinetic simulation (1.20d_i vs 1.15d_i), and is much larger than the highly peaked Hall MHD current sheet (0.18d_i). The peak current density is still higher than in the kinetic simulation due to the electron dynamics in the ten-moment model, which causes a more intense electron layer due to the isotropisation of the pressure tensor. Here the electron closure captures the off-diagonal terms of the pressure tensor which balance the reconnection electric field in the electron diffusion region self-consistently, but produces weaker electron pressure anisotropy than models which focus on particle trapping effects;^31–33 this does not affect the results of this paper as the reconnection rates are sensitive to the ion dynamics for this geometry.^20,21

The structure of the ion diffusion region is shown in Fig. 3, which contains the decomposition of the ion momentum equation on a cut along z = 0 at t = t_A. In Ref. 20 it was shown that the influence of the ion kinetic physics extends to a broader 2–3 d_i ion diffusion region, where the divergence of the pressure tensor balances the non-ideal electric field. This result is in contrast to the Hall MHD result, where the ion inertia balances the reconnection electric field below d_i scales.²⁰ 

While both local and nonlocal ten-moment models could capture the wider ion diffusion region, it was found that keeping the relaxation parameter k_i = 1/d_e, which was used in the Harris sheet study,¹³ caused larger Hall MHD-like reconnection rates to be observed and a narrower ion diffusion region to form. Only by adjusting k_i in the ion closure to 1/(3d_i), allowing the pressure tensor to deviate from isotropy at larger scales, could the kinetic reconnection rates be reproduced for small systems.²¹ 

In contrast, with the nonlocal closure, the wider ion diffusion region develops self-consistently without the need to fine tune the relaxation. While the reconnection rate does not match the kinetic result exactly, the larger region in which the divergence of the ion pressure tensor balances the non-ideal electric field is reproduced, with the associated reduction of the reconnection rate as compared to that of Hall MHD or the local ten-moment model with larger k_i.²¹ 

Additionally, due to the wavenumber-dependent damping, the model using the nonlocal closure shows good agreement in the ion velocity with the kinetic results, apart from the absence of particle-in-cell noise. This is because the new closure does not exhibit some of the spurious structures which have been observed when using the local closure. These are most evident along the line x = 0, where the ion velocity is reduced, and just upstream of the reconnection region in Fig. 4. The structures are present with the local closure due to the lower k_i used to obtain the wider diffusion region, which does not account for the faster damping at shorter wavelengths, unlike the nonlocal closure, which approximates this effect.

Further improvements when using the nonlocal closure are illustrated by considering a common metric used in reconnection, the agyrotropy $AØ$⁠, which measures the departure of the distribution function from the cylindrical symmetry.³⁴ Figure 5 shows the ion agyrotropy $AØi$ observed close to the central current sheet at t = t_A. In the kinetic simulation, there is strong agyrotropy in both a 2–3d_i thick region upstream of the current sheet and the exhaust region and it was shown in Ref. 20 that this corresponds to the region where ion motion is characterised by meandering orbits and ion pressure tensor effects contribute to momentum balance.

With the new nonlocal closure shown in the second panel of Fig. 5, the structure of the agyrotropy both up- and downstream of the x-point is reproduced. This is a qualitative improvement on the local model, which shows two distinct upstream regions and small agyrotropy in the exhaust region. The improvement is due to the nonlocal heat flux, which reduces the discrepancies in the diagonal components of the ion pressure tensor observed when only using the local model.¹² However, the values of the agyrotropy using the nonlocal closure are still lower than the kinetic results, as the fluid model cannot capture all the details of the non-gyrotropic distribution function caused by the meandering motion.

To summarise, we have implemented a nonlocal ion closure for the two-fluid, ten-moment equations that captures linear phase mixing in a fluid treatment. This provides heat flux at the x-point which is necessary for reconnection to take place in the collisionless ten moment model.¹² We demonstrate that using the new nonlocal closure gives improved agreement with kinetic results compared to previously implemented systems for moment closures and Hall MHD,^13,21 by allowing the ion pressure tensor to deviate from isotropy without introducing a free relaxation parameter. This sets the structure of the ion diffusion region self-consistently, and results in much improved agreement highlighted by the ion velocity profile and agyrotropy. However, the strong system-size dependence of the average reconnection rate observed in kinetic and hybrid simulations is not completely reproduced.²⁰ The scaling is due in part to the formation of secondary islands in the larger fluid systems, and also indicates that while the ion diffusion region physics is reproduced, there is still additional ion kinetic physics, such as the ion meandering motion and its effects on the distribution function,²⁰ that contributes to setting the reconnection rate in this system. While this system is less sensitive to the electron dynamics,^20,21 the nonlocal closure is also able to capture the off diagonal terms of the electron pressure tensor in the electron diffusion region. Future work could explore the extension to the strongly magnetised regime in a similar manner,^15,35 which could allow the treatment of guide-field reconnection and application to large systems, in which the development of strong electron pressure anisotropy is also an important process.^32,33

This work was supported by NSF Grant Nos. AGS-0962698 and AGS-1338944, DOE Contract No. DE-AC02-09CH11466 and NASA Grant No. NNH13AW51I. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and Trillian, a Cray XE6m-200 supercomputer at the UNH supported by the NSF MRI program under Grant No. PHY-1229408. We thank G. Hammett for fruitful discussions.