Multiphysics simulations of collisionless plasmas

Collisionless plasmas, mostly present in astrophysical and space environments, often require a kinetic treatment as given by the Vlasov equation. Unfortunately, the six-dimensional Vlasov equation can only be solved on very small parts of the considered spatial domain. However, in some cases, e.g. magnetic reconnection, it is sufficient to solve the Vlasov equation in a localized domain and solve the remaining domain by appropriate fluid models. In this paper, we describe a hierarchical treatment of collisionless plasmas in the following way. On the finest level of description, the Vlasov equation is solved both for ions and electrons. The next courser description treats electrons with a 10-moment fluid model incorporating a simplified treatment of Landau damping. At the boundary between the electron kinetic and fluid region, the central question is how the fluid moments influence the electron distribution function. On the next coarser level of description the ions are treated by an 10-moment fluid model as well. It may turn out that in some spatial regions far away from the reconnection zone the temperature tensor in the 10-moment description is nearly isotopic. In this case it is even possible to switch to a 5-moment description. This change can be done separately for ions and electrons. To test this multiphysics approach, we apply this full physics-adaptive simulations to the Geospace Environmental Modeling (GEM) challenge of magnetic reconnection.


INTRODUCTION
One of the most important challenges in astrophysical, space and fusion plasmas is the treatment of different spatial and temporal scales and the correct physical description on each of these different scales.
In order to give a rough estimate for different plasma systems, let us first consider the warm ionized phase (diffuse ionized hydrogen) in the interstellar medium. Here, the smallest relevant kinetic scales are in the order of magnitude of kilometres, while the global scale of the system is about 10 13 km. In the heliosphere the scales are altogether smaller (kinetic scales about 2 km, system scale about 10 8 km), but the ratio of global to kinetic scales is still astronomical in the truest sense. The situation is similar in fusion plasmas: the electron skin depth is about 5 · 10 −4 m and the vessel measures about 10 meters. In all these cases, it is not possible to carry out simulations which represent all scales with the finest level (kinetic equations) of the physical description. Most of these plasmas can be considered as collisionless, since collision times are orders of magnitude larger than time scales relevant for the dynamical evolution of the plasma. Such plasmas can be modelled with the kinetic Vlasov equation. Nevertheless, kinetic models are inherently computationally expensive, so that large-scale simulations of typical phenomena, as for example magnetic reconnection or collisionless shocks, are hardly feasible and only possible in localized regions of interest. As an alternative, much cheaper fluid models can be considered, but they lack the expressiveness and some physics of full kinetic models, even though some of the effects may be included. Simple treatments and modelling of Landau damping in the same context were proposed and analyzed in [1,2,3]. These studies Lautenbach et al.

