High-fidelity kinetic modeling of instabilities and gyromotion physics in nonuniform low-beta plasmas

Collisionless low-β plasmas, where β is the ratio of thermal pressure to magnetic pressure, can have complex behavior and transport properties. In particular, when gyroradii are comparable to other scale lengths, dynamics are not well-captured by fluid-model descriptions and may not satisfy the low-frequency assumptions of the gyrokinetic theory, thereby necessitating a fully kinetic treatment. Examples of such configurations include magnetically insulated transmission lines in pulsed power experiments,^1–3 magnetically insulated diodes,^4–7 low-density regions of magnetic confinement experiments,^8–10 and space plasmas.¹¹ Finite-gyromotion plasmas are known to have unique properties including cyclotron resonances, microinstabilities, non-Maxwellian distribution functions, charge separation, modified macroinstability growth rates, and polarization currents—aspects of which have been analyzed theoretically.^12–14 Since kinetic theory analysis is often intractable, kinetic simulations that accurately capture these effects are critical to expanding our understanding of low-beta collisionless plasmas.

To model collisionless low-beta plasmas with high fidelity, kinetic solvers have to: contend with the computational cost of discretizing a high-dimensional phase space, accurately evolve plasma dynamics over many gyroperiods, and facilitate means to systematically study different configurations. The first aspect has been partially addressed through advances in supercomputing capabilities,^15,16 which have facilitated calculations with a large number of degrees of freedom. Progress has been made on the second aspect through the development and deployment of continuum kinetic solvers,^17,18 which directly discretize the governing equations of kinetic theory and which can take advantage of high-order numerical methods that provide enhanced solution accuracy.^19–26 Inclusion of magnetic field physics has extended the utility and applicability of these methods.^3,22,24,26 Unlike traditional particle-in-cell (PIC) methods, continuum kinetic solvers are not subject to noise and have robust convergence properties, which allows them to maintain accuracy over long periods of simulated time.²² The last aspect, i.e., the systematic study of various plasma configurations, has—until recently—remained a challenge because few kinetic equilibria exist for magnetized plasmas, which has made it difficult to isolate and parse different physical processes in kinetic simulations.²⁷ The recent development of a method to construct customizable two-species kinetic equilibria for magnetized plasmas has overcome this hurdle and has thereby provided a means to systematically probe configurations that were previously inaccessible.^3,28

Here we demonstrate that the combination of a high-order continuum kinetic solver and initialization of customizable kinetic equilibria significantly expands the scope of physics that can be studied in kinetic simulations. The methodology opens the way to the targeted exploration of finite-gyromotion plasmas in a wide variety of nonuniform configurations. It offers a means to track the evolution of distribution functions, as well as their moments, and associated derivatives to a high degree of accuracy, thereby enabling detailed noise-free analysis. This paper presents the advantages of the methodology and its utility in application to the study of gradient-driven instabilities. In particular, we explore the kinetic properties of shear-driven Kelvin–Helmholtz (KH) instabilities and the diamagnetic-drift-driven lower hybrid drift (LHD) instabilities in plasmas where ion gyroradii are comparable to gradient scale lengths. The contents of the paper are organized as follows. Section II describes the methodology, which is based on a fourth-order accurate conservative finite-volume kinetic solver^22,23 and a method for constructing customizable kinetic equilibria.²⁸ Section III explores the kinetic properties of KH instabilities, with a focus on the pressure tensor and the shear stress tensor, which are used to assess the applicability of gyroviscous fluid theory. Section IV explores transport and energy partitioning properties associated with the LHD instability in the context of a low-beta collisionless plasma. Section V presents concluding remarks.

We consider a low-beta collisionless plasma immersed in a magnetic field. Assuming that the magnetic fields generated by plasma currents are negligible and that the plasma is uniform along the magnetic field, then all dynamics will occur in the perpendicular plane and the plasma can be treated as electrostatic and magnetostatic.²⁹ The evolution of such a plasma is described by the nondimensionalized Vlasov–Poisson equation system,

where $fs(x,v,t)$ is the probability distribution function for species s, M_s is the ratio of particle mass to the proton mass m_p, Z_s is the ratio of particle charge to the magnitude of the electron charge e, Ω_p = eB₀/m_p is the characteristic proton cyclotron frequency defined in terms of the characteristic magnetic field B₀, and ω_pp is the characteristic proton plasma frequency. In Eq. (1), the Vlasov equation is expressed in conservation-law form. The electric field is defined as $E=−∇ϕ$⁠, where $ϕ$ is the electrostatic potential. The species number density n_s is obtained from the zeroth velocity moment of the associated distribution function,

In Eqs. (1)–(3), time is normalized to the proton plasma frequency $ωpp$ and velocity is normalized to the proton Alfvén speed v_A. In effect, the characteristic scale length is $L0=vA/ωpp$⁠, the characteristic temperature is $mpvA2$⁠, and the potential is normalized to $mpvA2/e$⁠. For simplicity, we limit our consideration to a hydrogen plasma composed of two species, protons (which we will also refer to as ions) and electrons. Thus, $Zi=−Ze=1$ and M_i = 1. The normalized characteristic species cyclotron frequency is $Ωs=(ZsΩp)/(Msωpp)$⁠. For the computational studies presented here, we take the nondimensional magnetic field to be uniform $B=ẑ$ and fixed in time and we consider dynamics in the (x, y) plane.

The approach used to initialize and solve the Vlasov–Poisson equation system dictates the accuracy and physical fidelity of kinetic simulations. Importantly, it also sets the scope for the physics that can be captured in simulations. Here the methodology is based on a high-order accurate continuum kinetic Vlasov–Poisson solver and a systematic means to initialize customizable kinetic equilibria. While we consider electrostatic and magnetostatic systems, the approach extends directly to a Vlasov–Maxwell treatment and Vlasov–Maxwell equilibria.

