Kinetic effects play a critical role in many problems in plasma physics. Correct description of kinetic processes typically requires solution of the Vlasov-Maxwell equations. Unfortunately, kinetic physics brings into the picture very small spatial (Debye length, electron gyroradius) and temporal (plasma frequency, electron gyrofrequency) scales, which could be several orders of magnitude smaller than system scales for most problems of interest. For example, in the Earth's magnetosphere the microscopic physics is key to magnetospheric dynamics during geomagnetic storms and substorms (Reeves et al., 2013; Jordanova et al., 2017). In the magnetosphere, the Debye length is of the order of ~ 1 − 100 m, while the system scales are ~10⁹ m. Similar or even more extreme separation of scales is typical of the solar wind and the solar corona, where kinetic physics is tied for example to mechanisms of energy release by magnetic reconnection and turbulence, impacting such fundamental problems as as the origin of coronal heating or solar wind heating and acceleration (e.g., Aschwanden, 2005; Bruno and Carbone, 2005). Despite remarkable recent progress in the development of efficient numerical methods for the solution of the Vlasov-Maxwell equations (for instance Chen et al., 2011; Markidis and Lapenta, 2011; Juno et al., 2018 and references therein), no existing fully-kinetic method appears capable of bridging the gap between microscopic and system scales given computational resources available now or in the foreseeable future.

Recognizing the inability of fully-kinetic methods to simulate large-scale systems, alternative approaches are being sought that focus on the so-called fluid-kinetic coupling. One line of research, appropriate in situation where kinetic effects are expected to be localized in space, seeks to embed a kinetic method in selected regions of a large-scale domain described by a fluid model (Sugiyama and Kusano, 2007; Kolobov and Arslanbekov, 2012; Daldorff et al., 2014; Tóth et al., 2016). For example, Particle-In-Cell (PIC) kinetic solvers have been successfully coupled to magnetohydrodynamics (MHD) large-scale solvers (Sugiyama and Kusano, 2007; Daldorff et al., 2014; Tóth et al., 2016). Another approach seeks extended fluid models by including high-order moments in the fluid hierarchy and seeking approximate closures capturing the physics of interest. Examples include the 5- and 10-moments methods discussed in Wang et al. (2015, 2018). Many other ideas on addressing fluid-kinetic coupling have been proposed in the literature (Markidis et al., 2014; Ho et al., unpublished¹). Furthermore, many well-established methods for solving Vlasov-Maxwell system, such as gyrokinetics (Brizard and Hahm, 2007), analytically remove certain spatial and temporal scales while retaining kinetic effects associated with larger/slower scales and thus may be suitable for certain class of problems involving fluid-kinetic coupling. However, despite the progress made, all of the existing methods have significant shortcomings and a comprehensive approach capable of obtaining accurate solutions in all parameters regimes remains elusive.

A new method for the solution of the Vlasov-Maxwell equations that could resolve the issues outlined above was recently proposed by Delzanno (2015). The method is based on a spectral discretization in velocity space and leads to a hierarchy of partial differential equations for the coefficients of the spectral expansion, akin to the moments of the fluid hierarchy. By a suitable choice of the spectral basis (in this case the asymmetrically-weighted Hermite functions), the low order moments of the expansion give a fluid description of the system. The kinetic physics is retained by adding more moments to the expansion. Thus, the fluid-kinetic coupling is an intrinsic feature of the method. One can immediately recognize that the method unifies and encompasses the fluid-kinetic approaches that have been described above. On the one hand, it can be seen as a reduced method that is not limited to any specific number of moments, but rather can be applied with as many moments as necessary for convergence, limited only by available computational resources. Furthermore, the number of moments could vary in the computational domain (even adaptively), in some analogy with kinetic-MHD-coupling approaches mentioned above. Significant advantage of the expansion method is that the transition between the “fluid” and “kinetic” regions can be as smooth as necessary. In all, spectral methods might offer on optimal way to perform large-scale simulations that include kinetic physics accurately.

