{'fraction_non_alphanumeric': 0.042153580527656706, 'fraction_numerical': 0.016107636914913778, 'mean_word_length': 4.2888641425389755, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 multi-dimensional electromagnetic problems with multiple ion species, which features global mass, momentum and energy conservation. The scheme includes sub-cycling and orbit averaging of 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.


Introduction
Hybrid plasma models, in which some component of the plasma is treated kinetically and a fluid description is used for the remainder, are promising candidates to model multi-scale plasma phenomena [1,2,3,4,5]. Hybrid models typically feature higher fidelity than reduced fluid descriptions, while remaining more computationally feasible than solving full kinetic models. A variety of hybrid models are discussed in the literature, which differ mainly by the choice of the kinetic component and the calculation of the electric field. Here, the term will apply to the quasi-neutral model with full-orbit kinetic ions and fluid electrons. Such a model is valid for non-relativistic plasmas where the light wave is eliminated, and the electric field is calculated with Ohm's law rather than solving the full Maxwell equations. The kinetic equation for the ions is solved using the Particle-In-Cell (PIC) approach [6,7] where a large number of discrete macro-particles trace out the evolution of the distribution function using the method of characteristics. This approach avoids the need to solve for the distribution function on a 6D (3D-3V) grid and can be highly optimized to run on modern computer architectures with threading and vectorization, but suffers from discrete particle noise that scales as √ N with the finite number of macro-particles N. There is significant literature on the development of numerical algorithms based on the hybrid PIC ion and fluid electron model, driven by the applications of space weather modeling and magnetic fusion research, see Refs. [1,2] for reviews. The majority of existing codes employ explicit time-stepping, and much of the development has focused on key algorithmic issues related to accuracy and stability of the time integration schemes. One issue arises from the use of leapfrog time advance where the particle positions and momenta are staggered in time, while statically evaluating Ohm's law for the electric field using these as sources. This results in an interdependence in which the particle velocity update and the calculation of the electric field depend implicitly on each other [1,8]. To make this explicit, some form of the predictor-corrector method [9,10,11], moment method [12,13], or extrapolation [14,15,8] has been used previously. A further issue with explicit methods concerns the stable propagation of the dispersive whistler waves that are present in the hybrid model. In the absence of dissipation, it has been shown that two-step explicit schemes such as the second-order Runge-Kutta method are unconditionally unstable to whistler waves [16,11]. For stability in higher-order explicit methods, it is necessary to resolve the quadratic CFL-limited time-step associated with dispersive Whistler waves, ∆t ∝ n(∆x) 2 /B, where ∆t is the time-step, ∆x the grid spacing, n the density, and B include the options to use plasma resistivity to break the electron frozen-in condition, and also electron viscosity that can provide dissipation in a narrow band in k-space to target only the smallest scales. For energy conservation in the presence of dissipative terms in Ohm's law, however, it is necessary to use a separate electron pressure evolution equation with corresponding heating terms.
The continuum hybrid description that is used here is given by the following equations (γ − 1) −1 ∂ t p e + ∇ · (u e p e ) + p e ∇ · u e = H e − ∇ · q e , where f s is the distribution function for ion species s, v is the velocity space co-ordinate, E * is the electric field (E) minus the collisional friction (see below), B is the magnetic field, n = n e = 1 e s q s f s d 3 v = 1 e s q s n s is the quasineutral number density, u i = s q s f s vd 3 v /en = s q s n s u s /en is the ion current carrying drift velocity, p e is the electron pressure, j = ∇ × B/µ 0 is the current density, u e = u i − j/ne is the bulk electron velocity, ← → Π e = −µ e ∇u e is the electron viscous tensor, q e = −κ e ∇(p e /n) is the electron heat flux, and H e = η j 2 − ← → Π e : ∇u e is the Ohmic and electron viscous heating. The coefficients are the charge and mass of species s, q s and m s respectively, the ratio of specific heats γ, the resistivity η, the electron viscosity µ e , and the electron heat conductivity κ e .

Conservation properties in the continuum 2.1.1. Mass conservation and quasi-neutrality
Taking the zeroth moment of the Vlasov equation (1) gives the local mass conservation equation for ions of species s ∂ t n s + ∇ · (n s u s ) = 0.
Due to quasi-neutrality there is no need to solve a separate electron fluid mass conservation equation. Instead, this property requires the ambipolarity condition to be satisfied for consistency. Multiplying Eq. (5) by the ion charge q s , taking the sum over species s, and using the definitions of n, u i , u e , and j, gives

Momentum conservation
The ion momentum density equation for species s is found from the first moment of the Vlasov equation (1).
where ← → P s = m s vv f s d 3 v is the pressure tensor for species s. Taking the sum over species s gives the total ion momentum Using the definition of E * in Eq. (3), and the identity (∇ × B) × B = ∇ · [BB − (B 2 /2)I] gives the local momentum density conservation equation Here, the reason for using the frictionless electric field E * becomes clear: using the total electric field would give a frictional source of total momentum. The reason for this is the ad-hoc but widely used [17,2,1,16,20] choice to include resistive friction (due to ion-electron collisions) in Ohm's law (3) without an explicit collision operator in the kinetic ion equation (1). We follow this choice here, and leave consideration of self-consistent ion-electron collisions for future work.
The total momentum s m s n s u s dV is constant in a domain V only for suitable boundary conditions, where the flux terms vanish. Examples include periodic boundaries, but exclude open boundaries with a net momentum flux, or conducting walls.

Energy conservation
Taking the contracted second moment of the ion Vlasov equation gives an equation for κ s = m s v 2 f s d 3 v/2, the kinetic energy density for ions of species s where Q s = m s v 2 v f s d 3 v/2 is the energy flux. Taking the sum over species s, and using the definition of u i gives An equation for the magnetic energy density B 2 /2µ 0 (Poynting's theorem) can be found from the Faraday equation (2) which is the divergence of the Poynting flux and the magnetic dissipation respectively. Note the total electric field is used here. Taking the sum of the ion kinetic energy in Eq. (10), the magnetic energy in Eq. (11) and the electron internal energy in Eq. (4) gives the total energy Again, the total energy [ s κ s + p e /(γ − 1) + B 2 /2µ 0 ]dV is conserved in a domain V only for suitable boundary conditions, such as periodic or conducting walls (with zero viscous stress and heat-flux losses).

Discrete model
Our discrete momentum-energy conserving scheme is constructed on the key attributes of a cell-centered discretization in space, and an implicit midpoint method in time. For this choice, special consideration must be taken in the discretization of advection operators to handle problems with strong flows. This will be discussed further below and in the numerical examples. Derivatives are approximated by central differences, which are defined as These discretizations naturally preserve many of the continuum properties of the discrete operators, such as (∇ · (∇ × N) i jk ) i jk = 0, e.g. [41]. For quantities χ that are known at integer time-step, the definition at the half-integer time-step is given as χ n+1/2 = 1 2 (χ n+1 + χ n ). For the remainder of the paper, a vector potential formulation is adopted using the Weyl gauge (φ = 0), such that B = ∇ × A and E = −∂ t A. This choice is made for ease of implementation into an existing fluid framework, see below, but we note the conservation properties derived below and also the solenoidal condition (∇ · B) i jk = 0 hold for a magnetic field formulation discretized with the central differences defined above, e.g. [41]. For clarity, we also leave the numerical description of the electron viscosity/hyper-resistivity to Appendix A.
The discrete Ohm's law and electron pressure equations are given by p n+1 e,i jk − p n e,i jk ∆t + γ ∇ · u e p e n+1/2 i jk Here, j n+1/2 i jk are given below. For clarity, the conservation properties will be firstly demonstrated using the same time-step for the particle equations of motion as for the electromagnetic field and electron pressure equations. This is then generalized to include particle sub-cycling with orbit-averaged moments below. Also, from hereon, the species index s is incorporated into the macro-particle index p without loss of generality. The equations of motion for the macro-particle p with charge q p and mass m p are given by where x n p , v n p is the known position and velocity of particle p at time t = n∆t, and B n+1/2 p , E * , n+1/2 p are the magnetic and frictionless electric fields interpolated to the particle position at t = (n + 1/2)∆t as where S (x i − x p ) is, for now, an arbitrary order shape function for the particle-grid interpolations [6].
where ∆V = ∆x∆y∆z is the cell volume.

Discrete mass and quasi-neutrality
As is typical for particle methods, the particle mass, m p , and charge, q p , are constant along the characteristic particle trajectories, and are thus locally conserved in the Lagrangian sense. To conserve mass and charge globally it is then sufficient to conserve the total number of particles, which can be done for suitable boundary conditions in the absence of particle sources and sinks.
In the fully kinetic Vlasov-Maxwell system it is important to conserve charge locally on the spatial mesh, because small truncation errors in the local charge conservation equation can accumulate to give significant violations in Gauss' law (see Refs. [42,43] for discussions). For the quasi-neutral hybrid model, Gauss' law has been removed from the continuum equations, and instead the discrete form of Eq. (6) needs to be satisfied locally at the locations where the density moment is collected. This is satisfied for our discrete model, from the definition of u n+1/2 e,i jk and the discrete divergence and curl operators that give (∇ · j) n+1/2 i jk = 0. It is possible to gather density and momentum such that they satisfy a local mass conservation equation on the spatial grid using techniques developed for the full Vlasov-Maxwell model, e.g. Refs. [42,44,22,25]. Appendix B describes a generalization of the method of Refs. [22,25] for three spatial dimensions, where the density and momentum measures are defined at different locations on a staggered grid. We have chosen not to directly use such measures here, as we find they do not give momentum conservation for the hybrid model. However, as described in Section 4.2, we use the same particle sub-stepping scheme as Refs. [22,25] where the ions strictly deposit contributions to the density and momentum at all of the cells they pass through in a given time-step. The cell-centered measures we use lie within truncation error from the locally-mass-conserving measures, and can not increase unboundedly due to local mass conservation in the Lagrangian sense. Since global mass conservation also holds, these truncation errors cancel when integrated over the whole domain. Although we typically use quadratic spline shape functions (B.9,B.10), the derivations below hold for arbitrary order shape functions and, as such, we drop the subscript.

Discrete momentum conservation
We define the total momentum for all particles (of all species) at time t = n∆t by P n = p m p v n p . The discrete rate of change in total momentum is then given by  (21), and using Ohm's law from Eq. (14). For suitable boundaries, e.g. periodic, we can use the discrete integration by parts (telescoping the sums). That is, for discrete scalar fields χ i jk , ξ i jk , i jk (δ x χ) i jk ξ i jk = − i jk χ i jk (δ x ξ) i jk . Using this property, and noting that the sum over gradients vanishes in periodic domains, gives from the definitions of the discrete divergence and curl operators given above.
In the above, it has been shown that the total momentum of the system is conserved in the sense that it has no contribution from spatial or temporal truncation error. In practice, there are two other sources of numerical error in an implicit algorithm: numerical round-off error due to finite precision arithmetic, and non-linear convergence error associated with an iterative implicit scheme. It will be shown below that the largest of these errors depends on the precise implementation within the non-linear solver.

Discrete energy conservation
The total ion kinetic energy at t = n∆t is defined as K n i = p m p (v n p ) 2 /2. Then the discrete rate of change of K i is given by using Eqs. (17), (19), (21), and Ohm's law from Eq. (14). We define the sum of magnetic energy at cell centers and integer time-step as W n The discrete rate of change of total magnetic energy is then by discrete integration by parts, summing fluxes to zero, and using Eq. (13).
Since we neglect electron inertia, the sum of electron kinetic energy (electron thermal energy) is defined as K n e = i jk ∆V p n e,i jk /(γ − 1). To demonstrate discrete total energy conservation, it is sufficient to show that the total energy at the n + 1 time-step is equal to that of the nth time-step, ie.
The discrete change in total energy between the n and n + 1 time-step is given by the total ion kinetic energy, the magnetic energy, and the electron pressure equation (15), as by cancellation of terms, and using the definition of u n+1/2 e,i jk . The property i jk ∇ · u e p e n+1/2 i jk = 0 has also been assumed here, which is satisfied for all conservative definitions of this advective flux term. We have not yet fixed the definition of this term, to allow for the flexible use of monotonicity preserving schemes in problems with strong flows, see below. In this derivation, the implicit midpoint scheme allows for the cancellation of field and particle contributions to the total energy change without temporal truncation error. In practice, energy is conserved down to the non-linear tolerance of the iterative method, as will be discussed below.

Non-linear solver strategy
The temporal and spatial discretization described above yields the set of algebraic equations (13)(14)(15)(16)(17)(18)(19)(20)(21), which can be written in the form F(x n+1 ) = 0 for the sought solution vector at the new time-step x n+1 = (x n+1 p , v n+1 p , A n+1 i jk , p n+1 e,i jk ). The system contains a number of non-linear implicit couplings, and to solve it in a scalable manner requires the use of iterative methods. Here we use the Jacobian-free Newton-Krylov (JFNK) method to converge the residual F to within a specified non-linear tolerance t . Our general implementation has been well documented in Refs. [45,46], but we give an overview here to motivate the discussion of our specific implementation below. The Newton linearization gives an equation for the solution at non-linear iteration k + 1.

∂F ∂x
is an absolute tolerance to avoid trying to converge below round-off (N is the total number of degrees of freedom), r is the relative Newton convergence tolerance, and F(x 0 ) is the initial residual. At every Newton iteration Eq. (23) is solved using the flexible GMRES Krylov subspace method [47]. Constructing and storing the Jacobian matrix J = ∂F/∂x can be extremely memory intensive and it unnecessary in practice. Instead, the Gateux derivative is used to calculate the needed matrix-vector products as ∂F ∂x for a small and finite value of . As in previous JFNK implementations, we use an inexact Newton's method where the convergence tolerance of the Krylov method varies for each Newton iteration depending on the size of the residual, see Refs. [46,22] for details. A key advantage of Krylov methods is that they can be preconditioned by solving the alternate systems ( . When (P k ) −1 ≈ (J k ) −1 the convergence properties of the Krylov iterations can be dramatically improved, while the solution of the Jacobian system is unchanged upon convergence. We leave the discussion of the preconditioning strategy for a future publication, and take the preconditioning matrix to be the identity matrix P k = I in this paper. For the present discussion, we will use an electromagnetic field formulation taking the total solution vector to be In practice it is not necessary, and is indeed inefficient, to place the whole solution vector in the Newton residual as F(x n+1 ). Non-linear elimination (e.g. Ref. [22]) can be used to reduce the number of dependent variables in the residual from F(X 1 , X 2 ) to G(X 1 ) for new residual G when the dependent variable X 2 can be written explicitly as X 2 = f (X 1 ). For the present set of equations there are two notable choices in the form of the residual, which differ in the convergence properties for momentum conservation and in algorithmic efficiency. The first is to form the Newton residual as G(v n+1 p , B n+1 i jk , p n+1 e,i jk ), such that (v k+1 p , B k+1 i jk , p k+1 e,i jk ) is found from the (k + 1)th iteration of (23). The new position can be then directly calculated within the function evaluation for the next Newton iteration as x k+1 i jk , p k+1 e,i jk ). The primary advantage of this implementation is that the scatter and gather operations that are performed either side of the static evaluation of E k+1 i jk are done with the same particle position x k+1 p , in the same manner as in derivation (22). This ensures momentum conservation at each Newton iteration to round-off, and thus gives a momentum error that is independent of the non-linear convergence tolerance. However, a major drawback of this method is due to the presence of the particle quantities v n+1 p in the residual, which can be extremely memory intensive for typical problems where the number of particles greatly exceeds the number of grid points.
The second option is to solve for the Newton residual H(E n+1/2 i jk , p n+1 e,i jk ) = 0 and solve for B n+1 in the function evaluation. The major advantage of this method is the reduction in memory usage, and flexibility in the method for the particle push. Namely, a different time-step can be used to integrate the particle orbits (sub-cycling, see below), and the particle push can be optimized to run efficiently on the latest many-cores architectures using threading and vectorization. Due to the implicit coupling v n+1 p (x n+1 p (v n+1 p )), which arises from the interpolation of the fields to the particle position, an iterative method must be used to non-linearly converge the particles at each Newton iteration (see below). The momentum conservation is ensured down to the maximum of the particle orbit integration error and the outer non-linear convergence tolerance, which is demonstrated numerically in Appendix C. In this paper, we use the latter method modified for a potential formulation, i.e. with a residual H( A n+1 i jk , p n+1 e,i jk ). This allows us to easily implement the electromagnetic field and fluid equations in an existing Hall-MHD algorithm that uses a potential formulation.

Sub-cycled and orbit-averaged particle push
Taking the same time-step ∆t for both the particle orbit integration and the field advance can be inaccurate, e.g. Ref. [22]. For accuracy in the particle integration, it is necessary to resolve sufficiently the gyro-frequency Ω ci ∝ B that can vary strongly across the domain, depending upon the local magnetic field strength. However, solving for the electromagnetic fields at every time-step is unnecessary for the modeling of low-frequency phenomena where ω << Ω ci . To address this, we use sub-stepping in the particle orbit integration. The sub-step for each particle ∆τ ν p is calculated using a local error estimator using the same method as Ref. [25]. Here the macro time-step for the field solve ∆t relates to the particle sub-steps as ∆t = It is assumed that the electromagnetic fields are static over the particle sub-steps, such that variations in the particle force occur only due to the change in particle position and velocity. Thus, the particle equations of motion are given by where To close the system of equations, the moments are computed via orbit-averaging [48,22,21], where a contribution to the moment is made at each particle sub-step. A significant advantage of this approach is in the reduction of numerical noise in the gathered moments since a phase-space position for each macro-particle is sampled multiple times within ∆t. Our key consideration here is to preserve momentum and energy conservation for the orbit averaged formulation. It is straightforward to show that the above momentum-energy conservation holds when the moments are gathered as Note that the plasma density n n+1/2 i jk is also computed using orbit averaging here. In the implicit energy conserving Vlasov-Darwin algorithm [22], the density is calculated without orbit averaging at integer time-steps since it is only used in the charge conservation diagnostic. If the density is not orbit-averaged in the present case, then discrete momentum conservation is lost. The sub-cycled particle update equations are solved with a Picard-iterated implicit Boris method [22,25,49]. The positions and momenta at Picard iteration k are denoted by x k p , v k p , where x k=0 p = x ν p and v k=0 p = v ν p are the known values at the previous sub-step. The following steps are then iterated: where a particle Picard tolerance of p = 10 −13 is typically employed. When this is satisfied the particle positions and momenta are updated for the next sub-step as Here, α k = ∆τ k p q p /2m p , and dl x , dl y and dl z are the distances in each direction to the closest cell boundary in the direction of motion. These are used to prevent particles from crossing cell boundaries within a sub-step, and so contribute to the density and velocity moments within every cell that they cross.

Conservative smoothing
To reduce discrete particle noise, a user-defined number of binomial smoothing passes can be applied to the field and moment quantities. For a cell-centered quantity Q i jk , one pass of the smoothing operator is defined as SM(Q i jk ) = SM k (SM j (SM i (Q i jk ))) where the smoothing in each direction is defined as For suitable boundary conditions, e.g. periodic, the smoothing operator has the property i A i jk SM i (B i jk ) = i SM i (A i jk )B i jk . It can be shown that global momentum and energy conservation is unaffected by such smoothing when the moments are gathered as and the fields are scattered as This property is demonstrated numerically in Section 6.

Advection
The conservative advection term in the electron pressure equation (15), ∇ · u e p e n+1/2 i jk , gives no contribution to the global energy conservation error for suitably chosen boundary conditions. This allows some flexibility in the choice of its discretization. Here we allow a number of options, which are inherited from an existing MHD fluid framework [41]. For linear problems we typically use the symmetric ZIP differencing [50,41], which satisfies the discrete chain rule property and does not suffer from anti-diffusive truncation errors [41]. However, this can be non-monotonic, particularly for non-linear problems with strong flows such as the reconnection problem discussed below. For such problems we instead use the Sharp and Monotonic Algorithm for Realistic Transport (SMART) algorithm [51].

Normalization
The discrete model is implemented in normalized form using natural ion units. For a typical magnetic field strength B 0 , density n 0 , and ion mass m 0 , velocities are normalized by the Alfvén speed v 0 = v A0 = B 0 / √ m 0 µ 0 n 0 , times by the inverse cyclotron period t 0 = Ω −1 c0 = m 0 /eB 0 , and lengths by the ion skin depth L 0 = d i0 = v A0 /Ω c0 . The different species of kinetic ions are distinguished by their normalized charge Z s = q s /e, mass M s = m s /m 0 , and density N s = n s0 /n 0 . If they are initialized with a Maxwellian distribution function with temperature T s0 , we define the thermal velocity as v T s0 = √ T s0 /m s and plasma-beta β s0 /2 = M s N s (v T s /v A0 ) 2 = 2µ 0 n s0 T s0 /B 2 0 . The electron sound speed is defined as C s = √ T e0 /m 0 and the ratio of ion to electron temperatures is τ s = T s0 /T e0 .

Numerical test of the finite-grid instability for a cold drifting beam
A major limitation of particle-in-cell algorithms to model the full Vlasov-Maxwell system is the need to resolve the electron Debye length to prevent unstable numerical heating. Such heating is typically demonstrated for the problem of a cold electron beam moving through a uniform periodic mesh at constant speed with a charge neutralizing ion background [6,32]. The system is physically stable, however, spatial aliasing errors from particle-mesh interpolations can lead to a purely numerical instability that causes violation of momentum or energy conservation (or both), depending on the numerical scheme used [32,33,35]. Finite-grid instabilities have also been studied using an explicit momentum conserving algorithm for the quasi-neutral hybrid model. Ref. [36] presents a linear analysis of the hybrid finite grid instability with a discrete spatial grid for various shape functions, including nearest grid point (NGP) and linear interpolation/cloud-in-cell (CIC). Both the distribution function and the time advance is assumed to be continuous. It was found that a cold but finite ion temperature (0 < T i /T e < 1) can make the beam unstable, where the precise instability threshold depends on the order of particle-mesh interpolation used. For example, instability occurs for T i /T e < 1/12 for NGP and T i /T e < 1/25 for CIC. Finite drift is further destabilizing and leads to growth rates γ in the range 0.1 < γ∆x/C s < 1 for T i /T e < 20 using NGP interpolation.
Determining whether a PIC algorithm that simultaneously preserves both momentum and energy is absolutely stable with respect to finite-grid instabilities is a question left for future study, since it requires a dedicated analysis. Here, we present numerical evidence in support of this proposition using a particularly challenging case of a very cold (T i /T e = 1/600), drifting (u 0 /C s = 0.1) beam with NGP interpolation and with no smoothing. We compare the implicit conservative scheme with two other commonly used non-conservative schemes. The first scheme is a simple explicit leapfrog scheme that conserves momentum to numerical round-off, see Appendix D.1 and Ref. [36]. The second algorithm is an explicit predictor-corrector scheme, see Appendix D.2, that uses an additional particle push to calculate the time-advanced electric field, but does not conserve momentum or energy due to truncation error. The problem set-up is chosen to be the same as for Figures 3 and 4 of Ref. [36]; using a 1D-1V periodic domain of length 40∆x, 50 particles per cell, and time-step ∆t = 0.022∆x/C s . Figure 1 shows the numerical results for the cold-beam problem for the three algorithms. The top panel shows the particle phase-space distribution (x p , v xp ) after 400 time-steps (t = 8.8∆x/C s ) for the momentum-conserving leapfrog algorithm (red), the explicit predictor-corrector scheme (green), and the new implicit conservative scheme (blue). The phase space distributions for the first two (explicit) algorithms show characteristic fluctuations of the finite-grid instability, while the implicit algorithm appears stable. The effects of these fluctuations on the momentum and energy conservation are shown in the second and third panels in Fig. 1. In the second panel, the momentum error, (P n x − P 0 x ), remains of the order of 10 −14 for the two schemes that feature discrete momentum conservation, while there are significant errors in momentum conservation (≈ 0.03) for the predictor-corrector scheme. The total energy grows by 20 − 30% for both of the explicit schemes, in contrast to the new implicit scheme that conserves energy down to ≈ 10 −13 (for relative non-linear tolerance r = 10 −14 ). We note that discrete global energy conservation in the implicit conservative scheme does not preclude conversion of energy between the ion and electron kinetic energies. Indeed, spatial aliasing errors are present in all particle-in-cell schemes with a spatial grid, and can cause small variations in the energies of each ion particle. It is thus instructive to also consider the ion-to-electron temperature ratio, which is one measure of this energy partition between the kinetic and fluid components. The bottom panel shows T i /T e for the three schemes on a logarithmic axis. Here, the two explicit schemes exhibit linear growth which saturates at T i /T e = 0.8 for the momentum-conserving leapfrog-based algorithm and T i /T e = 0.85 for the predictor-corrector scheme. The saturation level differs by a small amount from the value of T i /T e = 0.625 mentioned in Ref. [36], and  we have verified that these differences are due to the use of an adiabatic, rather than isothermal, equation of state and our choice of discretization of the electric field (the explicit momentum-conserving algorithm in Appendix D.1 is equivalent to Ref. [36] in the linear regime). There is a small decrease (< 1%) in the ratio of T i /T e in the implicit conservative scheme which is non-monotonic over the timescale plotted. However, we see no sign of the unstable growth of T i /T e due to the finite grid instability for 10 6 time-steps (t = 22000∆x/C s ), which is as long as we have run the simulation.

Numerical Verification
In the remainder of the paper, we present numerical examples that are intended to verify the correctness of our implementation, and demonstrate the conservation properties as described above. For iterative implicit methods, such as the one discussed here, the cost of the algorithm strongly depends on the effectiveness of the preconditioner used. This is particularly the case for implicit particle methods, where the number of particle pushes is proportional to the number of function evaluations, and the cost of applying a grid-based preconditioner is small compared to that of each particle push. A thorough discussion of the preconditioning strategy and the cost to reach a given non-linear convergence is left for a separate publication. Here, we will typically converge the non-linear tolerance close to numerical round-off error in the following examples, to demonstrate the absence of any spatial or temporal truncation errors in the momentum and energy conservation. The choice of the macro time-step ∆t is not limited by stability considerations, but we typically choose it small to reduce stiffness in the absence of a preconditioner. In practice, the conservation properties are only ensured down to the chosen non-linear tolerance, which is demonstrated in Appendix C. We do not use heat conductivity for any of the problems (κ e = 0), and the resistive and hyper-resistive dissipation is set to zero in all examples except for the non-linear reconnection problem.

Landau damped ion acoustic wave (1D-1V electrostatic)
The first numerical example we consider is the 1D-1V electrostatic ion acoustic wave (IAW). Due to ion kinetic effects, the wave is damped via ion Landau damping. We initialize the problem with a single-ion species Maxwellian plasma with ion charge Z i = 1, mass M i = 1, thermal velocity v T i0 = √ 1/3, and quasi-neutral density N i = 1 + 0.01 sin (kx) with wavenumber k = π/8 and domain length 2π/k. The electron fluid has γ = 5/3 and we use a temperature ratio of τ = 0.2, which gives an electron sound speed of C s = √ 5/3. The dispersion relation is given by Z (ξ) = 2τ, where Z is the plasma dispersion function and ξ = (ω + iγ)/k √ 2v T i0 is the normalized complex frequency. Figure 2a) shows the L 2 -norm of the density perturbation ||dn|| = ||N i − 1|| from two simulations with 96 uniform cells, and 50, 000 particles per cell that are sampled from the ion distribution function using a lowdiscrepancy Hammersley sequence (quasi-quiet start). The simulations are performed with a time-step of ∆t = 0.02 and are converged to a relative Newton tolerance of r = 10 −12 using ion sub-cycling and orbit averaging with a Picard tolerance of p = 10 −13 . This tight tolerance typically requires 6 non-linear iterations and 6 linear iterations per time-step. The green curve shows the result with two applications of binomial smoothing to the fields and moments as described above, whereas the blue curve is the same run without any smoothing. In both cases, good agreement is found with linear theory for the real frequency ω = 0.6827 (2.13 k √ 2v T i0 ) and damping rate γ = −0.05265 where the latter is indicated by the red line in Fig. 2a). The initial quasi-quiet start degrades over time due to the accumulation of numerical error, until the noise floor overwhelms the damped signal. Since there is significant separation between the grid resolution and the wavelength of the ion acoustic wave, binomial smoothing can be safely applied in this problem to reduce the impact of discrete particle noise. As shown in Fig. 2a), two passes of binomial smoothing significantly improves the solution at late time. Figure 2b) shows the momentum error (P n x − P 0 x ), the electron thermal energy K n e , the ion kinetic energy K n i , and the relative total energy error [(K i + K e ) n − (K i + K e ) 0 ]/(K i + K e ) 0 for the two cases. The chosen relative tolerance gives errors on the level of double precision round-off: 10 −15 − 10 −16 for the momentum error, and 10 −14 for the energy error. These momentum and energy errors are unaffected by binomial smoothing, as discussed above, although there are some minor differences in the energy partition (K n e vs K n i ) between the two runs.

