A fully implicit, conservative, non-linear, electromagnetic hybrid particle-ion/fluid-electron algorithm

The quasi-neutral hybrid model with kinetic ions and fluid electrons is a promising approach for bridging the inherent multi-scale nature of many problems in space and laboratory plasmas. Here, a novel, implicit, particle-in-cell based scheme for the hybrid model is derived for fully 3D electromagnetic problems with multiple ion species, which features global mass, momentum and energy conservation. The scheme includes sub-cycling and orbit-averaging for the ions, with cell-centered finite differences and implicit midpoint time advance. To reduce discrete particle noise, the scheme allows arbitrary-order shape functions for the particle-mesh interpolations and the application of conservative binomial smoothing. The algorithm is verified for a number of test problems to demonstrate the correctness of the implementation, the unique conservation properties, and the favorable stability properties of the new scheme. In particular, there is no indication of unstable growth of the finite-grid instability for a population of cold ions drifting through a uniform spatial mesh, in a set-up where several commonly used non-conservative schemes are highly unstable.

Related papers:

title: Hybrid Simulation Codes: Past, Present and Future—A Tutorial
abstract: Hybrid codes, in which the ions are treated kinetically and the electrons are assumed to be a massless fluid, have been widely used in space physics over the past two decades. These codes are used to model phenomena on ion inertial and gyro-radius scales, which fall between longer scales obtained by magnetohydrodynamic simulations and shorter scales attainable by full particle simulations. In this tutorial, the assump- tions and equations of the hybrid model are described along with some most commonly used numerical implementations. Modifications to include finite electron mass are also briefly discussed. Examples of results of two-dimensional hybrid simulations are used to illustrate the method, to indicate some of the tradeoffs that need to be addressed in a realistic calculation, and to demonstrate the utility of the technique for problems of contemporary interest. Some speculation about the future direction of space physics research using hybrid codes is also provided.

title: The link between shocks, turbulence, and magnetic reconnection in collisionless plasmas
abstract: Global hybrid (electron fluid, kinetic ions) and fully kinetic simulations of the magnetosphere have been used to show surprising interconnection between shocks, turbulence, and magnetic reconnection. In particular, collisionless shocks with their reflected ions that can get upstream before retransmission can generate previously unforeseen phenomena in the post shocked flows: (i) formation of reconnecting current sheets and magnetic islands with sizes up to tens of ion inertial length. (ii) Generation of large scale low frequency electromagnetic waves that are compressed and amplified as they cross the shock. These “wavefronts” maintain their integrity for tens of ion cyclotron times but eventually disrupt and dissipate their energy. (iii) Rippling of the shock front, which can in turn lead to formation of fast collimated jets extending to hundreds of ion inertial lengths downstream of the shock. The jets, which have high dynamical pressure, “stir” the downstream region, creating large scale disturbances ...

title: Role of Ion Kinetic Physics in the Interaction of Magnetic Flux Ropes.
abstract: To explain many natural magnetized plasma phenomena, it is crucial to understand how rates of collisionless magnetic reconnection scale in large magnetohydrodynamic (MHD) scale systems. Simulations of isolated current sheets conclude such rates are independent of system size and can be reproduced by the Hall-MHD model, but neglect sheet formation and coupling to MHD scales. Here, it is shown for the problem of flux-rope merging, which includes this formation and coupling, that the Hall-MHD model fails to reproduce the kinetic results. The minimum sufficient model must retain ion kinetic effects, which set the ion diffusion region geometry and give time-averaged rates that reduce significantly with system size, leading to different global evolution in large systems.

title: Low frequency fully kinetic simulation of the toroidal ion temperature gradient instability
abstract: A fully kinetic ion model is useful for the verification of gyrokinetic turbulence simulations in certain regimes, where the gyrokinetic model may break down due to the lack of small ordering parameters. However, for a fully kinetic model to be of value, it must first be able to accurately simulate low frequency drift-type instabilities typically well within the domain of gyrokinetics. Here, a fully kinetic ion model is formulated with weak gradient drive terms and applied to the toroidal ion-temperature-gradient (ITG) instability for the first time. Implementation in toroidal geometry is discussed, where orthogonal coordinates are used for particle dynamics, but field-line-following coordinates are used for the field equation allowing for high resolution of the field-aligned mode structure. Variational methods are formulated for integrating the equation of motion allowing for accuracy at a modest time-step size. Linear results are reported for both the slab and toroidal ITG instabilities. Good agreement ...

title: Computer simulation using particles
abstract: Computer experiments using particle models A one-dimensional plasma model The simulation program Time integration schemes The particle-mesh force calculation The solution of field equations Collisionless particle models Particle-particle/particle-mesh algorithms Plasma simulation Semiconductor device simulation Astrophysics Solids, liquids and phase changes Fourier transforms Fourier series and finite Fourier transforms Bibliography Index