Numerical methods that discretize the Vlasov equation directly were—until recently—thought to be infeasible for practical problems, due to the high computational cost of discretizing four-, five-, or six-dimensional position-velocity phase space.¹⁸ Recent advances in supercomputing capabilities, however, have enabled direct discretization of the Vlasov equation in four,^{3,24,25,30–34} five,²⁴ and six^35,36 dimensions, including for applications where both ion and electron dynamics are governed by the Vlasov equation.^{3,24,25,32,33,36} These advances have transformed the paradigm for kinetic simulations and have led to wider use of continuum kinetic methods, which represent the distribution function in terms of piecewise smooth functions (e.g., polynomials) on a grid in phase space.

Continuum kinetic methods have advantageous mathematical properties. Unlike traditionally used particle-in-cell (PIC) methods, continuum kinetic methods are not subject to noise and have robust convergence properties. This means that accuracy increases systematically with resolution,^22–26 such that one can specify the resolution necessary to achieve a given error tolerance. A comparison between continuum kinetic and PIC simulations is presented in Ref. 37. Continuum kinetic Vlasov solvers have also benefited from recent advances in numerical methods for solving hyperbolic partial differential equations.^19–26,38 For example, high-order finite-difference and finite-volume discretization methods provide enhanced accuracy with less degrees of freedom—a feature that is particularly advantageous for high-dimensional phase space problems.^22,25 Furthermore, the use of conservative discretization methods, which solve the conservation-law form of the Vlasov equation [see Eq. (1)], ensures exact conservation of the primary variable, i.e., the distribution function, which results in exact conservation of mass and charge. Continuum kinetic methods can also conserve total energy to a high degree of accuracy,^{23–26,38–41} further ensuring high-fidelity modeling of plasma dynamics.

Leveraging these advances in hardware and numerical methods, a conservative fourth-order accurate finite-volume discretization has been developed for solving the Vlasov–Poisson equation system and has been implemented in the Vlasov Continuum Kinetics (VCK) code.^3,22,23,28 The unsplit method solves the conservation-law form of the Vlasov equation, which ensures that the distribution function, and by extension, its zeroth velocity moment is conserved to machine precision.^22,23 The solver achieves fourth-order accuracy in space and time by using a fourth-order quadrature rule⁴² to evaluate fluxes and a fourth-order explicit Runge–Kutta method to advance the solution in time.^22,23 Like all high-order methods, the solver requires the inclusion of high-order dissipation to damp modes whose motion is poorly represented on the grid. This is achieved by using a fifth-order upwinded reconstruction of the distribution function when evaluating fluxes.^22,23 The solver has been benchmarked and its convergence properties have been tested in Cartesian²² and cylindrical²³ phase space coordinates. Unlike more standard second-order solvers, the fourth-order solver results in high-accuracy converged solutions even at modest resolutions.²² The solver facilitates smooth solutions and enables detailed analysis of distribution functions, their moments, and associated derivatives.^3,22,23 Since the solver is not subject to noise and can preserve complex kinetic equilibria,²⁸ it is highly effective at tracking the evolution of low-amplitude perturbations and thereby facilitates the detailed study of both linear and nonlinear kinetic physics.^3,22

In addition to an accurate Vlasov solver, high-fidelity kinetic simulations require initial conditions that encapsulate the physics of interest, which may include finite Larmor motion, nonuniform density/velocity profiles, non-Maxwellian distribution functions, etc. Equilibrium initial conditions, in particular, constitute a useful starting point for computational investigations because they facilitate isolation of specific physical processes in kinetic simulations and enable detailed comparisons between simulations and linear theory.^3,22 For nonuniform finite-gyromotion plasmas, in which Larmor orbits are comparable to gradient scale lengths, equilibrium initial conditions that satisfy the steady-state governing equations are difficult to construct. This problem was recognized in the context of magnetized shear layer studies,^27,43,44 which demonstrated that using fluid equilibria to initialize kinetic simulations leads to dynamics that modify the plasma state and that complicate the interpretation of simulation results.

Efforts to derive better approximations to kinetic equilibria have utilized fluid equilibria that account for first-order finite Larmor radius effects⁴⁵ or have relied on integrating single-particle trajectories to construct distribution functions.^46,47 These magnetized shear-layer equilibria are restricted to uniform-density quasineutral plasmas.^45–47 In special cases, nonuniform kinetic equilibria can be constructed using truncated Hermite series expansions⁴⁸ or by stipulating rigid body rotation.⁴⁹ A recently developed method, based on solving a nonlinear ordinary differential equation, has made it possible to construct exact customizable kinetic equilibria for nonuniform low-beta plasmas.²⁸ To ensure that the steady-state Vlasov–Poisson equation system is satisfied, the method expresses ion and electron distribution functions in terms of constants of motion and relies on solving a nonlinear Poisson equation. The approach has similarities to the construction of Bernstein–Greene–Kruskal (BGK) modes in unmagnetized plasmas.⁵⁰ The associated equilibria are flexible in that they admit customizable ion density and electrostatic potential profiles, finite gyromotion for ions and electrons, finite space charge, and ion Larmor orbits that are comparable to gradient scale lengths.^3,28 Such equilibria enable targeted studies of isolated physics in finite-gyromotion plasmas in a wide range of configurations and are well-suited for situations where the ion drift speed is less than or approximately equal to the ion thermal speed.²⁸ When combined with the fourth-order continuum kinetic solver, these equilibria facilitate “quiet start” simulations, wherein single or multiple modes can be perturbed and their linear and nonlinear dynamics can be studied over long periods of time.³ 

