A new hybrid code (CHIEF) implementing the inertial electron fluid equation without approximation 1footnote 11footnote 1 ©2018. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/

A new hybrid code (CHIEF) implementing the inertial electron fluid equation without approximation 111 ©2018. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/

P. A. Muñoz munozp@mps.mpg.de N. Jain P. Kilian J. Büchner Max-Planck-Institut für Sonnensystemforschung, D-37077 Göttingen, Germany Max-Planck/Princeton Center for Plasma Physics Centre for Space Research, North-West University, 2520 Potchefstroom, South Africa
July 17, 2019
Abstract

We present a new hybrid algorithm implemented in the code CHIEF (Code Hybrid with Inertial Electron Fluid) for simulations of electron-ion plasmas. The algorithm treats the ions kinetically, modeled by the Particle-in-Cell (PiC) method, and electrons as an inertial fluid, modeled by electron fluid equations without any of the approximations used in most of the other hybrid codes with an inertial electron fluid. This kind of code is appropriate to model a large variety of quasineutral plasma phenomena where the electron inertia and/or ion kinetic effects are relevant. We present here the governing equations of the model, how these are discretized and implemented numerically, as well as six test problems to validate our numerical approach. Our chosen test problems, where the electron inertia and ion kinetic effects play the essential role, are: 0) Excitation of parallel eigenmodes to check numerical convergence and stability, 1) parallel (to a background magnetic field) propagating electromagnetic waves, 2) perpendicular propagating electrostatic waves (ion Bernstein modes), 3) ion beam right-hand instability (resonant and non-resonant), 4) ion Landau damping, 5) ion firehose instability, and 6) 2D oblique ion firehose instability. Our results reproduce successfully the predictions of linear and non-linear theory for all these problems, validating our code. All properties of this hybrid code make it ideal to study multi-scale phenomena between electron and ion scales such as collisionless shocks, magnetic reconnection and kinetic plasma turbulence in the dissipation range above the electron scales.

keywords:
Plasma simulation; Hybrid methods; Particle-in-cell method;
journal: Computer Physics Communications\biboptions

numbers

1 Introduction

Fully kinetic simulations, either via Lagrangian Particle-in-Cell (PiC) [1, 2, 3] or Eulerian Vlasov [4, 5, 6] methods, have become the standard tools to study kinetic plasma phenomena. They are essential to model self-consistently the dissipation at small scales and the multiscale processes in a wide variety of scenarios, often related with turbulence. But these simulations codes tend to be computationally very expensive, which hinders their applicability to many practical problems.

On the other hand, fluid models are used to model phenomena either on large scales and long timescales, like Magnetohydrodynamics (MHD) codes, or even on small space and timescales under appropriate conditions, like electron-MHD (EMHD) codes [7, 8]. Even though they are computationally more convenient than kinetic models, they cannot model self-consistently dissipation, particle acceleration, resonant processes or phenomena related with deviations from the thermal equilibrium or quasineutrality.

Hybrid codes bridge the gap between those different approaches and spatial/time scales. A large number of problems require kinetic effects of ions and fluid descriptions of electrons. This is especially true at frequencies lower than the ion cyclotron or spatial scales much larger than the ion skin depth , which are computationally expensive to simulate with fully-kinetic codes. Problems of this class are related to planetary magnetospheres and their associated collisionless shocks, collisionless magnetic reconnection, the acceleration of heavy ions by shocks and reconnection, instabilities driven by ion temperature anisotropy and turbulence in the kinetic dissipation range in the solar wind, among many other ones. In many of these scenarios, the gradient scales of different physical quantities can be very small, which requires to consider an electron fluid with finite mass.

In order to study these problems, we present a hybrid code using kinetic ions modeled via PiC methods, and a massive electron fluid using a modified EMHD model. We call it CHIEF: Code Hybrid with Inertial Electron Fluid. By not considering electron kinetic effects and going to the radiation-free limit, we can bypass the stringent stability conditions of a fully kinetic approach, making it computationally cheaper. Our code follows the standard quasineutral approach used in kinetic hybrid codes, i.e., the electron and ion densities are always equal, eliminating thus the radiation term in the Maxwell’s equation (Ampère’s) associated with the displacement current, light wave propagation (and the corresponding CFL condition of standard fully-kinetic codes) and other phenomena involving charge separation such as Langmuir waves. On the other hand, and different from other hybrid codes with inertial electrons (see, e.g., the review Ref. [9] or the textbook Ref. [10] and references therein), we consider without approximations all the electron inertial terms in the generalized Ohm’s law. Note that the inclusion of electron inertia reduces the numerical instability (at short wavelengths) known for hybrid codes with massless electrons. This instability is due to the increase of the phase speed of whistler waves without bound, since in that limit, where is the frequency and the wave number. All these features make our code ideal to deal with fast evolving phenomena at electron scales, which can break the frozen-in condition of ideal MHD allowing processes like magnetic reconnection.