Proton Cyclotron Anisotropy Instability 1D-3D electromagnetic
The Proton Cyclotron Anisotropy Instability (PCAI) is a fully electromagnetic instability driven by anisotropy in the proton temperature with respect to the background magnetic field with T p⊥ /T p > 1. The mode has finite real frequency with maximum growth at k × B 0 = 0 for wave vector k and background field B 0 . We simulate the problem in 1D-3V, where both the wave vector k and the background field B 0 are parallel to the x-axis. The first case considered is a pure hydrogen plasma with Z p = M p = N p = 1 drawn from a uniform bi-Maxwellian distribution with β p = 1, and T p⊥ /T p = 3. The second case uses the same parameters but is comprised of 20% doubly ionized helium ions (alpha particles), with Z α = 2, M α = 4, T α /T p = 2, T α⊥ /T α = T p⊥ /T p , N α /N p = 0.2 and Z α N α + Z p N p = 1. In both cases, the electron fluid is initialized with temperature ratio τ p = 1 and γ = 5/3. The simulation domain has a length of 10.5, using 64 cells with 50, 000 protons and 20, 000 helium ions per cell. The time-step is ∆t = 0.01 and relative Newton tolerance r = 10 −12 , which requires about 5 non-linear and 6 linear iterations per time-step. We use a quasi-quiet start and apply two passes of binomial smoothing. Figure 3 shows the L 2 -norm of the transverse bulk ion momentum ||n i u iy || vs time. For the pure hydrogen plasma (green) there is good agreement between the simulation and the theoretical linear growth rate (red, γ = 0.162 [52,53]) over approximately six decades of growth in the perturbation amplitude. Studies of pressure anisotropy driven instabilities in the magnetosheath have found that the linear growth rate of the PCAI can reduce in the presence of doubly ionized helium ions. Theory has shown that the reduction in the growth rate can lead to a change from cyclotron to mirror-like fluctuations in some regions of the Earth's magnetosheath [54,55], in agreement with in-situ spacecraft data, e.g. Ref. [56]. For the parameters used here, the linear growth rate is reduced to γ = 0.0785 (black fit) and there is good agreement with the fully non-linear simulation that includes the two ion species (blue). For a relative Newton tolerance of r = 10 −12 used here, the momentum and relative total energy errors (not shown) are on the order of 10 −14 and 10 −12 , respectively, for both of the non-linear simulations.

