A δf PIC method with forward–backward Lagrangian reconstructions

In order to reduce the statistical noise in numerical simulations of kinetic plasma problems, particle-in-cell (PIC) methods often follow a so-called $δf$ approach,^1–3 which consists of decomposing the transported density into two parts, a bulk density f₀ given by an analytical formula and a variation $δf$ represented with numerical particles. In Ref. 4, this approach was revisited as a variance reduction method in the scope of Monte Carlo algorithms, with f₀ playing the role of a control variate, and since then, several techniques have been devised to improve the reduction of statistical error, in particular for gyrokinetic simulation models^5–8 and collisional models.^9–11 

In many practical problems, the plasma either remains close to an equilibrium state^5,12 or evolves as a small perturbation of some analytically known background,^13–15 which can be used as a bulk density. In some cases, however, the plasma evolves in an unpredictable way and f₀ needs to be updated by a self-consistent algorithm to better follow the total density.

Typical examples are problems where full-f simulations are needed. In the physics of tokamak plasmas, one instance is the modeling of E × B staircases,¹⁶ which are long-lived patterns of quasiregular step-like profiles that develop spontaneously in turbulent plasmas. As they slowly move in the radial direction over large time scales, the separation assumption between an analytical background and fluctuations is no longer valid after some time and simulations need to involve full-f methods such as the semi-Lagrangian scheme used in the GYSELA code.¹⁷ Another example is the modeling of the tokamak edge region where the plasma density may strongly deviate from local Maxwellian distributions with large and intermittent fluctuations. This has motivated the development of various Eulerian full-f schemes such as those presented in Refs. 18 and 19, where large plasma blob structures can be seen propagating from the core region toward the tokamak edge.

If one desires to model such problems with a $δf$ PIC method, it is, thus, necessary to allow for general updates of the bulk density over time. An interesting approach in this direction was proposed in Ref. 20: It consisted of projecting the particle density $δf$ on a coarse spline basis and adding the resulting smoothed distribution to the bulk density.

In this article, we consider a variant of this approach where the bulk density is updated using a semi-Lagrangian approach based on the forward–backward Lagrangian (FBL) method.²¹ 

Inspired by Ref. 22, the core of the FBL method is to compute backward trajectories on arbitrary nodes by local inversions of the particle trajectories: As these describe the forward transport flow in phase space and are naturally provided by the PIC code, their local inversion allows to perform semi-Lagrangian updates of a smooth density represented on a coarse grid. In particular, the novelty of our approach is that it does not primarily rely on an accurate particle approximation of the density itself but rather rely on an accurate description of the particle trajectories. As these are, in general, much less noisy than the phase space density, we believe that this new paradigm can lead to efficient low-noise particle methods.

The outline is as follows: In Sec. II, we present our general ansatz for the discrete density, which may be seen as a hybrid discretization between particle and semi-Lagrangian density representations. In Sec. III, we recall the key steps of electrostatic particle-in-cell approximations, and in Sec. IV, we describe the $δf$ PIC method with FBL remappings of the bulk density. The proposed method is summarized in Sec. V, and in Sec. VI, we present a series of numerical results involving two-stream instability (TSI) test cases to illustrate the enhanced denoising properties of our approach. We conclude in Sec. VII with a summary of the proposed method, a discussion on its novelty, and two perspectives for future research.

Our ansatz for the general density at a discrete time tⁿ is a sum of two terms

(we shall often use $:=$ to highlight a definition), where the first one will be seen as the bulk density, a smooth approximation to the full solution, and the second one as the fine scale variations. Following the general principle of $δf$ methods, we require $f*n$ to have a simple expression that is easy to evaluate at arbitrary positions in a general d-dimensional phase space, and we represent the variation $δfn$ as an unstructured collection of numerical particles with coordinates $zkn$ and weights $δwkn$⁠,