Note that our main approximation is the choice of constant temperature as equation of state for the electrons. Other equations of state such as the polytropic one, where the electron pressure has an explicit dependence on other variables such as the density, can be straightforwardly implemented. On the other hand, numerical solutions to equations for the evolution of the electron pressure require additional numerical considerations, like adding a numerical viscosity, see Ref. [11]. That is especially true, in particular, when non-gyrotropic contributions are considered in the electron pressure tensor, since it makes necessary to include a heuristic isotropization term to take into account pitch angle scattering due to electron temperature anisotropy driven instabilities (see Ref. [12]).

In order to make the advantages of our code clear, let us briefly discuss the different approaches that previous hybrid codes have taken to implement the relevant electron fluid equations and their coupling with the Maxwell’s equations. Even since the early days of the development of (PiC-)hybrid codes, electron inertia has been considered to deal with fast phenomena close to the electron cyclotron frequency, as demonstrated by an early 1D electrostatic code [13]. A 2.5D electromagnetic hybrid code was developed by Ref. [14], allowing some electron mass effects. Those codes [13, 14] implemented the quasineutral approximation by using the radiation-free (Darwin) limit of the Maxwell’s equations. Correspondingly, they used a Helmholtz decomposition by solving separately the longitudinal and transverse parts of the electromagnetic potentials and current density. Longitudinal/transverse refers to the curl-free/divergence-free part of the corresponding fields, and not to physical directions in space. The electron inertia was considered in the equation for the transverse electron current density. An evolution equation for the electron temperature including terms up to the heat flux was considered. This method involved solving several elliptic equations for the fields (with the form of generalized Poisson’s equations). Hybrid codes based on potentials, however, have been rarely used since then, except 1D cases which use the electrostatic potential. One example is the 1D code developed in Ref. [15], which modeled a highly collisional plasma typical of pressure discharges with a complex chemistry. It consists of fully-kinetic electrons (using Monte-Carlo methods to take into account particle collisions) and several species of heavy ions modeled with a fluid approach.

Much later, Ref. [16] developed a (2.5D) hybrid code algorithm to simulate the entire Earth’s magnetosphere. To do so, it included a cold ion fluid component to take into account the ionospheric plasma, as well as curvilinear coordinates. In this approach, the electric field is determined from the electron momentum equation with electron inertia, while the electron velocity from the Ampère’s law and the magnetic field from the Faraday’s law. No electron pressure term was considered in the electron momentum equation. The electron inertia term was only considered as a correction in the form of an electron polarization drift in the (implicit) equation for the magnetic field, while neglected in all the other equations. This has the advantage, however, of making the entire code explicit. This code was applied to the substorm onset and other magnetospheric scenarios in Ref. [17].

After that, the standard approach in inertial hybrid codes has been to solve equations for generalized electromagnetic fields (see Sec. 5.7 in Ref. [10], Refs. [11, 18]). In this approach, the electromagnetic fields are obtained from a generalized magnetic and electric fields which satisfy a Faraday’s-like equation. Those 2.5D hybrid codes have been used to study collisionless magnetic reconnection [19, 20, 21]. The equations for the generalized fields and were solved by using a predictor-corrector scheme (which requires a staggered grid) or a trapezoidal leapfrog algorithm. As discussed in detail later, the electromagnetic fields are obtained from the generalized fields under different approximations. Spatio-temporal density variations at electron scales and the time derivative of the ion current were neglected in the expressions for and , respectively [18]. Although the contribution of electron inertia via the convective electron acceleration term was considered in calculating the electric field , the time-derivative of the electron inertia was neglected [18]. Some authors neglected even the convective electron acceleration term [11]. On the other hand, Ref. [11] used an evolution equation for a scalar electron pressure, while Ref. [18] included the full electron pressure tensor to take into account the non-gyrotropic effects that were shown to have an important role in balancing the reconnection electric field under some conditions. None of the effects related with an equation for the evolution of the electron pressure is considered in our code, yet.

A slightly different method to obtain the electric field from consists in writing implicitly the equation for the generalized electric field (instead of the magnetic field), combining the Ohm’s law with the other Maxwell’s equations, resulting in an elliptic equation for (see Sec. 5.2.4 in Ref. [10]). They pointed out that one of the main advantages of this approach is the possibility to simulate regions with very small density, something that the previous methods cannot handle easily due to the division by this quantity in several of the Ohm’s law terms. This kind of hybrid codes was applied mainly for simulations of collisionless shocks.

In hybrid codes, ions can be modeled not only via PiC methods, but also via Eulerian Vlasov methods. Ref. [22] developed a (2.5D) hybrid Vlasov code incorporating electron inertia. This code couples the ion equation of motion represented by the Vlasov equation, solved by standard splitting methods, with the Maxwell’s equations via the current advance method (CAM) [23]. The CAM method uses the ion momentum equation (obtained by calculating the first-order momenta of the Vlasov equation) including a tensor pressure, in order to advance the ion current and thus to calculate the electric field. The corresponding generalized Ohm’s law is a Helmholtz equation for the generalized electric field written in terms of the ion fluid velocity and the ion tensor pressure. It does not involve time derivatives, because of the quasineutral approximation and Faraday’s law. The solution of this equation is simplified by assuming negligible density perturbations in the expression for the electric field. The equation of state for the electrons is assumed to be isothermal. The code by Ref. [22] has been applied to several problems related to solar wind turbulence (see, e.g., Refs. [24, 25, 26]).