Non-linear magnetic reconnection
The final numerical example is a 2D-3V simulation of magnetic reconnection. Parameters are chosen to reproduce the GEM challenge reconnection problem [57], but with two current layers and double periodic boundary conditions. The simulation domain is a square box of dimensions Lx = 25.6, Ly = 25.6 with 384 × 384 cells and an average of 320 particles/cell. The time-step is ∆t = 0.01 and we use a relative Newton tolerance r = 10 −13 . A single ion species with Z i = M i = 1 is initialized from an isotropic Maxwellian distribution function with β i0 = 5/6 (the total ion+electron β = 1). The initial density is non-uniform with function The simulation is initiated with a large sinusoidal perturbation with amplitude 0.1B 0 . The electron fluid is initialized with γ = 5/3 and a temperature ratio of τ = 5. In contrast to the previous numerical examples, finite dissipation is included here to break the frozen in condition and allow reconnection. We use a uniform resistivity η = 5 × 10 −4 and hyper-resistivity η H = 2 × 10 −4 to break the frozen-in condition and allow reconnection. Figure 4 shows the results and conservation properties from the simulation. The top left panel shows the current density j z and contours of the vector potential A z at t = 20 during the non-linear phase. As is typical for the hybrid model, e.g. Ref. [8], the out-of-plane current j z is localized around the magnetic X-point. We note that an anisotropic closure for the electron pressure can capture elongated current sheets within the hybrid model [38], which are more similar to those found in full kinetic simulations, particularly in the low-β regime. Associated with these magnetic and current structures are strong outflows, which can lead to unphysical oscillations when using ZIP differencing to treat the conservative advection term in the electron pressure equation. As discussed in Section 4.4, we find that using the SMART algorithm for this advection term preserves monotonicity and removes such oscillations. The bottom left panel in Fig. 4 shows the reconnected flux over time. Both the reconnection rate, (∂ t A z )| X ≈ 0.25 measured at the Xpoint, and the saturation level agree with those reported in Ref. [57] for a single current layer set-up. The momentum and energy conservation results are shown in the right panel. In this case the momentum error in each component direction is on the order of 10 −10 − 10 −11 and the relative total energy error is of the order 10 −10 . Due to the significant stiffness from the whistler wave dispersion and the high-order hyper-resistive dissipation, the simulation typically requires 12 nonlinear and 200 linear iterations per time-step to converge to this tight tolerance. This motivates the development of an optimal preconditioning strategy for the hybrid model, which we will pursue in future work.

