Splitting methods for Fourier spectral discretizations of the strongly magnetized Vlasov-Poisson and the Vlasov-Maxwell system

by   Jakob Ameres, et al.

Fourier spectral discretizations belong to the most straightforward methods for solving the unmagnetized Vlasov--Poisson system in low dimensions. In this article, this highly accurate approach is extended two the four-dimensional magnetized Vlasov--Poisson system with new splitting methods suited for strong magnetic fields. Consequently, a comparison to the asymptotic fluid model is provided at the example of a turbulent Kelvin--Helmholtz instability. For the three dimensional electromagnetic Vlasov--Maxwell system different novel charge conserving implementations of a Hamiltonian splitting are discussed and simulation results of the Weibel streaming instability are presented.



There are no comments yet.


page 7

page 12

page 13

page 14

page 15

page 16

page 24

page 26


Poisson Integrators based on splitting method for Poisson systems

We propose Poisson integrators for the numerical integration of separabl...

Hamiltonian Particle-in-Cell methods for Vlasov-Poisson equations

In this paper, Particle-in-Cell algorithms for the Vlasov-Poisson system...

Splitting integrators for stochastic Lie–Poisson systems

We study stochastic Poisson integrators for a class of stochastic Poisso...

A fast spectral method for electrostatics in doubly-periodic slit channels

We develop a fast method for computing the electrostatic energy and forc...

Efficient 6D Vlasov simulation using the dynamical low-rank framework Ensign

Running kinetic simulations using grid-based methods is extremely expens...

Stability and conservation properties of Hermite-based approximations of the Vlasov-Poisson system

Spectral approximation based on Hermite-Fourier expansion of the Vlasov-...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The Vlasov equation can be discretized following a Lagrangian or Eulerian approach. Lagrangian particle methods such as Particle in Cell have been dominant for a long time because they share the characteristics with actual physical particles and are easy to implement and parallelize [birdsall2004plasma]. Although their convergence rate is strongly limited by the Monte-Carlo approach and they based on moving large amounts of data in memory (the particles) it is possible to yield excellent conservation properties [kraus2016gempic]

. Semi-Lagrangian methods still use the particles for transporting the distribution function but yield higher convergence rates with an intermediate interpolation step using an Eulerian grid, yet their conservative form remains expensive 

[crouseilles2010conservative, crouseilles2014new]. There exists a variety of Eulerian Vlasov–Poisson solvers [filbet2001conservative] where lately geometric methods gained popularity [kraus2013variational]

. One of the simplest Eulerian solvers are pseudo-spectral solvers. They, of course, suffer from the curse of dimensionality but not on the computational level here, since the FFTW library is well optimized, see fig. 


. Constant coefficient advection in a periodic domain can be solved exactly in Fourier space. In all cases treated here, there is a Hamiltonian splitting available yielding constant advection possible. Fourier spectral solvers for the Vlasov equation, that employ also a Fourier transform in velocity space date back to 

[joyce1971numerical, izrar1989integration]. Such Fourier-Fourier solver were further developed for higher dimensions [eliasson2001outflow, eliasson2003numerical] and also extended to the Vlasov–Maxwell equation [eliasson2003numerical, eliasson2010numerical],[crouseilles2015hamiltonian]. It is also possible to use the Fourier basis as the interpolator underlying a Semi-Lagrangian scheme [FATONE2019349]. For Vlasov–Poisson it has been shown that Fourier filtering can be used to suppress the recurrence phenomenon [einkemmer2014strategy] or filter filamentations [klimas1994splitting]. For Vlasov–Poisson the Hamiltonian splitting has also been known [watanabe2005vlasov], but for Maxwell, none of these splitting methods is of geometric origin.
It should be mentioned that for the velocity space discretization also Chebyshev and Hermite polynomials have been used [joyce1971numerical, vencels2016spectralplasmasolver]. There the discretization by low degree Hermite polynomials provides an elegant way to approximate a fluid model on the numerical level.
A priori structure should be conserved for long terms and e.g. energy conservation is just a consequence but not the goal itself. Fourier spectral methods do not conserve positivity of the distribution function. In this context, we neglect the question on positivity conserving schemes although for other forms of discretizations there have been improvements in that direction [filbet2001conservative, filbet2003comparison, arber2002critical].
We begin by recalling the mechanisms for the Fourier spectral discretization of the two-dimensional Vlasov–Poisson and Vlasov–Ampère systems. In the second part, the four-dimensional magnetized Vlasov–Poisson system is obtained by the introduction of an external homogeneous magnetic field. There the Fourier spectral counterparts of known exponential splitting methods [knapp2015splitting] are presented and their performance is investigated under a stronger magnetic field with the use of the Kelvin Helmholtz instability. In the third part, we turn to electromagnetic physics by the means of the three dimensional Vlasov–Maxwell system, where the methods based on a Hamiltonian splitting are discussed at various test-cases following [kraus2016gempic]. The implementation in MATLAB used for the numerical examples can be found in a repository [FourierSpectralVlasovGITLAB].

