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 . 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
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].

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)
(1) |
and the Poisson equation
(2) |
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
(3) | ||||
(4) | ||||
(5) |
where the wave vectors are
-
Advection in
(6) -
Advection in and Poisson solve
(7) Here we solve the Poisson equation with constant background (for ), but other fields are also possible.
(8)
For the splitting we consider the time to be one time step.
-
Advection in in spatially transformed space
Considering to be a fixed parameter the constant coefficient advection yields an ODE for each Fourier coefficient(9) which can be solved exactly over this splitting step:
(10) -
Advection in in velocity transformed space
(11) Note that in this step the advection in cancels out under the velocity integral.
(12) Therefore, the electric field can be obtained in the spatially transformed space before or at the end of the split step.
(13)
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
(14) |
This leaves us with the following splitting:
(15) | |||
(16) |
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
(17) |
In spatially Fourier transformed space eqn. (17) can be solved by inserting the solution of the constant coefficient advection given in (10) as follows:
(18) |
2 Magnetized Vlasov–Poisson (2d2v)
The magnetized Vlasov equation reads
(19) |
which, reduced to two dimensions for and the magnetic field , reads
(20) |
The canonical Hamiltonian splitting for the magnetized Vlasov–Poisson system reads
(21) | |||
(22) | |||
(23) |
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.
(24) | |||
(25) | |||
(26) |
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),
(27) | ||||
(28) |
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
(29) | |||
(30) | |||
(31) |
which leads us to the distribution counterpart
(32) | |||
(33) | |||
(34) |
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
(35) |
the solution to eqn. (31) reads
(36) |
In eqn. (34) the spatial position is only a parameter such that the solution to (34) by the methods of characteristics reads
(37) |
This corresponds to a rotation in the velocity plane for each position. In [paeth1986fast] the two-dimensional rotation matrix is decomposed into three shears:
(38) |
Note that the two shears
(39) |
Both shears, and correspond merely to a single dimensional advection and can be calculated in Fourier space using one dimensional transforms:
(40) | ||||
(41) |
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
(42) |
which correspond exactly to the second order Strang splitting of eqn. (34), where the Lie steps read
(43) |
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 | |||
![]() |
![]() |
![]() |
![]() |
shearing | |||
![]() |
![]() |
![]() |
![]() |
Strang splitting and reordering for multiples of | |||
![]() |
![]() |
![]() |
![]() |
shearing and reordering for multiples of | |||
![]() |
![]() |
![]() |
![]() |
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:
(44) | |||
(45) |
Note the following properties of the rotation matrix:
(46) |
(47) |
The exact solution of the characteristics in eqn. (44) reads then
(48) | ||||
(49) | ||||
(50) | ||||
(51) | ||||
(52) |
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
(53) |
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
(54) | ||||
(55) | ||||
(56) | ||||
(57) | ||||
(58) | ||||
(59) |
The discrete counterpart of (55) reads then:
(60) |
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:
(61) |
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
(62) |
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
(63) |
such that the Ampère update for , following eqn. (18), reads
(64) |
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
(65) |
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
(66) |
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
(67) |
coupled to the same fields stemming from the Poisson equation
(69) | ||||
(70) | ||||
(71) |
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
(72) |
and on the kinetic time scale
(73) |
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 | |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
rel. energy error | |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
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. 5, 7, 8 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 | |
![]() |
![]() |
(weak) | (strong) |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
(weak) | (strong) |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
(weak) | (strong) |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
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 .
(74) |
For a density , the two components of the electric field and the magnetic field the reduced Vlasov equation is given in eqn. (75).
(75) |
Dropping the species index yields the corresponding characteristics in eqn. (76).
(76) |
The time dependent Maxwell equations reduce then to a system of three equations (77).
(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, .
(78) |
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 ,
Comments
There are no comments yet.