Summary
A novel implicit, non-linear, particle-in-cell scheme has been formulated for the electromagnetic quasi-neutral hybrid model featuring global conservation of mass, momentum, and energy. The scheme uses the implicit-midpoint method for time advance, with field and fluid quantities defined at cell centers and using central differences for the numerical derivatives. Momentum and energy conservation is ensured by exact cancellation between particle moment quantities with terms in Ohm's law and the electron pressure equation that that exist on the spatial grid. These conservation properties hold for arbitrary-order shape functions (provided that the same shape function is used throughout), and are also preserved when sub-cycling the particle advance using a Picard-iterated implicit Boris push, collecting orbit-averaged moments, and applying conservative smoothing. The discrete equations are implemented within an iterative Jacobian-free Newton Krylov solver, and we have compared two possible forms of the Newton residual in terms of the convergence of the momentum conservation error, and the potential for optimization of the particle advance routines. The favorable stability properties of the algorithm are demonstrated for the problem of a cold beam moving through a uniform periodic mesh [36]. There is no indication of the finite grid instability for the new conservative scheme, for a challenging set of parameters in which two commonly used non-conservative schemes are unstable. Future work will explore these stability properties in more detail, comparing with linear theory for the instability that includes the effects of a spatial grid.
The correctness of implementation, and the conservation properties, are then demonstrated for a number of test problems including the Landau damping of the ion acoustic wave, the Proton Cyclotron Anisotropy Instability (PCAI), and non-linear magnetic reconnection. For the non-linear problem, it was found that oscillations due to sharp gradients and strong flows were eliminated by the use of the SMART scheme [51] for the conservative advection term in the electron pressure equation. Additional future work will explore the possibilities of discrete particle noise reduction via deltaf methods [58] or remapping [59], and the implementation of an optimal preconditioner for the electron model.
In a typical symmetric set-up, where the location of the magnetic X-point coincides with the stagnation point for the ion and electron flows, some non-ideal term is required to balance the reconnection electric field. Using only resistivity is often inadequate, because it is unable to support sub-layers in the electron current which thin over time [60,61]. Here, we include the flexibility to use higher-order hyper-resistive friction. This is able to support the electron layer even when reconnection becomes quasi-steady [61], and also has an added benefit that it can set a dissipation scale for dispersive whistler waves that can otherwise lead to grid-scale noise in non-linear simulations, e.g. Ref. [62].
We treat the hyper-resistivity as a friction term, which is added to the right hand side of Eq. (13) but is not included in E * , n+1/2 i jk . We discretize the term using a compact stencil as where the second derivative is defined as (d 2 x j) n+1/2 = ( j n+1/2 i+1 jk − 2 j n+1/2 i jk + j n+1/2 i−1 jk )/∆x 2 . The term has no effect on momentum conservation, but, since it causes dissipation of magnetic energy, it is necessary to add a corresponding hyper-resistive viscous heating term in the electron pressure equation for total energy conservation. In discrete form, global energy is conserved for suitable boundary conditions when the heating term is defined as where the spatial derivatives are defined as (d x j α ) n+1/2 i+1/2 jk = ( j α ) i+1 jk − ( j α ) i jk n+1/2 /∆x. The conservation of global energy is demonstrated numerically for the reconnection problem in Section 6.3. , p n+1 e,i jk ) and with particle orbit Picard-iteration tolerances p = 10 −4 (red x), 10 −8 (green ·), and 10 −12 (blue ).
Appendix C. Non-linear convergence of momentum and energy errors Figure C.5 shows the convergence properties of the momentum and energy errors for the two different implementations of the discrete model, discussed in Section 4.1, for the ion acoustic wave (IAW) problem. The problem set-up is the same is described in Section 6.1. Since there is zero magnetic field for the IAW problem, the two choices for the Newton residual are G(v n+1 p , p n+1 e,i jk ) (black +) and H(E n+1/2 i jk , p n+1 e,i jk ). For the latter, the particles are advanced with the Picard-iterated implicit Boris push as described in Section 4.2, but, for a fair comparison between the two methods, we fix the number of particle sub-steps N νp = 1. For this implementation, the results are plotted for a particle orbit convergence tolerance of p = 10 −4 (red ×), p = 10 −8 (green ·), and p = 10 −12 (blue ). The left panel shows how the momentum error, max|(P n − P 0 ) ·x|, reduces with the L 2 -norm of the residual in each implementation. For the first implementation with residual G(v n+1 p , p n+1 e,i jk ), the momentum is conserved to numerical round-off error regardless of the size of the residual norm. However, as discussed above, such an implementation is impractical for multi-scale problems. For the second implementation with residual H(E n+1/2 i jk , p n+1 e,i jk ) the momentum error depends on both the size of the residual and the particle convergence tolerance. For a loose particle tolerance of p = 10 −4 , errors from the particle advance become dominant for ||residual|| < 10 −10 , where the momentum error becomes independent of residual norm. For tighter particle convergence tolerances of 10 −8 and 10 −12 , the particle error is negligible and the momentum error reduces with the residual norm down to the level of numerical round-off. The relative total energy error, max[(K i + K e ) n − (K i + K e ) 0 ]/(K i + K e ) 0 , is shown in the right panel for the same simulations. In contrast to the momentum error, there is strong dependence of the energy error on the non-linear convergence for both implementations. For the loose particle tolerance p = 10 −4 , errors from the particles can dominate the energy error for ||residual|| < 10 −6 , but for the tighter tolerances the energy error can be converged down to numerical round-off. Note that small differences in the residual norm between implementations are due to different sized discrete jumps with each non-linear iteration, and should not be taken to be significant. We fix the particle convergence tolerance to p = 10 −13 to ensure that the errors are controlled by the outer Newton iteration. This typically requires about 4 Picard iterations per particle sub-step at every function evaluation.