title: On magnetic reconnection regimes and associated three‐dimensional asymmetries: Hybrid, Hall‐less hybrid, and Hall‐MHD simulations
abstract: Magnetic reconnection in a plane current sheet is investigated in both two and three dimensions, using three different types of simulation codes, Hall MHD, hybrid (electron fluid, kinetic ions), and a new code called Hall-less hybrid. The latter code, which is similar to the hybrid code but has the Hall term removed, enables us to clarify the differences between kinetic ion and Hall MHD approaches. The major findings of this research are (1) Sweet-Parker regime of reconnection cannot be maintained and does not reach a steady state in a kinetic plasma for physically interesting parameter regimes. (2) Fast asymptotic reconnection rate of 0.15VA0B0 is obtained both in the hybrid and Hall-less hybrid simulations with outflow boundaries. VA0 and B0 are the Alfven velocity and magnetic field strength in the upstream region. This finding has two immediate implications. First, ion kinetics are sufficient to lead to fast reconnection even in the absence of the Hall term, and second, explanation of fast reconnection in terms of quadratic disperion of whistlers needs to be reconsidered, as whistlers are dispersionless in Hall-less hybrid limit. (3) Unlike in MHD, diffusion region is different in size that the region of localized resistivity. (4) While both Hall and hybrid codes show that reconnection is inherently asymmetric in three dimensions, there are differences in the nature of the asymmetry. In Hall MHD we show that the X-line grows in the direction of the electron drift, propagating as a (reconnection) wave because the current is carried by electrons, although the wave direction can change in the presence of a substantial ion flow. However, in the hybrid simulations here, as is the case for typical conditions at the magnetopause and magnetotail, ions carry the bulk of the current, and the observed asymmetry is found to be due to ion flow and not a wave motion unless the extent of finite resistivity in the third dimension is very thin, comparable to the current sheet thickness. Thus aside from scenarios where electrons are the dominant current carriers, such as very thin current sheets that are on electron scales, we do not expect the reconnection wave to form. This result is relevant to the magnetotail where dawn-dusk asymmetries are observed in the motions of auroral brightenings and surges, as well as in the statistical location of pressure decreases, flows, and magnetic signatures associated with the near-Earth neutral line and early plasmoids. We attribute the observed dawn-dusk asymmetries to ion flows. One interesting question left for future work is the possibility that reconnection waves may form in thin electron-scale current layers that are sometimes observed embedded within a thicker sheet in the magnetotail.

title: Quasineutral hybrid simulation of macroscopic plasma phenomena
abstract: A method for solving the quasineutral hybrid plasma equations in two dimensions is presented, using full ion dynamics and inertialess electrons. The method uses a predictor-corrector field solver and is extended to allow plasma-vacuum interfaces of arbitrary shape. A three-region method for treating the plasma-vacuum interfaces makes possible the simulation of slowly evolving phenomena over time scales much longer than the ion cyclotron period. The algorithm is applied to the study of rotational instabilities in theta pinch Vlasov equilibria.

title: Electromagnetic ion beam instabilities: Comparison of one- and two-dimensional simulations
abstract: Numerical simulations of electromagnetic instabilities, driven by a cool tenuous ion beam propagating along an ambient magnetic field, have been conducted in one and two spatial dimensions. The calculations employ particle ions, fluid electrons, and a predictor corrector scheme for solving the electromagnetic field in two dimensions that is described in some detail. While the principal features of the one-dimensional calculations (which reproduce previous work) are retained, the two-dimensional simulations show some reduction of the overall level of the magnetic field fluctuations. Enhancement of the heating of the beam ions at the expense of the core ions also occurs in the case where the beam density is sufficiently large that the right-hand nonresonant instability dominates. Implications of the results for modeling of the ion foreshock and quasi-parallel shocks are discussed.

title: Pegasus: A new hybrid-kinetic particle-in-cell code for astrophysical plasma dynamics
abstract: We describe Pegasus, a new hybrid-kinetic particle-in-cell code tailored for the study of astrophysical plasma dynamics. The code incorporates an energy-conserving particle integrator into a stable, second-order-accurate, three-stage predictor-predictor-corrector integration algorithm. The constrained transport method is used to enforce the divergence-free constraint on the magnetic field. A @df scheme is included to facilitate a reduced-noise study of systems in which only small departures from an initial distribution function are anticipated. The effects of rotation and shear are implemented through the shearing-sheet formalism with orbital advection. These algorithms are embedded within an architecture similar to that used in the popular astrophysical magnetohydrodynamics code Athena, one that is modular, well-documented, easy to use, and efficiently parallelized for use on thousands of processors. We present a series of tests in one, two, and three spatial dimensions that demonstrate the fidelity and versatility of the code.

title: Magnetic field and density fluctuations at perpendicular supercritical collisionless shocks
abstract: Perpendicular collisionless shocks are studied by means of two-dimensional hybrid (particle ion, fluid electron) simulations, with emphasis on the relaxation of the highly anisotropic ion distribution that arises primarily from reflection of some of the incident ions but also adiabatic compression of the other, directly transmitted ions and the related growth of low frequency electromagnetic waves. It is commonly assumed that the waves are due to Alfven ion cyclotron instability that propagate parallel to the ambient magnetic field (k/sub perpendicular/ = 0) and that the isotropization of the ions due to pitch angle scattering by the waves and the corresponding modification of the wave spectrum is quasi-linear. It is shown that this is indeed a reasonably good description downstream of the shock front, behind the magnetic overshoot. However, at the shock ramp there is a large discrepancy between the wavelengths measured in the simulations and those predicted by linear theory, and large density and magnetic field oscillations parallel to the ambient magnetic field are also seen. By comparing results for both high and low ion beta cases, it is shown that these effects can be understood in terms of obliquely propagating (k/sub perpendicular/not =0) modes, more likely due to themore » Alfven ion cyclotron instability instead of the (drift) mirror instability, although a more complete explanation awaits the derivation and analysis of an appropriate dispersion relation describing the growth and coupling of low-frequency modes in the inhomogeneous high beta environment of the shock. The observational consequences of these results and their application to improving nonlocal leakage models for the ion foreshock are also discussed.