Figure 1:

Fourier transforming a multidimensional array along one particular dimension yields a strided access pattern resulting in a slowdown. Timings are shown for forth- and back-transform in MATLAB (using FFTW) on a laptop. Although a slowdown is visible, it is not prohibitive for high dimensional spectral methods.

1.1 Vlasov–Poisson (1d1v)

We consider the one dimensional Vlasov equation (1)


and the Poisson equation


Here we Fourier transform in velocity and spatial space where denotes a transformation. For notational simplicity the transformed dimension is indicated by or in the argument. The spatial, velocity and fully Fourier transformed densities are defined as


where the wave vectors are

and for . Note that one can easily by a Fourier forth and back-transform switch between those three representations on a discrete level. We split the integration in three parts in , where the Vlasov steps can be integrated exactly in Fourier space.

  1. Advection in

  2. Advection in and Poisson solve


    Here we solve the Poisson equation with constant background (for ), but other fields are also possible.


For the splitting we consider the time to be one time step.

  1. Advection in in spatially transformed space
    Considering to be a fixed parameter the constant coefficient advection yields an ODE for each Fourier coefficient


    which can be solved exactly over this splitting step:

  2. Advection in in velocity transformed space


    Note that in this step the advection in cancels out under the velocity integral.


    Therefore, the electric field can be obtained in the spatially transformed space before or at the end of the split step.


The Lie steps can be composed by symmetric composition, see [watanabe2005vlasov]. The symplectic Runge Kutta scheme from Forest and Ruth [forest1989fourth] also works as it is just shifted by a half step and, therefore, adjoint symplectic for the Eulerian discretization.

1.2 Vlasov–Ampère (1d1v)

For the Vlasov–Ampère formulation the Poisson equation needs to be solved only once at , such that the electric field evolves in time by the Ampère equation


This leaves us with the following splitting:


The second split step (16) is now missing the Poisson equation but can be solved as before, whereas the second one (16) incorporates now the Ampère equation. It can be integrated exactly, since the solution to the constant coefficient advection is known to be which can be inserted into the Ampère equation reading


In spatially Fourier transformed space eqn. (17) can be solved by inserting the solution of the constant coefficient advection given in (10) as follows:


2 Magnetized Vlasov–Poisson (2d2v)

The magnetized Vlasov equation reads


which, reduced to two dimensions for and the magnetic field , reads


The canonical Hamiltonian splitting for the magnetized Vlasov–Poisson system reads


but has the disadvantage that spatial Fourier transform in and requires a convolution between and . We avoid this by separating the advection in each velocity component.


The Poisson equation in (24) and (25) has precisely the same solution for both split steps since the charge density actually stays constant over the advection and therefore, needs to be only solved once. If we take a look at the characteristics corresponding to (24)-(26),


we realize that the circular gyromotion for a strong magnetic field is not described very well, since it is split along each dimension. So we desire a spectral counterpart to more robust methods for strong magnetic fields like the exponential Boris algorithm [hongqin2013, knapp2015splitting].

2.1 Exponential splitting

The characteristics of the splitting underlying the exponential Boris algorithm reads


which leads us to the distribution counterpart


The exponential Boris scheme, along with many other integrators, use the fact that eqn. (31) can be solved exactly. With the two-dimensional rotation matrix


the solution to eqn. (31) reads


In eqn. (34) the spatial position is only a parameter such that the solution to (34) by the methods of characteristics reads


This corresponds to a rotation in the velocity plane for each position. In [paeth1986fast] the two-dimensional rotation matrix is decomposed into three shears:


Note that the two shears


Both shears, and correspond merely to a single dimensional advection and can be calculated in Fourier space using one dimensional transforms:


This is a commonly known method for image rotation by the discrete Fourier transform [larkin1997fast]. Contrary to splitting this rotation into two sub steps as in (24)-(25) the rotation by shearing is independent of the relation between time step and magnitude of the potentially strong magnetic field. Using the Taylor expansion for small , by approximating and we obtain the shears