In this paper, we describe a 3D parallel simulation code SpectralPlasmaSolver (SPS), which implements the method proposed in Delzanno (2015). The first application of the method to the problem of the three-dimensional turbulent cascade in the solar wind is demonstrated and the results are in good overall agreement with corresponding particle-in-cell simulations. The paper is organized as follows. In section 2 we briefly summarize the Hermite spectral method and its numerical formulation. In section 3, we discuss some aspects of the implementation, focusing on the parallelization for distributed-memory architectures. In section 4, linear properties of the methods and the simulations of decaying turbulence are discussed. Finally, a brief summary is presented in section 5.

We consider Vlasov-Maxwell equations for the evolution of the distribution function f_s(x, v, t) of a collisionless, magnetized plasma consisting of electrons (labeled by subscript s = e) and singly-charged ions (s = i) in a three-dimensional Cartesian domain in physical space x. The three-dimensional velocity space is represented by v. The Vlasov equation reads:

where t is time, ∇_{x, (v)} represents the gradient operator in physical (velocity) space and E (B) is the electric (magnetic) field. Equation (1) is written according to the following normalization: time is normalized to inverse of the electron plasma frequency ω_pe computed with a reference plasma density n₀, ωpe=e2n0ε0me, where e is the positive elementary charge, m_s is the mass of species s, ε₀ is the permittivity of vacuum, t → ω_pet; velocities are normalized to the speed of light c, v→vc; spatial coordinates are normalized to the electron inertial length d_e = c/ω_pe, x→xde; the magnetic field is normalized to a reference magnetic field B₀, B→BB0; the electric field is normalized to cB₀, E→EcB0; and the distribution function is normalized to n0c3, fs→fsc3n0. Furthermore, in Equation (1) q_s is the charge of each plasma species (qse=±1) and Ωcs=eB0ms is the corresponding cyclotron frequency defined to be always positive. Equation (1) is coupled to Maxwell's equations by the electromagnetic field:

We assume periodic boundary conditions in physical space and consider domain [0, L_x] × [0, L_y] × [0, L_z]. The distribution function in velocity space is assumed to go to zero at infinity.

The Vlasov-Maxwell Equations (1)–(3) are solved with a spectral discretization in both physical and velocity space. In particular, we expand the distribution function as

with N_{x, y, z} (N_{n, m, p}) the number of modes in physical (velocity) space along each spatial direction. Here index k refers to a set of three indices k = {k_x, k_y, k_z}.

The electric field is expanded as

and similarly for the magnetic field B. We adopt a Fourier decomposition in physical space

where δ = x, y, z. The decomposition in velocity space is obtained by asymmetrically-weighted (AW) Hermite basis functions, Ψn(x)=(π2nn!)-1/2Hn(x)e-x2, where H_n is the Hermite polynomial of the n-th order (Holloway, 1996). The argument of the basis functions is defined as ξδs=(vδ-uδs)/αδs, where uδs and αδs are constant shift and scaling parameters. The method could also be formulated to allow spatial and/or temporal variations in uδs and αδs.

By substituting expansion (4) into Equations (1)–(3) and using the orthogonality relations for both Fourier and Hermite functions, we can rewrite the Vlasov-Maxwell equations as a system of non-linearly coupled ordinary differential equations (ODE) for the coefficients of the expansion. These are (Delzanno, 2015):

for the Vlasov equation, and

for Maxwell's equations (written in index notation with summation over repeated indices implied). Here ε_βγδ is the Levi-Civita tensor and (a, b, c) = (1, 0, 0), (0, 1, 0), (0, 0, 1) for δ = x, y, z, respectively. Equations (7) are defined for n ∈ [0, N_n − 1], m ∈ [0, N_m − 1] and p ∈ [0, N_p − 1] and the convolution operator is defined as

where H is an arbitrary function. Note that Equation (7) couples a basis function (mode) indexed by n, m, p with successive modes of the hierarchy and it is necessary to specify how the system of equations is closed: We use Cn,m,pk,s=0 for n ≥ N_n, m ≥ N_m, p ≥ N_p. In order to address potential problems arising from the filamentation typical of collisionless plasmas, we add a collisional term on the right hand side of Eq. (7)

The collisional operator given by Equation (10) is not meant to model any specific physical dissipation mechanism, such as Coulomb collisions. Rather, it is added to damp higher-order Hermite modes. This motivates the specific functional form of the operator, together with the requirement that it does not act on the first three modes of the Hermite expansion, which have connections to density, momentum and energy of the system (Delzanno, 2015; Camporeale et al., 2016). Other forms of the operator could be used to model Coulomb collisions or to provide a physically-inspired closure (Loureiro et al., 2016).

