1. Introduction
Space and astrophysical plasmas are characterised by the extreme separation between system scales, where large-scale drivers act, and kinetic scales, where wave/particle interaction processes result in particle heating, acceleration and energy dissipation. The interaction between the system and kinetic scales is not trivial, as we highlight in the examples that follow: first, in the terrestrial magnetosphere, magnetotail reconnection (Zweibel & Yamada Reference Zweibel and Yamada2009; Yamada, Kulsrud & Ji Reference Yamada, Kulsrud and Ji2010) is driven by the global-scale Dungey cycle (Dungey Reference Dungey1965), triggered by magnetic field lines opening and reconnecting at the electron scales (Kuznetsova, Hesse & Winske Reference Kuznetsova, Hesse and Winske2001), and results in global reshaping of the magnetospheric field, as well as particle heating and acceleration (Angelopoulos et al. Reference Angelopoulos, Runov, Zhou, Turner, Kiehas, Li and Shinohara2013). Second, in the solar wind, the scales of energy injection and energy dissipation are separated by several decades in wavenumber and frequency (Frisch Reference Frisch1995; Bruno & Carbone Reference Bruno and Carbone2013). Intermediate scales are, however, of extreme relevance, since meso-scale processes lead to the formation of structures, such as elongated reconnecting current sheets, which fundamentally influence how energy flows across the scales (Sundkvist et al. Reference Sundkvist, Retinò, Vaivads and Bale2007). Furthermore, again in the solar wind, further characteristic length scales are introduced by plasma expansion in interplanetary space, which results in the variation in the radial direction of background plasma parameters and, hence, in the triggering of ion and electron instabilities contributing to the constraining of large-scale solar wind properties such as temperature anisotropies, plasma betas and heat flux (Matteini et al. Reference Matteini, Hellinger, Goldstein, Landi, Velli and Neugebauer2013; Štverák et al. Reference Štverák, Trávníček, Maksimovic, Marsch, Fazakerley and Scime2008; Micera et al. Reference Micera, Zhukov, López, Boella, Tenerani, Velli, Lapenta and Innocenti2021; Verscharen et al. Reference Verscharen2022; Micera et al. Reference Micera, Verscharen, Coburn and Innocenti2025).
For this reason, one of the main challenges in plasma modelling is to capture, in a single simulation, the non-trivial interaction between small (kinetic) and large (system) scales. Fully kinetic models, i.e. Vlasov (Cheng & Knorr Reference Cheng and Knorr1976) or particle in cell (PIC) (Hockney & Eastwood Reference Hockney and Eastwood1988) methods, model kinetic-scale processes self-consistently and without approximations to the state of the plasmas by discretising the Vlasov equations for particle evolution. They do so either directly in velocity space or by sampling the distribution function using computational particles each representing a number of point particles inhabiting nearby regions of phase space, respectively. However, their computational costs are often prohibitive for the study of physical systems of interest, whose sizes very often exceed characteristics kinetic scales by several orders of magnitude. Methods can be devised to reduce the computational costs of both Vlasov and PIC simulations, including semi-implicit (Hewett & Langdon Reference Hewett and Langdon1987; Vu & Brackbill Reference Vu and Brackbill1992; Markidis et al. Reference Markidis2010; Innocenti et al. Reference Innocenti, Johnson, Markidis, Amaya, Deca, Olshevsky and Lapenta2017; Lapenta, Gonzalez-Herrero & Boella Reference Lapenta, Gonzalez-Herrero and Boella2017; Kormann & Sonnendrücker Reference Kormann and Sonnendrücker2021) and fully implicit (Chen, Chacón & Barnes Reference Chen, Chacón and Barnes2011; Markidis & Lapenta Reference Markidis and Lapenta2011) discretisations, adaptive grids (Vay et al. Reference Vay2004; Fujimoto & Sydora Reference Fujimoto and Sydora2008; Innocenti et al. Reference Innocenti, Lapenta, Markidis, Beck and Vapirev2013), low-rank tensor compression (Kormann Reference Kormann2015; Einkemmer & Lubich Reference Einkemmer and Lubich2018; Allmann-Rahn et al., Reference Allmann-Rahn, Grauer and Kormann2022a ; Ye & Loureiro Reference Ye and Loureiro2024; Einkemmer et al. Reference Einkemmer, Kormann, Kusch, McClarren and Qiu2025), asymptotic-preserving schemes (Degond et al. Reference Degond, Deluzet, Navoret, Sun and Vignal2010; Ji et al. Reference Ji, Yang, Li, Wu, Jin and Xu2023) and ad hoc techniques addressing specific characteristic of the system such as expansion or contraction (Riquelme et al. Reference Riquelme, Quataert, Sharma and Spitkovsky2012; Ahmadi, Germaschewski & Raeder Reference Ahmadi, Germaschewski and Raeder2017; Innocenti, Tenerani & Velli Reference Innocenti, Tenerani and Velli2019). This notwithstanding, the computational cost of large-scale kinetic simulations requires that reduced models have to be employed.
At the opposite extreme of the computational cost versus physical accuracy, one finds magnetohydrodynamics (MHD). The MHD description makes a number of assumptions on the state of the plasma, including quasi-neutrality, massless electrons (the plasma is considered as a single fluid), only ideal electric field components and scalar pressure (Ledvina, Ma & Kallio Reference Ledvina, Ma and Kallio2008). As long as one is not interested in processes that require violating one of these assumptions, MHD models can be used for relatively cheap large-scale simulations including, e.g. global heliospheric simulations (Gombosi et al. Reference Gombosi, van der, Bart, Ward and Sokolov2018; Poedts et al. Reference Poedts2020), which are too expensive for practically any other model. If one aims for simulations that can tackle larger boxes with respect to fully kinetic simulations, but include more physical details than MHD simulations, two types of possibilities present themselves. On the one hand, MHD and fully kinetic (most often of the PIC flavour) codes can be coupled, under the assumption that the MHD assumptions are broken only in a small fraction of the domain, where the kinetic description is needed. Coupled MHD/PIC codes are used in a number of environments (such as the terrestrial (Lapenta et al. Reference Lapenta, Ashour-Abdalla, Walker and El Alaoui2016; Chen et al. Reference Chen, Tóth, Cassak, Jia, Gombosi, Slavin, Markidis, Peng, Jordanova and Henderson2017) and planetary (Chen et al. Reference Chen, Tóth, Jia, Slavin, Sun, Markidis, Gombosi and Raines2019) magnetospheres and the solar corona (Haahr, Gudiksen & Nordlund Reference Haahr, Gudiksen and Nordlund2025)) where system-scale drivers result in kinetic-scale processes, often magnetic reconnection; see also Makwana, Keppens & Lapenta (Reference Makwana, Keppens and Lapenta2017). On the other hand, intermediate models between fully kinetic and MHD can be used, which require a reduced number of assumptions on the state of the plasma with respect to MHD. Table 1 in Ledvina et al. (Reference Ledvina, Ma and Kallio2008) details the assumptions made in multi-fluid, Hall-MHD and in hybrid models. Multi-fluid multi-moment models describe plasmas as composed of multiple fluids, and evolve a varying number of their velocity moments (density, current and energy in the 5-moment model; density, current and the pressure tensor in the 10-moment model), as we describe in detail in § 2.
The muphyII framework (Lautenbach & Grauer Reference Lautenbach and Grauer2018; Allmann-Rahn et al. Reference Allmann-Rahn, Lautenbach, Deisenhofer and Grauer2024) used in this work implements both a Vlasov and multi-fluid (5 and 10 moment) descriptions, which can be used either stand-alone or in coupled mode. In this regard, muphyII is a peculiarity among plasma modelling frameworks, since one can use both reduced descriptions and code coupling to extend domain size and physical duration of plasma simulations.
In this paper we benchmark 10-moment simulations of plasma turbulence against their fully kinetic, Vlasov counterparts. We focus on the 10-moment model because plasma turbulence is a typical multi-scale plasma process, where one wants to simulate both very large domains (hence, fluid models are preferable to fully kinetic models, due to their lower computational cost), as well as the scales and mechanism of energy dissipation via wave–particle interactions. This transfer process, fundamentally a kinetic phenomenon, plays a crucial role in determining the shape of energy spectra, particularly at scales intermediate between the ion and electron scales. For this reason, using a high number of velocity moments and a good choice of closure of the moment hierarchy are of particular importance. Landau-fluid closures, first introduced by Hammett & Perkins (Reference Hammett and Perkins1990), approximate linear kinetic effects within a fluid representation, particularly for parallel heat transport, by introducing specific relations between higher- and lower-order moments of the distribution function (see also Hunana et al. Reference Hunana, Tenerani, Zank, Goldstein, Webb, Khomenko, Collados, Cally, Adhikari and Velli2019). They are derived in regimes where linear approximation holds, meaning the velocity distribution function is only slightly perturbed with respect to the Maxwellian distribution. However, in this work, we show that this approximation can still be used far from equilibrium, by adjusting closure parameters.
In the local (as opposed to non-local, Ng et al. Reference Ng, Hakim, Bhattacharjee, Stanier and Daughton2017) three-dimensional generalisation of the Hammett–Perkins closure in weakly magnetised plasmas, a free parameter
$k_0$
emerges that represents a characteristic wavenumber related to the maximum damping rate (Wang et al. Reference Wang, Hakim, Bhattacharjee and Germaschewski2015). Allmann-Rahn, Trost & Grauer (Reference Allmann-Rahn, Trost and Grauer2018) and Ng et al. (Reference Ng, Hakim, Wang and Bhattacharjee2020) have shown that 10-moment simulations closed with different flavours of Landau-fluid closures reproduce well fully kinetic magnetic reconnection results. Ten-moment models with Landau-fluid closures have also delivered favourable results in reproducing pressure-gradient- and pressure-anisotropy-driven instabilities, such as the lower hybrid drift instability (Ng et al. Reference Ng, Hakim, Wang and Bhattacharjee2020; Allmann-Rahn et al. Reference Allmann-Rahn, Lautenbach, Grauer and Sydora2021) and the firehose instability (Walters et al. Reference Walters, Klein, Lichko, Juno and TenBarge2024). It is often the case that the free parameter
$k_0$
that delivers the best comparison with fully kinetic simulations is related to characteristic scales of the plasma species (Ng et al. Reference Ng, Hakim, Wang and Bhattacharjee2020).
The goal of our 10-moment turbulence simulations is not to perfectly reproduce the kinetic damping and energy transfer process, which would require non-local closures and possibly a higher number of moments, but rather to identify if a closure parameter exists, that can reproduce sufficiently well the global effects of energy dissipation at kinetic scales, i.e. the slope and shape of the energy spectra. We will do so via a rather brute force approach: after identifying a suitable Landau-fluid closure, i.e. the one described and validated in Allmann-Rahn et al. (Reference Allmann-Rahn, Lautenbach, Grauer and Sydora2021), we will run multiple 10-moment turbulence simulations with varying values of the
$k_0$
parameters. We will then identify as ‘optimal’
$k_0$
parameters, for both ions and electrons, the ones which deliver turbulent spectra which better compare with the Vlasov results.
We proceed in two distinct phases: first, we establish an optimal electron closure parameter by comparing hybrid simulations (kinetic ions, fluid electrons) against fully kinetic reference simulations, focusing on matching energy spectra across the broadest possible range of scales. Second, we apply this optimised electron closure in full fluid simulations while varying the ion closure parameter to determine its optimal value. We compare and comment simulation results across different models (Vlasov, hybrid, two fluid) and with different parameter choices.
This paper is organised as follows. First, in § 2, we describe the so-called closure problem for 10-moment models. In § 3, we describe the three types of problems that we tackle in this paper: Landau damping, Kelvin–Helmholtz instability and, finally, decaying turbulence simulations. While our final aim is decaying turbulence simulations, we first address ‘simpler’ problems such as Landau damping and Kelvin–Helmholtz simulations because they allow us to tackle, in isolation, processes that will occur, together, in turbulence simulations. Using the Landau damping test case, we show that an optimal value of the
$k_0$
parameter can be identified, allowing a fluid description to capture this intrinsically kinetic process. Naturally, this comes in reduced form, leading to discrepancies compared with fully kinetic simulations. With the Kelvin–Helmholtz simulations we demonstrate that 10-moment models can satisfactorily reproduce kinetic results for velocity shear instabilities, which are expected to develop in the decaying turbulence simulation. Parameter optimisation and model comparison for our three test cases of interest are presented in § 4. Discussion and conclusions follow in § 5.
2. Theoretical framework
The behaviour of collisionless plasmas is governed by the evolution of the distribution function
$f_s(\boldsymbol x, \boldsymbol v, t)$
for each species
$s$
, as described by the Vlasov equation
\begin{equation} \partial _t f_s +\boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla} _x f_s + \frac {q_s}{m_s} (\boldsymbol E + \boldsymbol v \times \boldsymbol B) \boldsymbol{\cdot }\boldsymbol{\nabla} _v f_s = 0, \end{equation}
where
$\boldsymbol E$
and
$\boldsymbol B$
are the electric and magnetic field, respectively,
$q_s$
and
$m_s$
are the charge and mass of species
$s$
.
$x$
is space and
$v$
is velocity. This equation couples to Maxwell’s equations, which determine the self-consistent electromagnetic fields
\begin{equation} \partial _t \boldsymbol B = -\boldsymbol{\nabla }\times \boldsymbol E,\quad \partial _t \boldsymbol E = c^2 \left (\boldsymbol{\nabla }\times \boldsymbol B - \mu _0 \boldsymbol J\right )\!,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol B = 0,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol E = \frac {\rho _c}{\varepsilon _0}, \end{equation}
where
$\rho _c$
represents the charge density,
$\boldsymbol J$
the current density and
$c$
,
$\mu _0$
and
$\varepsilon _0$
are the speed of light, vacuum permeability and permittivity, respectively.
The so-called 10-moment model (Hakim Reference Hakim2008; Johnson & Rossmanith Reference Johnson and Rossmanith2010) is a suitable compromise between reduced computational cost and expressivity of the description. The 10-moment model includes the first ten fluid moments which are obtained by integrating over velocity space: the particle density
$n_s = \int f_x(\boldsymbol x, \boldsymbol v, t) \, \text{d} \boldsymbol v$
, mean velocity
$\boldsymbol u_s = \int \boldsymbol v f(\boldsymbol x,\boldsymbol v, t) \, \text{d} \boldsymbol v/n_s$
and the momentum flux density
$\mathcal{P}_{s} = m_s \int \boldsymbol v \otimes \boldsymbol v f_s(\boldsymbol x, \boldsymbol v, t) \, \text{d} \boldsymbol v$
. The equations are then obtained from integration of the Vlasov equation as follows:
\begin{align} \frac {\partial n_{s}}{\partial t} &+ \boldsymbol{\nabla }\boldsymbol{\cdot }(n_{s} \boldsymbol{u}_{s}) = 0, \end{align}
\begin{align} m_{s} \frac {\partial (n_{s} \boldsymbol{u}_{s}) }{\partial t} &- n_{s} q_{s} (\boldsymbol{E} + \boldsymbol{u}_{s} \times \boldsymbol{B}) + \boldsymbol{\nabla }\boldsymbol{\cdot }\mathcal{P}_{s} = 0, \end{align}
\begin{align} \frac {\partial \mathcal{P}_{s}}{\partial t} &- q_{s} \left(n_{s} \text{sym}[\boldsymbol{u}_s \boldsymbol{E}] + \frac {1}{m_{s}} \text{sym}[\mathcal{P}_{s} \times \boldsymbol{B}]\right) + \boldsymbol{\nabla }\boldsymbol{\cdot }\mathcal{Q}_{s} = 0. \end{align}
Here, we denote the generalised tensor product by
$(\mathrm{a}\times \boldsymbol b)_{ijk} = \epsilon ^{klm}a_{ijl}b_m$
, the standard outer product by
$\otimes$
and
$\textrm {sym}[\boldsymbol{\cdot }]$
represents the symmetrisation operator
$\textrm {sym}[\mathrm{a}]_{ijk} = a_{ijk}+a_{jki}+a_{kij}$
. The last equations involves the next higher moment, the third moment
$\mathcal{Q}_{s} = \int \boldsymbol v \otimes \boldsymbol v \otimes \boldsymbol v f(\boldsymbol x, \boldsymbol v, t) \, \text{d} \boldsymbol v$
, which is related to the heat flux
$\mathrm{Q}_{i j k}$
as
\begin{equation} \mathrm{Q}_{s}=\mathcal{Q}_{s}-\operatorname {sym}[\boldsymbol{u} \otimes \mathcal{P}]_{s}+2 m_s n_s \boldsymbol{u}_{s} \otimes \boldsymbol{u}_{s} \otimes \boldsymbol{u}_{s} . \end{equation}
The critical step in developing a practical fluid model lies in the closure of the moment hierarchy.
The simplified 5-moment model is obtained when assuming isotropic pressure and zero heat flux. In this case, the tensorial pressure reduces to the energy density
$\mathcal{E}_s = \frac {m_s}{2} \int |\boldsymbol{v}|^2 f_s(\boldsymbol x , \boldsymbol v, t) \, \text{d}\boldsymbol v$
and the fluid equations are given by
\begin{align} \partial _t n_s &+ \boldsymbol{\nabla }\boldsymbol{\cdot }(n_{s} \boldsymbol u_s) = 0, \end{align}
\begin{align} m_{s} \partial _t (n_s \boldsymbol u_s) &- n_s q_s (\boldsymbol E + \boldsymbol u_s \times \boldsymbol B) + \frac {1}{3}\boldsymbol{\nabla }\left ( 2 \mathcal{E}_s - m_s n_s |\boldsymbol{u}_s|^2 \right ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( m_s n_s \boldsymbol{u}_s \otimes \boldsymbol{u}_s \right )= 0, \end{align}
\begin{align} \partial _t \mathcal{E}_s &- q_s n_s \boldsymbol{u}_s \boldsymbol{\cdot }\boldsymbol{E} + \frac {1}{3} \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \boldsymbol{u}_s (5\mathcal{E}_s - m_sn_s |\boldsymbol{u}_s|^2\right ) = 0. \end{align}
This 5-moment model represents the simplest fluid description and completely neglects kinetic effects. In the case of a 10-moment model, one needs to find a suitable formulation for the divergence of the heat flux tensor in (2.6).
We use here a local (Sharma et al. Reference Sharma, Hammett, Quataert and Stone2006) closure that receives information from the temperature gradient
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }Q_s = -\frac {\chi }{k_{s,0}}\,n_s\,v_{\text{th},s}\,{\nabla} ^2\,\mathrm{T}_s. \end{equation}
With respect to other closures based directly on pressure tensors, the temperature-gradient approach offers better numerical stability during adaptive coupling between kinetic and fluid regions (Allmann-Rahn et al. Reference Allmann-Rahn, Trost and Grauer2018, Reference Allmann-Rahn, Lautenbach and Grauer2022b
), making it particularly suitable for multi-scale simulations with localised kinetic effects. It introduces the parameter
$1/k_{s,0}$
, which has units of length and can be interpreted as a characteristic scale for heat transport processes. The dimensionless factor
$\chi$
is typically set to
$\sqrt {8/\pi }$
based on the linear theory of Landau damping (Hammett & Perkins Reference Hammett and Perkins1990), while
$k_{s,0}$
remains a free parameter that must be calibrated to match the desired kinetic response. Throughout this paper, unless otherwise specified, we normalise
$k_0$
to the skin depth of the respective particle species, i.e. electron inertial length
$d_e = c/\omega _{pe}$
or ion inertial length
$d_i = c/\omega _{pi}$
, where
$\omega _{pe}$
and
$\omega _{pi}$
are the electron and ion plasma frequencies based on the initial uniform plasma density and background magnetic field. The optimal value of
$k_0$
depends on the dominant physical processes in the plasma regime under consideration.
It is interesting to examine the relation of the 10-moment model with the 5-moment and (by extension) MHD limits. Unlike some fluid models that can formally reduce to either under appropriate limiting conditions, our Landau-fluid closure for the 10-moment model cannot inherently recover the 5-moment limit for any choice of the free parameter
$k_{s,0}$
. This occurs because the parameter creates a mathematical tension between two incompatible conditions: as
$k_{s,0}$
decreases (or
$1/k_{s,0}$
increases), the heat flux term in (2.10) becomes stronger, effectively enhancing thermal diffusion and driving the distribution function more aggressively toward an isotropic and thus Maxwellian state. However, the 5-moment model assumes zero heat flux rather than enhanced heat transport. Conversely, as
$k_{s,0}$
increases (or
$1/k_{s,0}$
decreases), the heat flux contribution diminishes, potentially allowing non-isotropic distributions to persist, another deviation from the 5-moment model, which fundamentally assumes local thermodynamic equilibrium. This inherent tension underscores that our closure represents an intermediate regime between full kinetic physics and 5-moment models, capturing important kinetic effects while retaining computational efficiency, but without the ability to formally reduce to either limiting case through parameter adjustment alone.
3. Simulation set-up
We study three test cases, Landau damping, the Kelvin–Helmholtz instability and a decaying turbulence test case. We will study the influence of the choice of
$k_0$
in the 10-moment approximation. The distinct models employed in this study are:
-
(a) fully kinetic Vlasov simulation of both species: high fidelity reference;
-
(b) hybrid model with fully kinetic ions, 10-moment electrons: next model in the hierarchy to separate the influence of the moment model approximation on the two species;
-
(c) 10-moment model for both species;
-
(d) 10-moment ions, 5-moment electrons;
-
(e) 5-moment models for both species (with isotropic adiabatic closure).
Table 1 describes which model is used for the three simulation test cases described in the following sections.
Models employed across different simulation cases. Each cell indicates whether a particular plasma model (column) was used for a specific simulation set-up (row).
The numerical implementation uses the muphyII framework (Allmann-Rahn et al. Reference Allmann-Rahn, Lautenbach, Deisenhofer and Grauer2024) with efficient algorithms chosen for their conservation properties and numerical stability. For the Vlasov equation, we utilise the positive flux-conservative scheme developed by Filbet, Sonnendrücker & Bertrand (Reference Filbet, Sonnendrücker and Bertrand2001), enhanced with the backsubstitution method (Schmitz & Grauer Reference Schmitz and Grauer2006) and a moment-matching algorithm that significantly improves energy conservation (Allmann-Rahn et al. Reference Allmann-Rahn, Lautenbach and Grauer2022b ). The electromagnetic fields evolve according to Maxwell’s equations through the finite difference time domain method on a staggered grid, following the classical approach of Yee (Reference Yee1966).
The fluid models employ a central weighted essentially non-oscillatory reconstruction scheme (Kurganov & Levy Reference Kurganov and Levy2000), combined with third-order Runge–Kutta time integration (Shu & Osher Reference Shu and Osher1988). To manage the difference in temporal scales between the electron and ion dynamics, we implement electron subcycling within each ion time step, improving computational efficiency while maintaining numerical stability.
3.1. Landau damping
To investigate electron Landau damping, we initialise a one-dimensional periodic system with a small-amplitude electron density perturbation. The simulation domain spans
$L_x = 4\pi \lambda _{D,e}$
where
$\lambda _{D,e}$
is the electron Debye length, chosen to accommodate a single wavelength of the perturbation. The initial condition consists of a Maxwellian electron velocity distribution with density modulation
\begin{equation} f_e(x,v,t=0) = n_0\sqrt {\frac {m_e}{2\pi \,T_{0,e}}}\left((1 + \alpha \cos [k_px]\right )\exp \left [-m_e\frac {v^2}{2T_{0,e}}\right ]\!, \end{equation}
where
$\alpha = 0.01$
ensures operation in the linear regime, and
$k_p = 2\pi /L_x$
is the wavenumber of the perturbation. The background ion population remains uniform and stationary, providing a neutralising background. Space is resolved with
$256$
cells, the ion velocity grid spans over
$\pm 6\sqrt {2} \,\lambda _{D,e}\omega _{p,e}$
, with
$\omega _{p,e}$
the electron plasma frequency, and is resolved with
$256$
cells. We use a realistic mass ratio
$m_i/m_e=1836$
and the temperature ratio is set to
$T_i/T_e=1$
.
3.2. Kelvin–Helmholtz instability
The Kelvin–Helmholtz instability is initiated through a velocity shear layer in both ions and electrons in a doubly periodic domain of size
$L_x \times L_y = 18 L_0 \times 27 L_0= 900 \; d_e \, \times \, 135 \, d_e$
, where
$L_0$
is the half-width of the initial shear layer. The shear layer extends along the
$y$
direction and is located at the centre of the box along
$x$
. At the boundary layer, together with the velocity shear, we initialise the gradient of density and out-of-plane magnetic field following hyperbolic tangent profiles, maintaining pressure balance across the shear layer. To trigger the instability we impose a small perturbation on the ion velocity with
$N=10$
modes according to
\begin{align} &\delta \boldsymbol u = \hat {\boldsymbol e}_z \times \boldsymbol{\nabla }\varPsi \quad \text{with}\quad \varPsi = \frac {\epsilon }{N} \exp \left [-\left (\frac {x-x_0}{L_0}\right )^2\right ] \sum _{m=1}^N \frac {1}{m} \cos \left [2\pi m \frac {y}{L_y} + \phi _m\right ]\!,\nonumber \\[-2pt]& \end{align}
where
$\epsilon$
is the perturbation amplitude and
$\phi _m$
the individual phase amplitude. We use
$2400 \times 1600$
cells in
$x$
and
$y$
, respectively, velocity is resolved with
$48^3$
cells and spans
$\pm 4 \,v_S$
for ions and
$ \pm 9.2 \,v_S$
for electrons, with
$v_S$
the shear velocity, the mass ratio between ions and electrons is
$m_i/m_e = 25$
, the ratio between the initial ion and electron temperature is
$T_{0i}/T_{0e} = 5$
and
$\omega _{p,e}/\omega _{c,e} = 4$
, with
$\omega _{c,e}$
being the electron cyclotron frequency.
3.3. Decaying turbulence
To initiate decaying two-dimensional turbulence with low velocity–magnetic field correlation, we adopt a Biskamp–Welter vortex configuration (Biskamp & Welter Reference Biskamp and Welter1989). The initial magnetic and velocity fields are given by
\begin{align} \boldsymbol B &= \frac 13 \delta B\ \hat {\boldsymbol{e}}_z\times \nabla\left (\cos \left [(2x+2.3)\frac {2\pi }{L}\right ]+\cos \left [(y+4.7)\frac {2\pi }{L}\right ]\right ) + B_0\hat {\boldsymbol{e}}_z, \end{align}
\begin{align} \boldsymbol u_s &= \delta u\ \hat {\boldsymbol{e}}_z\times \nabla\left (\cos \left [ (x+1.4)\frac {2\pi }{L}\right ]+\cos \left [(y+0.5)\frac {2\pi }{L}\right ]\right ) + u_{z,s}\hat {\boldsymbol{e}}_z. \end{align}
The out-of-plane electron velocity
$u_{z,e}$
is determined by the MHD current condition
$\boldsymbol J = \boldsymbol{\nabla }\times {\boldsymbol B}/\mu _0$
, while maintaining
$u_{z,i}=0$
for ions. The initial electrostatic field follows from ideal Ohm’s law,
$\boldsymbol E_0 = \boldsymbol u_i \times B_0\hat {\boldsymbol{e}}_z$
, with corresponding electron density perturbation
$n_e = n_0 - \boldsymbol E_0 \varepsilon _0$
.
Our parameter choices follow the benchmark case A1 from Grošelj et al. (Reference Grošelj, Cerri, Navarro, Willmott, Told, Loureiro, Califano and Jenko2017). In particular, we use equal electron and ion temperatures
$T_e/T_i=1$
, a low plasma beta of
$\beta _s=0.1$
, reduced mass ratio
$m_e/m_i = 1/100$
, reduced vacuum speed of light
$c=18.174 v_{\text A}$
and a moderate initial fluctuation amplitude of
$\varepsilon = \delta B/B_0 = \delta u/v_{\text A} = 0.3$
.
The periodic computational box for the Biskamp–Welter turbulence set-up extends to
$L \times L = 8\pi d_i \times 8\pi d_i$
. This domain is resolved with
$1024 \times 1024$
grid points in configuration space, providing sufficient resolution to capture structures down to sub-ion scales while maintaining reasonable computational costs. Velocity-space spans
$\pm 16\,v_A$
for electrons and
$\pm 1.6\,v_A$
for ions in each direction, resolved with
$20^3$
cells. This relatively low resolution is made possible by the moment-fitting procedure described in Allmann-Rahn et al. (Reference Allmann-Rahn, Lautenbach and Grauer2022b
), which allows us to capture the relevant physics accurately even with a coarse grid. See the accompanying discussion of figure 4 in Allmann-Rahn et al. (Reference Allmann-Rahn, Lautenbach and Grauer2022b
), which addresses the velocity-space resolution.
4. Parameter optimisation and model comparison
In this section, we will compare several reduced 10-moment simulations with different values of the
$k_0$
parameter in (2.10) with their fully kinetic Vlasov counterparts. In each case, we will identify the optimal
$k_0$
for that specific problem, and strive to relate it to characteristic plasma parameters, for the sake of generalisation to different simulation set-ups. Each subsection addresses a different aim: § 4.1 examines whether there exists at all an optimal value of
$k_0$
that best characterises the damping behaviour of the heat flux. Then, we proceed to investigate if the closure is able to reproduce the specific plasma dynamics, i.e. plasma mixing processes (§ 4.2). Finally, we tackle our problem of interest, decaying turbulence, § 4.3.
4.1. Landau damping
We begin our investigation with the electron Landau damping, a fundamental kinetic process that represents the archetypal wave–particle interaction mechanism in collisionless plasmas (Nicholson Reference Nicholson1983). This phenomenon presents an ideal starting point for our investigation for several reasons. First, as a linear process with well-established theoretical understanding, it provides a clean benchmark against which we can evaluate fluid closures designed specifically to capture such kinetic effects. Second, the one-dimensional nature of the problem isolates the physics of interest without the complications of multi-dimensional coupling. Third, and most importantly, Landau damping represents the primary mechanism through which electromagnetic fluctuations transfer energy to particles in turbulent collisionless plasmas – precisely the physics our optimised closure needs to capture correctly.
Landau damping in hybrid (Vlasov ion, fluid electron) simulations as a function of
$k_0 \lambda _{D,e}$
. The damping rate
$\gamma$
is plotted in blue on the left axis, the half-period of the electric field oscillations in red on the right axis. The dashed lines at the top and bottom are results from the fully kinetic benchmark simulation. A parabolic fit of the damping rate over the highest values (dotted blue) gives the closest approach to the benchmark value at
$k_0=0.11\,\lambda _{D,e}^{-1}$
with a damping rate of
$0.22\,\omega _{p,e}^{-1}$
.
We conduct a systematic series of electron Landau damping simulations with a fluid treatment of electrons while varying the closure parameter
$k_{0,e}$
. The ions are simulated kinetically.
The resultant damping rate
$\gamma$
and the half-period
$\tau$
of the electric field oscillations (time between peaks in the absolute value of the electric field) are presented in figure 1, on the left (blue) and right (red) axis, respectively, as a function of the
$k_0$
value. Results from the benchmark Vlasov simulations are depicted as a dashed line.
The observed behaviour of the damping rate
$\gamma$
and half-period
$\tau$
as functions of
$k_0$
can be understood by considering the heat flux closure’s physical interpretation. At very small
$k_0$
values, the heat flux becomes very large, causing rapid thermalisation that overwhelms the wave–particle resonance responsible for Landau damping. Conversely, at large
$k_0$
values, the heat flux contribution becomes negligible, essentially decoupling the higher-order moments from the lower-order dynamics and preventing effective energy transfer between waves and particles. The optimal value represents a balance where the closure-induced dissipation most accurately mimics the kinetic resonance condition that drives Landau damping in the fully kinetic description.
Among the various
$k_0$
values tested, we observe that
$k_0 \lambda _{D,e} \approx 0.11$
produces a damping rate that most closely approximates the kinetic value. The optimal value
$k_0 \lambda _{D,e} \approx 0.11$
corresponds to a characteristic heat transport length
$1/k_0 \approx 9 \lambda _{D,e}$
, which is of the order of the wavelength of the damped wave (
$\lambda = 4\pi \lambda _{D,e}$
). The oscillation frequency shows less sensitivity to parameter variation but still exhibits reasonable correspondence near this optimal value. We select this optimal value based primarily on the damping rate, as it directly reflects the efficiency of energy transfer between fields and particles, a critical process for accurately modelling turbulent cascades in more complex scenarios.
It is worth emphasising that the wave–particle energy transfer mechanism examined in this simple Landau damping set-up forms an important building block of dissipation in collisionless plasma turbulence, particularly for the turbulence scenarios we study here where wave–particle resonant interactions play a dominant role. In turbulent plasmas, electromagnetic fluctuations cascade to smaller scales until they can efficiently transfer energy to particles through resonant interactions, precisely the process we are optimising our closure to represent. While the turbulent case involves a spectrum of waves rather than a single mode, and occurs in a multi-dimensional setting with complex nonlinear interactions, the basic physics of field-to-particle energy transfer captured by Landau-type closures remains relevant in these turbulent regimes. Therefore, this one-dimensional calibration provides a crucial foundation for our subsequent optimisation in more complex turbulent scenarios. We remark, however, that wave/particle resonant interactions result in heat transport also in non-turbulent systems (Komarov et al. Reference Komarov, Schekochihin, Churazov and Spitkovsky2018; Yerger et al. Reference Yerger, Kunz, Bott and Spitkovsky2025) and that other mechanisms (e.g. magnetic reconnection in current sheets (Sundkvist et al. Reference Sundkvist, Retinò, Vaivads and Bale2007; Dahlin, Drake & Swisdak Reference Dahlin, Drake and Swisdak2017), diffusive shock acceleration at shock fronts (Bell Reference Bell1978), stochastic heating (McChesney, Stern & Bellan Reference McChesney, Stern and Bellan1987)) transfer energy to particles in different plasma contexts and involve additional physics beyond the wave–particle resonances captured here.
4.2. Kelvin–Helmholtz instability
Having established a preliminary closure parameter from the Landau damping case, we next examine the Kelvin–Helmholtz instability (KHI), which introduces critical multi-dimensional and plasma mixing dynamics absent in the previous test. The KHI represents an intermediate step between the simple linear physics of Landau damping and the full complexity of turbulent cascades, combining large-scale fluid-like behaviour with the emergence of smaller-scale structures where kinetic effects become important.
This test case serves multiple purposes in our closure study. First, it allows us to evaluate how well our Landau damping-optimised closure translates to multi-dimensional settings with more complex energy dissipation processes, possibly operating at multiple scales. Second, the natural development of vortices and secondary instabilities during KHI evolution tests the closure’s ability to handle strongly non-equilibrium scenarios. Third, the KHI set-up provides an excellent framework for comparing our hierarchy of models before proceeding to fully developed turbulence. The KHI has already been used as a test case for model comparison in Henri et al. (Reference Henri2013), where MHD, two-fluid 5-moment and fully kinetic PIC simulations of KHI are compared. It is observed there that the different models compare well in the linear stage of the instability, when vortices develop at larger scales (tens of ion skin depth), while differences emerge in the nonlinear stage, when vortices roll up, and the smaller-scale dynamics, eventually reaching kinetic scales, appear. In this section, we employ all five models described in our methodology, see figure 2.
The quantity we depict in figure 2 is the out-of-plane magnetic field
$B_z$
, normalised to the initial field strength
$B_0$
and at time
$t = 200\,\omega _{p,e}^{-1}$
. We use
$k_{0,e}\,d_e = 100$
and
$k_{0,i}\,d_i = 10$
for the electron and ion closures in the 10-moment models for the respective species. We notice, however, that the specific values of the parameters
$k_{0,e} d_e$
and
$k_{0,i} d_i$
do not affect simulation evolution significantly: the crucial ingredient seems to be the possibility to develop agyrotropies and anisotropies (i.e. evolution of the full pressure tensor), rather than the specific value of the closure parameter, see discussion below.
Value of
$B_z/ B_{0}$
as a function of space in KHI simulations with different models, at
$t = 200 \omega _{p,e}^{-1}$
: fully kinetic Vlasov simulation (panel a); hybrid simulation, with Vlasov ions and fluid, 10-moment electrons (panel b); two-fluid 10-moment simulation (panel c); two-fluid simulation, with 10-moment ions and 5-moment electrons (panel d); two-fluid, 5-moment simulation (panel e). The colour indicates the strength of the out-of-plane magnetic field
$B_z$
normalised to the initial field strength
$B_0$
. Horizontal and vertical axes show position normalised to the half-width of the initial shear layer
$L_0$
.
By comparing the Vlasov and hybrid simulations, panels (a) and (b) in figure 2, we preliminarily observe that the electron model does not affect vortex position and shape: this is consistent with the fact that the KHI is a shearing instability and the fluid velocity in a plasma is essentially the ion velocity. This observation is also supported by the fact that also panels (c) and (d) (10-moment or 5-moment electrons, 10-moment ions) are essentially identical: the electron dynamics does not affect the growth phase and early nonlinear evolution of the KHI. The model we use to simulate ions instead affects the internal shape of the vortex, as evident comparing panels (b) and (c), i.e. Vlasov or 10-moment ions: in the presence of 10-moment ions, higher wavenumber (smaller-scale) oscillations appear inside the larger-scale vortices. When describing ions with a 5-moment model as well, panel (e), we observe the development of even smaller-scale vortices due to secondary instabilities at the boundary layers: since heat flux dissipation is absent in the 5-moment model (see discussion in last paragraph of § 2), anisotropic ions are needed to introduce a high-wavenumber cutoff in the modes unstable to KHI. Our results are consistent with previous work. Cerri et al. (Reference Cerri, Henri, Califano, Sarto, Daniele and Pegoraro2013) argument that a sheared velocity field introduces anisotropies in the pressure tensor on very fast temporal scales, and that at least the parallel and perpendicular pressure components have to be evolved independently. They call their model, which evolves 6 equations per plasma species (1 for density, 3 for the current, 1 for the parallel and 1 for the perpendicular pressure component), the extended two-fluid model. Del Sarto et al. (Reference Del Sarto, Pegoraro and Califano2016) investigated initial condition configurations for the simulation of the KHI. They found that velocity shears introduce not only anisotropies, but also agyrotropies in the velocity distribution function (VDF), hence the full pressure tensor has to be evolved. Electrons are treated as a massless fluid, supporting our results that the electron dynamics is not essential during the early stages of KHI evolution, although the electron dynamics can become important in later nonlinear stages, see e.g. Goodwill et al. (Reference Goodwill, Adhikari, Li, Pucci, Yang, Guo and Matthaeus2025). This model comparison shows that accurate simulations of the plasma mixing dynamics in magnetised plasmas requires a 10-moment description, to avoid introducing unphysical wavenumbers in the simulation and to reproduce faithfully the anisotropies and agyrotropies in the VDF introduced by shear flows. This is particularly important in turbulence simulations, addressed in § 4.3, where we want an accurate description of spectra across a wide range of wavenumbers and in the presence of shearing motions. A 10-moment model for the ions, with appropriate closure, appears sufficient to reproduce the main aspects of plasma mixing at a boundary.
4.3. Decaying turbulence
Building upon the insights gained from our Landau damping and Kelvin–Helmholtz investigations, we now address the central challenge of this work: identifying the closure parameters that best reproduces key characteristics of decaying turbulence simulations. Decaying turbulence provides an ideal testbed for this work as it naturally generates a broad spectrum of fluctuations across scales, from MHD-like behaviour at large scales to kinetic physics at small scales. This multi-scale energy cascade, and particularly the transition region between fluid and kinetic regimes, is precisely where our Landau-fluid closure must perform accurately. Our choice of
$k_{0,e}^o$
and
$k_{0,i}^o$
will be guided by the attempt to match at best the spectra from the reference Vlasov simulations, and to reproduce the spatial structure of electron and ion current, as well as the divergence of the electron and ion heat fluxes in the hybrid and 10-moment simulations.
We employ the Biskamp–Welter set-up described in § 3, which initialises low correlation velocity and magnetic fluctuations that rapidly develop into complex turbulent structures through nonlinear interactions. Unlike forced turbulence, this decaying configuration allows us to observe how different models handle the natural evolution of the cascade, with particular attention to the transition region from fluid-like behaviour (above the ion skin depth) to the complex dynamics that emerges at kinetic scales.
We proceed through distinct phases: first, we evaluate varying electron closure parameters
$k_{0,e}$
in hybrid simulations against a fully kinetic reference simulation to establish the optimal
$k_{0,e}^o$
for the electron dynamics, where the superscript stands for ‘optimal’. Subsequently, we perform 10-moment fluid simulations maintaining the established
$k_{0,e}^o$
for electrons while varying the ion parameter
$k_{0,i}$
to determine the optimal
$k_{0,i}^o$
for ion dynamics. To identify the most suitable time range for spectral analysis, we examine the temporal evolution of the enstrophy-type quantities
$\int |\boldsymbol{\nabla }\times \boldsymbol{u_{cm}}|^2 d^3x$
and
$\int |\boldsymbol{\nabla }\times \boldsymbol{B}|^2 d^3x$
for the centre-of-mass velocity
$\boldsymbol{u}_{cm}= {(m_e \boldsymbol{u}_e + {m_i} \boldsymbol{u}_i)}/{(m_i+m_e)}$
and magnetic field
$\boldsymbol B$
in the Vlasov simulation, in figure 3. Starting from the smooth Biskamp–Welter initial condition, strong, isolated current sheets form with nearly exponential growth (Friedel, Grauer & Marliani Reference Friedel, Grauer and Marliani1997) up to time
$t \approx 35 \varOmega _i^{-1}$
. After this initial phase, the current sheets become unstable and turbulence develops.
Enstrophy evolution as a function of time in the fully kinetic Vlasov simulation for the centre-of-mass velocity and the magnetic field.
List of the turbulence runs examined in this work. For each run, we provide a run ID (BW stands for Biskamp–Welter), the model type (kinetic, hybrid or two-fluid 10-moment) and, when needed, the
$k_0$
value used for the heat flux closure, normalised to the respective skin depth. The
$k_0$
value that compares best with the kinetic benchmark is marked in bold.
We select the interval
$t \in [50, 100]\,\varOmega _i^{-1}$
for our spectral analysis, corresponding to the early decay phase where nonlinear interactions remain significant while avoiding the initial transient dynamics. This interval captures the fully developed turbulent state where the differences between kinetic and fluid models become most apparent, providing an ideal window for comparing spectral characteristics. The runs we examine are summarised in table 2. For each run, we report the run ID and the model used. Run K is the fully kinetic run of reference. The runs marked as H are hybrid runs, where ions are simulated with the Vlasov model and electron with the 10-moment model and the closure described in (2.10). The
$k_{0,e}$
we use in each run is listed in the corresponding row. The ‘optimal’
$k_{0,e}^o$
, chosen as described in the rest of the section, is marked in bold. The F10 runs employ the 10-moment model for both electrons and ions. We close the electrons using the optimal
$k_{0,e}^o$
selected in the hybrid runs. The different
$k_{0,i}$
we use are listed in the corresponding lines. The ‘optimal’
$k_{0,i}^o$
is again marked in bold.
We start by examining the hybrid simulations. In figure 4, we depict the electron velocity, ion velocity, magnetic field and electric field spectra, as a function of the wavenumber
$k\;d_i$
, for the Vlasov, BW-K, simulation and for the hybrid simulations, run with the different
$k_{0,e} d_e$
. The spectra are compensated by multiplying them by
$k^{po}$
where
$k^{po}=1.5$
for the velocity spectra and
$k^{po}=2$
for the field spectra. In all cases, we observe a remarkable similarity between the hybrid spectra, notwithstanding the different
$k_{0,e}$
closure, and the Vlasov reference.
Compensated spectra for the hybrid simulations with varying
$k_{0,e}$
, compared with the Vlasov reference: (a) electron velocity, (b) ion velocity, (c) magnetic field, (d) electric field spectra. The different hybrid simulations are marked in colours. The Vlasov reference is depicted in dashed line.
In figure 5 we plot several fields in the
$xy$
plane at time
$t=100\,\varOmega _{c,i}^{-1}$
. The spectral norm of the closure expression for each point in space is calculated as
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s\|(x,y) = \sigma _{\textrm {max}}((\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s)(x,y))$
, where
$\sigma _{\textrm {max}}()$
returns the largest singular value of the argument. First, in the first row, panel (a), we depict the magnitude of the electron current density
${\boldsymbol{j}}_e= n_e {\boldsymbol{u}}_e$
in the reference Vlasov simulation, to highlight the evolution of the turbulent patterns in the simulations. In panel (f),
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e\|$
is calculated directly from the velocity distribution function in the Vlasov simulations: this is the ‘ground truth’ we strive to reproduce by selecting
$k_{0,e}$
appropriately. In panel (g), as a comparison, we depict the absolute value of the heat flux calculated from the moments of the Vlasov simulation using (2.10) and a
$k_{0,e}$
value,
$k_{0,e} d_e= 201$
, obtained by comparison with the results in panel (b). Comparing panel (f) and panel (g), we are reassured that (2.10) is a good choice for our closure. We expect similar, but not identical, results when plotting the closure obtained from (2.10) for the hybrid simulations (panels (h)–(j)), because the lower moments used to calculate the closure are affected by the closure itself. In panels (h)–(j) we depict the closure obtained from (2.10), with the
$k_0 d_e$
value indicated in the caption. In the centre row of figure 5, we rescale the individual plots to ad hoc colour scales: the aim is to highlight how the turbulent patterns change with the different
$k_{0,e}$
. In the top row, panels (b)–(e) we use for all plots the same colour range as the kinetic benchmark (f): we want to identify which
$k_{0,e}$
value delivers magnitude of the heat flux closer to the kinetic benchmark. As one could expect from (2.10), higher
$k_0 d_e$
values result in lower
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e\|$
magnitudes. After analysing figures 4 and 5, we choose
$k_{0,e} d_e= 200$
as our best fit. We, however, remark that both spectra and turbulent patterns do not change significantly with varying
$k_{0,e}$
in the hybrid simulations. In figure 5(k–o) we depict the electron
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e\| \boldsymbol{\cdot }|j_e|$
, with
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e\|$
in the different panels calculated as in panels (f)–(j). By comparing the different panels with corresponding plots for the electron current (see panel (a) for the electron current density in the Vlasov simulation), we see that the peaks of
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e\|$
are well aligned with the areas of maximum current density. This suggests that it is particularly critical to correctly capture heat flux closures in correspondence of the current density peaks, if one intends to reproduce fundamental turbulence properties such as intermittency (Wan et al. Reference Wan, Matthaeus, Karimabadi, Roytershteyn, Shay, Wu, Daughton, Loring and Chapman2012).
Comparison of two-dimensional colour maps at time
$t=100\,\varOmega _{c,i}^{-1}$
showing: panel (a):
$|j_e|$
from the Vlasov reference simulations. Panel (f):
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e \|$
calculated from the distribution function of the Vlasov simulation; panel (g): calculated with (2.10) using moments from the Vlasov simulations. By matching amplitude with (f), we a posteriori determine
$k_{{0,e}}\,d_e= 201$
;
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e \|$
from hybrid simulations using (c, h)
$k_0 d_e= 20$
, (d, i)
$k_0 d_e= 200$
and (e, j)
$k_0 d_e= 2000$
. In (g)–(j), we adjust the colour scale in every panel. (b)–(e) Show the same data with colour scales set to the range of the kinetic ‘ground truth’ (f). (k)–(o) depict
$\| \boldsymbol{\nabla }\boldsymbol{\cdot }Q _e\| \boldsymbol{\cdot }|j|$
, with
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_e \|$
calculated as in the row above.
After having identified
$k_{0,e}^o d_e= 200$
as the optimal closure parameter for our hybrid simulations, we repeat the same exercise with the two-fluid 10-moment simulation series in table 2. As for the hybrid simulations, we plot in figure 6 the compensated electron velocity, ion velocity, magnetic field and electric field spectra as a function of the wavenumber. We depict results for the Vlasov, K, simulation and for the 10-moment simulations, run with
$k_{0,e}= k_{0,e}^o$
for the electrons and the different
$k_{0,i}$
. for the ions. We observe that, in this case, the choice of the closure parameter affects quite significantly the different spectra, with the exception of the magnetic field spectra.
Compensated spectra for the fluid simulations with
$k_{0,e} d_e = 200$
and with varying
$k_{0,i}$
, compared with the vlasov reference: (a) electron velocity, (b) ion velocity, (c) magnetic field, (d) electric field spectra. The different two-fluid simulations are marked in colours. The Vlasov reference is depicted in dashed line.
Comparison of ion heat flux divergence
$\| \boldsymbol{\nabla }\boldsymbol{\cdot }Q _i\|$
at time
$t=100\,\varOmega _{c,i}^{-1}$
for 10-moment simulations with
$k_{0,e} d_e = 200$
and varying
$k_{0,i}$
: the ion closures use (c, h, m)
$k_0 d_i= 2$
, (d, i, n)
$k_0 d_i= 20$
and (e, j, o)
$k_0 d_i= 200$
. Panel compositions are the same as in figure 5.
Comparison of electron heat flux divergence
$\| \boldsymbol{\nabla }\boldsymbol{\cdot }Q _e\|$
at time
$t=100\,\varOmega _{c,i}^{-1}$
for 10-moment simulations with
$k_{0,e} d_e = 200$
and varying
$k_{0,i}$
: the ion closures use (c, h, m)
$k_0 d_i= 2$
, (d, i, n)
$k_0 d_i= 20$
and (e, j, o)
$k_0 d_i= 200$
. Panel compositions are the same as in figure 5.
Figures 7 and 8 highlight the effects of closures applied to the ion species to both ion (figure 7) and electron (figure 8) quantities. Note that we only vary the ion closure, the variation in figure 8 shows the effect of the ion closure on the electron species. First, in panel (a) we depict the electron and ion current density respectively. We observe that, as expected, electron current structures develop at smaller scales than ions. Panels (f) and (g) depict
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s\|$
from Vlasov simulations, with
$s$
the particle species. In panel (f) we depict the ‘true’
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s\|$
, calculated from the distribution function, in panel (g) (2.10), calculated using the moments from the Vlasov simulations. Quite surprisingly, we observe in figure 7 that the ion heat flux divergence spectral norm
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_i\|$
exhibit smaller scale signatures than the ion current density. In panels (c)–(e) and (h)–(j), we depict the
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s\|$
closure, calculated following (2.10) from the two-fluid 10-moment simulations reported in the respective column titles. The
$k_{0,i} d_i$
for the ions is
$k_{0,i} d_i=2, 20, 200$
, while for the electron we keep it fixed at
$k_{0,e}^o d_e= 200$
. As before in figure 5, we use in the first row the same colour scale as the ‘ground truth’, individual colour scales in the second row. Comparing the fluid ion closures with the Vlasov ‘ground truth’ we observe, this time rather starkly, the effects that the choice of the
$k_{0,i}$
parameters has on the formation of patterns at different scales: the higher the
$k_{0,i} d_i$
(i.e. the lower the
$\boldsymbol{\nabla }\boldsymbol{\cdot }Q_i$
), the more small-scale signatures emerge in the fluid closures. Comparing the electron closures from the fluid simulations with the Vlasov ‘ground truth’ and the respective ion closure, we notice the effect of large-to-small-scale coupling: the high magnitude signatures in the electron closures resemble more closely the respective ion signatures rather than the Vlasov electron ‘ground truth’. This is particularly evident in the last column,
$k_{0,i} d_i= 200$
. The rather high value of
$k_{0,i} d_i$
results in unphysical small-scale structures in the ion closure: the same small-scale structures are visible in the electron closure in figure 8 as well. Furthermore, both the
$k_{0,i} d_i= 2$
and
$k_{0,i} d_i= 20$
electron closures (figures 8
c, d) exhibit a specific ‘double-arcade’ pattern at
$10 \lesssim x/ d_i \lesssim 20$
;
$y/ d_i \lesssim 2.5$
. The same pattern appears in the respective ion closures, figure 7, but not in the ‘ground truth’ Vlasov electron heat flux, figure 8(f), nor in the depiction of (2.10) calculated with moments from the Vlasov simulations in figure 8(g): it must have been driven by the selected ion closure.
In figures 7 and 8(k–o) we depict
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s\| \boldsymbol{\cdot }|j_s|$
, with
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s\|$
calculated as in figure 5 and
$s=i,\:e$
, respectively. Also here, we observe that the peaks of
$\|\boldsymbol{\nabla }\boldsymbol{\cdot }Q_s\|$
are well aligned with those of
$ |j_s|$
.
Analysing together figures 6, 7 and 8 we choose
$k_{0,i} d_i= 20$
as the optimal closure parameters for the ions. This value corresponds to a characteristic heat transport length scale of
$1/k_{0,i} \approx 0.05\, d_i$
, which lies in the kinetic range where ion dissipation processes become important.
Energy evolution for the Vlasov simulation (first column), ‘best’ hybrid simulation, with
$k_{0,e}^0 d_e = 200$
(second column) and ‘best’ two-fluid 10-moment simulation, with
$k_{0,e}^0 d_e = 200$
and
$k_{0,i}^0 d_i = 20$
(third column). Top row: total energy and each component is normalised to their respective initial value. Bottom row: total energy and each component is normalised to initial value of total energy.
In figure 9 we depict the energy evolution as a function of time for three simulations: the Vlasov simulation (first column), the ‘best’ hybrid simulation, with
$k_{0,e}^0 d_e = 200$
(second column) and the ‘best’ two-fluid 10-moment simulation, with
$k_{0,e}^0 d_e = 200$
and
$k_{0,i}^0 d_i = 20$
(third column). In the first row, the total energy and the different energy components (electric, magnetic, electron kinetic and thermal, ion kinetic and thermal) are normalised to their initial respective value; in the second row, to the total initial energy. Comparing the three simulations, we observe excellent energy conservation in all cases (minimal energy decrease is observed in the two-fluid simulation), and also good agreement between the evolution of the different energy components between the Vlasov, hybrid and two-fluid cases. In the first row of figure 9 we observe some difference between the electron kinetic and thermal energy evolution in the Vlasov simulation, on the one hand, and hybrid and two-fluid simulations on the other: at the end of the Vlasov simulation, the electrons exhibit lower kinetic energy than in the hybrid and two-fluid simulations, but are hotter. We observe, however, in the second row that both the electron thermal and, in particular, the electron kinetic energy components are rather low, in percentage, with respect to the total energy content of the simulation. This is the reason why this different distribution of the electron energy seems to impact minimally on the Vlasov and hybrid spectra, which are remarkably similar (see figure 4). On the other hand, comparing the Vlasov, hybrid and two-fluid energy evolution, we can conclude that the difference in spectra between the two-fluid and Vlasov cases in figure 6 is not related to different energy distribution among the different energy components, but rather to different energy distribution between different scales: in figure 6 we see that higher
$k_{0,i} d_i$
values result in higher energy at high wavenumber for the electron and ion velocity spectra, and for the electric field spectra. In the energy plot, however, no significant difference between the Vlasov and two-fluid simulations is visible for the electric field energy and for electron and ion kinetic energy.
5. Discussion and conclusions
This work is motivated by one of the most challenging aspects of computational plasma physics, and computational space physics in particular: since plasmas are inherently multi-scale systems, simulations need to accommodate, together, the large box sizes and temporal durations needed to capture large-scale, slow processes and the small-scale, fast, kinetic processes that critically influence energy dynamics and hence even large-scale evolution. We investigate here two-fluid 10-moment models as a compromise between fully kinetic and MHD models, which have the respective drawbacks of being either extremely computationally expensive or inaccurate in the modelling of kinetic processes. Ten-moment methods evolve velocity moments up to the pressure tensor included, and therefore need a closure for the divergence of the heat flux tensor. The choice of the closure is essential in reproducing the effects at large scales of kinetic-scale energy transport processes, at least for some of the test cases we examine.
In this paper we have focused on decaying turbulence, a problem where, on the one hand, one needs large separation in wavenumber between the scale of energy injection and that of energy dissipation and, on the other hand, energy dissipation at kinetic scales has to be properly modelled for its effects in determining energy spectra at intermediate and kinetic scales. We have chosen a specific local Landau-fluid closure, (2.10), and we have investigated if specific choices of the closure parameters
$k_{0,s}$
for both electrons and ions can deliver spectra and turbulent structures which compare well against reference fully kinetic simulations. Before addressing the decaying turbulence problem, we have tuned our methodology on ‘simpler’ problems, i.e. Landau damping and KHI, where we could address separately the role of the heat flux closure in wave damping and velocity shear instability, respectively.
The results of our work are the following: first, our investigations on the Landau damping problem shows that 10-moment models with our specific Landau-fluid closure and an appropriately chosen closure parameter
$k_{0,s}$
reproduce wave damping into the plasmas (damping rate, period of the electric field oscillation) with sufficient accuracy for our intended usage in large-scale simulations. This is not a surprise, since Landau-fluid closures have been designed to reproduce exactly this process. Second: our Kelvin–Helmholtz simulations show that a ‘minimum’ model exists for the simulation of velocity shear-driven instabilities such as the KHI. This ‘minimum’ model is a 10-moment model for the ions, 5-moment model for the electrons, that allows for the formation of ion anisotropies and agyrotropies triggered by the presence of a velocity shear. Within this model, a wide range of
$k_{0,i}$
closure parameters for the ions deliver comparable results. In contrast, when a 5-moment model is used for the ions as well, high-wavenumber spurious oscillations are observed in the simulations. This result is of extreme interest for a wide range of physical problems, from turbulence simulations to global magnetospheric simulations, where velocity shear instabilities occur. Third, our decaying turbulence simulations highlight the role of the ions in energy transport from fluid-like to kinetic scales: while in 10-moment hybrid simulations (fluid electrons, kinetic ions) energy spectra for both fields and particles are well reproduced notwithstanding the choice of the electron closure parameters, the choice of the closure parameter for the ions is more critical. When testing multiple
$k_{0,i}$
closure parameters for the ions, we aim at matching against the reference Vlasov simulations both the energy spectra and the magnitude and structures in the divergence of the heat flux calculated in the 10-moment simulations from lower moments with (2.10). The
$k_{0,i}$
value we consider ‘optimal’ after comparison with the Vlasov reference,
$k_{0,i}\,d_i= 20$
, satisfies both these requirement. Lower and higher
$k_{0,i}$
values tend to respectively over- and under-estimate heat flux divergence, as one could expect from the closure we use. Interestingly, higher
$k_{0,i}$
values also results in the appearance of high-wavenumber spurious oscillations in the calculated heat flux divergence for both the ions and the electrons, showing that the ions are essential in transmitting energies from the fluid-like to the kinetic scales. This conclusion is further strengthened by the observation that structures in the electron closure from 10-moment simulations appear to map directly to similar structures in the ion closures, which are absent in the Vlasov ‘ground truth’. We remark at this stage that, while for some of the test cases we have benchmarked our 10-moment results against Vlasov simulations using quantitative metrics (e.g. damping rate and oscillation period for the Landau damping case, spectral slope for the turbulence simulations), additional quantitative metrics, such as spectral slope fits, cross-correlation coefficients and normalised error norms could be employed in future or similar work to further assess the degree of agreement between the two sets of simulations and make benchmarking activities more quantitative.
Further considerations can be derived from our work. First, we have demonstrated the robustness of Landau-fluid closures with respect to some of the assumptions under which they are derived. Landau-fluid closures are derived under the assumption that the plasma is close to LTE. However, they still perform quite well in our Kelvin–Helmholtz and decaying turbulence test cases (which, as with most of the plasma processes of interest, are far from LTE) as long as appropriate closure parameters
$k_{0,s}$
are found. Second, the closure parameters
$k_{0,s}$
that we have identified as ‘optimal’ in all our test cases not only reproduce well the large-scale consequences of energy dissipation, but also are in the range of the spatial scales where we expect dissipation processes to become of relevance for the different plasma species (Ng et al. Reference Ng, Hakim, Wang and Bhattacharjee2020). This confirms that we can interpret the inverse of the
$k_{0,s}$
parameter as a characteristic length for heat transport and gives a hint for how the parameter should be chosen. Third, we have highlighted the importance of appropriately approximating heat flux closures especially in regions of high current densities, which are associated with intermittent and dissipative behaviour. This hints to the possibility of using clustering methods in simulations (see e.g. Donaghy & Germaschewski Reference Donaghy and Germaschewski2023; Köhne et al. Reference Köhne, Boella and Innocenti2023) to identify dissipative regions, and then design heat flow closures, whereby the training process is conditioned on regions with high dissipation.
We acknowledge a number of limitations in our current investigation, which will constitute the basis for future work. First, in this work we have considered one single term, of analytical origin, as our closure. In recent years, a number of works have focused on extracting closure equations for reduced models (in our case, the 10-moment model) from simulated data obtained using more expressive models (e.g. fully kinetic models), considered as the ‘ground truth’. The specific strategy for extracting closures from the data varies: regression methods such as the sparse identification of nonlinear dynamics method, trained on the more expressive model, may be used to extract a closure in symbolic form (Alves & Fiuza Reference Alves and Fiuza2022; Donaghy & Germaschewski Reference Donaghy and Germaschewski2023). Alternatively, neural methods, again trained on the more expressive model, may be used to predict non-symbolic closures on the basis of quantities available in the reduced-model simulations as well (Cheng et al. Reference Cheng, Fu, Wang, Dong, Jin, Jiang, Ma, Qin and Liu2023; Huang, Dong & Wang Reference Huang, Dong and Wang2025; Miloshevich et al. Reference Miloshevich, Vranckx, de Oliveira, Nathan, Dazzi, Arrò and Lapenta2025). Especially this latter approach may prove more robust in tackling spatially highly non-uniform or intermittent problems, where applying the same closure term with the same
$k_{0,s}$
coefficient over the entire simulations may give sub-optimal results. Validation of such strategies are left as future work. Second, we have demonstrated that optimal closure parameters exist under very specific plasma conditions, without verifying how well these parameters perform when plasma or simulation parameters change. Examples of the former are plasma betas, parallel to perpendicular temperature ratio; of the latter, ion to electron mass ratio, box size. Proving that the same closure parameters can be used under different plasma conditions or that effective closure parameters can be derived as a function of plasma conditions is the focus of ongoing investigation. Finally, and perhaps most important: another assumption behind the Landau-fluid closure we used is that transport is isotropic. This is the case in the simulations we have presented here, but very much not in the space physics systems we ultimately target, such as the solar wind, where parallel and perpendicular transports are quite different. In future work, we intend to verify if this anisotropy can be incorporated into 10-moment models while retaining our current closure, (2.10), perhaps using different closure parameters in the parallel and perpendicular directions, or if different closure equations have to be used in the different directions.
Acknowledgments
The authors acknowledge support from the German Research Foundation DFG within the Collaborative Research Centre SFB1491 and projects 497938371. We also gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) and SuperMUC-NG at the Leibniz Supercomputing Centre. Computations were conducted on JUWELS/JUWELS-booster, SuperMUC-NG and on the Galileo cluster of the SFB1491.
Editor Thierry Passot thanks the referees for their advice in evaluating this article.
Decleration of interests
The authors report no conflict of interest.
Open access