which correspond exactly to the second order Strang splitting of eqn. (34), where the Lie steps read


This also explains why there is no visible difference between the Strang splitting and the exact shears in the third and fourth row in fig. 2. The rotation of an image multiples of can be implemented exactly by permutation involving transposing and flipping arrays, hence we can restrict the rotation in Fourier space on as recommended in [larkin1997fast] and also suggested by fig. 2.

Strang splitting
Strang splitting and reordering for multiples of
shearing and reordering for multiples of
Figure 2: Rotating an asymmetric two-dimensional Gaussian (Maxwellian) by Strang splitting and shearing, and with reordering respectively. Shearing in Fourier space for more than leads to heavy distortions (second row). Array rotations by multiples of can be implemented exactly by reordering such that shearing in Fourier space is only necessary for which leads to much better results (third and fourth row).

Another option is to use cubic B-spline interpolation for image rotation which is e.g. provided by imrotate in MATLAB. This corresponds to a backward Semi-Lagrangian discretization.
In the special case of a homogeneous magnetic field we can consider the splitting underlying Scovel’s method:


Note the following properties of the rotation matrix:


The exact solution of the characteristics in eqn. (44) reads then


If we suppose to be constant then the advection in velocity space (52) is independent of eqn. (48). Hence it is straightforward so solve (48) first and (52) thereafter. By spatial Fourier transform the advection in can be integrated exactly, which yields


Since only one dimensional Fourier transforms are used and the entire problem is two-dimensional for each

, aliasing can be suppressed by zero padding at small costs compared to padding the entire distribution function. The

rotation in (52) is the same as in eqn. (52) and hence can be discretized as before e.g. by Fourier transform in with rotation by shearing. For symmetric composition the adjoint method is needed, which means the rotation in has to be applied before the rotation and advection in . To account for the fact, that the rotation (52) is applied first, we rewrite eqn. (52) and (48) into:


The discrete counterpart of (55) reads then:


Since applying the adjoint after a forward time step with a negative time step () corresponds exactly to the identity map we obtain a symmetric method by combining the Scovel and the adjoint Scovel. For the discrete rotation by shearing in this is obvious, since . By combining (53) and (60) the in time symmetry is also easily verified:


Note that in the case of a constant homogeneous magnetic field the symmetrically composed Scovel coincides with the splitting presented in [EINKEMMER20142144] and the symmetric methods (28) and (29) in [knapp2015splitting].

2.2 Extension to Ampère

In the case of a homogeneous Maxwellian background one can solve the Ampère instead of the Poisson equation in order to obtain an update on the fields. In Fourier space this reads


Now any split step containing an advection in has to update the electric field according to eqn. (62). For the exponential Boris the only relevant split step is


such that the Ampère update for , following eqn. (18), reads


By this technique Gauss’ law is satisfied at any time. We recall that the electric field is obtained from the Poisson equation at any time as


For Scovel’s method it is not as straightforward, such that this shall be treated elsewhere.

2.3 Kelvin Helmholtz Instability

We consider a two-dimensional periodic domain with the lengths and the initial condition for the electrons


along with a constant ion background in the Poisson equation. In case of a strong magnetic field the dynamics of the fully kinetic model described by the Vlasov equation is very well approximated by the corresponding fluid model, a scaled version of the vorticity equation


coupled to the same fields stemming from the Poisson equation


The detailed scalings and techniques can be found in [golse1999vlasov, frenod2000long, filbet2016asymptotically, filbet2017asymptotically]. Here eqn. (67) is set on the kinetic time scale, but by introducing the fluid time scale with and a perturbed fluid density, the same as in [shoucri1981two] by removing the constant background . This allows us to use the results on the linear stability of the Kelvin Helmholtz instability derived in [shoucri1981two]. Depending on the wave number the growth rate on the fluid time scale in a periodic domain is by using a Taylor expansion approximated as


and on the kinetic time scale