More recently, Ref. [27] developed a (3D) hybrid code with electron inertia and kinetic ions modeled via the method, by means of a second-order accurate semi-implicit algorithm using a current closure scheme. The method consists in evolving only the variations of the ion distribution function, which are assumed to be small compared to a given background, and thus to reduce the numerical noise predominant in conventional full- methods. The electron inertial effects are included in the generalized Ohm’s law for the electric field, and are calculated by using the Ampère’s law and the ion momentum equation. This approach was first proposed by Ref. [28] for the 1D gyrofluid electron equations, being later applied to study the propagation of dispersive Alfvén waves in the Earth’s magnetosphere-ionosphere region and the associated mechanisms of electron acceleration [29, 30], and the Alfvén dynamics in an Io-Jupiter flux tube [31]. The electric field is advanced first by decomposing the Ohm’s law into a perturbed and an equilibrium part. The fields are solved with finite differences in the inhomogeneous direction and with spectral methods in the homogeneous direction. The equilibrium part of the Ohm’s law is cast in the form of a linear equation and solved by matrix inversion methods, while the perturbed part is solved iteratively. After that, the magnetic field is obtained explicitly by using Faraday’s law and finally, the electromagnetic fields can be used to advance the (perturbed) particle distribution. The closure for the electron pressure was chosen to be an isothermal equation of state. This code was benchmarked against the (2D) resistive tearing instability, among other (1D) problems.

Finally, one of the more recent hybrid codes with electron inertia was published in Ref. [32]. The special feature of this 1D code is that can handle low density and vacuum regions without becoming unstable. In order to do that, a correction term is introduced to the electric instead of the magnetic field, in addition to a variable (ion-to-electron) mass ratio which reduces the phase speed of whistler waves to satisfy the CFL stability condition everywhere. This is done not only because of the divisions by density in the equation for the electric field, but also because the Alfvén speed is inversely proportional to (the square root of) the density, implying a more stringent condition in the choice of the time step and the relevant CFL condition (associated to the electron Alfvén speed). The approach by Ref. [32] consisted in not using the generalized electromagnetic fields as in Refs. [11, 18, 10], but instead solving a modified form of the Ohm’s law for the electric field (without divisions by density) and then obtaining the magnetic field through the Faraday’s law. The first equation is solved by matrix methods while the second one by a direct iterative algorithm. The electron inertia effects are considered in both equations, however, with some approximations. In the vacuum limit, and following the idea of a previous inertia-less hybrid code [33], the equation for the electric field becomes the Laplace’s equation, being easily solvable without numerical singularities as in other approaches due to the division by density. The electron pressure is considered to follow an adiabatic equation of state, with a corresponding evolution equation in case of a finite resistivity (constant otherwise). Finally, the variable mass ratio procedure is done in practice by varying locally and temporally the electron mass, adjusting it to always satisfy the CFL condition based on the electron Alfvén speed for given timestep and grid cell sizes. This works well as long as the scales of interest are not too close to the electron inertial scales.

The rest of this paper is organized as follows. In Sec. 2 we present a detailed description of the simulation model used by our code. In Sec 3 we describe the numerical implementation of our simulation model. In Sec. 4 we show a detailed comparison of simulations results of one numerical and six physical test problems with the predictions of linear (and non-linear) theory: 0) excitation of parallel eigenmodes (to check numerical convergence and stability), 1) parallel (to a background magnetic field) propagating electromagnetic waves, 2) perpendicular propagating electrostatic waves (ion Bernstein modes), 3) ion beam right-hand instability (resonant and non-resonant), 4) ion Landau damping, 5) ion firehose instability and 6) 2D oblique ion firehose instability. These test problems cover different aspects of ion kinetic and electron inertial effects, essential to understand phenomena in a wide variety of laboratory, space and astrophysical plasmas. All the results match within the expected error range. We summarize our results in Sec. 5.

2 Simulation model

We treat ions as Lagrangian macro-particles modeled via the PiC method and electrons as a (Eulerian) fluid with finite mass and temperature. Each ion macro-particle (index “”) represents a (large) number of physical ions (index “”) according to , where , are the macro-ions’ position, velocity and is called the shape function (see details in Sec. 3). The (physical) ion distribution function is obtained as the summation over all ion macro-particles . In the collisionless plasmas to be considered, this macro-ion distribution function evolves according to a Vlasov-like equation governing the full ion distribution function (see, e.g., Ref. [34]):

(1)

Here, and are the charge (with the fundamental charge) and mass of a physical ion, respectively. By taking the first order momenta in and of Eq. (1) (i.e., multiplying and integrating the Vlasov equation by and , respectively. See details in, e.g., Ref. [34]), one obtains the equations of motion of the macro-ions (setting , i.e., singly charged ions, for simplicity):

(2)
(3)

