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 kilometers, while the global scale of the system is about 10¹³ km. In the heliosphere the scales are altogether smaller (kinetic scales about 2 km, system scale about 10⁸ 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⁻⁴ 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 modeled 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 modeling of Landau damping in the same context were proposed and analyzed in Wang et al. [1], Ng et al. [2], and Allmann-Rahn et al. [3]. These studies 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 Sharma et al. [7]. An excellent overview is given in Chust et al. [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 [10–14]). In the context of plasma physics Sugiyama et al. [15, 16], Markidis et al. [17] and Daldorff et al. [18] show ways to combine PIC and MHD fluid models, Kolobov and Arslanbekov [19] describe the transition from neutral gas models to models of weakly ionized plasmas. In Rieke et al. [20] a coupling between the Vlasov equation and a 5-moment model is presented. In recent simulations by Tóth et al. [21, 22] and Makwana et al. [23] not only the multiphysics but also the multiscale aspects have been addressed by utilizing adaptive mesh refinement (AMR) in the coupling strategy.

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 [20], 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 modeled 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 Wang et al. [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 generalized 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 as follows: 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 [24] and perform comparisons to pure kinetic and pure fluid simulations.

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 summarize these sets of equations.

We apply the described models and the multiphysics coupling strategy to the Geospace Environmental Modeling (GEM) Magnetic Reconnection Challenge [24], which has been set up to study 2d reconnection. It is a Harris sheet equilibrium with magnetic field

density

magnetic flux perturbation

initially uniform electron and ion temperatures of ratio T_i/T_e = 5, and mass ratio m_i/m_e = 25. Moreover, B₀ = 1, n₀ = 1, ψ₀ = 0.1, λ = 0.5 in the standard normalized units. While the original domain size is L_x = 25.6d_i, L_y = 12.8d_i, we use the slightly smaller L_x = 8πd_i, L_y = 4πd_i with ion inertial length d_i, as we have depicted in Table 2. 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=π64di in 2d physical space and a uniform resolution of 32 cells in each direction in 3d velocity space. To reduce the computational costs, we use the reduced speed of light of 20 times the Alfvén speed. The numerical setup is depicted in Table 2.

In the simulation we use the u_z,s-based criterion with a thresholds depicted in Table 3. The criterion is reassessed every 0.1ωc,i-1 (inverse ion-gyrofrequencies) for every subregion.

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. A video showing the dynamical evolution depicted in Figure 3 is available in the Supplementary Material section.

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.

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.

Alone the presence of electron and ion kinetic regions in the very center of reconnection zone and 10-moment fluid regions around it seems to ensure the characteristic behavior of the full kinetic reconnection scenario.

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 [40, 41]. 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 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.

The datasets generated for this study can be found on the data server https://vlasov.tp1.ruhr-uni-bochum.de/data/paper-frontiers-2018.

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.