Equations (7) and (8) can be cast in matrix form as

where C represents the vector of Cn,m,pkx,ky,kz,s coefficients and F contains the electric and magnetic field. Furthermore, L_1,2,3 represent the linear part of Equations (7) and (8) (advection in the case of Vlasov equations) while N(C,F) represents the nonlinear part with the convolutions.

System (11) is discretized in time with a fully-implicit Crank-Nicolson scheme (Crank and Nicolson, 1947). In residual form, we have

where Δt is the time step, and superscript θ labels time C(t^θ) = C(θΔt) = C^θ, and C^θ+1/2 = (C^θ+1+C^θ)/2 (and similarly for F), X=[CF] and ℝ=[R1R2].

The resulting numerical method formulated above and implemented in the SpectralPlasmaSolver (SPS) has several important properties. The main motivation for the choice of AW-Hermite spectral discretization of velocity space is the fact that the low order modes of the expansion correspond to the classic fluid description of the plasma (with a particular closure), while the kinetic physics is retained by adding more modes to the expansion. In other words, the fluid-kinetic coupling is an intrinsic feature of SPS. The Fourier discretization ensures that the constraints of Maxwell's equation (i.e., ∇·B = 0 and Poisson's equation) are automatically satisfied at all times if they are at t = 0. The fully-implicit Crank-Nicolson discretization improves the numerical stability of the method and ensures the exact conservation of total mass, momentum and energy in a finite time step. More details can be found in Delzanno (2015).

The implementation of SPS described here extends 2D Fortran 90 code for shared-memory architectures that was described in Vencels et al. (2016). The new code implements full 3D algorithm described in section 2 and is parallelized for distributed memory architectures using Message Passing Interface (MPI) calls. Similarly to the previous implementation, the new code utilizes Jacobian-free Newton-Krylov solver provided by PETSc library (Balay et al., 1997, 2017). Since PETSc supports parallel computations with MPI, minimal modifications were required to parallelize the solver.

Computation of convolutions is based on a parallel Fast Fourier Transforms (FFTs). The code utilizes 1D and 2D domain decomposition of the solution space, such that each domain owns a portion of k-space for all fields and for all coefficients C. The parallel FFTs are implemented using a modified version of 2DECOMP&FFT library (Li and Laizet, 2010). Our modifications to 2DECOMP&FFT ensure that padding with zeros, necessary for performing convolutions, is performed in an efficient manner and zeros are never communicated between processors. In three dimensions, this results in significant reduction in overall computational cost of the 3D transform relative to directly padding the 3D array. FFTW library (Frigo and Johnson, 2005) or Intel MKL² (via its FFTW interface) can be used as the backend to compute FFTs.

Figure 1 shows examples of strong (top panel) and weak (bottom panel) scalings obtained on Blue Waters supercomputer. The strong scaling was obtained for a system with 255³ Fourier modes and N_h = 4 Hermite modes in each direction. The weak scaling study was performed by proportionally scaling up the system size and the MPI rank count for a system with 63³ Fourier modes and N_h = 4 Hermite modes per direction. These tests were performed on XE nodes of Blue Waters supercomputer and used 8 MPI ranks per node. GCC compilers version 4.9.3 and FFTW library version 3.3.4 were used to compile the code. Each Blue Waters XE node is equipped with 2 AMD Interlagos model 6276 CPU processors. Significant variability in timings for a given simulation size was observed in some tests, presumably due to how a particular job was placed with respect to the topology of the communication links. The results shown in Figure 1 represent the lowest value. The break point observed in both strong and weak scalings is most likely related to the fact that simulation time is limited by communication time at high rank count. This is further confirmed by results shown in Figure 2, which demonstrates strong scaling of SPS on Skylake nodes of Electra supercomputer at NASA Ames research center. Each Skylake node is equipped with two Xeon Gold 6148 processors. Intel compilers and MKL library version 110302 were utilized in these tests. Furthermore, mpiprof tool³ was used to collect separate timings for computation and MPI communication. The scaling shown here was obtained for a system with 255³ spatial modes, 4³ Hermite functions, and initial conditions corresponding to Maxwellian plasma. It is apparent that saturation of the overall scaling is related to increasing relative cost of the communication with increasing number of processes.

