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

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)

\partial_{t} f (x, v, t) + v \partial_{x} f (x, v, t) + \frac{q}{m} (E (x, t) + E_{e x t} (x, t)) \partial_{v} f (x, v, t) = 0

(1)

and the Poisson equation

- \partial_{x x} Φ (x, t) = 1 + q \int_{- \infty}^{\infty} f (x, v, t) d v, E (x, t) = - \partial_{x} Φ (x, t)

(2)

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

$^f (k_{x}, v, t)$	$= \frac{1}{L} \int_{0}^{L} f (x, v, t) exp - i x k_{x} d x,$	(3)
$^f (x, k_{v}, t)$		(4)
$^f (k_{x}, k_{v}, t)$		(5)

where the wave vectors are

k_{x} = n \frac{2 π}{L}

and

k_{v} = \frac{2 π}{v_{max} - v_{min}}

for

n \in Z

. 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

τ [0, t]

, where the Vlasov steps can be integrated exactly in Fourier space.

Advection in $x$

$\partial_{t} f (x, v, t) + v \partial_{x} f (x, v, t) = 0$ (6)
Advection in $v$ and Poisson solve

$\partial_{t} f (x, v, t) + \frac{q}{m} (E (x, 0) + E_{e x t} (x, 0)) \partial_{v} f (x, v, t) = 0$ (7)

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

$\partial_{x} E (x, t) = 1 + q \int f (x, v, t) d v$ (8)

For the splitting we consider the time $[0, t]$ to be one time step.

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

$\partial_{t}^f (k_{x}, v, t) = - v i k_{x}^f (k_{x}, v, t),$ (9)

which can be solved exactly over this splitting step:

(10)
Advection in $v$ in velocity transformed space

(11)

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

$\begin{matrix} \int_{R} f (x, v, t) d v & = \int_{R} f (x, v + t \frac{q}{m} [E (x, 0) + E_{e x t} (x, 0)], 0) d v = \int_{R} f (x, v, 0) d v \end{matrix}$ (12)

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