Multiphysics simulations of collisionless plasmas
were based on the closure introduced by Hammett and Perkins [4] and successive work in this direction [5,6]. An extension providing heat fluxes in the parallel and perpendicular directions (with respect to the magnetic field) was presented in [7]. An excellent overview is given in Chust and Belmont [8].
Fortunately, many relevant problems like magnetic reconnection or collisionless shocks exhibit a rather clear separation of scales and regimes such that an adaptive approach is promising and might combine the best of the two worlds: cheap models where they are sufficient and detailed models where they are necessary and interesting. The idea of coupling different physical models is not new and has been applied in different physical contexts. Schulze et al. [9] couple kinetic Monte-Carlo and continuum models in the context of epitaxial growth. Considerable efforts have been made to couple kinetic Boltzmann descriptions with fluid models (see e.g. [10,11,12,13,14]). In the context of plasma physics Sugiyama and Kusano [15], Markidis et al. [16] and Daldorf et al. [17] show ways to combine PIC and MHD fluid models, and Kolobov and Arslanbekov [18] describe the transition from neutral gas models to models of weakly ionized plasmas.
We take a slightly different route in solving the Vlasov equation on the finest relevant scales and then adaptively use less and less detailed fluid models outside the kinetic region. In this way we have some control where to use which kind of physical model at the expense of dealing with a substantially more complicated computational infrastructure.
Our group has developed and is continuously developing and improving methods and codes that are capable of combining kinetic and fluid models during runtime [19], making it possible to consider problems of the type mentioned above at much lower expenses than before.
A sketch of this hierarchy is depicted in figure 1. In the inner zone, both ions and electrons are treated kinetically and solved with the Vlasov equation. Adjacent to this zone, ions are still modelled with the Vlasov equation but electrons are described with a 10-moment fluid model. On the next coarser level of description, the ions are also described by a 10-moment fluid model. To ease the transition from the kinetic to the 10-moment fluid description we apply the Landau closure developed in [1] in the fluid description. It may turn out that in some spatial regions outside the reconnection zone the temperature tensor in the 10-moment description is nearly isotopic. In this case it is even possible to switch to a 5-moment description. This change can be done separately for ions and electron. In future studies we will also try to include the coupling of the 5-moment model to magnetohydrodynamic (MHD) models (with generalised Ohms law) which would represent the last step in this hierarchy.
With this multiphysics strategy, these codes can be applied to problem sizes that are otherwise impossible to reach with kinetic simulations and the understanding of the impact of small scale phenomena on the dynamics on global scales is in reach.
The outline of the paper is the following: first we briefly describe all the plasma models and the necessary numerical schemes (Vlasov equation, 10-and 5-moment fluid equations, Maxwell's equations, the coupling procedure, the Landau fluid closure). We will then study the Geospace Environmental Modeling (GEM) reconnection setup [20] and perform comparisons to pure kinetic and pure fluid simulations.

PLASMA MODELS
The plasma models that we have to consider are: i) the Vlasov equation, ii) Maxwell's equations and iii) the 10-and 5-moment fluid equations. We will briefly summarise these sets of equations.  Figure 1. Oversimplified sketch of a multiphysics approach for tail reconnection

Vlasov equation
Collisionless plasmas on the finest level of description are governed by the Vlasov equation where f s (x, v, t) denotes the phase-space density, q s and m s the particle charge and mass for species s ∈ {e, i} (electrons and ions). The electric and magnetic fields E and B are given by Maxwell's equations: with speed of light c and electric constant ε 0 . Maxwell's equations depend on charge and current densities ρ and j, which are obtained from the phase-space densities f s (x, v, t): Vlasov equation (1) and Maxwell's equations (2) form a closed set of equations and constitute the most fundamental description of a collisionless plasma.

Two-species fluid equations
Fluid descriptions can be obtained from the Vlasov equation (1) by taking moments of the phase-space density f s , Here, v n denotes the n-fold tensor product of v with itself, v 0 := 1. Typically, only the first few moments are considered since a Gaussian distribution f s (x, v, t) is exactly represented by the moments µ 0,s , µ 1,s , µ 2,s (and all other moments equaling zero).
We will subsequently the describe the 10-and 5-moment equations. Consider the lowest moments up to µ 3,s : particle density: energy density tensor: heat flux tensor: n s ,û s , E s are evolved by the following equations, obtained from the Vlasov equation (1): Naturally, these equations are not closed. Designing appropriate fluid closures have a long history. An excellent overview is given in Chust and Belmont [8]. In order to mimic kinetic Landau damping effects, several closures have been developed (see [5,6]), all based on the early Hammett and Perkins [4] model.
Wang et al. [1] suggested a heat flux closure which approximates a spectrum of wave numbers by one single wave number k 0 . Following this idea, Allmann-Rahn et al. [3] developed an improved model that is able to correctly describe the kinetic scaling of average reconnection rate (λ/d i ) −0.73 as a function of the distance between the islands' O-points λ and where d i denotes the ion skin depth (see figure 9 in [3]).
In the simulations used in this paper the original k 0 -closure from Wang et al. [1] is used. It approximates the divergence of the heat flux with an expression that forces an anisotropic pressure tensor to a more isotropic one. More precisely, the expression reads with the pressure tensor P s = E s − m s n sûsûs , the scalar pressure p s = 1 3 tr P s and thermal velocity v th .
The parameter k 0 is choosen on the order of the inverse Debye length. Together with (6), this constitutes a closed set of ten fluid equations.
In future simulations it planed to switch to the improved model introduced in Allmann-Rahn et al. [3].
As an alternative to the 10-moment description, an even simpler 5-moment description can be introduced where the energy density tensor and heatflux tensor are replaced by the scalar energy density E s = 1 2 tr E s and vector heat flux Q s . The scalar energy density evolves in time according to: Together with the assumption of adiabaticity, ∇ · Q s ≡ 0 , (6a, 6b, 8) form the set of five moment equations.
Completely analogous to the case of the Vlasov equation, the source terms ρ and j in Maxwell's equations (2) are formed from the particle densities n s (see equation (5a)). Actually, as will be described in section 3.3, only the current density j is needed to propagate Maxwell's equations.