where and are the first order momenta in and of , respectively. These equations resemble Newton’s equation of motion for the macro-particles, with the difference that and are calculated integrating the full electromagnetic fields over the shape function (see Eqs. (20)-(21)). The code actually solves the relativistic version of these equations, obtained by replacing with the four-velocity , with the relativistic gamma factor. But relativistic effects can be safely neglected for the regime appropriate for a quasineutral hybrid code (, with the speed of light). The ion density , ion bulk velocity and the ion current density are obtained as the first two order velocity momenta of the ion distribution function , i.e., and (see details of the implementation in Sec. 3). With a similar procedure it is possible to obtain higher order momenta such as the ion pressure tensor (associated with the ion temperature).

For fluid electrons, we consider the electron momentum equation with finite electron inertia and resistive effects, i.e., a generalized Ohm’s law for the electric field:

(4)

Here, is the collisional resistivity with the collision frequency (between electrons and ions), is the total current density, is the electron fluid velocity and is the component of the electron pressure tensor. Note that is a summation index, in such a way that represents the component.

The electric and magnetic fields are related to the plasma current via the Faraday’s law and Ampère’s law without displacement current:

(5)
(6)

Note that neglecting the displacement current in Eq. (6) is equivalent to set the dielectric constant , preventing light (and Langmuir) wave propagation and assuming quasi-neutrality for this parameter regime, with the electron/ion density, respectively. This condition implies the restriction to frequencies much smaller than the electron plasma frequency and length-scales much larger than the Debye length . The Poisson’s equation is automatically satisfied because of the quasi-neutral approximation and boundary conditions, while is also fulfilled if the initial conditions do so [9].

Combining the electron momentum Eq. (4) and Faraday’s law Eq. (5), we can eliminate the electric field and obtain an evolution equation for the magnetic field, which has the form of a generalized continuity equation,

(7)

where

(8)

is the generalized vorticity. For simplicity, we have assumed a scalar electron pressure given by the isothermal equation of state,

(9)

with the Boltzmann constant and the electron temperature, constant in time. Note that if the temperature and density are spatially uniform, the second term in Eq. (7) vanishes.

In conventional hybrid codes with electron inertia [11, 18, 10], the electric field is obtained from Eq. (4) while the magnetic field is obtained from Eqs. (7) and (8). This is done, however, by neglecting some terms. Substituting for from Eq. (6) and neglecting terms proportional to and and electron-scale spatial variations of the density, the l.h.s. (left hand side) of Eq. (7) can be written entirely in terms of the time derivative of the magnetic field, i.e., as , where is the electron inertial length. The magnetic field is then obtained by solving an elliptic equation for . Neglecting and in Eq. (7) is justified for a large mass ratio . However, these approximations are not necessarily valid in all the situations of interest and need to be checked for each case, especially for simulations using an artificially low mass ratio . The electric field is calculated directly from the electron momentum Eq. (4) by neglecting , which is clearly inconsistent with keeping this term in Eq. (7).

In our hybrid simulation model, we solve Eqs. (2)-(9) without making any of these approximations. Eq. (7) is solved for the generalized vorticity and then, the magnetic field is calculated from an elliptic equation obtained by combining Eq. (8) with Ampère’s law Eq. (6),

(10)

Note that in Eq. (10) we do not neglect the density variations at electron scales, and explicitly calculate the curl of the ion bulk velocity. Next, the electric field is calculated from the electron momentum Eq. (4) (in the form of a generalized Ohm’s law) by explicitly evaluating the time derivative term . For this purpose, we assume that ions quantities do not change much in a single time step, which is a small fraction of the electron gyroperiod. Thus, for a single time step () of the simulation, is obtained at and by advancing Eq. (7) in two steps of length from to . This allows the calculation of at by a central finite difference scheme. In the second of the two time steps, the code reuses the same ion quantities and used in the first time step (with index ).

2.1 Importance of the full electron inertial terms compared to other hybrid algorithms

In general, any hybrid code with electron inertia should be considered in scenarios with length scales on the order of the electron skin depth and frequencies on the order of . But in order to demonstrate the importance of the additional electron inertial terms of the generalized Ohm’s law, let us combine the relevant equations of our hybrid plasma model (Eqs. (3) and (4)-(6)), in a way similar to the typically used in other hybrid codes, without electron pressure and without resistivity for simplicity:

(11)
(12)
(13)
(14)

where in the last equation we have dropped the indices related to the coarse-graining of particles and electromagnetic fields by the particle shape functions (index “p”). Note that certain terms in the set of the equations (11)-(14) have been multiplied by parameters () whose values can be chosen to be either unity or zero. The full set of inertial terms equations (used in our code) can be obtained by setting (). The first and third terms in Eq. (12) are combined to give . The second term in Eq. (12) comes from one of the terms resulting from the expansion of the fourth term in the Ohm’s law (Eq. (4)). Since this term has been combined with the term in Eq. (12) (first parenthesis), we have simultaneously subtracted with the second term in the square brackets, so that for we recover the Ohm’s law in its usual form (Eq. (4)). We compare explicitly our equations with the approximated set of hybrid plasma equations simulated by Ref. [18], which can be obtained by setting and , and those simulated by [11] by setting ().