The equilibrium construction method, which is described in detail in Ref. 28, is summarized here. Let $gi(x)$ be the desired ion density profile and $ϕ*(x)$ be the desired electrostatic potential profile. The number density and electric-field profiles from a uniform-temperature two-fluid theory equilibrium can be used as a basis for the choice of g_i and $ϕ*$⁠. Using g_i and $ϕ*$⁠, we can construct auxiliary ion and electron distribution functions,

where for a uniform magnetic field, $psy=x+vy/Ωs$ is the scaled canonical momentum for species s, and the nondimensionalized species cyclotron frequency Ω_s can be positive or negative depending on the species charge Z_s. The term $niaux(X)$ in Eq. (5) is the zeroth velocity moment [see Eq. (3)] of the distribution function $fiaux$⁠, evaluated at position argument X. For general choices of g_i and $ϕ*$⁠, the velocity moment of $fiaux$ cannot be evaluated analytically, in which case the analytic function $niaux(X)$ can be obtained by computing the velocity integral in Eq. (3) numerically on a discrete set of points and interpolating the result.²⁸ When $ϕ*=0, niaux(X)$ can be approximated by $gi(X)$⁠, which results in a quasineutral equilibrium. On a discrete set of points in the x domain, a nonlinear Poisson equation is constructed,

where $nsaux$ is the zeroth velocity moment of $fsaux$⁠. Given boundary conditions, Eq. (6) can be solved numerically for the equilibrium potential $ϕ$ by discretizing the Laplacian operator to fourth-order accuracy and using the Newton method with $ϕ=ϕ*$ as the initial guess.²⁸ Notably, as magnetization increases, $ϕ$ approaches $ϕ*$⁠. The exact equilibrium distribution function for each species is then given by

Provided plasma β $≪$ 1, an electromagnetic equilibrium, in which the magnetic field is nonuniform and Ampère's law is satisfied, can be obtained by iterating through the above procedure. In the first iteration, the magnetic field is taken to be uniform and $fs(1)$ is evaluated as described above. In the second iteration, the magnetic field $Bz(1)(x)$ is obtained from Ampère's law with species currents computed from $fs(1)$ such that the scaled canonical momentum in Eqs. (4) and (5) in the second iteration is $psy=(ZsAy(x)+Msvy)/(Ms$Ω_s), where $Ay(x)$ is the y-component of the magnetic vector potential associated with $Bz(1)(x)$⁠. Additional iterations may be necessary, depending on the value of plasma β.

To initialize continuum kinetic simulations, Eq. (7) is evaluated on a grid in $(x,y,vx,vy)$ phase space. Simply evaluating Eq. (7) at cell centers results in a second-order accurate initialization. Since initial conditions should have at least the same accuracy as the spatial discretization of the Vlasov solver (otherwise there is virtually no benefit to using a high-order method), a fourth-order quadrature rule based on discrete second derivatives is used to ensure a fourth-order accurate finite-volume initialization.²² For discretizations based on the finite element method, which typically requires analytic functions for initialization, Eq. (7) can be evaluated analytically so long as the discrete $ϕ$ that solves Eq. (6) can be evaluated analytically. An analytic expression for $ϕ$ can be facilitated by performing an interpolation of the discrete $ϕ$⁠, much like a polyharmonic spline interpolation can be used to evaluate $niaux(X)$⁠.²⁸ PIC simulations can also be initialized using distribution functions like those in Eq. (7), but require careful particle loading that accounts for the gyrophase.^51–53 

Together, the fourth-order accurate continuum kinetic Vlasov–Poisson solver and the systematic method for constructing equilibrium initial conditions provide an accurate and flexible means to study the linear and nonlinear physics of gradient-driven instabilities, while fully capturing gyromotion physics. We proceed to apply this methodology to explore two gradient-driven instabilities in 4D $(x,y,vx,vy)$ phase space.

The KH instability is a classical fluid instability that arises when the fluid velocity profile has a gradient and an inflection point. In classical fluid theory and in single-fluid theory for a highly magnetized plasma, the growth rate of the instability is proportional to the velocity shear in the vicinity of the inflection point. In plasmas that have ion gyroradii that are comparable to gradient scale lengths, the KH instability is strongly influenced by kinetic physics. A recent study employed the methodology described here to investigate kinetic and two-fluid physics in KH instabilities.³ The study showed that finite gyromotion sets a lower bound on gradient scale lengths, has a stabilizing effect on the KH instability, can suppress the development of secondary instabilities, and leads to ion distribution functions that deviate significantly from a Maxwellian. The study also demonstrated the importance of polarization drift in KH instabilities and its effect on growth rates and mass transport.³ Kinetic physics of KH instabilities is also explored in Refs. 27, 33, 43, and 53–55. Here the focus is not on the details of KH instability physics, but rather on the utility of the methodology described in Sec. II and how it facilitates extracting detailed information from kinetic simulations of magnetized plasmas.

Equilibrium initial conditions for KH instability simulations are given by Eq. (7) for both electrons and ions, with desired ion number density and potential profiles set to