Vlasov equation
In order to circumvent the complexity that could arise from the high dimensionality of the phase space, the Vlasov equation is split into five one-dimensional problems using Strang splitting [21]. These onedimensional advection problems are solved with a third order semi-Lagrangian flux-conservative scheme introduced by Filbet et al. [22]. In order to minimize the error due to the Strang splitting when calculating the backward characteristics needed in the semi-Lagrangian method, the cascade interpolation [23] is combined with the Boris step [24] to form the backsubstitution method introduced in Schmitz and Grauer [25]. Details of this procedure can be found in [26].
The code is fully parallelized using the message passing interface (MPI) [27] where the Vlasov part is solved in parallel on distributed graphics cards using CUDA programming tools [28].

Two-species fluid equations
Both the 10-moment and the 5-moment two-species fluid models are all discretized with the same numerical methods.
For the discretization in space, we use the CWENO scheme introduced by Kurganov and Levy [29], an easy to implement third order finite-volume scheme which is a perfect compromise between sharp shock resolution and high-order approximation in spatially smooth regions. A third order strong-stabilitypreserving Runge-Kutta scheme [30] is employed for the time integration.

Maxwell's equation
The electromagnetic fields are positioned on a staggered Yee grid [31] in order to maintain the divergence free condition for the magnetic field: ∇ · B = 0. Equations (2c, 2d) are evolved through the FDTD method presented in Taflove and Brodwin [32]. Here, only the current density j enters as a source term. Since the speed of light exceeds all other speeds found in the plasma by far, subcycling is used in order to resolve lightwaves while keeping the global timestep as large as possible. In addition, the speed of light is artificially reduced to 20 times the Alfvén speed.

Adaptive Coupling
The coupling strategy is the most important and at the same time the most critical part of the multiphysics simulations. The coupling strategy involves two separate problems: first, providing the correct boundary conditions at interfaces between different physical models and second, designing criteria to decide which model can be used in which part of the computational domain in an adaptive way.
We start with discussing the strategy for obtaining boundary conditions at the interfaces. Providing boundary conditions for the fluid part at the kinetic/fluid interface is rather straightforward: the fluid boundary conditions are obtained by taking the necessary moments of the phase-space densities f s at the interface. Providing boundary conditions for the phase-space densities f s from the fluid description is far less trivial and described in detail in Rieke et al. [19]. In short the procedure can be summarised as follows: we first extrapolate the phase-space density f s to the boundary region. Next, we adjust the extrapolated phase-space density f s such that the moments equal the moments from the 10-moment fluid description. In this way, we only manipulate the phase-space density f s rather "smoothly" with minimal changes and do not force f s to a Gaussian shape. The coupling between the 10-moment and 5-moment fluid regions is done in a very natural way. The boundary conditions for the 5-moment description is simply obtained by calculating the energy scalar from the trace of the energy density tensor. In the other direction, the energy tensor has to be constructed from the energy scalar by assuming a diagonal shape at the 10-/5-moment interface.
The criteria to decide which of the available models shall be used in which subregion of the domain is a highly non-trivial issue. Presently, our strategy is still in a phase of proof of concept and further work has to be invested. For the case of magnetic reconnection, we implemented heuristic criteria based on the current density j z since it is a good indicator for regions of high reconnection. In order to allow for a finer detachment of electrons and ions, we actually use the velocity u z,s = n sûz,s for electrons and ions separately. This is reasonable as the current density is given by j z = s q s u z,s (compare (3b)). It should be stated clearly that a criterion based on the current density is rather heuristic and far from universally applicable approach and further work has to be invested on robust criteria.
In addition, once a criterion is considered satisfactory for the context, thresholds have to be defined that mirror a good trade-off between the need to use a higher-information model for a correct representation and the opportunity to save computational resources with a lower-information model. Up to know we can only state that this is based on educated guesses.