Here, $φε$ is a smooth shape function of scale ε, and $z=(z1,…,zd)$ is a phase-space coordinate. In typical problems where the transported density slightly deviates from an initial profile, the bulk density is often set to this initial value, $f*n=f0$⁠, or to some analytical equilibrium.^5,12 In this note, we investigate an alternate approach where the bulk density is represented as an arbitrary collection of B-splines on a coarse grid with the mesh size $h*>ε$⁠, in the spirit of Ref. 20. This yields an expression that is formally similar to that of $δfn$⁠,

where $φ*(·−jh*)$ is now the coarse B-spline shape centered on a general d-dimensional grid node $jh*$ (⁠$Z$ is the set of integers), and $w*,jn$ is its weight. In practice, only a finite number $N*$ of such nodes is used, and as outlined in the introduction, we will periodically update the coefficients of this spline bulk density using forward–backward Lagrangian (FBL) reconstructions, which involves a relatively small set of passive markers designed to track the characteristic flow in phase space.

The main numerical parameters are as follows:

• N_r is the remapping period, i.e., the number of time steps between two updates of the bulk density. The limit value $Nr=∞$ corresponds to a frozen bulk density, namely, $f*n=f*0$ for every time step n.

• $N*$ is the number of coarse splines used to represent the bulk density in the computational domain $Ω⊂Rd$⁠, and it is on the order of $Vol(Ω)(h*)−d$⁠. The limit value of $N*=0$ (an empty grid) corresponds to a “full-f” particle approximation.

• N_p is the number of numerical particles describing the fine scale structures. The limit value of N_p = 0 corresponds to a semi-Lagrangian ansatz²³ where the full density is represented on a structured grid and can be updated in time with an ad hoc scheme such as the FBL method described in Ref. 21.

For simplicity, we consider spline shape functions for both the bulk density and the fine scale variations. Specifically, we set

with a reference shape function $φ$ defined as a centered cardinal B-spline of degree p,

involving standard univariate B-splines defined recursively by

In the sequel, it will be convenient to denote an arbitrary collection of weighted splines as

where $W=(wk)k=1…N, Z=(zk)k=1…N$⁠. With this convention, the two components of our general ansatz (1) read

where $Z*:=(jh*)j∈Zd$ are the nodes of the spline grid.

Our method may be described for general non-linear transport problems of the form

where $z∈Rd$ is the phase-space variable, and $U[f]$ is the generalized velocity field associated with the solution f.

The characteristic trajectories associated with Eq. (5) are the curves $Z(t)=Z(t;s,z)∈Rd$⁠, solution to the ordinary differential equations (ODEs)

for arbitrary $s,t∈[0,T]$ and $z∈Rd$⁠, see, e.g., Ref. 24. The (forward) characteristic flow between two times $tn=nΔt$ and $tn+1=(n+1)Δt$ is then defined as

and the inverse mapping $Bexn,n+1:=(Fexn,n+1)−1$ is the backward flow

Using the backward flow, we can write the solution to Eq. (5) over the time interval $[tn,tn+1]$ as

Here, we may restrict ourselves to divergence-free velocity fields: $divzU[f]=0$⁠. The characteristic flows are then measure-preserving, and the transport is conservative.

A simple example is the periodic Vlasov–Poisson equation in a two-dimensional phase space, i.e., $z=(x,v)$ with a periodic space coordinate $x∈[0,L]$ and velocity $v∈R$⁠,

for $t≥0, (x,v)∈[0,L]×R$ with a normalized initial density $∫0L∫Rf0(x,v) dv dx=1$ and a periodic electric field $E=E[f]$ defined by

This corresponds to an electrostatic, normalized (⁠$ε0=qe=m=1$⁠) periodic electron plasma in 1D with constant neutralizing background ion density, so that $∫0Lρ dx=0$⁠. Here, the generalized velocity field is

Due to the non-linear nature of this transport equation, the characteristic flow has no explicit expression.

Particle approximations represent the transported density $fn(z)≈f(tn,z)$ as a sum of numerical particles of the form

with weights initially set to

where g⁰ is the sampling distribution of the initial markers $Z0=(zk0)k=1,…,Np$⁠, see, e.g., Ref. 11. As the problem is conservative, the weights are kept constant in time, $Wn+1=Wn$⁠, and the markers are pushed forward