Here d = 0.05 is shear layer half-width, $E0=2.0×10−2$ is the value of the electric field in the high-velocity region, b = −20 is the inverse scale length that characterizes the density gradient to the right of the shear layer, and $T=6.25×10−4$ is the nominal temperature. The mass ratio is $Mi/Me=25$ with M_i = 1. The initial ion Larmor radius is nominally $rLi=1ΩiTMi$⁠, but dips slightly in the middle of the shear layer, where the temperature—when computed from the velocity moments of the distribution function—is anisotropic.³ This particular choice for g_i and $ϕ*$ facilitates constructing $E×B$ shear layers resembling those that have been analyzed using linear theory.³ Specifically, g_i given by Eq. (8) produces an ion number density profile that is uniform to the left of the shear layer and decays exponentially to the right of the shear layer; and the combination of g_i and $ϕ*$ yields an ion velocity profile that closely resembles a hyperbolic tangent, wherein $E×B$ drift is the dominant drift. The ion number density and ion velocity profiles are shown in Fig. 1 for two different levels of magnetization: low magnetization with $Ωp/ωpp=1.0$ and high magnetization with $Ωp/ωpp=2.0$⁠. These two equilibria correspond to parameter cases A3 and A4 in Ref. 3, which includes additional details regarding initialization. Plasma β, evaluated in the region of maximum density, is $2.5×10−3$ for both magnetization cases. The phase space simulation domain is $[−Lx/2,Lx/2]×[−Ly/2,Ly/2]×[−vmax,s,vmax,s]×[−vmax,s,vmax,s]$⁠, where $Lx=1,Ly=π/4,vmax,e=1.0,vmax,i=0.2$⁠. The domain is periodic in the y direction and has reflecting-wall boundaries in the x-direction. To ensure conservation of the distribution function to machine precision, the velocity domain is set to have zero-flux boundary conditions. This condition does not affect the physics of the solution, provided that the value of the distribution function is sufficiently close to zero at velocity boundaries—as is the case here.^22,23 At x boundaries, Dirichlet conditions are specified for the potential $ϕ$⁠, with values set to $ϕ*$ boundary values. The grid resolution is $Nx×Ny×Nvx×Nvy=256×64×96×96$ and time step size is $Δt=4×10−4$ for both magnetization cases. In the absence of perturbations, the time-dependent Vlasov–Poisson solver used here, at the specified resolution, can maintain shear-layer equilibria to within about 0.04% of the initial condition for ions and within about 3% of the initial condition for electrons, and these deviations decrease with increasing resolution.²⁸ 

A small-amplitude perturbation of a shear layer equilibrium results in the growth of the KH instability, with distinguishable linear and nonlinear phases of evolution. Figure 2 shows the evolution of the perturbation, as measured by the magnitude of the x-directed ion velocity $|uix|$ evaluated in the center of the domain. During the linear phase, $|uix|$ grows exponentially and has a characteristic frequency. Eventually, instability growth saturates and a nonlinear phase ensues. Since continuum kinetic simulations are noise-free, they capture perturbation amplitude growth over many decades, which means that growth rates and oscillation frequencies can be easily extracted and quantitatively compared to linear theory, as demonstrated in Ref. 3. Importantly, continuum kinetic solvers—unlike many PIC solvers—can preserve complex magnetized-plasma equilibria,²⁸ which means that a perturbation can be tailored to excite a single mode or multiple modes. In continuum kinetic simulations, fast-growing and/or slow-growing modes can be studied in isolation, whereas in PIC simulations it is more difficult to isolate modes due to noise. The ability to preserve kinetic equilibria and perturb specific modes allows for targeted comparisons with linear theory and with simulations based on fluid models and thereby can elucidate the limits of applicability of a given model.³ 

As an example of analysis facilitated by continuum kinetic simulations and the differentiable solutions they provide, we investigate the properties of the pressure tensor in KH instabilities. A nondiagonal pressure tensor is a known feature of magnetized plasmas and is characterized in terms of a stress tensor Π_ij. For each species, the stress tensor is defined as

where P_ij is the total pressure tensor, P₀ is the scalar isotropic pressure, δ_ij is the Kronecker delta, and i, j are coordinate indices. The scalar pressure $P0=tr(Pij)/D$⁠, where $tr$ denotes trace and D is the number of velocity dimensions. By definition and by terminology commonly used in stress analysis, Π_ij is the stress deviator tensor for P_ij. The total pressure tensor can be evaluated from kinetic simulations using velocity moments of the distribution function, such that

where number density n is given by the zeroth moment of the distribution function [Eq. (3)] and drift velocity u_j is obtained from the zeroth and first velocity moments, such that

In Eqs. (10)–(12) and in the proceeding analysis, we omit the subscript s, denoting species, and consider only the stress tensor computed from the ion distribution function. Ions—due to their larger mass—are the drivers of the KH instability, have larger Larmor orbits, and are therefore more strongly influenced by kinetic processes. For a two-dimensional plasma in the (x, y) plane, the ion scalar pressure is $P0=12(Pxx+Pyy)$⁠. In the limit where ion gyroradii are much smaller than gradient scale lengths, the stress tensor can be approximated by the Braginskii form^56,57 of the gyroviscous stress tensor, which we will denote by $ΠijG$⁠. See Refs. 58 and 59 for a discussion of various forms of gyroviscous stress and associated assumptions. In 2D, the components of the gyroviscous stress are given by^56,57

In kinetic simulations, drift velocities u_x, u_y are evaluated using Eq. (12). For convenience, we use tilde (⁠$̃$⁠) to denote the stress tensor scaled by the scalar pressure, such that

To understand the role of finite Larmor radius effects and the degree to which they can be captured by reduced models, it is useful to compare the kinetic stress tensor Π_ij to the gyroviscous stress tensor $ΠijG$⁠. A direct comparison between components of Π_ij and $ΠijG$ is presented in Ref. 54 for KH instability evolution in a nonequilibrium shear layer configuration with $Ωp/ωpp=0.05$⁠. The study used a second-order accurate continuum kinetic Vlasov–Maxwell solver and showed that for configurations where $(∇×u)·B<0$⁠, the relative difference $(Π̃ij−Π̃ijG)$ for each component is 1%–3% in the linear and early nonlinear stage of the instability. While the stress tensor is straightforward to evaluate, each component of the stress tensor can be difficult to interpret because during the nonlinear stage of the KH instability, local features like the ion velocity gradient are not generally aligned with the x or y coordinates. A metric for the stress that encapsulates the magnitude of the local shear, independent of the orientation of the local features, can provide a clearer picture of local stress. As such, we use the eigenvalues of the stress tensor to quantify the scalar magnitude or intensity of the local stress and denote this magnitude symbolically by $⟦·⟧$⁠. In two dimensions, the stress magnitude is