title: Current Advance Method and Cyclic Leapfrog for 2D Multispecies Hybrid Plasma Simulations
abstract: CAM-CL (current advance method and cyclic leapfrog) is a new algorithm for hybrid plasma simulations. In common with existing methods, its physical basis is a "hybrid" plasma model which treats the ions as particles and the electrons as a massless fluid. CAM-CL is distinguished from previous 2D hybrid algorithms by four main features: (1) Multiple ion species may be treated with only a single computational pass through the particle data: this is achieved without extrapolation of the electric field in time. The particles are advanced by a leapfrog procedure which requires the electric field to be a half time-step ahead of the particle velocities. The electric field depends on the ionic current density and hence the particle velocities. In order to avoid a time-consuming "pre-push" of the velocities, CAM advances the ionic current density a half-step with an appropriate equation of motion. This is similar in concept to the moment method (D. Winske and K. B. Quest,  (1988), Appendix A), except for the next two features: (2) CAM advances the ionic current density, whereas the moment method advances the fluid velocity. Consequently, multiple ion species may be easily treated (3) A free-streaming ionic current density is collected (velocities are collected at positions a half time-step ahead). An equation of motion is then applied, in which the advective term and the ionic stress tensor in the moment method are not needed, since transport effects are included in the free-streaming current. (4) CL is a leapfrog scheme for advancing the magnetic field, an adaptation of the modified midpoint method described by W. H. Press et al. (Numerical Recipes ( 1986)). It is stable and allows sub-stepping of the magnetic field (the magnetic field time-step may be different to the particle time-step). A two-dimensional version of the algorithm has been tested on a quiet plasma, MHD wave propagation, and ion beam instabilities, the results of which are discussed.

title: Re‐forming supercritical quasi‐parallel shocks: 1. One‐ and two‐dimensional simulations
abstract: Hybrid (particle ion, massless fluid electron) simulations are carried out in one and two spatial dimensions to show that the supercritical quasi-parallel shock is not steady. Instead of a well-defined narrow transition region (shock ramp) between upstream and downstream states that remains at a fixed position in the flow, the ramp periodically steepens, broadens, and then re-forms upstream of its former position. During part of the re-formation cycle, the number of ions backstreaming from the shock increases to a significant fraction (up to 40%) of the incoming ions. This dense ion beam then couples to the incident ions to form a heated region, the leading edge of which becomes the new shock ramp. Two-dimensional simulations demonstrate that this process is not a consequence of restricting the calculations to one spatial dimension and show that the re-formation cycle is not in phase along the shock surface, although the surface itself remains fairly laminar. Generally, the downstream wave turbulence is comparable in all three magnetic components (unlike one dimensional where B{sub n} = const) and shows no tendency to be aligned in any particular direction. Both the two-dimensional simulations as well as the one-dimensional runs indicate that the downstream plasma is amore » complex mixture of regions of merged incident and reflected ion beams. The simulations are consistent with recent observational studies at the quasi-parallel bow shock showing the presence of magnetic pulsations, cold ion beams, and a downstream consisting of alternating regions of low-density hot ions and higher density cooler ions.« less

title: A second-order semi-implicit δfδf method for hybrid simulation
abstract: A second-order accurate semi-implicit Lorentz force ions, fluid electrons δf hybrid model has been developed using a “current closure” scheme. The model assumes quasi-neutrality and is fully electromagnetic. The implicit field solver improves numerical accuracy by separating the equilibrium terms in the presence of small perturbations. The equilibrium part of the generalized Ohm’s law is solved by direct matrix inversion along the direction of gradients for every Fourier mode in the other two directions, while the nonlinear part is solved iteratively. The simulation has been benchmarked on Alfven waves, ion sound waves and whistler waves against analytical dispersion relation in a slab. In particular, the first-order and second-order schemes are compared by studying the numerical damping of whistler waves. The full evolution of the resistive tearing mode using the Harris sheet equilibrium is also investigated. The linear growth rate and mode structure are compared with the resistive MHD theory. Important tearing mode nonlinear phenomena such as the Rutherford regime and saturation are demonstrated. We also presented systematic study of Rutherford growth rates and saturation island width, which is consistent with previous MHD studies.