using some approximation to the forward flow in Eq. (6), which takes as parameter the numerical solution at time tⁿ,

For the Vlasov–Poisson equation [Eq. (9)], a standard numerical flow is given by a leap-frog (Strang splitting) scheme

where $°$ is the usual function composition operator associated with a grid with M points and the step size $Δx=L/M$⁠. The first split flow reads

and the second one takes the form

Here, we have set

and the electric field is computed from the particles with a discrete Poisson solver such as the one described just below, which takes as input the particle weights and positions,

To clarify how the particles are coupled with the field in Eq. (18), we recall the precise form of a basic Finite Difference (FD) Poisson solver for a general collection of weighted particles, i.e.,

We note that for periodic systems spectral solvers involving fast Fourier transforms are often used, see, e.g., Ref. 25, Sec. 8. In both cases, the charge deposition algorithm consists of computing grid values for the charge density

for $i=1,…M$⁠. A discrete Poisson equation is then solved. With a FD scheme, this reads

for $i=1,…,M$⁠. The electric field (19) is then defined using the same spline shape functions as the particles but a different normalization and scaled with the FD grid as follows:

In the usual case where $ε=Δx$⁠, we write $EΔxpic=EΔx,Δxpic$⁠.

We now describe the main steps of a $δf$ method involving a decomposition of the form [Eq. (1)]. Denoting by N_r the remapping period, we first observe that the bulk part is frozen between two remapping steps. This yields, in particular,

The $δf$ part is then evolved as follows.

A typical $δf$ time step consists of (i) pushing the markers forward to follow the characteristic curves as in Eq. (13)

and (ii) updating the weights. Indeed, since the numerical particles now represent a variation $δfn=fn−f*n$ with frozen $f*n$ on the time step, their weights no longer carry a conserved density and they must be evolved in time. In the direct$δf$ method (see, e.g., Ref. 11, Sec. 3), they are set as

where

correspond to the transported density and the markers density probability, respectively, evaluated at the particle positions. Since both values are constant along exact trajectories, a reasonable choice is

Up to these constant values, we see that the $δf$ weights only depend on the bulk density, which is characterized by the weights $W*n=W*m$⁠, at the marker positions $Zn$⁠. Hence, we may summarize the weighting scheme in Eqs. (23)–(25) as

In the case of a $δf$ ansatz, the leap-frog flow reads

with split flows defined similarly as in Eqs. (15)–(16) and a field solver adapted to the ansatz from Eq. (1), namely,

with

Since both parts of the density are formally written as weighted collections of spline shapes, we may use a standard PIC solver for each parts. The field induced by the bulk density $f*m$ is

and the field corresponding to the variation $δfn+12=Φε[δWn+12,Zn+12]$⁠, with scaling $ε=Δx$⁠, reads

see (17), (26), and (23). The resulting field solver corresponding to (29) may then be written as

Since the field $E*m$ is constant between two remappings, regular time steps only involve a deposition of $δf$ particles to update the $Eδf$ part. Note that even if $|δfn|≪|f*m|$⁠, there is no reason why $Eδf$ should be small compared to $E*m$⁠. In particular, we have $E*m≈0$ when $f*m$ is close to an equilibrium density.

For the remappings, namely, the updates of the bulk density $f*m↦f*n$ with $n=m+Nr$⁠, we resort to forward–backward Lagrangian (FBL) reconstructions.²¹ As mentioned above, this involves a collection of auxiliary markers that have been reset on the Cartesian grid at the previous remapping step,

and pushed forward similarly as the classical ones in Eq. (22),

in order to track the forward characteristic flow. The FBL method then performs a semi-Lagrangian step

where $A*$ is a spline interpolation or quasi-interpolation operator on the $h*$ grid, and $Tfbl[Z̃m]:f↦f°Bfblm,n$ is a transport operator approximating the exact one in Eq. (8), based on the FBL backward flow

As a result, we obtain a numerical scheme to update the coefficients of the bulk density, which we summarize as

where n is a multiple of N_r and $m=n−Nr$⁠.

In Sec. IV D, we give some details about the FBL algorithm.