where $Π1$ and $Π2$ are the eigenvalues of the stress tensor. The argument to the square root in Eq. (16) is the second invariant of the stress tensor. It is worth noting that tensor invariants are also used to characterize turbulence structures in fluid simulations^60,61 and have been used to quantify gyrotropy in kinetic simulations.⁶² As defined in Eq. (16), $⟦Πij⟧$ is proportional to the von Mises stress⁶³ in structural mechanics and has been used in continuum mechanics, e.g., to quantify stress in fluids confined in narrow channels.⁶⁴ Note that $⟦Π̃ij⟧=⟦Πij⟧/P0$⁠.

KH instability evolution for low magnetization, where $Ωp/ωpp=1$ and initially $rLi/d=0.5$⁠, is shown in Fig. 3. The figure shows the evolution of the scalar pressure P₀, the magnitude of the stress tensor relative to the scalar pressure $⟦Π̃ij⟧$⁠, the relative difference between the kinetic stress and the gyroviscous stress $⟦Π̃ij−Π̃ijG⟧$⁠, and the ratio of shear magnitude ν_s to the ion cyclotron frequency Ω_i. The velocity shear magnitude ν_s is a local quantity and we define it as

The argument to the square root in Eq. (17) is the second invariant of the deviator tensor for the tensor $12(∇u+(∇u)T)$⁠, which is closely related to $ΠijG$⁠. As the instability evolves into the nonlinear stage, gradients increase in localized regions, leading to larger velocity shear and thereby larger shear stress. At the same time, Larmor radii evolve in response to local heating/cooling, which occurs due to the fact that the plasma dynamics are not adiabatic (i.e., while the magnetic field is constant, the magnetic moment is not conserved). The lack of adiabaticity can be explained by the presence of resonances, which occur when local velocity shear takes on values close to the ion cyclotron frequency. Ion Larmor orbits set a lower bound on the local gradient scale length, which means that gyromotion acts to limit the local shear. When gradient scale lengths become comparable to the local Larmor orbit size, scalar pressure features become smeared out—an effect that is analogous, but not equivalent, to the effect of collision-induced viscosity. Regions with significant smearing are those that exhibit large shear and a local increase in temperature. The kinetic stress, as measured by the value of $⟦Π̃ij⟧$⁠, grows from 12% at the initial time to over 42% during the nonlinear stage. At the initial time, the relative difference between the kinetic stress and the gyroviscous stress is about 1%; however, as the instability develops, the gyroviscous stress becomes a poor approximation for the kinetic stress, such that the relative difference between the two stresses eventually exceeds 20%. The largest discrepancies between Π_ij and $ΠijG$ roughly coincide with regions where $νs/Ωi$ is peaked. Local values of $⟦Π̃ij−Π̃ijG⟧$ become significant, i.e., above 15%, when $νs/Ωi≳0.5$⁠. This happens around the same time and location that the ratio of ion Larmor orbit size to shear scale length approaches unity.

As shown in Fig. 4, the same general trends are observed for the case of high magnetization, where $Ωp/ωpp=2$ and initially $rLi/d=0.25$⁠. At initial time $⟦Π̃ij⟧$ is less than 4% and eventually exceeds 50%, indicating that the system supports larger stress late in time when the initial ion Larmor orbits are smaller. In the high-magnetization case, the plasma is more adiabatic since $max(νs/Ωi)$ at initial time is about a factor of three smaller than in the case of low magnetization. This means that at higher magnetization, resonance effects take more flow-through periods to develop and to become significant. Figure 4 also shows that throughout the initial nonlinear stage of the instability evolution, gyroviscous stress, and kinetic stress are approximately equal with a relative difference of less than 5%. However, later in the nonlinear evolution at t = 595.2, the relative difference exceeds 35%, indicating that strong magnetization does not obviate the need for kinetic simulations. These results support the general notion that when the ratio $νs/Ωi$ exceeds a threshold value—here about 0.5 (analogously when the local shear scale length approaches the ion Larmor orbit scale), the gyroviscous stress tensor becomes a poor approximation for the kinetic stress tensor. In the case of low magnetization, $νs/Ωi≳0.5$ in the early nonlinear stage of the evolution, whereas in the case of high magnetization $νs/Ωi$ locally exceeds a value of 0.5 only in the later stages of eddy roll-up. Regions with large values of $⟦Π̃ij−Π̃ijG⟧$ coincide with regions where gyroradii are comparable to the shear scale length and constitute a subset of the region where distribution functions deviate significantly away from a Maxwellian—see Ref. 3. These findings are also consistent with previous kinetic KH instability studies,⁵⁴ which showed that for an initial condition where vorticity is antiparallel to the magnetic field and $rLi/d=0.25$ (parameters that match our high-magnetization case), the relative discrepancy between gyroviscous stress and the kinetic stress is approximately 3% in the early evolution of the instability. Here we demonstrate that the discrepancy between $Π̃ij$ and $Π̃ijG$ during the nonlinear stage can actually be over 20% for low magnetization and late-time high-magnetization scenarios.

Kinetic equilibria are integral to being able to isolate the KH instability in simulations, and study it in the absence of dynamic sheaths, Langmuir waves, lower-hybrid waves, and upper-hybrid waves, which can be excited when a kinetic equilibrium is not initialized.²⁸ Equilibria are also what enable the study of the detailed structure of the pressure tensor in situations where finite Larmor radius effects play a significant role. The use of a fourth-order continuum kinetic solver facilitates the long-time accurate tracking of kinetic distribution functions and their moments over the course of many ion gyroperiods. The smooth solutions allow for the calculation of velocity moments and derivatives of variables to high-order accuracy and facilitate comparisons—such as the comparison between the kinetic stress tensor and the gyroviscous stress tensor—that shed light on the limitations of first-order approximations to kinetic physics. The KH instability evolves in such a way that first-order gyroviscous theory eventually becomes inapplicable. Specifically, the shear scale does not saturate until it has decreased to a value that is close to the ion Larmor orbit scale, at which point the assumptions of gyroviscous theory are invalid.