$^E (k_{x}, 0) = q \frac{1}{i k_{x}} \int^f (k_{x}, v, 0) d v, for k_{x} \neq 0$ (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 $t = 0$ , such that the electric field evolves in time by the Ampère equation

\partial_{t} E (x, t) = - j (x, t) = - \sum s \frac{q}{m} \int v f_{s} (x, v) d v .

(14)

This leaves us with the following splitting:

	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + v \partial_{x} f_{s} (x, v, t) = 0 \partial_{t} E (x, t) = - \sum_{s} \frac{q}{m} \int v f_{s} (x, v, t) d v \end{matrix}$		(15)
	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + \frac{q}{m} (E (x, t) + E_{e x t} (x, t)) \partial_{v} f_{s} (x, v, t) = 0 \end{matrix}$		(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 $f_{s} (x, v, t) = f_{s} (x - t v, 0)$ which can be inserted into the Ampère equation reading

\begin{matrix} E (x, t) & = E (x, 0) - \sum s \frac{q}{m} \int_{0}^{t} \int v f_{s} (x, v, t) d v d τ = E (x, 0) - \sum s \frac{q}{m} \int_{0}^{t} \int v f_{s} (x - τ v, v, 0) d v d τ . \end{matrix}

(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

\partial_{t} f_{s} (x, v, t) + v \cdot \nabla_{x} f_{s} (x, v, t) + \frac{q}{m} [E (x, t) + v \times B (x, t)] \cdot \nabla_{v} f (x, v, t) = 0

(19)

which, reduced to two dimensions for $x = (x_{1}, x_{2}), v = (v_{1}, v_{2})$ and the magnetic field $B (x, t) = (0, 0, B_{3} (x, t))$ , reads

\partial_{t} f_{s} (x, v, t) + v \cdot \nabla_{x} f_{s} (x, v, t) + + \frac{q}{m} [(E_{1} (x, t) + v_{2} B_{3} (x, t)) \partial_{v_{1}} f (x, v, t) + (E_{2} (x, t) - v_{1} B_{3} (x, t)) \partial_{v_{2}} f (x, v, t)] = 0.

(20)

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

	$H_{E} {\begin{matrix} \partial_{t} f_{s} (x, v, t) + \frac{q}{m} E (x, t) \cdot \nabla_{v} f_{s} (x, v, t) = 0 E (x, t) = - \nabla Φ (x, t), - Δ Φ (x, t) = \sum_{s} q_{s} \int f_{s} (x, v, t) d v, \end{matrix}$		(21)
	$H_{p_{1}} {\begin{matrix} \partial_{t} f_{s} (x, v, t) + v_{1} \partial_{x_{1}} f_{s} (x, v, t) - \frac{q}{m} v_{1} B_{3} (x, t) \partial_{v_{2}} f (x, v, t) = 0, \end{matrix}$		(22)
	$H_{p_{2}} {\begin{matrix} \partial_{t} f_{s} (x, v, t) + v_{2} \partial_{x_{2}} f_{s} (x, v, t) + \frac{q}{m} v_{2} B_{3} (x, t) \partial_{v_{1}} f (x, v, t) = 0, \end{matrix}$		(23)

but has the disadvantage that spatial Fourier transform in $H_{p_{1}}$ and $H_{p_{1}}$ requires a convolution between $B_{3}$ and $f$ . We avoid this by separating the advection in each velocity component.

	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + \frac{q}{m} [E_{1} (x, t) + v_{2} B_{3} (x, t)] \partial_{v_{1}} f_{s} (x, v, t) = 0 E_{1} (x, t) = - \partial_{x_{1}} Φ (x, t), - Δ Φ (x, t) = \sum_{s} q_{s} \int f_{s} (x, v, t) d v, \end{matrix}$		(24)
	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + \frac{q}{m} [E_{2} (x, t) - v_{1} B_{3} (x, t)] \partial_{v_{2}} f_{s} (x, v, t) = 0 E_{2} (x, t) = - \partial_{x_{1}} Φ (x, t), - Δ Φ (x, t) = \sum_{s} q_{s} \int f_{s} (x, v, t) d v, \end{matrix}$		(25)
	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + v \cdot \nabla_{x} f_{s} (x, v, t) = 0 \end{matrix}$		(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),

	$⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} _{1} (t) & = \frac{q}{m} [E_{1} (X (t), t) + V_{2} (t) B_{3} (X (t), t)]_{2} (t) & = 0 ˙ X (t) & = 0 \end{matrix}$			(27)
	$⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} _{1} (t) & = 0_{2} (t) & = \frac{q}{m} [E_{2} (X (t), t) - V_{1} (t) B_{3} (X (t), t)] ˙ X (t) & = 0 \end{matrix}$	$⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} _{1} (t) & = 0_{2} (t) & = 0 ˙ X (t) & = V (t) \end{matrix}$		(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

	${\begin{matrix} ˙ V (t) & = 0 ˙ X (t) & = V (t) \end{matrix}$		(29)
	${\begin{matrix} ˙ V (t) & = \frac{q}{m} E (X (t), t) ˙ X (t) & = 0 \end{matrix}$		(30)
	$⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ \begin{matrix} _{1} (t) & = \frac{q}{m} V_{2} (t) B_{3} (X (t), t)_{2} (t) & = - \frac{q}{m} V_{1} (t) B_{3} (X (t), t) ˙ X (t) & = 0 \end{matrix}$		(31)

which leads us to the distribution counterpart

	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + \frac{q}{m} E (x, t) \cdot \nabla_{v} f_{s} (x, v, t) = 0, \end{matrix}$		(32)
	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + v \nabla_{x} f_{s} (x, v, t) = 0, \end{matrix}$		(33)
	${\begin{matrix} \partial_{t} f_{s} (x, v, t) + \frac{q}{m} [v_{2} B_{3} (x, t) \partial_{v_{1}} f (x, v, t) - v_{1} B_{3} (x, t) \partial_{v_{2}} f (x, v, t)] = 0 . \end{matrix}$		(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

R (θ) = (\begin{matrix} cos (θ) & - sin (θ) sin (θ) & cos (θ) \end{matrix}),

(35)

the solution to eqn. (31) reads

V (t) = R (- θ) V (0) for θ = t \frac{q}{m} B_{3} (X (0), 0) .

(36)

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

f_{s} (x, v, t) = f_{s} (x, R (- θ) \cdot v, 0), θ = t \frac{q}{m} B_{3} (x, 0) .

(37)

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

R(θ)=(cos(θ)−sin(θ)sin(θ)cos(θ))=(1−tan(\nicefracθ2)01)(10sin(θ)1)(1−tan(\nicefracθ2)01)

(38)

Note that the two shears

S_{1} (α) = (\begin{matrix} 1 & α 0 & 1 \end{matrix}) and S_{2} (β) = (\begin{matrix} 1 & 0 β & 1 \end{matrix})

(39)

Both shears, $S_{1}$ and $S_{2}$ correspond merely to a single dimensional advection and can be calculated in Fourier space using one dimensional transforms:

	$F_{v_{1}} [f (x, S_{1} (α) \cdot v)]$			(40)
	$F_{v_{2}} [f (x, S_{2} (β) \cdot v)]$			(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 $sin (x) \approx x$ and $tan (x) \approx x$ we obtain the shears

R(−θ)=(cos(θ)sin(θ)−sin(θ)cos(θ))≈(1\nicefracθ201)(10−θ1)(1\nicefracθ201),

(42)

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

{\partial_{t} f_{s} + \frac{q}{m} v_{2} B_{3} \partial_{v_{1}} f = 0 and% {\partial_{t} f_{s} - \frac{q}{m} v_{1} B_{3} \partial_{v_{2}} f = 0.

(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 $\nicefracπ2$ can be implemented exactly by permutation involving transposing and flipping arrays, hence we can restrict the rotation in Fourier space on $θ \in [- \frac{π}{4}, \frac{π}{4}]$ as recommended in [larkin1997fast] and also suggested by fig. 2.

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 $90^{\circ}$ leads to heavy distortions (second row). Array rotations by multiples of $90^{\circ}$ can be implemented exactly by reordering such that shearing in Fourier space is only necessary for $- 45^{\circ} \leq θ \leq 45^{\circ}$ 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:

	${\begin{matrix} ˙ V (t) & = \frac{q}{m} E (X (t), t) ˙ X (t) & = 0 \end{matrix}$		(44)
	$⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ \begin{matrix} _{1} (t) & = \frac{q}{m} V_{2} (t) B_{3}_{2} (t) & = - \frac{q}{m} V_{1} (t) B_{3} ˙ X (t) & = V (t) \end{matrix}$		(45)

Note the following properties of the rotation matrix:

ddθR(θ−\nicefracπ2)=R(θ) and R(θ−\nicefracπ2)=(01−10)

(46)

R (θ) = (\begin{matrix} cos (θ) & - sin (θ) sin (θ) & cos (θ) \end{matrix}), R (θ)^{- 1} = (\begin{matrix} cos (θ) & sin (θ) - sin (θ) & cos (θ) \end{matrix}) = R (- θ) .

(47)

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

$X (t)$	$= X (0) + \int_{0}^{t} ˙ X (τ) d τ = X (0) + \int_{0}^{t} V (τ) d τ$	(48)
	$= X (0) + \int_{0}^{t} R (- τ \frac{q}{m} B_{3} (X (0), 0)) V (0) d τ$	(49)
	$=X(0)+mqsB3[R(−θ+\nicefracπ2)−R(0)]V(0)$	(50)
	$= X (0) + \frac{m}{q_{s} B_{3}} (\begin{matrix} sin (θ) & 1 - cos (θ) cos (θ) - 1 & sin (θ) \end{matrix}) (\begin{matrix} V_{1} (0) V_{2} (0) \end{matrix}),$	(51)
$V (t)$	$= R (- θ) V (0), for θ = t \frac{q}{m} B_{3} (X (0), 0) .$	(52)

If we suppose $B_{3}$ 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 $x$ can be integrated exactly, which yields

(53)

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

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

v \times B

rotation in (52) is the same as in eqn. (52) and hence can be discretized as before e.g. by Fourier transform in

v

with rotation by shearing. For symmetric composition the adjoint method is needed, which means the rotation in

v

has to be applied before the rotation and advection in

x

. To account for the fact, that the rotation (52) is applied first, we rewrite eqn. (52) and (48) into:

$V (0)$	$= R (- θ)^{- 1} V (t) = R (θ) V (t), for θ = t \frac{q}{m} B_{3} (X (0), 0),$	(54)
$X (t)$	$= X (0) + \int_{0}^{t} V (τ) d τ$	(55)
	$= X (0) + \int_{0}^{t} R (- τ \frac{q}{m} B_{3} (X (0), 0)) V (0)      = R (θ) V (t) d τ$	(56)
	$=X(0)+mqsB3[R(−θ+\nicefracπ2)−R(0)]R(θ)V(t)$	(57)
	$= X (0) + \frac{m}{q_{s} B_{3}} (\begin{matrix} sin (θ) & 1 - cos (θ) cos (θ) - 1 & sin (θ) \end{matrix}) (\begin{matrix} cos (θ) & - sin (θ) sin (θ) & cos (θ) \end{matrix}) (\begin{matrix} V_{1} (0) V_{2} (0) \end{matrix})$	(58)
	$= X (0) + \frac{m}{q_{s} B_{3}} (\begin{matrix} sin (θ) & cos (θ) - 1 1 - cos (θ) & sin (θ) \end{matrix}) (\begin{matrix} V_{1} (0) V_{2} (0) \end{matrix}) .$	(59)

The discrete counterpart of (55) reads then:

(60)

Since applying the adjoint after a forward time step with a negative time step ( $φ (Δ t) φ^{*} (- Δ t) = I d$ ) 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 $v$ this is obvious, since $R (θ) R (- θ) = I$ . 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

{^E}_{j} (k_{x}, t) = {^E}_{j} (k_{x}, 0) + \int_{0}^{t} \iint v_{j}^f (k_{x}, v, τ) d (v_{1}, v_{2}) d τ, j = 1, 2.

(62)

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

{\begin{matrix} ˙ V (t) & = 0 ˙ X (t) & = V (t) \end{matrix}

(63)

such that the Ampère update for $j = 1, 2$ , 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

{^E}_{j} (k_{x}, t) = \frac{i k_{x_{j}}}{k_{x} \cdot k_{x}} q \int^f (k_{x}, v, t) d v .

(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 $Ld=\nicefrac2πkd, d=1,2$ and the initial condition for the electrons $(q_{e} = - 1, m_{e} = 1)$

f_{e} (x_{1}, x_{2}, v_{1}, v_{2}, t = 0) = (1 + sin (k_{2} x_{2}) + ν cos (k_{1} x_{1})) \frac{1}{2 π} exp (\frac{v_{1}^{2} + v_{2}^{2}}{2})

(66)

along with a constant ion background $ρ_{i} = 1, q_{i} = 1$ in the Poisson equation. In case of a strong magnetic field $B_{3}$ 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

\frac{\partial_{t} ρ_{e} (x, t)}{∥ B_{3} ∥} + q_{i} [E_{1} (x, t) \partial_{x_{2}} ρ_{e} (x, t) - E_{2} (x, t) \partial_{x_{2}} ρ_{e} (x, t)] = 0

(67)

coupled to the same fields stemming from the Poisson equation

$- Δ Φ (x, t)$	$= ρ_{i} + ρ_{e} (x, t),$	(69)
$E (x, t)$	$= - \nabla Φ (x, t),$	(70)
$ρ_{e} (x, t)$	$= q_{e} \int f_{e} (x, v, t) d v .$	(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 $t^{ψ} = \frac{t}{∥ B_{3} ∥}$ and a perturbed fluid density, the same as in [shoucri1981two] by removing the constant background $ψ (x, t^{ψ}) = [ρ_{i} + ρ_{e} (x, ∥ B_{3} ∥ t^{ψ})]$ . This allows us to use the results on the linear stability of the Kelvin Helmholtz instability derived in [shoucri1981two]. Depending on the wave number $k$ the growth rate on the fluid time scale in a periodic domain is by using a Taylor expansion approximated as

ω^{ψ} = i 2 (1 - k) k

(72)

and on the kinetic time scale

ω = \frac{i 2 (1 - k) k}{∥ B_{3} ∥} .

(73)

In [shoucri1981two] the neutrally stable mode in the periodic domain was found in agreement to eqn. (72) at $k = 1$ . Therefore, we install a linearly stable mode in the second dimension by $k_{2} = 1$ and excite a linearly unstable mode in the first one $k_{1} = 0.4$ with small amplitude $ν = 0.015$ in order to observe a Kelvin-Helmholtz instability with growth rate $ω = i 0.480$ .
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 $\frac{1}{2} | | E_{1} | |^{2}$ . The growth rate is only known for the limit $| | B_{3} | | \to \infty$ , hence we do not expect agreement for small $| | B_{3} | |$ . But fig. 4 shows that we approach the fluid model with increasing $| | B_{3} | |$ . For $| | B_{3} | | = 1$ all integrators show the same performance, but in the case of the strong magnetic field Scovel’s splitting is clearly better. For $| | B_{3} | | = 16$ the standard and exponential Boris splitting fail entirely whereas Scovel’s splitting remains unaffected. For $| | B_{3} | | = 32$ 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.

Figure 3: Electrostatic energy of the unstable mode and relative energy error in the Kelvin Helmholtz instability for increasing magnetic field strength $| | B | |$ . (See continuation)

Figure 4: Electrostatic energy of the unstable mode and relative energy error in the Kelvin Helmholtz instability for increasing magnetic field strength $| | B | |$ , $N_{x} = 32, N_{v} = 32, Δ t^{ψ} = 0.01$ and $t_{max}^{ψ} = 40$ . 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. 5, 7, 8 and 6. The entire system is driven by the unstable mode in $x 1$ , which deteriorates the stable mode in $x 2$ leading to an energy loss in $E_{2}$ 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.

Figure 5: Energies in the Kelvin Helmholtz instability under weak and strong magnetic field for $N_{x} = N_{v} = 32, Δ t^{ψ} = 0.01, t_{max}^{ψ} = 100$ with Scovel’s splitting. In both cases the stable mode $E_{2}$ looses energy with in the transition to the turbulent phase, where the system is driven by the unstabe mode $E_{1}$ . Although the dynamics are highly nonlinear the energy error remains constant.

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.

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.

Figure 8: Projection of the phase space onto the $(x_{1}, v_{1})$ 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 $z$ -direction is denoted by $B_{3}$ .

x = x_{1}, v = (v_{1}, v_{2}), E = (E_{1}, E_{2}), B = B_{3}

(74)

For a density $f (x, v_{1}, v_{2}, t)$ , the two components of the electric field $E_{1} (x, t), E_{2} (x, t)$ and the magnetic field $B (x, t)$ the reduced Vlasov equation is given in eqn. (75).

\partial_{t} f_{s} + v_{1} \partial_{x} f_{s} + \frac{q}{e} \frac{m_{e}}{m} [E_{1} \partial_{v_{1}} f_{s} + E_{2} \partial_{v_{2}} f_{s} + B (v_{2} \partial_{v_{1}} f_{s} - v_{1} \partial_{v_{2}} f_{s})] = 0

(75)

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

\begin{matrix} \frac{d}{d t} V_{1} (t) & = \frac{q}{e} \frac{m_{e}}{m} [E_{1} (X_{s} (t), t) + V_{2} (t) B (X (t), t)] \frac{d}{d t} V_{2} (t) & = \frac{q}{e} \frac{m_{e}}{m} [E_{2} (X_{s} (t), t) - V_{1} (t) B (X (t), t)] \frac{d}{d t} X (t) & = V_{1} (t) \end{matrix}

(76)

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

\begin{matrix} \partial_{t} E_{1} (x, t) & = - \sum s \frac{q}{e} \int v_{1} f_{s} (x, v_{1}, v_{2}, t) d v \partial_{t} E_{2} (x, t) & = - {(\frac{c}{v_{t h, e}})}^{2} \partial_{x} B (x, t) - \sum s \frac{q}{e} \int v_{2} f_{s} (x, v_{1}, v_{2}, t) d v, \partial_{t} B (x, t) & = - \partial_{x} E_{2} (x, t) \end{matrix}

(77)

At the initialization for $t = 0$ the Poisson eqn. (78) needs to be solved in order to obtain the first component $E_{1}$ of the electric field. The second component is always initialized as zero, $E_{2} (x, 0) = 0$ .

- \partial_{x x} Φ (x, t) = \sum s \frac{q}{e} \int_{R^{d}} f_{s} (x, v_{1}, v_{2}, t) d v, E_{1} (x, t) = - \partial_{x} Φ (x, t)

(78)

Here we chose $\frac{c}{v_{t h, e}} = 1$ and consider only the electrons $f = f_{e}, q = e m = m_{e}$ , with a constant ion background $f_{i} = exp (- \frac{v_{1}^{2} + v_{2}^{2}}{2})$ . 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 $f (x, v_{1}, v_{2}, t)$ denote the plasma density and $^f$ the Fourier transform. Since there are six different combinations of transforms $^f$ denotes a transformation, where the transformed dimension is indicated as before by $k_{x}$ , $k_{v_{1}}$ or $k_{v_{2}}$ in the argument. That means $^f (k_{x}, v_{1}, k_{v_{2}})$ denotes the Fourier transform of $f$ in $x$ and $v_{2}$ . We begin by treating the Hamiltonian splitting for time integration from $0$ to $t$ .

Kinetic energy $(d = 1)$ , $H_{p_{1}} = \frac{1}{2} ∭ v_{1}^{2} f (x, v, t) d x d v_{1} d v_{2}$

$\begin{matrix} \partial_{t} f (x, v_{1}, v_{2}, t) + v_{1} \partial_{x} f (x, v_{1}, v_{2}, t) - \frac{q}{m} B_{3} (x, t) v_{1} \partial_{v_{2}} f (x, v_{1}, v_{2}, t) & = 0 \partial_{t} B_{3} (x, t) & = 0 \partial_{t} E_{1} (x, t) = - \end{matrix}$

Strang splitting $(120^{\circ} \leq θ \leq 180^{\circ})$

shearing $(120^{\circ} \leq θ \leq 180^{\circ})$

Strang splitting $(- 45^{\circ} \leq θ \leq 45^{\circ})$ and reordering for multiples of $90^{\circ}$

shearing $(- 45^{\circ} \leq θ \leq 45^{\circ})$ and reordering for multiples of $90^{\circ}$

$\frac{1}{2} \| \| E_{1} \| \|^{2}$	rel. energy error
$\| \| B_{3} \| \| = 1$

$\| \| B_{3} \| \| = 2$

$\| \| B_{3} \| \| = 4$

$\| \| B_{3} \| \| = 1$ (weak)	$\| \| B_{3} \| \| = 32$ (strong)
electrostatic energies

kinetic energies

relative energy error

$Δ t = 0.01, t_{max} = 100$	$Δ t = 0.32, t_{max} = 3200$

$\| \| B_{3} \| \| = 1$ (weak)	$\| \| B_{3} \| \| = 32$ (strong)
$t^{ψ} = 18$

$t^{ψ} = 58$

$t^{ψ} = 98$

$\| \| B_{3} \| \| = 1$ (weak)	$\| \| B_{3} \| \| = 32$ (strong)
$t^{ψ} = 18$

$t^{ψ} = 58$

$t^{ψ} = 98$

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

Authors

Poisson Integrators based on splitting method for Poisson systems

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

Splitting integrators for stochastic Lie–Poisson systems

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

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

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

1 Introduction

1.1 Vlasov–Poisson (1d1v)

1.2 Vlasov–Ampère (1d1v)

2 Magnetized Vlasov–Poisson (2d2v)

2.1 Exponential splitting

2.2 Extension to Ampère

2.3 Kelvin Helmholtz Instability

3 Vlasov–Maxwell (1d2v)

$\frac{1}{2} \| \| E_{1} \| \|^{2}$	rel. energy error
$\| \| B_{3} \| \| = 8$

$\| \| B_{3} \| \| = 16$

$\| \| B_{3} \| \| = 32$

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

Authors

Related Research

Poisson Integrators based on splitting method for Poisson systems

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

Splitting integrators for stochastic Lie–Poisson systems

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

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

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

This week in AI

1 Introduction

1.1 Vlasov–Poisson (1d1v)

1.2 Vlasov–Ampère (1d1v)

2 Magnetized Vlasov–Poisson (2d2v)

2.1 Exponential splitting

2.2 Extension to Ampère

2.3 Kelvin Helmholtz Instability

3 Vlasov–Maxwell (1d2v)