Now, let us estimate the orders of magnitude of each of the commonly neglected terms related to the electron inertia in Eqs. (11)-(14).

In order to estimate the importance of the terms with , it is necessary to expand the second term in square brackets by using vector identities and the Ampère’s law. In Ref. [18], the full terms with are retained, while in Ref. [11] only part of the aforementioned expanded second term is retained, neglecting in particular . Note that , as part of the electric field is always canceled after substituting it in the curl of the Faraday’s law Eq. (11). But it has to be considered in the ion equation of motion Eq. (14), which is not done in Ref. [11]. Then, approximating the scale of variation of the electron fluid velocity as , we can compare with the combined form of the first and third term of Eq. (12), ,

(15)

Therefore, the full term with becomes important whenever the length scale variation of the fluid velocity (shear) is on the order of magnitude of the electron skin depth and/or the electron fluid velocity is on the order of the electron Alfvén speed (which is very common in magnetic reconnection, even at ion scales).

The term with in Eq. (12) can be compared with the first term in the same equation. This yields:

(16)

i.e., it is not correct to neglect this term whenever the length scale of the magnetic fields is similar to the electron skin depth: . This term was considered by Ref. [18] but not in Ref. [11].

The term with in Eq. (12) can be also compared with the combined form of the first and third terms of the same equation, approximating , a typical fluctuation frequency:

(17)

where we have assumed that the ratio between scales (at least) as the inverse of the square root of the mass ratio, as in many typical plasma speeds (e.g: Alfvén speed and the electron Alfvén speed). This term, in general, will be small, unless the frequencies are greater than the electron cyclotron frequency: and the mass ratio is small, as is often the case in many plasma simulations where an artificial mass ratio is used for computational convenience. Note that this term was ignored in both Refs. [18, 11].

The term with in Eq. (13) can be compared with the second term in that equation. Since , the gradient , with the typical length scale of variation of the density. Approximating , the scale length of variation of magnetic field, the estimation yields:

(18)

This term is relevant wherever the magnetic length scale gradient is similar or greater than the density gradient: (and both terms are, of course, on the order of ). Such a case might take place during collisionless guide field reconnection, in the cavities of the low density separatrix [35]. Note that this term was ignored in both Refs. [18, 11].

Finally, the term with in Eq. (14) can be compared with the last term in the same equation. The ratio between those two terms yields:

(19)

The first factor is usually negligibly small, but the last one can be very large. So, in case the multiplication of both terms leads to a factor of order 1, the term with should not be neglected in case of steep gradients of magnetic field, i.e.: . Note that is the microscopic velocity which can be considered as part of the ion fluid velocity when the ions are cold and the first-order fluid momenta are taken from Eq. (14). Note that this term was ignored in both Refs. [11, 18].

3 Numerical implementation

The simulation model discussed in Sec. 2 is numerically implemented by combining the Particle-in-Cell (PiC) code ACRONYM [36] 222 http://plasma.nerd2nerd.org/ , with an electron-magnetohydrodynamics (EMHD) code [8, 37]. Both codes have been tested independently and used in the past for simulations of magnetic reconnection [38, 39], instabilities [40, 41], particle acceleration via instabilities [42], CME-driven shocks [43], wave coupling [44], resonant wave-particle interaction [45], etc. Fig. 1 shows a scheme with the numerical implementation and coupling of the two codes. This choice of algorithm was the fastest and safer way (less prone to numerical errors) to couple both codes without further modification and taking advantage of their good properties. But additional schemes and algorithms to solve those equations can be attempted in the future to check whether if they can provide an advantage to the existing ones.

Figure 1: Scheme of the equations of our hybrid simulation code developed by coupling a PiC and an EMHD code. “Sx” stands for step number x (explained in the text).

On the PiC side of the code, the electromagnetic fields are assumed to be known at the staggered grid called the Yee lattice [46]. They are, however, provided by the EMHD side of the code on a grid where all quantities are defined at cell centers except the magnetic field which is defined at the edges of each cell. Therefore, we employ an interpolation procedure between the two grids schematized in Fig. 2, in order to align correctly those quantities in the Yee lattice. Once the electromagnetic fields are interpolated from the EMHD to the PiC cell, the following steps (shown in Fig. 1) are taken to move the ions:

Figure 2: Scheme of the allocations of different physical quantities in the grids associated to the PiC (Yee lattice) and EMHD parts of the hybrid code. For simplicity, only the 2D case is shown. The indices of the grid points are indicated in the vertices of the cell ( for the direction and for the direction).

3.1 Particle mover (PiC part)