The combination of a noise-free high-order-accurate Vlasov solver and the ability to initialize exact equilibria provides a versatile and comprehensive means of investigating other gradient-driven instabilities that are modified by finite gyromotion. Here, we examine the lower-hybrid drift (LHD) instability^65–68 in a low-beta plasma wherein both ions and electrons have finite Larmor radii, with ion Larmor orbit size comparable to the gradient scale length. While the instability exists for a range of wavenumbers and can be influenced by electromagnetic physics,^67,69 we focus on the fastest-growing electrostatic mode of the instability for which $k⊥rLe$ is of order unity and $k·B=0$⁠.^66–68 

The LHD instability is typically driven by cross-field diamagnetic current in the presence of an inhomogeneity, e.g., nonuniform density. The instability converts drift energy into other forms, such that the relative drift between ions and electrons decreases as the instability evolves. The nonlinear state is characterized by fluctuations in density and electric field. The instability is of interest because it is a collisionless process that gives rise to anomalous resistivity⁶⁶—a phenomenon that is observed in a variety of magnetic confinement experiments and that is not predicted by the single-fluid magnetohydrodynamic model. The saturation mechanism for the instability has long been a subject of theoretical and computational investigations and remains an active area of research.^66,70–72 Here the focus is on how the methodology described in Sec. II broadens the scope of LHD instability physics that can be accessed through kinetic simulations.

The LHD instability is often studied in the context of a Harris current sheet,⁷³ which is a high-beta charge-neutral kinetic equilibrium in which density and magnetic field vary spatially and drift velocity is uniform for both ions and electrons. This configuration has been studied extensively via theoretical analysis^69,74,75 and kinetic simulations.^{69,72,76–78} While the Harris equilibrium has a convenient analytic form, it is not conducive to the study of low-beta plasmas or to the customization of density and/or velocity profiles. Leveraging the versatility of the methodology described in Sec. II for the study of the LHD instability, we consider a low-beta equilibrium in which the plasma transitions from a high-density region with constant number density n_A to a low-density region with a constant number density n_B. Such a configuration facilitates the study of LHD instabilities separate from kink instabilities, which are seeded when current-sheet configurations are perturbed.^69,72,76,77

While we can stipulate a variety of different electrostatic potential profiles, for simplicity, we consider a quasineutral equilibrium and thereby set $ϕ*=0$⁠. The equilibrium is constructed with desired ion number density,

where x₀ is a constant and the transition between n_A and n_B is a monotonically decreasing portion of a Gaussian on the subdomain $x∈[0,x0]$⁠. We could have also chosen g_i to be a hyperbolic-tangent profile, or any other profile that facilitates a smooth transition between two regions of different density and proceeded to successfully construct a kinetic equilibrium. The advantage of the piecewise form given in Eq. (18) is that it results in $fiaux$ [see Eq. (4)] whose velocity moments can be analytically evaluated, which makes it conducive to theoretical analysis. The choice of g_i also serves to highlight that the methodology of Sec. II allows kinetic equilibria to be specified in terms of different kinds of functions—including those that are piecewise defined. This means equilibria can be customized and tailored to a variety of different applications. Substituting Eq. (18) into Eq. (4) and evaluating Eq. (5) with $niaux(X)=gi(X)$ yields the auxiliary distribution function, which—like g_i—is piecewise defined, such that for each species,

where p₀ is a constant and $Hs=(Ms2πT)exp (−Msvx22T−Ms2T[(psy−x)Ωs]2)$ is the energy-dependent portion of the distribution function, with velocity v_y expressed in terms of the scaled canonical momentum p_sy. Evaluating $nsaux$ by taking the zeroth velocity moment of Eq. (19), substituting the result into Eq. (6), setting fixed values for ${T,p0,nA,nB,Ωp/ωpp}$⁠, numerically solving Eq. (6) for the equilibrium potential $ϕ$ subject to a homogeneous Neumann boundary condition at the left x-boundary and a homogeneous Dirichlet boundary condition at the right x-boundary, and substituting $ϕ$ into Eq. (7), yields the equilibrium distribution function for each species. The equilibrium distribution functions for ions and electrons are non-Maxwellian, the associated drift velocities are nonuniform, and gyromotion effects are fully encapsulated.

To set the equilibrium, we choose scalar parameters $T=6.25×10−4, p0=0.5$⁠, n_A = 1, $nB=exp (−5p02), Ωp/ωpp=0.1$⁠. On account of computational cost considerations, and as is common in kinetic simulations, an artificial mass ratio of $Mi/Me=25$ is used. For this set of parameters, plasma β is $2.5×10−3$ in the high-density region. Since the equilibrium potential $ϕ$ is nonuniform, the equilibrium plasma is not charge-neutral—it is, however, very close to charge-neutral since the maximum deviation between ion and electron number densities, i.e., $max(|ni−ne|/ni)$⁠, is about 2% and the mean deviation is 0.02%. The simulation parameters are chosen so that the peak relative drift velocity $Vd=uiy−uey$ is greater than the ion thermal speed v_Ti. For the chosen parameters, $max(|Vd|/vTi)≈1.43$⁠. This ensures that the growth rate⁶⁸ is large enough so that the initial stage of the instability can be observed over the course of several ion gyroperiods. The phase space simulation domain is defined as $[−1,1.5]×[0,Ly]×[−vmax,s,vmax,s]×∈[−vmax,s,vmax,s]$⁠, where $Ly=2π/15,vmax,i=0.23, vmax,e=1.0$⁠. The domain is periodic in the y direction and reflecting-wall boundary conditions are applied at x domain boundaries. Like the simulations described in Sec. III, the velocity domain is set to have zero-flux boundary conditions. In the simulation, the boundary conditions for the electrostatic potential are $∂ϕ/∂x=0$ at x = −1 and $ϕ=0$ at x = 1.5, consistent with the boundary conditions used to solve Eq. (6). The boundaries are situated at least one ion Larmor orbit away from the density variation to ensure that boundary conditions have a negligible effect on the dynamics. In the equilibrium, the temperature—computed from moments of the species distribution function—is nonuniform and consequently, the ion Larmor radius varies throughout the domain such that $rLi∈[0.236,0.267]$ with a mean value of 0.250 and the electron Larmor radius is $rLe∈[0.0497,0.0511]$ with a mean value of 0.0500.