title: A non-staggered, conservative, V×B=0' finite-volume scheme for 3D implicit extended magnetohydrodynamics in curvilinear geometries
abstract: In the development of spatial difference schemes for magnetohydrodynamics (MHD), the preservation of continuum properties such as conservation of mass, momentum, and energy, as well as required electromagnetic constraints (∇⋅B→=∇⋅j→=0, where j→=∇×B→ is the electrical current), is desirable to preserve numerical accuracy. Moreover, simplicity of the scheme is also a desirable feature, particularly when an implicit implementation is considered (the focus of this paper). We propose here a finite-volume, cell-centered (non-staggered) scheme for the extended MHD formulation that: (1) is suitable for implicit implementations in arbitrary curvilinear geometries, (2) is conservative, (3) preserves both the magnetic field and the electrical current solenoidal to machine precision, and (4) is linearly and nonlinearly stable in the absence of numerical and physical dissipation. Crucial to the viability of the scheme is the use of a clever interpolation scheme (ZIP [Hirt,  (1968) 339–355]), the proper treatment of boundary conditions in curvilinear geometry, and a novel treatment of geometric source terms in the momentum equation that ensures their exact cancellation in the absence of pressure forces.

title: Rigorous charge conservation for local electromagnetic field solvers
abstract: In this paper we present a rigorous method for finding the electric field in two-and-one-half and three dimensional electromagnetic plasma simulations without resorting to computationally expensive transforms. A finite grid interpretation of the divergence equation ∇· J = -ϖρ/ϖt is offered which allows the current density and thus new local electric and magnetic field strengths to be determined directly from knowledge of charge motion.

title: Charge conservation in electromagnetic PIC codes; spectral comparison of Boris/DADI and Langdon-Marder methods
abstract: Non-charge conserving current collection algorithms for electromagnetic PIC plasma simulations may cause errors in Gauss's law. These errors arise from violations of the charge continuity equation, ∇ · J = −ϵπ/ϵt, which is turn cause errors in the irrotational part of E. 
 
Two techniques for reducing these errors are examined and compared: a modified Marder correction which corrects electric fields locally and primarily affects short wavelengths, and a Boris divergence correction, which solves Poisson's equation to correct the electric fields so that Gauss's law is enforced globally. The effect of each method on the spectrum of the error is examined. Computational efficiency and accuracy of the two techniques are compared: neither method is clearly superior. 
 
Cases examined include corrections in electromagnetic relativistic beam simulations, and a hot thermal plasma. In addition, the spectral comparison provides insight into the behavior of the schemes applied.

title: Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor
abstract: The new method of local electric current density assignment in the Particle-in-Cell code in Cartesian geometry is presented. The method is valid for an arbitrary quasi-particle form-factor assuming that quasi-particle trajectory over time step is a straight line. The method allows one to implement the PIC code without solving Poisson equation. The presented formula for the current density associated with the motion of a single quasi-particle is the unique linear combination of form-factor differences in consistency with the discrete continuity equation. The computation scheme is demonstrated in 2D and 3D.

title: An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm
abstract: This paper discusses a novel fully implicit formulation for a one-dimensional electrostatic particle-in-cell (PIC) plasma simulation approach. Unlike earlier implicit electrostatic PIC approaches (which are based on a linearized Vlasov–Poisson formulation), ours is based on a nonlinearly converged Vlasov–Ampere (VA) model. By iterating particles and fields to a tight nonlinear convergence tolerance, the approach features superior stability and accuracy properties, avoiding most of the accuracy pitfalls in earlier implicit PIC implementations. In particular, the formulation is stable against temporal (Courant–Friedrichs–Lewy) and spatial (aliasing) instabilities. It is charge- and energy-conserving to numerical round-off for arbitrary implicit time steps (unlike the earlier “energy-conserving” explicit PIC formulation, which only conserves energy in the limit of arbitrarily small time steps). While momentum is not exactly conserved, errors are kept small by an adaptive particle sub-stepping orbit integrator, which is instrumental to prevent particle tunneling (a deleterious effect for long-term accuracy). The VA model is orbit-averaged along particle orbits to enforce an energy conservation theorem with particle sub-stepping. As a result, very large time steps, constrained only by the dynamical time scale of interest, are possible without accuracy loss. Algorithmically, the approach features a Jacobian-free Newton–Krylov solver. A main development in this study is the nonlinear elimination of the new-time particle variables (positions and velocities). Such nonlinear elimination, which we term particle enslavement, results in a nonlinear formulation with memory requirements comparable to those of a fluid computation, and affords us substantial freedom in regards to the particle orbit integrator. Numerical examples are presented that demonstrate the advertised properties of the scheme. In particular, long-time ion acoustic wave simulations show that numerical accuracy does not degrade even with very large implicit time steps, and that significant CPU gains are possible.