The FBL scheme in Eq. (34) performs a semi-Lagrangian approximation with a backward flow derived from the position of the auxiliary markers $Z̃m$ in Eq. (33), as follows:

• For any index $j∈Zd$⁠, let $z̃im$ with $i=i(j)$ be an FBL marker close to the node $jh*$⁠;

• Compute the quadratic local backward flow associated with $z̃in$⁠, namely,

  $Bim,n(z)=zi*+Din(z−z̃in)+12(z−z̃in)TQin(z−z̃in),$

  (37)

  where $zi*=ih*$ and $Din, Qin$ correspond, respectively, to the Jacobian matrix and Hessian tensor of the backward flow at $z=z̃in$⁠, computed from the positions of the logical neighbors of $z̃in$ (for more details see Ref. 21);

• Reconstruct the transported density at node $zj*=jh*$ using this approximated flow

  $f*,jn:=f*m(Bim,n(zj*)).$

  (38)

• Compute the new weights $W*n$ using a spline approximation operator $A*$⁠.

For the operator $A*$⁠, we can use a standard interpolation scheme or even a local quasi-interpolation scheme, see, e.g., Ref. 26. Both compute high-order B-spline approximations, but the latter only involves local pointwise evaluations of the target function. The resulting approximation takes the form

with

with coefficients $ai:=ai1⋯aid$ on a local stencil of size s_p, as given by the symmetric method (⁠$ai=a−i$⁠) in Ref. 26, Sec. 6.

Gathering the above steps, we may write the solution at time tⁿ as

where we remind that $Z*=(jh*)j∈Zd$ are the coarse grid nodes and $m=m(n)$ is the last remapping step preceding n, see Eq. (21). The proposed scheme can then be summarized by the flow chart in Fig. 1 and reads as follows.

• Compute the weights $W*0$ of the initial bulk density $f*0=Φh*[W*0,Z*]$ with

  $f*0:=A*f0;$

• Use a spline approximation scheme such as the one in Eq. (39);

• Draw the $δf$ markers according to some initial sampling distribution

  $Z0∼g0;$
• Initialize the auxiliary FBL markers on the $h*$ grid (32),

  $Z̃0:=Z*.$

For $n=0,…,Nt−1$⁠, let $m:=⌊nNr⌋Nr$ denote the previous remapping step as in Eq. (21), and do

• If $n=m>0$⁠, remap

  • (i)
    Update the bulk density weights (coefficients) with the FBL method described in Sec. IV D

    $W*m:=Remap[Z̃m](W*m̂).$

    (41)
  • (ii)

    Reset the FBL markers with (32), $Z̃m:=Z*$⁠.