The PiC side of the hybrid code, as done in all fully-kinetic PiC codes, including ACRONYM, updates ion velocities and positions in the given electric and magnetic fields and known at the time index .

  • Step 1: The electromagnetic fields are interpolated from the grid to the macro-particle positions with a weighting given by the shape function :

    (20)
    (21)

    where the subscript indicates the location of each macro-particle. In the code, the interpolations between the EMHD and PiC grid and from the PiC grid to macro-particles’ position are combined in a single operation. The deposition of the particle quantities density and current density is not done to the grid points as in the usual Yee lattice, but to the grid centers as required by the EMHD solver. This is possible by centering the particle shape functions around the grid center, and considering that there is no need to calculate these quantities in the PiC cell, since they are used as input to advance the electromagnetic fields in the EMHD side of the code, which requires only quantities living in the EMHD cell. This combination of interpolations has the advantage of less calculations that could introduce numerical errors otherwise, in addition to use the same identical scheme as the Step 3 below by construction.

  • Step 2: Then, the ions can be moved via the second-order accurate leap-frog algorithm, which means that at each timestep , the ion velocities are advanced from the half-timestep to the half-timestep and the ion positions from the timestep to the timestep by using the discretized version of Eqs. (2) and (3) (see, e.g., Sec. 4.3 of Ref. [1]):

    (22)
    (23)

    Note that the first equation is explicit while the second implicit. We advance the ion velocity via the Boris method [47], which involves a rotation of (i.e., its magnitude is kept constant) making Eq. (23) also explicit.

  • Step 3: Because of the staggering of the position and velocity updates, the sources of the electromagnetic field, the ion number density and current density , are computed at the time index and , respectively. This deposition is done via an interpolation scheme using the same shape function as for the electromagnetic field interpolation from the grid to the macro-particles position (Eqs. (20)-(21)):

    (24)
    (25)

    The code has available, at compilation time, an option to smooth both quantities by performing a binomial filter and, if required, a compensation pass [48]. Note that, different from ACRONYM and other state-of-the-art fully-kinetic PiC codes, we do not need to use the charge-preserving Esirkepov scheme [49] to deposit the current onto the grid, due to the quasineutrality condition of the hybrid approach333 Other hybrid codes have used this method due to other reasons. For example, the hybrid code by Ref. [17] implemented a charge preserving scheme [50] because it was computationally convenient for his curvilinear coordinate system. . On the other hand, the reason to use the same shape function in Eqs. (24) and (20) is to avoid self-forces and preserve the global momentum (see, e.g., Secs. 8.5-6 in Ref. [1], Sec 5.3.3 in Ref [2] or Ref. [51] for fully-kinetic PiC codes, or Sec. 4.5.2 in Ref. [10] for hybrid-PiC codes). Our hybrid code has available, as inherited from ACRONYM, the following shape functions: NGP (Nearest Grid Point), CIC (Cloud in Cell), TSC (Triangular Shaped Cloud), PQS (Piecewise Quadratic Spline) and other ones, useful under some specific conditions. By default, we use the second order shape function TSC to reduce the intrinsic PiC shot noise and the associated numerical heating without affecting the computational performance too much [52, 53, 41]. Once the sources of the electromagnetic fields in Eqs. (24)-(25) are known at the Yee lattice, they are passed and interpolated to the EMHD grid following the inverse procedure as in Fig. 2 (similarly as Step 1, the operation is actually combined with the deposition in Eqs. (24)-(25)). Then, the EMHD part of the code updates the electromagnetic fields by using the electron fluid equations and Maxwell’s equations described later. The updated electric and magnetic fields are then passed back to the PiC side and the full cycle is repeated.

3.2 Field updater (EMHD part)