title: A multi-dimensional, energy- and charge-conserving, nonlinearly implicit, electromagnetic Vlasov-Darwin particle-in-cell algorithm
abstract: For decades, the Vlasov–Darwin model has been recognized to be attractive for particle-in-cell (PIC) kinetic plasma simulations in non-radiative electromagnetic regimes, to avoid radiative noise issues and gain computational efficiency. However, the Darwin model results in an elliptic set of field equations that renders conventional explicit time integration unconditionally unstable. Here, we explore a fully implicit PIC algorithm for the Vlasov–Darwin model in multiple dimensions, which overcomes many difficulties of traditional semi-implicit Darwin PIC algorithms. The finite-difference scheme for Darwin field equations and particle equations of motion is space–time-centered, employing particle sub-cycling and orbit-averaging. The algorithm conserves total energy, local charge, canonical-momentum in the ignorable direction, and preserves the Coulomb gauge exactly. An asymptotically well-posed fluid preconditioner allows efficient use of large cell sizes, which are determined by accuracy considerations, not stability, and can be orders of magnitude larger than required in a standard explicit electromagnetic PIC simulation. We demonstrate the accuracy and efficiency properties of the algorithm with various numerical experiments in 2D–3V.

title: A global method of solving the electron-field equations in a zero-inertia-electron-hybrid plasma simulation code
abstract: A formulation of the electron momentum equation and Maxwell's field equations suitable for global solution in an r-z hybrid plasma simulation code has been derived. The assumption of zero electron inertia is made in the electron momentum equation and Maxwell's equations are used in the radiation-free or Darwin limit. These techniques make explicit use of the axisymmetric properties of the model to decouple the components of the model equations. Equations to self-consistently advance the electron temperature are not presently included in this scheme. The model equations which result from these considerations are two coupled, nonlinear, second order partial differential equations. These two equations are integrated in time by a noniterative ADI procedure along with the explicit particle-in-cell ion time integration procedure. The resulting nearly implicit electron-field algorithm treats wide variations in the local signal velocity without instability; this consideration is most important since pure vacuum regions are allowed. The global nature of the solution requires boundary conditions only on the boundaries of the simulation region; arbitrary intermixing of plasma vacuum regions requires only the simple detection of low density cells and does not require monitoring of plasma vacuum interfaces.

title: A robust method for handling low density regions in hybrid simulations for collisionless plasmas
abstract: A robust method to handle vacuum and near vacuum regions in hybrid simulations for space and astrophysical plasmas is presented. The conventional hybrid simulation model dealing with kinetic ions and a massless charge-neutralizing electron fluid is known to be susceptible to numerical instability due to divergence of the whistler-mode wave dispersion, as well as division-by-density operation in regions of low density. Consequently, a pure vacuum region is not allowed to exist in the simulation domain unless some ad hoc technique is used. To resolve this difficulty, an alternative way to introduce finite electron inertia effect is proposed. Contrary to the conventional method, the proposed one introduces a correction to the electric field rather than the magnetic field. It is shown that the generalized Ohm's law correctly reduces to Laplace's equation in a vacuum which therefore does not involve any numerical problems. In addition, a variable ion-to-electron mass ratio is introduced to reduce the phase velocity of high frequency whistler waves at low density regions so that the stability condition is always satisfied. It is demonstrated that the proposed model is able to handle near vacuum regions generated as a result of nonlinear self-consistent development of the system, as well as pure vacuum regions set up at the initial condition, without losing the advantages of the standard hybrid code.

title: An implicit δf particle-in-cell method with sub-cycling and orbit averaging for Lorentz ions
abstract: A second order implicit ?f Lorentz ion hybrid model with sub-cycling and orbit averaging has been developed to study low-frequency, quasi-neutral plasmas. Models using the full Lorentz force equations of motion for ions may be useful for verifying gyrokinetic ion simulation models in applications where higher order terms may be important. In the presence of a strong external magnetic field, previous Lorentz ion models are limited to simulating very short time scales due to the small time step required for resolving the ion gyromotion. Here, we use a simplified model for ion Landau damped ion acoustic waves in a uniform magnetic field as a test bed for developing efficient time stepping methods to be used with the Lorentz ion hybrid model. A detailed linear analysis of the model is derived to validate simulations and to examine the significance of ion Bernstein waves in the Lorentz ion model. Linear analysis of a gyrokinetic ion model is also performed, and excellent agreement with the dispersion results from the Lorentz ion model is demonstrated for the ion acoustic wave. The sub-cycling/orbit averaging algorithm is shown to produce accurate finite-Larmor-radius effects using large macro-time steps sizes, and numerical damping of high frequency fluctuations can be achieved by formulating the field model in terms of the perturbed flux density. Furthermore, a CPU-GPU implementation of the sub-cycling/orbit averaging is presented and is shown to achieve a significant speedup over an equivalent serial code.

title: Jacobian-free Newton-Krylov methods: a survey of approaches and applications
abstract: Jacobian-free Newton-Krylov (JFNK) methods are synergistic combinations of Newton-type methods for superlinearly convergent solution of nonlinear equations and Krylov subspace methods for solving the Newton correction equations. The link between the two methods is the Jacobian-vector product, which may be probed approximately without forming and storing the elements of the true Jacobian, through a variety of means. Various approximations to the Jacobian matrix may still be required for preconditioning the resulting Krylov iteration. As with Krylov methods for linear problems, successful application of the JFNK method to any given problem is dependent on adequate preconditioning. JFNK has potential for application throughout problems governed by nonlinear partial differential equations and integro-differential equations. In this survey paper, we place JFNK in context with other nonlinear solution algorithms for both boundary value problems (BVPs) and initial value problems (IVPs). We provide an overview of the mechanics of JFNK and attempt to illustrate the wide variety of preconditioning options available. It is emphasized that JFNK can be wrapped (as an accelerator) around another nonlinear fixed point method (interpreted as a preconditioning process, potentially with significant code reuse). The aim of this paper is not to trace fully the evolution of JFNK, nor to provide proofs of accuracy or optimal convergence for all of the constituent methods, but rather to present the reader with a perspective on how JFNK may be applicable to applications of interest and to provide sources of further practical information.