In [shoucri1981two] the neutrally stable mode in the periodic domain was found in agreement to eqn. (72) at . Therefore, we install a linearly stable mode in the second dimension by and excite a linearly unstable mode in the first one with small amplitude in order to observe a Kelvin-Helmholtz instability with growth rate .
By rescaling to the fluid time scale we are able to compare the three different schemes from weak to strong magnetic field whilst holding the actual number of time steps constant. We focus on the unstable mode and there fore the electrostatic energy in the first dimension . The growth rate is only known for the limit , hence we do not expect agreement for small . But fig. 4 shows that we approach the fluid model with increasing . For all integrators show the same performance, but in the case of the strong magnetic field Scovel’s splitting is clearly better. For the standard and exponential Boris splitting fail entirely whereas Scovel’s splitting remains unaffected. For the exponential Boris reverts back to a lower frequency, which is a known effect from integrating particle trajectories [patacchini2008explicit], hence the steeper growth rate. This demonstrates that it is worthwhile to actually include the spatial rotation into the numerics.

rel. energy error
Figure 3: Electrostatic energy of the unstable mode and relative energy error in the Kelvin Helmholtz instability for increasing magnetic field strength . (See continuation)
rel. energy error
Figure 4: Electrostatic energy of the unstable mode and relative energy error in the Kelvin Helmholtz instability for increasing magnetic field strength , and . The actual number of time steps stays constant, which allows a comparison between the standard Strang splitting for the rotation.

In the following we extend our investigation into the nonlinear phase using the superior Scovel method in case of a weak and strong field in figures. 578 and 6. The entire system is driven by the unstable mode in , which deteriorates the stable mode in leading to an energy loss in especially in the nonlinear phase. This behavior is more pronounced in seems to be present for the stronger magnetic field. While the kinetic energy remains almost constant the only difference is the frequency of the oscillation which is directly linked to the gyro-frequency. Most importantly, fig. 5 shows that the energy error remains despite the nonlinear dynamics constant over long time. It is also slightly higher in the case of a strong field. This can be also seen in fig. 6, where the turbulence is much more pronounced for the strong field. The finer mode structure can also be seen in fig. 8 explaining the higher electrostatic energy. Fig. 7 indicates already that kinetic effects are only present in the weak case and a fluid model based on a Maxwellian velocity distribution is a fairly good approximation under a strong magnetic field.

(weak) (strong)
electrostatic energies
kinetic energies
relative energy error
Figure 5: Energies in the Kelvin Helmholtz instability under weak and strong magnetic field for with Scovel’s splitting. In both cases the stable mode looses energy with in the transition to the turbulent phase, where the system is driven by the unstabe mode . Although the dynamics are highly nonlinear the energy error remains constant.
(weak) (strong)
Figure 6: Projection of the phase space onto the spatial plane. For the weak case the initial mode structure is less pronounced as there is less confinement. The lack of confinement appears to introduce diffusion like effects. Nevertheless turbulence evolves in both cases at later times.
(weak) (strong)
Figure 7: Projection of the phase space onto the velocity plane. Due to the lack of confinement by the weak magnetic field kinetic effects smear out the Maxwellian. For the strong field the distribution resembles a sharp Maxwellian, such that a fluid model based on precisely that assumption describes the dynamics very well and kinetic effects are not essential.
(weak) (strong)
Figure 8: Projection of the phase space onto the plane in order to observe the kinetic structure of the unstable mode. For the strong magnetic field a pronounced turbulence in fig. 6 is observed. Here this leads to much finer mode structure in the reduced phase space for the strong case compared to the weak one.

3 Vlasov–Maxwell (1d2v)

We consider a reduction of the full six-dimensional Vlasov–Maxwell model onto one spatial and two velocity components. Elimination of the second and third spatial component, leaves us with two components of the electric field and one component of the magnetic field. Here the single magnetic component in -direction is denoted by .


For a density , the two components of the electric field and the magnetic field the reduced Vlasov equation is given in eqn. (75).


Dropping the species index yields the corresponding characteristics in eqn. (76).


The time dependent Maxwell equations reduce then to a system of three equations (77).


At the initialization for the Poisson eqn. (78) needs to be solved in order to obtain the first component of the electric field. The second component is always initialized as zero, .


Here we chose and consider only the electrons , with a constant ion background . The Hamiltonian splitting was already discussed extensively for Lagrangian particles [kraus2016gempic], nevertheless, it is also possible to derive the same method for a spectral discretization. For a different, but incorrect [qin2015comment], splitting this has already been done in [crouseilles2015hamiltonian]. Here we use the correct Hamiltonian splitting from [he2015hamiltonian]. Let denote the plasma density and the Fourier transform. Since there are six different combinations of transforms denotes a transformation, where the transformed dimension is indicated as before by , or in the argument. That means denotes the Fourier transform of in and . We begin by treating the Hamiltonian splitting for time integration from to .

  • Kinetic energy ,