After step 3 in Fig. 1, just before entering into the EMHD part of the code, the input quantities are known at the following time indices: , , and . The ion density is also calculated at the same index, , by taking the average value of and .

  • Step 4: The quasineutrality condition is applied to get the electron density .

  • Step 5: With and a specified value of the electron temperature (which is constant and does not evolve in time), the code computes the electron pressure via the equation of state, Eq. (9).

  • Step 6 (first half-loop): In order to update the electric and magnetic fields, the generalized vorticity in Eq. (7) is advanced from the timestep to by using the flux-corrected transport algorithm [54] of the LCPFCT package 444 http://www.nrl.navy.mil/lcp/LCPFCT , which is specially convenient to resolve steep gradients, one of the main applications of our code. The full update to the next timestep () requires and , which are not known yet. Because of this, the field solver first advances the generalized vorticity by half a timestep (with the input quantities , , , , ) to estimate , , and , by solving the discretized version of Eq. (7):

    (26)

    A multi-dimensional time operator splitting technique (see, e.g., Sec. 5.3.1 of Ref. [10] or Sec. 20.3.3 of Ref. [55]) is used to solve this equation by decomposing it into six sub-equations, grouped in three sets (one for each direction of integration , , ) composed of terms containing the spatial derivative along the integration direction and a time derivative of the generalized vorticity. Then, step 6 to step 8 are performed for each direction of integration.

  • Step 7 (first half-loop): Next, the solution for of the discretized version of the elliptic Eq. (10) is obtained,

    (27)

    Its solution allows to obtain from . This elliptic partial differential equation (PDE) is solved by the subroutines provided by the MUDPACK 5.0 libraries 555 https://www2.cisl.ucar.edu/resources/legacy/mudpack , which use the multigrid iteration method [56, 57]. In particular, the form of the equation is appropriate to be solved via the “mud2/3” subroutine, which is a second-order difference approximation for non-separable PDEs on a rectangular/cubic domain (for 2D/3D cases, respectively).

  • Step 8 (first half-loop): The next step consists in obtaining the components of from the discretized form of Ampère’s law Eq. (6):

    (28)
  • Step 9: Steps 6-8, second half-loop: At this point of time, the code knows all relevant quantities at the time index , i.e., , and . Then, all these quantities are advanced to the time index , by repeating the steps 6 to 8 in Fig. 1. Explicitly, the equations to be solved in this second half-loop are:

    (29)
    (30)
    (31)

    With this, we have the quantities , and . Note that we have used the ion quantities ( and ) at the time index , same as for the first half-loop. This is justified by the heavy mass of the ions compared to the electrons.

    The two half-loops described before (steps 6 to 8) are repeated for the integration along each one of the spatial directions and (in the 3D case), as required by the multi-dimensional time operator splitting technique in step 6, thus advancing in time all the components of the quantities used by field solver.

  • Steps 10: Steps 6-9, full-loop, second time: Before calculating the electric field, we repeat the above procedure (Steps 6-9) by advancing the equations from time step to but also using the ion quantities at the time index (i.e., assuming that they are temporally fixed). This gives us , and .

  • Step 11: Now, the electric field at the time step is calculated from the discretized form of the generalized Ohm’s law Eq. (4):

    (32)

    We use the value of provided by the equation of state Eq. (9) in the same way as before (steps 4 and 5 in Fig. 1). Note that we have assumed and .

    Finally, the values of and are passed to the PiC side of the hybrid code, which now provides and . This ends a timestep and the cycle is repeated for the following one.

The algorithm used in the code is second order accurate in space, as a result of the previous discretization of the equations. A measurement of this error, with ion effects included, is hindered because the PiC shot noise affects any measure of error below the noise level, and even the initial conditions will show some fluctuations at this level. But without ion effects (only fluid quantities) a comparison of the error with analytical solutions can be performed with higher accuracy showing the second order convergence in space. This is proven numerically in the first of our test problems (see Sec. 4.1).

3.3 Boundary conditions

The EMHD Maxwell solver of the code relies on specific libraries to solve some of the equations. Both of them, LCPFCT and MUDPACK, allows to specify periodic boundary conditions or a given value or its derivative to the input variables. That easiness is the reason because periodic and reflecting/PEC (perfect electric conductor) are currently implemented in the code. Other more complicated boundary conditions such as absorbing or open (already implemented in the ACRONYM PiC-code) would require a modification of the EMHD Maxwell solver.

3.4 Numerical stability

The time advancing part of our hybrid code is solved by the flux corrected transport method, and mostly their properties control the stability of the overall algorithm. A local von Neumann stability analysis of this scheme results in an expression for the amplification factor , whose specific form varies for the specific coefficients of diffusion and anti-diffusion used in the algorithm (see Refs. [58, 59, 60]). The stability condition (strictly for , but it can also be applied to due to the second order discretization in space of the elliptic equation (10)) can then be calculated as , and most of the different versions lead to the same local CFL condition: , where is a typical (electron) speed in the system. As explained below, and based on the maximum speed of the whistler wave for this system, can be taken in many cases as the electron Alfvén speed. Of course, some specific simulated plasmas can develop high electron fluid velocities, in which case would correspond to the maximum of those electron velocities in the system. This (linear and local) requirement is, however, not strict: the flux corrected transport algorithm can be stable with larger timesteps, with an upper bound that is usually found empirically. This is due to the non-linear flux correction (suppressing spurious short-wavelength maxima/minima caused by numerical instabilities) and some non-local properties of the algorithm (since it uses a stencil of more than 7 points in space). For a specific hydrodynamic case, this scheme was found (by means of numerical tests) to be stable for CFL numbers  [58].

A basic plasma wave mode present in all magnetized plasmas is the whistler wave. That is why it is convenient to use it as a basis for the calculation of a CFL condition and choice of the time step. Without ion effects and assuming a constant density, the system of equations (4)-(6) (in , and ) simplifies significantly. Linearizing around an equilibrium with uniform magnetic field, the resulting whistler dispersion relation becomes [61],

(33)

where is the electron skin depth. From here it can be clearly seen that in the limit without electron inertia (as in most of the hybrid codes), the denominator in the r.h.s. becomes 1 and the frequency increases quadratically with the wavenumber , implying an increasing phase speed without bounds for short enough wavelengths. Therefore, those codes will become unstable unless a very small time step is chosen, or other (numerical) diffusive term is added in their algorithm. For a finite , the frequency of the whistler branch increases in a similar way proportional to for relatively small frequencies . But for higher frequencies, reaches an asymptote at with a zero phase speed. This makes the algorithm automatically stable against whistler wave propagation, as long as the maximum wave speed ( reached at ) is considered by the CFL condition mentioned above. Another important constraint is that the ion gyration should be well resolved, leading to . Finally, it is worth to mention that in cases where all the important phenomena take place uniquely at ion time and spatial scales, electron inertia can be switched off in the code, with the consequent possibility to choose a larger timestep, on the order of the (ion) Alfvén speed . But note that, in this case, all the electron physics and frequencies above are severely modified.