These results point to a number of further enhancements in the code that would significantly increase parallelism and performance: (i) implementing decomposition in the Hermite coefficient space in addition to the k-space; (ii) overlapping communication and computation by using non-blocking collective MPI operations; (iii) performing transforms for several variables at a time.

Solar wind is believed to represent one of the best accessible examples of astrophysical large-scale plasma turbulence (Bruno and Carbone, 2005). The spectrum of magnetic fluctuations in the solar wind, as measured in-situ by spacecraft, exhibits power law behavior over many orders of magnitude in scales. A large range of scales, where fluctuations have predominantly Alfvenic character, exhibits many features expected of the inertial range of turbulence. This is followed by steepening to another power law range at proton kinetic scales kd_p ~ 1 (where d_p is the proton inertial length). The behavior at electron scales is less well characterized, but generally power law spectra are observed to extend to scales of the order of electron gyroradius, with subsequent steepening to either another power law or an exponential range reported in the literature (Sahraoui et al., 2010; Alexandrova et al., 2012).

In a classical paradigm of turbulence, the energy is transported by nonlinear interactions with no loss through the inertial range and is dissipated at small scales. In a weakly collisional plasma, such as solar wind, dissipation of the turbulent fluctuations is due to collective wave-particle interactions. These interactions can deposit the energy into various plasma components, such as thermal or suprathermal electrons, protons, and heavy ion species. Since turbulence has been proposed to play a significant role in the energy balance of solar wind, solar corona, and other space and astrophysical systems, understanding of the details of energy dissipation is a question of both significant theoretical interest and of great practical importance.

The problem of correctly describing energy dissipation in weakly collisional turbulence can be viewed as a problem of understanding fluid-kinetic coupling. Indeed, the fluctuations in the inertial range are widely thought to be well-described by MHD approximation. At the same time, adequate description of the dynamics at sub-proton scales and especially of the energy dissipation mechanisms requires electron and ion kinetic effects. The large-scale fluctuations drive dynamics at sub-proton scales, while the small-scale fluctuations may influence the large scales e.g., by supplying dissipation of energy and momentum, breaking topological constraints through magnetic reconnection, or placing constraints on values of temperature anisotropy enforced by instabilities. Significant progress has been achieved in numerical modeling of turbulence in weakly collisional plasmas (e.g., Howes, 2017), but the complexity of the problem drives interest in the methods capable of simulating sufficiently large range of scales, while describing kinetic effects with sufficient fidelity. In this section we describe an example of applying SPS to this problem. As the first step, we focus on the ability of the method to adequately describe kinetic dissipation and the small-scale dynamics.

In this paper we summarized recent development of a highly parallelized version of the code SPS. SPS is a spectral code that implements dual Fourier-Hermite transform to obtain numerical solution of the Vlasov-Maxwell system. The method is particularly suitable for describing fluid-kinetic coupling since a direct analogy could be drawn between traditional fluid hierarchy and the equations for the evolution of the coefficients of expansion in the case of asymmetrically-weighted Hermite basis. In contrast to the traditional multi-moment models, the number of coefficients in the expansion is limited only by available computational resources, which opens up intriguing possibilities for building models capable of adaptively increasing the number of coefficients in selected regions of physical or wavenumber space.

As an example, we have discussed application of SPS to the problem of turbulence in weakly collisional plasmas. A model of plasma dynamics capable of capturing essential features of turbulence at electron scales could be obtained with only a few Hermite modes, as was demonstrated by considering linear properties of Kinetic Alfvén waves. Fully nonlinear simulations of decaying plasma turbulence were performed and successfully compared against corresponding PIC simulations. While some differences in the decay of the imposed perturbation of ion flow were observed between PIC and SPS, SPS simulations successfully reproduce most of the statistical properties of magnetic and electron fluctuations observed in PIC simulations at kρ_e ≲ 1. Significant differences in the spectra are observed at scales approaching Debye scale and below. In practice, this range of scales might not be significant for turbulence investigations, since most of cascading energy is expected to dissipate before it reaches Debye scales.

The datasets for this study can be obtained by sending request to VR.

VR and GLD closely collaborated on development and validation of the version of SPS discussed here and preparation of the manuscript.