The equilibrium is perturbed by multiplying the electron distribution function by a factor $(1+ϵ)$⁠, where

The exponential term in Eq. (20) ensures that the perturbation is localized near the region with maximum drift velocity, which helps to isolate the physics of interest and to avoid exciting additional dynamics.³ The wavenumbers $k1=15$ and $k2=2k1$ are chosen so that $k1≤k*≤k2$⁠, where $k*≈rLe−1$ is the wavenumber associated with maximum growth for the LHD instability, according to local-approximation linear theory.^66,68 Perturbing with wavenumbers k₁ and k₂ thus ensures that the fastest-growing mode can be captured in the simulation domain. The phase space grid resolution is $Nx×Ny×Nvx×Nvy=144×64×128×128$ cells for both ion and electron distribution functions. The simulation time step size is $Δt=3×10−3$⁠.

The initial stage of the LHD instability exhibits exponential growth of the seeded perturbation, with the electric field and potential amplitude growing in time. The structure for the electrostatic potential is shown in Fig. 5—both during the linear phase at time t = 144 and the nonlinear phase at t = 1920. Instability feature scales along the y direction are close to the electron gyroscale, with $kyrLe=0.75$⁠, as is consistent with local linear theory predictions for fast-growing modes.^66,68 The frequency associated with the temporal evolution of the domain-averaged electron momentum is $ω=0.42$⁠, which is close to the nominal nondimensional lower-hybrid frequency $ωLHωpp=((Ωi/ωpi)2+1(ωpe/Ωe)2+1)1/2$⁠, which has a value of 0.45 in the high-density region. Note that the assumption $(ωpe/Ωe)2≫1$ is not invoked in this definition of lower-hybrid frequency since it does not apply to the plasma under consideration.

The ion and electron distribution functions associated with the LHD instability evolve to have complex structures in phase space. The ion distribution function, in particular, exhibits particle-wave resonance features such that the distribution function strongly deforms in the vicinity of the phase velocity v_y = |v_ph| = |ω/k_y| = 0.028. As a caveat, this definition of phase velocity is applicable in the linear local approximation limit, which is never fully satisfied in the simulation, where frequencies exhibit spatiotemporal evolution, making it difficult to define a single characteristic phase velocity. Due to the magnetic field, the distribution function rotates in the (v_x, v_y) plane, and consequently resonance-induced features can be observed at positive and negative velocities. The ion distribution function $fi(x,vy)$ at fixed coordinates $(y=0.213,vx=1.8×10−3)$ and at different snapshots in time is shown in Fig. 6. The v_y coordinates corresponding to $±vph$ are shown in white dashed lines. Significant changes in the distribution function topology that lead to the creation of local extrema occur shortly after the end of the linear phase, around $t∈[200,400]$⁠, when the electric-field energy peaks. Such features tend to form near the phase velocity and where the magnitude of $∂f/∂vy$ is relatively large. Since continuum kinetic solvers, like the fourth-order finite-volume solver used here, can capture complex distribution function topologies, they provide a valuable means to develop and assess fluid models that encapsulate velocity moments beyond the first three, such as those explored in Refs. 79–81.

Over the course of the LHD instability, the relative drift velocity V_d decreases and the initial density profile becomes more diffuse. This process is shown in Fig. 7, where the y-averaged species number density $⟨ns〉$ and the y-averaged drift velocity to ion thermal speed ratio $⟨Vd/vTi⟩$ are plotted as a function of x position at different snapshots in time. The outward spread of the instability from the region with the largest initial density gradient is consistent with localized LHD instability behavior in current sheets.^69,76,77 However, unlike in current sheet configurations, wherein the ion density profile is almost unaffected by the instability,^77,78 here the ion and electron densities evolve together. Figure 7 shows that the instability leads to cross-field transport of plasma from the high-density region to the low-density region. The collisional analog to this process is the diffusive cross-field motion of the plasma due to resistive drag between ions and electrons. Since the LHD instability is a collisionless process, it is not characterized by smooth diffusion—instead the number density profile $⟨ns⟩$ changes in a stair-step fashion and exhibits x-direction spatial features that are approximately the size of an ion Larmor orbit. The evolution of the LHD instability slowly approaches a saturated nonlinear state, characterized by fixed-amplitude fluctuations.

Through the instability, a portion of the drift energy is converted into other forms of energy. The electric-field energy W_E, species thermal energy W_sT, species x-directed drift energy W_sx, and species y-directed drift energy W_sy are defined as

In Eqs. (22)–(24), n_s is evaluated using Eq. (3) and u_sj is evaluated using Eq. (12). Summing over species, the total thermal, x-directed drift, and y-directed drift energies are