Code performance
The described numerical codes, the adaptive coupling procedures and the parallelisation framework based on space-filling curves [33] is build in our framework called muphy. This framework has been developed over the last 10 years. It is written in C++/CUDA, runs partly on GPUs and partly on CPUs and employs MPI for parallelization.
Scaling runs have been performed on the JURECA supercomputer at the FZ Jülich, Germany [34], on a fully kinetic Whistler-wave setting [17]. Scaling results are excellent as shown in table 1 and figure 2. Note that the number of GPUs was only restricted by the actual configuration of JURECA.

RESULTS
We apply the described models and the multiphysics coupling strategy to the Geospace Environmental Modeling (GEM) challenge [20]. The domain size is chosen as 4π d i in xand as 2π d i in y-direction, where d i denotes ion skin depth. The symmetry properties of the GEM problem make it sufficient to calculate only one quarter of the spatial domain. We use a uniform cell-width of dx = dy = π 128 d i in 2d physical space and a uniform resolution of 32 cells in each direction in 3d velocity space. To reduce the computational costs, we apply the common values for the reduced mass ratio m i me = 25 and reduced speed of light of 20 times the Alfvén speed. The numerical setup is depicted in table 2.  Table 3. Thresholds for the u z -criterion in units of Alfvén-speed v A . They are evaluated for every subregion separately. The five-moment model is used iff neither the kinetic nor the ten-moment thresholds are met kinetic iff . . . ten-moment iff not kinetic and . . .
In figure 3 the fields j z , u z,e and u z,i are shown together with the areas depicting the different physical models for different times of the simulation. From these figures one can deduce that substantial saving in simulation time can be achieved since the performance gain is approximately proportional to the ratio of the computational domain to the area where the Vlasov equation is solved.
In figure 4, a comparison of different simulations are shown and compared to the multiphysics run. Depicted is the current density j z at the time of the highest reconnection rate. The fully kinetic Vlasov simulation agrees rather well with the multiphysics simulation. The overall agreement is substantially better than the results obtained from the 10-and 5-moment simulations. However, differences especially in the precise values of the absolute maxima of the current density, are visible. Whether this is an effect of the model selection criteria has to be tested in further investigations.   Red areas are solved with the kinetic solver, blue areas with the ten-moment and yellow areas with the five-moment fluid solver. Note that the initial models at time t = 0 have been prescribed.  As a more quantitative comparison, the reconnecting flux for the multiphysics simulation is plotted in figure 5 together with purely kinetic and fluid runs. There are a number of things to observe from the plot: While the reconnecting flux of the fluid runs does not saturate within the simulation time, the kinetic and the multiphysics runs saturate. In addition, they both saturate at the same level and thus capture essentially the same small scale physics which is not possible with the fluid models.

Multiphysics simulations of collisionless plasmas
Alone the presence of electron and ion kinetic regions in the very center of reconnection zone and 10moment fluid regions around it seems to ensure the characteristic behaviour of the full kinetic reconnection scenario.

DISCUSSION
We showed that the proposed multiphysics coupling hierarchy can give excellent results even when only a small part of computational domain near the reconnection zone is captured with a kinetic model. However, still many questions and challenges remain and it is clear that the present simulations are only on the level of a proof of concept. Most important is the issue of designing robust physics refinement criteria and their thresholds. First attempts based on the heat flux are under investigation. In addition, the multiphysics coupling strategy should be formulated as an asymptotic preserving scheme [35,36]. The coupling of the 10-and 5-moment models is already in this state when incorporating the effect of Landau damping [1,3]. Presently, we are also reformulating the coupling between the Vlasov and the 10-moment model. For this, we formulate the kinetic description as an adaptive (in time and space) δf method and Lautenbach et al. multi-physics simulation fully kinetic fully two-fluid with ten moments fully two-fluid with five moments Figure 5. Reconnection fluxes for the coupled run and some uncoupled comparison runs with a single solver. Crosses mark the point of highest reconnection rate throughout the respective run ease the transition to the fluid description as an asymptotic preserving scheme. Finally, the multiphysics hierarchy should not stop at the level of the 5-moment fluid description. Work to couple the 5-moment model to MHD is in progress.