title: An optimal, parallel, fully implicit Newton–Krylov solver for three-dimensional viscoresistive magnetohydrodynamicsa)
abstract: The conceptual development and implementation of a scalable nonlinear solver for the time-dependent three-dimensional (3D) compressible resistive magnetohydrodynamics model (MHD) is discussed. The approach is based on Jacobian-free Newton–Krylov technology, preconditioned with multigrid methods for algorithmic scalability. The key to the approach is the reformulation of the hyperbolic MHD system into a parabolic one, which is amenable to multigrid techniques. Such reformulation (parabolization) aims to render the modified system block diagonally dominant (unlike the original MHD system, which is diagonally submissive for implicit time steps Δt larger than the explicit Courant–Friedrichs–Lewy time step). The algorithm has been tested on a variety of 2D and 3D configurations which demonstrate its excellent algorithmic scalability properties. In particular, it is shown that, serially, the CPU time for a given simulation scales linearly with the number of unknowns, and that large CPU gains (∼30) are attainabl...

title: Orbit-averaged implicit particle codes
abstract: The merging of orbit-averaged particle code techniques with recently developed implicit methods to perform numerically stable and accurate particle simulations are reported. Implicitness and orbit averaging can extend the applicability of particle codes to the simulation of long time-scale plasma physics phenomena by relaxing time-step and statistical constraints. Difference equations for an electrostatic model are presented, and analyses of the numerical stability of each scheme are given. Simulation examples are presented for a one-dimensional electrostatic model. Schemes are constructed that are stable at large-time step, require fewer particles, and, hence, reduce input-output and memory requirements. Orbit averaging, however, in the unmagnetized electrostatic models tested so far is not as successful as in cases where there is a magnetic field. Methods are suggested in which orbit averaging should achieve more significant improvements in code efficiency.

title: Heuristic stability theory for finite-difference equations☆
abstract: A simple method is proposed for investigating the computational stability of finite-difference equations. The technique is especially powerful because of its applicability to nonlinear equations with variable coefficients. The method, which is based on an examination of certain kinds of truncation errors, is illustrated by applying it to a simple linear difference equation. Then it is used to explain the origin of instabilities observed in calculations of one- and two-dimensional fluid flows.

title: Curvature‐compensated convective transport: SMART, A new boundedness‐ preserving transport algorithm
abstract: The paper describes a new approach to approximating the convection term found in typical steady-state transport equations. A polynomial-based discretization scheme is constructed around a technique called ‘curvature compensation’; the resultant curvature-compensated convective transport approximation is essentially third-order accurate in regions of the solution domain where the concept of order is meaningful. In addition, in linear scalar transport problems it preserves the boundedness of solutions. Sharp changes in gradient in the dependent variable are handled particularly well. But above all, the scheme, when used in conjunction with an ADI pentadiagonal solver, is easy to implement with relatively low computational cost, representing an effective algorithm for the simulation of multi-dimensional fluid flows. Two linear test problems, for the case of transport by pure convection, are employed in order to assess the merit of the method.