For the electrostatic and magnetostatic simulation explored here, the total energy is defined as $Wtot=WE+WT+WX+WY$⁠. Note that since the magnetic field is constant in time and space, the magnetic field energy is not an energy source/sink and is omitted in the evaluation of total energy. By contrast, in high-beta current sheet configurations, the magnetic field energy can change significantly in response to the LHD instability and can be converted into thermal energy.⁷⁶ In computational studies of the LHD instability, PIC simulations are able to preserve total energy to within about 1% of the initial value.^69,72,76 By contrast, the VCK solver at the specified resolution is able to preserve total energy to within 0.002% of the initial value, which facilitates tracking and quantifying energy partitioning to a high degree of accuracy. Note that the degree to which energy is conserved in the finite-volume Vlasov solver depends on the resolution used—higher resolution leads to better energy conservation, consistent with the convergence properties of the algorithm.²³ 

The partitioning of energy for the LHD instability simulation is shown in Fig. 8, where for each type of energy W, the quantity $ΔW/Wtot$⁠, i.e., the net change $ΔW=W|t−W|t=0$ in the energy relative to the total energy W_tot, is plotted as a function of time. At initial time, $Wtot=8.557×10−4, WE=1.920×10−8, WT=8.483×10−4,WiT=4.204×10−4, WeT=4.278×10−4, WX=2.251×10−36$⁠, and $WY=7.424×10−6$⁠. Figure 8 shows that a portion of the initial drift energy is converted into electric-field energy and into fluctuations associated with x-directed motion. The bulk of the converted drift energy ultimately goes into thermal energy—specifically into ion thermal energy. The net change in the y-directed drift energy is primarily due to a decrease in W_iy. This partitioning is different from the energy partitioning observed in high-beta current-sheet configurations, in which the LHD instability leads primarily to magnetic field energy being converted into electron thermal energy.⁷⁶ The difference can be attributed to the fact that in the low-beta configuration studied here, the instability dynamics are electrostatic, whereas in current-sheet configurations electromagnetic physics can dictate dynamics.⁶⁹ The underlying physics that leads to a given partitioning is unknown and is left as a topic for future research. While the overall change in each type of energy is relatively small, i.e., less than 0.2% of the total energy, the instability's effect on plasma transport—shown in Fig. 7—is significant. As demonstrated, the noise-free Vlasov–Poisson solver and equilibrium initialization together offer a powerful means to explore the nonlinear saturated state of the LHD instability over long time intervals. The approach thereby provides advanced computational capability to probe cross-field anomalous transport physics.

A new methodology is deployed to facilitate high-fidelity modeling of kinetic physics in nonuniform low-beta plasmas, with a focus on configurations where ion gyroradii are comparable to gradient scale lengths. The approach combines a conservative fourth-order-accurate Vlasov–Poisson solver and a recently developed method for constructing exact customizable kinetic equilibria. The continuum kinetic solver provides enhanced-accuracy noise-free solutions, while customizable kinetic equilibria provide a means to isolate specific physical processes in kinetic simulations. The combined methodology is applied to simulate the shear-driven KH instability and the diamagnetic-drift-driven LHD instability in 4D $(x,y,vx,vy)$ phase space. The approach is shown to be effective for targeted physics investigations, for the detailed study of linear and nonlinear dynamics, and for capturing gyromotion physics in a variety of configurations.

The simulations are used to study the properties of the nondiagonal pressure tensor in KH instabilities where ion Larmor orbits are comparable to the shear length scale. During the nonlinear stage of the instability, when significant shear stresses develop, the stress tensor is shown to deviate substantially from the gyroviscous stress tensor. The deviation is quantified in terms of the local shear stress magnitude, which is related to the second invariant of the stress tensor. Local velocity shear relative to the ion cyclotron frequency is shown to be an important ratio, as is the evolving ratio of gradient scale to the ion Larmor orbit size, which tends to saturate near unity, thereby invalidating small Larmor orbit theory.

In the context of the LHD instability, the methodology is used to study a low-beta plasma configuration that smoothly transitions from a high-density to a low-density region. Such an equilibrium initialization is a stark departure from the well-studied high-beta Harris current sheet equilibrium and provides a means to investigate the LHD instability separate from the kink instability. The dynamics are shown to be consistent with local-approximation theoretical descriptions and are shown to lead to cross-field plasma transport. Distribution function features are captured in detail and particle-wave resonances are shown to influence the evolution of the ion distribution function. The ability of the solver to preserve the total energy within 0.002% of the initial condition allows for the detailed study of energy partitioning over the course of the instability. It is shown that the bulk of the converted drift energy goes into thermal energy.

Due to computational cost, instability simulations presented here were carried out at a fixed grid resolution. It is worth noting, however, that the resolution necessary to achieve a desired error tolerance is problem-dependent and can be determined through convergence studies. Rigorous convergence studies have been carried out in three phase space dimensions for the Vlasov–Poisson solver used here,^22,23 but are challenging for four-dimensional simulations since a factor of two increase in grid resolution in all directions would result in a factor of sixteen increase in the number of cells. Future continuum kinetic solver algorithms that leverage next-generation supercomputers are expected to overcome this challenge by enabling higher-resolution simulations in 4D, 5D, and 6D phase space. The ability to capture more degrees of freedom will facilitate rigorous convergence studies, detailed quantification of error, and will enable continuum kinetic simulations with larger dynamic range.

Through simulations of gradient-driven instabilities, it is demonstrated that a fourth-order continuum kinetic Vlasov–Poisson solver and the ability to initialize customizable kinetic equilibria together expand the scope of plasma dynamics that can be captured using kinetic simulations. The methodology opens the way for targeted high-accuracy computational studies of experimentally relevant plasma configurations—particularly those in which gyromotion plays an important role.

G.V.V. would like to thank U. Shumlak for helpful discussions. The kinetic simulations described herein used features of the Chombo library⁸² and the PETSc library⁸³ for parallelization. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52–07NA27344 and was supported in part by the LLNL-LDRD Program under Project No. 18-ERD-048.

This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.

The data that support the findings of this study are available from the corresponding author upon reasonable request.