• Push the markers forward

  ${Zn+1:=Fn(Zn),Z̃n+1:=Fn(Z̃n)$

  (42)

  using some discrete flow $Fn=FΔtδf[W*m,Zn]$ such as the leap-frog scheme in Eq. (27), which involves deposition and field solver steps as detailed in Secs. IV B and III E.

A semi-Lagrangian FBL variant of the above scheme can be obtained by (i) setting N_p = 0 and (ii) treating the FBL markers $Z̃n$ as “active” particles in the field solver.

Specifically, this amounts to discarding the $δf$ part (⁠$Zn=∅$⁠) so that the full density is represented on the structured spline grid on the successive remapping steps and to defining the electric field in Eqs. (27) and (28) by

as one would do with a standard full-f PIC scheme, see (16)–(18), using the weighted FBL markers as standard particles. Note that here the weights are those set at the last remapping step (41), i.e., $W*n+12=W*m$ with m as in (21).

In this section, we assess the proposed method on the periodic 1D1V Vlasov–Poisson system (9) and (10), where we compare it with other PIC and semi-Lagrangian schemes.

For the purpose of comparison, we use different values for the main parameters listed in Sec. II A, including some limiting values, which allow to emulate full-f and $δf$ PIC schemes, as well as a semi-Lagrangian FBL scheme. All the methods use cubic splines for particle shape functions.

As the solutions are typically supported in a region occupying roughly one half of the computational phase-space domain, for each grid resolution $Nx×Nv$⁠, we indicate as $N*eff≈N*/2$ with $N*=NxNy$⁠, the approximate number of active spline coefficients on the grid, corresponding to the “effective” grid resolution. The resulting methods are summarized in Table I.

Here, the active/passive/frozen status of FBL markers $Z̃n$ refers to their role in the scheme and, in particular, in the field solver that is involved in the push step in Eq. (42): Active markers (used in the FBL variant) carry a charge which is deposited as explained in Sec. V C, while passive markers (used in the proposed FBL-$δf$ scheme) do not, as described in Sec. IV B. Finally, frozen markers are passive and never pushed: In the absence of remappings, only their initial position on the grid is used to represent the bulk density $f*n=f*0$⁠.

We point out that in this study we do not attempt to assess the computational cost of various methods, as the testcases are very simple and our implementation is far from being optimal. In particular, one should keep in mind that full-f and basic $δf$ methods, even if noisier, may be quicker to run and easier to parallelize than their FBL counterparts. Efficient implementations of the latter are left for future works.

Our numerical experiments are based on a textbook 1D1V two-stream instability test case where the initial distribution is of the form

with a wave number of $k=12$ and a perturbation amplitude of $α=0.01$⁠. This solution is known to develop thin filaments in the phase space that are difficult to resolve numerically, see, e.g., Refs. 27 and 28. Here, we will simulate two-stream instabilities on the time interval $[0,100]$⁠. We shall consider two versions of the problem:

• the basic Two-Stream Instability (TSI) test case where the solution is evolved with the self-consistent Vlasov–Poisson system (9) and (10);

• a Driven Two-Stream Instability (DTSI) test case where an oscillating external field

  $Eext(t,x)=110cos (πt50)$

  (45)
  is added to the self-consistent field from Eq. (10). The effect of this external drive is to add a periodic perturbation to TSI trajectories, uniform in x and v,

  $(x,v)↦(x+250π2(1−cos (πt50)),y+5πsin (πt50)),$

  so that the solutions to both test cases coincide for t = 100.

The computational phase-space domain is set to

with $L=2πk$ and $vmax=10$⁠. Note that this value is significantly larger than what is usually seen in the literature, in order to correctly represent the variations in velocity induced by the oscillating drive. Finally, the initial sampling distribution for the particles is taken as $g0=(∫f0(z) dz)−1f0$ for the full-f PIC scheme and as

with σ = 2 for the $δf$ and FBL-$δf$ PIC methods.

In this section, we assess the accuracy and denoising capabilities of the basic $δf$ PIC and FBL-$δf$ PIC methods by plotting various densities in phase space. In Figs. 2–9, the axes show $x∈[0,L]$ and $v∈[−vmax,vmax]$ with L and $vmax$ as indicated above, and the same color range is used for the purpose of comparison.

In Figs. 2 and 3, we first plot the densities fⁿ corresponding to different methods listed in Table I for both the TSI and DTSI test cases, at the time $tn=100$ where the exact solutions coincide.

For both test cases, we observe very noisy (oscillating) profiles for the coarse full-f PIC solutions, which use about 20 000 particles, and smoother profiles for the fine PIC solutions, which use 16 times as many particles. FBL solutions are essentially free of spurious oscillations with many fine structures well resolved in the fine solutions compared to the coarse ones. We also note that the TSI and DTSI solutions are in very good agreement at $tn=100$ for these full-f numerical simulations.

Regarding the $δf$ solutions, we observe that in the TSI test case (Fig. 2) both the basic $δf$ PIC and the FBL-$δf$ PIC methods allow to strongly reduce the level of oscillations compared to the coarse PIC simulation—all three using approximatively 20 000 active particles. In the DTSI case, however (Fig. 3), we see a significant reduction in the accuracy of the basic $δf$ solution, where the global central vortex suffers from a visible deterioration and spurious streams of particles have appeared at high velocities. In contrast, the accuracy of the FBL-$δf$ PIC method in the DTSI case remains very comparable to that of the TSI case. In particular, we observe that the new method is able to capture several low-density filaments, i.e., fine regions with zero or very low density distinct from the central vortex. In the basic $δf$ simulation, these fine structures are poorly resolved, as we can see from Fig. 3 where the low-density filaments on the left of the central vortex are virtually merged with the latter. To better visualize the difference in the two methods, we next show in Figs. 4–9 a series of plots corresponding to snapshots of various densities at successive times $tn=0,12.5,…,87.5$ for the DTSI test case: In Figs. 4 and 5, we show the evolution of the $δf$ and FBL-$δf$ solutions, completed by that of a reference scheme (an FBL simulation using a 1024 × 1024 grid) in Fig. 6. In Fig. 7, we then show the evolution of the bulk density $f*n$ computed by the FBL-$δf$ PIC method, which we remind is periodically remapped on a 50 × 50 grid of cubic splines, see Table I. Finally, the evolution of the residual $δfn$ densities is shown in Fig. 8 for the basic $δf$ method and in Fig. 9 for the FBL-$δf$ method.

From these plots, we see that the density strongly deviates from its initial value due to the external drive, which clearly challenges the basic assumption of a standard $δf$ approach. Here, we have chosen a deviation so large that the bulk density cannot be compensated by a proper $δf$ weighting of the particles, leading to significant errors in the full fⁿ density visible in Fig. 4. As expected, the solution computed with the FBL-$δf$ PIC method in Fig. 5 does not show this erratic behavior. The main reason for this is visible in Fig. 7 where we see that the coarse bulk density computed by the FBL-$δf$ scheme properly follows the total density. As a result, the $δf$ particles can effectively represent small scale variations with respect to the main part of the plasma density. This is made evident in Figs. 8 and 9 where we see that the amplitude of the $δf$ particle density for the FBL-$δf$ scheme is significantly lower than that with a basic $δf$ method.

In Fig. 7, we also observe that the numerical bulk density is subject to a visible numerical diffusion. This is an expected consequence of the periodic coarse-grid interpolations involved in the semi-Lagrangian updates, and over moderate times, it should not pose a serious issue since the total density is described by the sum of the bulk and the perturbed particle density. Over long simulation times, however, this numerical diffusion will lead to bulk densities that do not approximate well the total f. In such cases, a proper re-initialization of the bulk density $f*$ will have to be devised, based on the available representation of the total density, and various methods do exist in this direction such as adaptive kernel density estimation algorithms^29–31 or adaptive sparse-grid methods.³² We note that such an algorithm is likely be more costly than a simple update of the coarse bulk density, but this should be mitigated by the fact that it will be applied on much larger time periods.

In Figs. 10 and 11, we next plot the snapshots of the deposited charge densities ρⁿ for the DTSI test case, corresponding to $δf$ and FBL-$δf$ solutions shown in Figs. 4 and 5. For a better assessment of the accuracy, we also plot in dashed line a reference charge density computed by a fine semi-Lagrangian scheme. The associated Eⁿ fields are plotted in Figs. 12 and 13.

These plots show that despite the velocity averaging, the high level of noise visible in the $δf$ phase-space densities is also present in the deposited charge densities, which results in significant errors in the self-consistent electric field.

We next assess the accuracy of the different methods in approximating the system momentum and energy, namely,

with $z=(x,v)∈Ω=[0,L]×[−vmax,vmax]$ as above. In Fig. 14, we plot the time evolution of these discrete quantities for the TSI test-case, where they are invariants of the exact solution (assuming f = 0 for $|v|>vmax$⁠). The first finding is that both the full-f PIC and the FBL scheme (for the coarse and fine meshes) preserve these quantities with very good accuracy, while visible errors accumulate for the two $δf$ methods. More precisely, we observe oscillating errors for both $δf$ schemes, with no clear winner: Although the oscillations are clearly of larger amplitude with the basic $δf$ scheme, a secular growth is visible on the numerical energy for the FBL-$δf$ method, which may be due to the numerical diffusion of the bulk density (visible in Fig. 7) caused by the periodic coarse updates as discussed above. We note, however, that for the TSI test-case both $δf$ schemes did also show similar accuracy in the density plots of Fig. 2, so that their similar performance in the conservation properties may not come as a big surprise.

For the DTSI test-case, the total momentum and energy are no longer system invariants due to the external drive. To compare the performance of the different methods, we, thus, plot in Fig. 15 the energy of the self-consistent electric field, for both the TSI (left) and DTSI (right) test-cases. Here, the first observation is that the results present a high level of variability for the different methods. However, not all methods perform equally. If we consider that in both test-cases, the fine FBL method (purple curve) can be used as reference, then we find that the curves of the FBL-$δf$ scheme (in red) are of reasonable accuracy. More precisely, they are less damped than those of the coarse FBL scheme (in green), and they show less oscillations than those of the coarse full-f PIC scheme (in blue). In the TSI test case, the curve of the basic $δf$ scheme (in orange) shows a similar behavior with slightly more oscillations, while in the DTSI test case, the oscillations are of very high amplitude (close to the actual level of the energy itself). As a result, we, thus, find that for the DTSI test case the new scheme also shows a superior accuracy in terms of energy approximation.

We conclude our numerical experiments by considering the empirical variances for the total particle number and current, as defined by Eqs. (4.80) and (4.81) in Ref. 8, namely,

and

respectively. We note that here the weight normalization differs from that of Ref. 8: For $δf$ and full-f simulations (where the weights are constant), we have

respectively, see (23) and (12). In particular, $σN=0$ for full-f runs using a sampling distribution $g0=f0/(∫f0(z) dz)$⁠.

Following the interpretation of the $δf$ method as a variance reduction technique,⁴ these respective quantities may be used as quantitative indicators of the statistical errors associated with the Monte Carlo approximation of the total density $∬f(t,x,v) dx dv$ and normalized current $∬f(t,x,v)v dx dv$⁠, see Ref. 8, Sec. 2.

In Fig. 16, we plot these quantities for the two test cases, and for the three coarse PIC schemes, namely, the coarse full-f PIC, the basic $δf$⁠, and the FBL-$δf$ scheme. For the TSI test case corresponding to the top plots, we observe that the $δf$ schemes indeed result in a reduction in the current variance σ_J, of a factor between two and three (the full-f density variance is zero as previously noticed). We also note a secular growth in the case of the FBL-$δf$ scheme: This is probably due to the slow diffusion of the bulk density caused by the periodic remappings and should be mitigated by proper re-initialization methods for the bulk density. We refer to the perspective section below for a discussion on possible improvements.

For the DTSI test case, we observe similar values for the current variances σ_J of the full-f FBL-$δf$ PIC simulations (with values around 20 and 10, respectively). In contrast, strong oscillations are visible in the variances of the basic $δf$ PIC scheme, which are closely correlated with the deviation in velocity already seen in Fig. 4, and the associated offset with the initial density. The amplitude of these oscillations (around 40) confirms our previous observations that a basic $δf$ approach is no longer justified when the density deviates significantly from its initial profile and further validates the proposed method in such regimes.

In this article, we have presented a $δf$ PIC method where the bulk density is represented by smooth B-splines on a coarse grid and periodically updated with a forward–backward Lagrangian (FBL) method. This method amounts in updating the spline coefficients of the bulk density with a backward semi-Lagrangian method, where the backward trajectories are obtained by local inversions of the forward trajectories provided by the PIC code. In these inversions, the local backward flow is represented as a quadratic mapping whose Jacobian and Hessian matrices are computed from the relative displacements of regular passive markers pushed forward by the PIC code. After describing the representation of the bulk and variation density in terms of weighted collections of shape functions with variable width and recalling the general form of full-f and $δf$ particle approximations to transport problems, we have provided a detailed description of a remapped $δf$ scheme where the bulk density updates are performed using the FBL method. Finally, we have studied the basic performance of the method on two textbook examples involving one-dimensional (i.e., 1D1V phase-space) electrostatic problems: a simple two-stream instability where the deviation from the initial density is mild and does not justify any update of the bulk density, and a driven one where an oscillating electric field induces strong displacements of the phase-space density in velocity. In the latter case, the basic equilibrium/deviation separation assumption is no longer valid and the $δf$ method fails as expected to correctly resolve the solution. In contrast, the FBL-$δf$ method shows good results and a measurable reduction of the statistical error.

Although we did not compare the numerical performance of our method with the direct projection scheme of Ref. 20, we believe that opting for semi-Lagrangian updates of the bulk density is a significant novelty which has the key feature that it does not rely on an accurate particle approximation of the total density, but rather on accurate particle trajectories. Two typical configurations may illustrate why this difference is meaningful: The first one is a regime where a localized plasma blob moves through an otherwise empty phase space. In a region crossed by the blob, the bulk density should ideally go from zero to some positive density and then go back to zero as the blob leaves the region. With a direct projection scheme, however, this cannot happen: As the $δf$ particles enter the region, the bulk density will be updated by their coarse projection and take a positive value, but as they leave, they will not be able to exactly compensate the smooth bulk density (and neither will their projection), because they have precisely moved away from it. As a result, the updated bulk density will retain a spurious residual of that fine discrepancy, which will remain there since the $δf$ particles have now left. Another problematic scenario is the opposite one, where both the density and the transport flow are smooth in phase space. If the bulk density grid is not too coarse, then one will always find some phase-space cells devoid of markers: On such cells, the projected $δf$ density will be zero by construction, leading to strong oscillations. Because the forward–backward Lagrangian approach relies on the particle trajectories and a smooth grid-based bulk density, it is exempt of both these pitfalls: For the case of a moving blob, this has already been seen on the DTSI test case where the external drive moves the plasma density over large portions of the phase space, with Fig. 7 showing a clean transport of the bulk density (e.g., between times $tn=0$ and 25). For the second one, the key point is that a smooth flow leads to accurate backward trajectories, according to the analysis of Ref. 21. On every grid node, the new value will, thus, correspond to that of the previous bulk density at a point where the smooth spline density should be also accurate, independent of the markers resolution.

We also point out that this approach not only allows to approximate general bulk densities, but it can be applied to other representations, as long as their evaluation is easily computable: Instead of B-splines, one could, for instance, use piecewise polynomials or a collection of Fourier modes. The semi-Lagrangian updates would then involve specific algorithms (such as interpolation or FFT) for updating the weights, based on nodal evaluations of the previous density at the feet of the backward characteristics computed by the FBL scheme.

Moving forward, an important direction of research is to extend the present method to more challenging problems such as higher-dimensional ones. We are, in particular, interested in gyrokinetic models where the transport phase space has typically four dimensions. In their current form, the various algorithms (namely, the FBL flow reconstruction and the spline quasi-interpolation) can be directly applied to a 4D phase space; however, preliminary numerical simulations have shown that typical transport flows in higher dimensions have an anisotropic smoothness which is not well captured by our isotropic algorithms. Our plan is, thus, to design extensions of the current approach that will be able to efficiently reconstruct high-dimensional characteristic flows with anisotropic smoothness.

Another improvement will be to design good re-initialization schemes for the bulk density that could be applied on larger time scales to counter the numerical diffusion inherent to the coarse semi-Lagrangian updates, as well as the secular growths in the energy error and in the variance, observed in Figs. 14 and 16. Adaptive approximation methods do exist which have already shown good results for particle denoising, such as adaptive wavelet or Fourier filters^29,30 kernel density estimation algorithms,³¹ which automatically select the optimal width for the particle shape functions, or sparse-grid adaptive methods,³² which have also proven effective for particle denoising in higher dimensions.

The authors thank Roman Hatzky for fruitful discussions on statistical error reductio, as well as the anonymous reviewers for their comments, which have improved the presentation of our results. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement (No. 633053). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

The authors have no conflicts to disclose.

Martin Campos Pinto: Conceptualization (lead); Methodology (equal); Software (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal). Merlin Pelz: Conceptualization (supporting); Methodology (equal); Software (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Pierre-Henri Tournier: Conceptualization (supporting); Methodology (equal); Software (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

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