title: Nonphysical noises and instabilities in plasma simulation due to a spatial grid
abstract: A series of plasma numerical simulation has been performed in order to understand the enhancement of nonphysical noises and instabilities due to the use of a spatial grid. Several different superparticle models including the Nearest Grid Point (NGP) model, Cloud-in-Cell (CIC) or Particle-in-Cell (PIC) models, Lewis energy conserving code, and the multipole expansion code have been examined for a Maxwellian plasma and a one beam plasma using a one-dimensional, one-specie (electron) plasma. An instability was observed for all of the models when the Debye length was too small compared with the grid size. When the Debye length is comparable to the grid size, no instabilities were observed. However, the enhancement of noises at high frequencies (ω > 3ω pe may not always be negligible- even for long wavelength modes for the NGP model. For the NGP and CIC, PIC models, the experimental results are in good agreement with Langdon's theory. It is observed that the dipole expansion model, which is the first-order approximation to the multipole expansion scheme, is similar to CIC, PIC models in many respects and appears to be the same order of approximation.

title: Finite grid instability and spectral fidelity of the electrostatic Particle-In-Cell algorithm
abstract: The origin of the Finite Grid Instability (FGI) is studied by resolving the dynamics in the 1D electrostatic Particle-In-Cell (PIC) model in the spectral domain at the single particle level and at the collective motion level. The spectral fidelity of the PIC model is contrasted with the underlying physical system or the gridless model. The systematic spectral phase and amplitude errors from the charge deposition and field interpolation are quantified for common particle shapes used in the PIC models. It is shown through such analysis and in simulations that the lack of spectral fidelity relative to the physical system due to the existence of aliased spatial modes is the major cause of the FGI in the PIC model.

title: Finite-grid instability in quasineutral hybrid simulations
abstract: Analysis and simulation of numerical instability caused by the spatial grid is presented for the case of quasineutral hybrid plasma simulation. The simplest case of particle ions advanced in the ambipolar electric field due to an isothermal, charge neutralizing electron fluid is examined. Stationary plasma is destabilized by finite ion temperature with the threshold temperature Ti characterized by the ratio ZTe/Ti , where Z is the charge state of the ions and Te is the electron temperature. A cold drifting plasma is found to be unstable with typical growth rate γ ≈ 0.4-0.6 Cs/Δx (Δx is the grid spacing and Cs is the ion acoustic speed); simulations show that unstable growth is saturated by particle trapping leading to heated ions. For nearest grid point weighting, saturation, and threshold for growth are observed to be Ti ≈ 0.6 ZTe while for linear particle weighting, Ti ≈ 0.1 ZTe. These results are relevant to hybrid particle simulations of laser generated plasmas, a regime where ZTe/Ti ⪢ 1 is encountered.

title: A linear dispersion relation for the hybrid kinetic-ion/fluid-electron model of plasma physics
abstract: A dispersion relation for a commonly used hybrid model of plasma physics is developed, which combines fully kinetic ions and a massless-electron fluid description. Although this model and variations of it have been used to describe plasma phenomena for about 40 years, to date there exists no general dispersion relation to describe the linear wave physics contained in the model. Previous efforts along these lines are extended here to retain arbitrary wave propagation angles, temperature anisotropy effects, as well as additional terms in the generalized Ohm's law which determines the electric field. A numerical solver for the dispersion relation is developed, and linear wave physics is benchmarked against solutions of a full Vlasov-Maxwell dispersion relation solver. This work opens the door to a more accurate interpretation of existing and future wave and turbulence simulations using this type of hybrid model.

title: Numerical simulation of nonoscillatory mirror waves at the Earth's magnetosheath
abstract: The generation of nonoscillatory mirror waves is studied using a one-dimensional periodic hybrid electromagnetic simulation. The ion dynamics are treated exactly; the electrons are approximated as a finite pressure, massless fluid. Compression of the flux tubes in the magnetosheath causes a large pressure anisotropy, and it has been proposed that this anisotropy drives a mirror instability. The mirror waves have been identified by large amplitude fluctuations of the magnetic field, anticorrelated with pressure fluctuations. The simulations are initiated in a homogeneous high beta (beta = 2.5) plasma with the ambient magnetic field at various angles to the simulation axis. It is found that ion cyclotron waves are also driven by the pressure anisotropy, in competition with the nonoscillatory mirror waves. Simulations indicate that in a pure ¹H+ plasma the much faster growing ion cyclotron waves absorb the free energy in the anisotropy to the extent that mirror waves should not be observed. Analysis of the dispersion relations of mirror waves and ion cyclotron waves in the multi-component plasma indicates that 4He2+ and 16O6+ ions in the solar wind should stabilize the ion cyclotron waves sufficiently that the mirror waves become the dominant instability.

title: Ion anisotropy instabilities in the magnetosheath
abstract: Recent observations in Earth's magnetosheath have delineated several different kinds of magnetic fluctuation spectra below the proton cyclotron frequency. This paper provides a theoretical interpretation for some of these observations, describing solutions of the linear Vlasov dispersion equation for fully electromagnetic instabilities for particle distributions which model those observed in the magnetosheath. The model assumes a high-[beta], homogeneous plasma with bi-Maxwellian, anisotropic (T[sub [perpendicular]p]/T[sub [parallel]p]) protons; a similarly anisotropic, tenuous, doubly ionized helium component; and isotropic cool electrons. This model also assumes a uniform magnetic field and no relative average drifts among the components. This model yields three growing modes: the proton cyclotron anisotropy, the helium cyclotron anisotropy, and the mirror instabilities. The regimes in [beta][sub p] versus T[sub [perpendicular]p]/T[sub [parallel]p] space which correspond to maximum growth rates of the mirror and proton cyclotron anisotropy instabilities are computed, including helium densities appropriate for the magnetosheath; the results show very good agreement with the observations of mirror-like and proton-cyclotron-like events. This agreement with observations implies that the transition between cyclotron and mirror fluctuation dominance is consistent with linear theory. 33 refs., 7 figs., 1 tab.

title: Observation of mirror waves downstream of a quasi-perpendicular shock
abstract: We have performed auto and cross spectral analysis of the electron density and vector magnetic field fluctuations recorded aboard ISEE 1 and 2 in the magnetosheath close the Earth's bow shock when in quasiperpendicular configuration. The analysis shows the characteristic signatures of mirror type modes downstream of the undershoot of the shock.

title: Geospace Environmental Modeling (GEM) magnetic reconnection challenge
abstract: The Geospace Environmental Modeling (GEM) Reconnection Challenge project is presented and the important results, which are presented in a series of companion papers, are summarized. Magnetic reconnection is studied in a simple Harris sheet configuration with a specified set of initial conditions, including a finite amplitude, magnetic island perturbation to trigger the dynamics. The evolution of the system is explored with a broad variety of codes, ranging from fully electromagnetic particle in cell (PIC) codes to conventional resistive magnetohydrodynamic (MHD) codes, and the results are compared. The goal is to identify the essential physics which is required to model collisionless magnetic reconnection. All models that include the Hall effect in the generalized Ohm's law produce essentially indistinguishable rates of reconnection, corresponding to nearly Alfvenic inflow velocities. Thus the rate of reconnection is insensitive to the specific mechanism which breaks the frozen-in condition, whether resistivity, electron inertia, or electron thermal motion. The reconnection rate in the conventional resistive MHD model, in contrast, is dramatically smaller unless a large localized or current dependent resistivity is used. The Hall term brings the dynamics of whistler waves into the system. The quadratic dispersion property of whistlers (higher phase speed at smaller spatial scales) is the key to understanding these results. The implications of these results for trying to model the global dynamics of the magnetosphere are discussed.

title: Hybrid simulations of magnetic reconnection with kinetic ions and fluid electron pressure anisotropy
abstract: We present the first hybrid simulations with kinetic ions and recently developed equations of state for the electron fluid appropriate for reconnection with a guide field. The equations of state account for the main anisotropy of the electron pressure tensor. Magnetic reconnection is studied in two systems, an initially force-free current sheet and a Harris sheet. The hybrid model with the equations of state is compared to two other models, hybrid simulations with isothermal electrons and fully kinetic simulations. Including the anisotropic equations of state in the hybrid model provides a better match to the fully kinetic model. In agreement with fully kinetic results, the main feature captured is the formation of an electron current sheet that extends several ion inertial lengths. This electron current sheet modifies the Hall magnetic field structure near the X-line, and it is not observed in the standard hybrid model with isotropic electrons. The saturated reconnection rate in this regime nevertheless ...

title: A fully nonlinear characteristic method for gyrokinetic simulation
abstract: A new scheme that evolves the perturbed part of the distribution function along a set of characteristics that solves the fully nonlinear gyrokinetic equations is presented. This low‐noise nonlinear characteristic method for particle simulation is an extension of the partially linear weighting scheme, and may be considered an improvement over existing δf methods. Some of the features of this new method include the ability to keep all nonlinearities, particularly those associated with the velocity space, the use of conventional particle loading techniques, and also the retention of the conservation properties of the original gyrokinetic system in the numerically converged limit. The new method is used to study a one‐dimensional drift wave model that isolates the parallel velocity nonlinearity. A mode coupling calculation for the saturation amplitude is given, which is in good agreement with the simulation results. Finally, the method is extended to the electromagnetic gyrokinetic equations in general geometry.

title: A Particle-in-cell Method with Adaptive Phase-space Remapping for Kinetic Plasmas
abstract: We present a new accurate and efficient particle-in-cell (PIC) method for computing the dynamics of one-dimensional kinetic plasmas. The method overcomes the numerical noise inherent in particle-based methods by periodically remapping the distribution function on a hierarchy of locally refined grids in phase space. Remapping on phase-space grids also provides an opportunity to integrate a collisional model and an associated grid-based solver. The positivity of the distribution function is enforced by redistributing excess phase-space density in a local neighborhood. We demonstrate the method on a number of standard plasma physics problems. It is shown that remapping significantly reduces the numerical noise and results in a more consistent second-order method than the standard PIC method. An error analysis is presented which is based on prior results of Cottet and Raviart's work [ 21 (1984), pp. 52-76].

title: Fast magnetic reconnection with large guide fields
abstract: In this letter, it is demonstrated using two-fluid simulations that low-β magnetic reconnection remains fast, regardless of the presence of fast dispersive waves, which have been previously suggested to play a critical role. To understand these results, a discrete model is constructed that offers scaling relationships for the reconnection rate and dissipation region (DR) thickness in terms of the upstream magnetic field and DR length. We verify these scalings numerically and show how the DR self-adjusts to process magnetic flux at the same rate that it is supplied to a larger region where two-fluid effects become important. The rate is therefore independent of the DR physics and is in good agreement with kinetic results.

title: A 2D high-ß Hall MHD implicit nonlinear solver
abstract: A nonlinear, fully implicit solver for a 2D high-β (incompressible) Hall magnetohydrodynamics (HMHD) model is proposed. The task in non-trivial because HMHD supports the whistler wave. This wave is dispersive (ω ˜ k2) and therefore results in diffusion-like numerical stability limits for explicit time integration methods. For HMHD, implicit approaches using time steps above the explicit numerical stability limits result in diagonally submissive Jacobian systems. Such systems are difficult to invert with iterative techniques. In this study, Jacobian-free Newton-Krylov iterative methods are employed for a fully implicit, nonlinear integration, and a semi-implicit (SI) preconditioner strategy, developed on the basis of a Schur complement analysis, is proposed. The SI preconditioner transforms the coupled hyperbolic whistler system into a fourth-order, parabolic, diagonally dominant PDE, amenable to iterative techniques. Efficiency and accuracy results are presented demonstrating that an efficient fully implicit implementation (i.e., faster than explicit methods) is indeed possible without sacrificing numerical accuracy.