4 Test problems

In the following we describe one numerical and six physical test problems used to validate our numerical method.

Test problem
(0) Excitation of parallel eigenmodes 1
(1) Parallel EM modes 1
(2) Ion Bernstein modes 0.1 1
(3) Ion beam R instability =0.1 10
(4) Ion Landau damping 2.0 [0.1-0.66]
(5) Ion firehose instability 1
(6) 2D oblique ion firehose instability 1
Table 1: Main physical parameters of the tests problems. The ion Landau problem is unmagnetized, but we still use a reference magnetic field for calculating and all the normalizations depending on it. We use two values for in the parallel EM modes, and a range of values for for the ion Landau damping test problem. Since there are different values of for the different components of the ion beam R instability, here we give instead .
Test problem ppc
(1) 1024 102.4 600 512
(2) 1920 63.6 600 64
(3) 1024 256 120 1024
(4) 256 42.24 60 4096
(5) 480 300 2400 2048
(6) 400 256
Table 2: Main numerical parameters of the test problems. is the number of grid points along the resolved direction, is the simulation box length, is the grid resolution, is the time-step, is the total simulation time and ppc are the number of particles per cell. We use grid points in the transverse direction for all 1-D runs, while the grid size for the 2D firehose test problem is explicitly indicated. The ion Landau problem is unmagnetized, but we still use a reference magnetic field for calculating and all the normalizations depending on it. The numerical parameters of the first numerical eigenmodes test problem (0) are not shown because all of them are varied.

4.1 Excitation of parallel eigenmodes

In order to test the accuracy of our algorithm, we carried out several sets of simulations based on the parallel propagating normal modes. In the framework of the cold two-fluid plasma model, the two basic parallel (to a background magnetic field) propagating electromagnetic waves are the L and R modes. Their name stand for their polarization: L/R for left/right-handed circularly polarized waves, sometimes called ion-cyclotron Alfvén and whistler waves, respectively. These waves have electric and magnetic field fluctuations perpendicular to their propagation direction . Their dispersion relations are given by (see, e.g., Eqs. 6.49-50 of Ref. [62] or Eq. 2.5 of Ref. [63]):

(34)

where the signs in the denominator of the first () and second term () of the right hand side correspond to L/R modes, respectively. Here, are the electron/ion plasma frequencies and are the electron/ion cyclotron frequencies. In order to compare our results with the models used in the classical hybrid codes with massless electrons, we use the approximation of Eq. (34) for low frequencies and :

(35)

where the signs correspond to the R/L modes (see also the Eq. 2.45 in Ref. [64]). Note that we have used the relations and , where is the Alfvén speed and is the ion skin depth. For higher frequencies around , the R mode is usually called whistler wave, with a dispersion relation given by Eq. 33.

The initial setup is the following for a parallel propagating wave along the direction with a background and constant magnetic field . A perturbation of the magnetic field in the form:

(36)
(37)

(for , and so only is initially chosen, not which is selected by the simulated plasma), where the signs are chosen to take into account the polarization of the left (L) and right (R, whistler) handed waves, respectively (and thus to have in each case forward propagating waves). The perturbation strength is chosen as , in order to be substantially above the numerical noise level, but low enough in order to avoid non-linear effects such as harmonics and parametric decay into other waves [65]. The density is kept constant, while other physical parameters are summarized in Tables 1. Boundary conditions are periodic and the tests are quasi-1D (mostly along one direction, averaging over 4 cells in the transverse direction).

The tests use a different set of numerical parameters in order to isolate the individual effects to be analyzed, classified as follows:

4.1.1 Convergence on of a whistler wave: Accuracy of the EMHD solver.

In order to isolate the effects of the EMHD solver from ion kinetic effects, we ran these tests without ion response, i.e., with density constant and the ion current always zero. Both ions and electrons are also cold .

We fit four wavelengths of a whistler wave (R mode) with in the simulation box (the wavelength is kept fixed and equal to ). We resolve each wavelength by a variable number of grid points in the range , with a corresponding grid cell size (). The timestep is chosen to fulfill the CFL condition for the electron Alfvén speed at the level , and so is variable between . Note that this restrictive condition is strictly not needed for this kind of relatively low frequency wave, since their phase speed does not approach to , but it is only chosen to keep consistency with later investigations (which develop high frequency waves).

We let the wave travel for four periods , as calculated from the theoretical dispersion relation (see Eq. (34)). First, we compare the wave shape at this point in time with the analytical solution Eq. 36, which has only one free parameter to adjust, the phase . This is because we are assuming that the wave does not damp out (so is constant), strictly true for all the R branch and the low frequency part () of the L mode (see Fig. 7), and also that the initial wave number is not modified (no mode conversion). The comparison of wave shapes is done with the so-called L1-norm,