<figure><img src="image_1.png" /></figure>

HAL
open science

# A binary filter inspired from the PIC sparse grid technique – Illustration on the XTOR-K code

T. Nicolas, V. Dubois, Q. Fang, Hinrich Lütjens

► To cite this version:

T. Nicolas, V. Dubois, Q. Fang, Hinrich Lütjens. A binary filter inspired from the PIC sparse grid technique – Illustration on the XTOR-K code. Computer Physics Communications, 2024, 299, pp.109151. 10.1016/j.cpc.2024.109151. hal-04512472

HAL Id: hal-04512472

https://hal.science/hal-04512472v1

Submitted on 14 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Computational Physics

# A binary filter inspired from the PIC sparse grid technique – Illustration on the XTOR-K code ⭐

T. Nicolas *, V. Dubois, Q. Fang, H. Lütjens

CPHT, CNRS, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

<figure><img src="image_3.png" /></figure>

Check for
updates

## ARTICLE INFO

**Keywords:**

Sparse grid

Particle in cell

Hybrid kinetic/fluid code

## ABSTRACT

It is known [1,2] that the sparse grid method for Particle-In-Cell (PIC) solvers acts as a filter to reduce the PIC noise. In this paper, a simple rule to discard or keep modes in Fourier space (a binary filter with values either 0 or 1) is derived using the sparse grid combination formula. Its relation to the standard sparse grid filter, which is characterized quantitatively, is explained. The relations between the sparse grid filters on grids of arbitrary levels are also outlined. Namely, in two (resp. three) dimensions and for bi-linear (resp. tri-linear) moment deposition, it is proven rigorously that the sparse grid filter, for a grid of size equal to an arbitrary power of two, can be expressed in terms of two (resp. three) unique real valued functions. The advantage of the binary filter over the standard sparse grid filter is the reduction of signal deformation introduced by the latter, for the same noise reduction capability. By applying the filter to moments of a marker distribution coming from the XTOR-K code, it appears the noise could be significantly reduced, with a moderate overhead in the moment deposition part of the algorithm.

## 1. Introduction

The Particle-in-cell (PIC) method is one of the basic tools of computational plasma physics. It is used in purely kinetic codes as well as hybrid kinetic/fluid codes. In a PIC code, markers representing chunks of the distribution function are evolved in continuous space-time according to an electro-magnetic field known on a grid. The moments of the marker distribution are periodically deposited on the grid in order to evolve said fields according to some subset of Maxwell's equations. The time evolution of charged particle in a magnetized plasma usually follows the Boris-Buneman algorithm [3,4]. PIC methods are robust, versatile and easily parallelized. In spite of this, the well known problem common to all PIC codes is the statistical noise associated with the limited number of markers per grid cell. The noise being inversely proportional to the square root of the latter, reducing the PIC noise by brute force increasing the total number of markers often leads to a prohibitive increase of computational complexity.

In the past, different strategies have been deployed in PIC codes in order to address this noise issue, such as the $\delta f$ technique [5,6], quiet start [7–9], and noise filtering, which can be spatial [2,6,10,11] or temporal [9,12,13]. In ref. [1], noise reduction is achieved by applying the

sparse grid combination technique [14,15] to the PIC method. The idea is to collect moments of the marker distribution (e.g. the density) on several grids that are refined in only one direction at a time, and sparse in the other directions. The result is that the number of markers per cell is significantly increased. The different results on the different grids can finally be interpolated on the desired fine grid and combined with the so-called sparse grid combination technique [14]. The resulting numerical solution is significantly denoised with a moderate increase in numerical complexity. Although the sparse grid combination technique had already been used in plasma physics [16–18], it was first suggested in reference [1] as a way to reduce the PIC noise. Following this seminal paper, the merits of the sparse grid PIC technique have been exploited in plasma applications pertaining to radio frequency discharges [19] and to the electron drift instability in Hall thrusters [20]. The mathematical properties of the sparse grid PIC technique and its parallel implementation have been extensively studied in refs. [21,22]. In these references, the numerical error in the sparse grid PIC method is rigorously estimated in terms of grid based error on the one hand, and particle based error on the other hand.

The present article can be seen as an alternative way to explain the action of the sparse grid combination technique in the context of the PIC

<footnote>* The review of this paper was arranged by Prof. David W. Walker.</footnote>

<footnote>* Corresponding author.
E-mail address: timothee.nicolas@polytechnique.edu (T. Nicolas).</footnote>
method. Namely, we regard the sparse grid technique as a filter, and focus on the action of the filter in Fourier space. Before going any further, let us note that the sparse grid PIC technique can be applied in two distinct ways, as explained in ref. [1]. In the first, the electromagnetic (EM) field is evaluated on the subgrids using the velocity moments collected on the subgrids. The sparse grid combination technique is then used on the EM field to reconstruct the fine grid field. In the second, the sparse grid combination technique is applied to the moment only to reconstruct the denoised moments on the fine grid, which is used to compute the fine grid EM field. The present paper is limited to the latter case. In this case, the sparse grid combination technique acts as a sort of filter, as has been fully realized in ref. [2]. Nonetheless, ref. [2] does not precisely characterize the filter, although estimates for the error caused by different combination techniques are precisely computed. By characterizing the filter, we mean specifying how the filter acts in Fourier space, that is, by what coefficient each mode is multiplied. It turns out to be very simple. The essential idea is that the moment deposition step convolves the marker distribution with the deposition shape function. In Fourier space, the different modes of the signal are then multiplied by the Fourier transform of the shape function. The shape of the filter results from the combination of these factors using the sparse grid combination technique. In this paper, we prove that in two (resp. three) dimensions, the shape of the filter can be inferred from two (resp. three) single real-valued functions, that are independent of the size of the grid. Even the truncated combination formulae of ref. [2] are strongly connected with these functions. More importantly, we derive a binary filter, consistent with the combination formula, which has almost the same action as the sparse grid filter, except that its values are either 0 (the mode is discarded) or 1 (the mode is kept). The merits of this filter over the standard one are emphasized.

An advantage of the point of view in terms of filters is that the very same denoising effect as the sparse grid combination technique can now be obtained by applying the filter to the fine grid moment. This means that the process of collecting the moment on $N_g$ different grids, with a numerical complexity scaling as $N_g N_p$, where $N_p$ is the number of markers, can be replaced with only one deposition step, followed by a forward and backward fast Fourier transform. The issue is that $N_g$ can be large, especially in three dimensions. Therefore, it can be expected that replacing the sparse grid combination technique with the sparse grid filter, which has exactly the same action, could be faster. Note that efficient parallel implementations of the fast Fourier transform, such as advertized in ref. [23], may be required to fully leverage the advantage of the point of view developed in the present article.

The rest of the manuscript is organized as follows. In section 2, the reduction factor of Fourier modes due to the convolution with the shape function is introduced after preliminary definitions. The tent function, corresponding to linear interpolation, is assumed. The filter is then defined using the sparse grid combination technique. In section 3, the properties of the filter are detailed. Most mathematical proofs of these properties are confined in the appendices. In section 4, the action of the sparse grid combination technique is reinterpreted in the case where moment deposition is done directly in Fourier space. We obtain a binary filter that has qualitatively the same properties as the sparse grid filter, but is much simpler to compute and more intuitive to understand. In section 5, the different filters are applied to distributions coming from the hybrid kinetic/fluid magnetohydrodynamic code XTOR-K, in order to analyze how the PIC noise in the radial direction can be reduced. The different filters are compared. In section 6, the case of truncated schemes of ref. [2] is analyzed. In section 7, it is shown that the results of the paper for the standard sparse grid filter easily extend to any order of the shape function, but with some restrictions, namely, the filter is positive definite only when the order of the shape function is odd. Finally, we conclude in section 8.

## 2. Signal reduction and analytical expression of the sparse grid filter in two and three dimensions

### 2.1. Preliminary definitions and notations: grid level, degrees of freedom and moment deposition

We start by giving the framework of our work. We will only use grids with a number of intervals equal to a power of 2 in all directions. The obvious advantage of such grids is that they can easily be divided in coarser grids, which naturally combines well with the sparse grid formalism. A grid that has $2^p$ intervals (hence, $2^p + 1$ points) will be said to be level $p$. In $d$ dimensions, the level is a $d$-tuple $(p_1, ..., p_d)$. If $p_1 = ... = p_d = p$, then we will speak of the grid of level $p$ without further specification. The coarsest grid considered will be level 1. The $d$-dimensional grid is always supposed to describe the $[0, 1]^d$ volume. The coordinates are denoted $X_1$ to $X_d$. To differentiate from this notation, the $j$th element of the grid for the coordinate $X_i$ is denoted $X_i^{(j)}$. The coordinates of the grid are sampled according to

$$X_i^{(j)} = j2^{-p}, \quad j \in [0, 2^p] \qquad (1)$$

that is, the interval size is $h = 2^{-p}$. We add the restriction that all signals considered have periodic boundary conditions. The restriction is not as stringent as it seems: in many applications of interest, one is interested only in what happens in a region interior to the grid, so that the perturbations can be assumed to vanish at the boundaries, making periodic boundary conditions suitable. With periodic boundary conditions, a one-dimensional grid of level $p$ has $2^p$ degrees of freedom (dof). According to the Nyquist theorem, the fastest signal that can be represented with such a sampling has a mode number of $k = k_{max} \equiv 2^{p-1}$, corresponding to a signal $s = \cos(2\pi 2^{p-1}(X_i^{(j)})_{j \in [0, 2^p]}) = ((-1)^j)_{j \in [0, 2^p]}$. The corresponding sine wave identically vanishes and, therefore, does not contribute any dof. The zero frequency mode also contains only one dof. Therefore, all modes between $k = 1$ and $k = k_{max} - 1$ contribute two dof (cosine and sine), while $k = 0$ and $k = k_{max}$ contribute only one dof. The total number of dofs is $N_{dof} = 2(k_{max} - 1) + 1 + 1 = 2^p$. Without surprise, it is found that the Fourier representation has exactly the same number of dofs as the underlying grid, and this statement is true in any dimension.

Finally, we assume the following for the moment deposition method of the PIC algorithm. Unless stated otherwise, it is assumed that the markers' shape function is derived, in any dimension, from the so-called tent function:

$$\tau_1(x) = \begin{cases} 1 - |x| & \text{if } |x| \le 1 \\ 0 & \text{if } |x| > 1 \end{cases} \qquad (2)$$

The tent function is the order one of a series of shape functions of any order, built by successive convolutions of the hat function, which is non zero and equal to one only for $-0.5 < x < 0.5$. In the following, we will denote the order of the deposition function by the letter $q$. The hat function and tent function correspond respectively to $q = 0$ and $q = 1$. The case of shape functions of arbitrary order $q$ is examined in section 7. For now, the tent function of equation (2) is assumed. The shape function is then defined by [1]

$$S(x) = \prod_{i=1}^{d} \frac{\tau_1(x_i/h_i)}{h_i}, \qquad (3)$$

where $h_i = 2^{-p_i}$ is the interval size in direction $i$. Thus, a marker located in phase space at $(x_\ell, v_\ell)$ contributes

$$\frac{N}{N_P} S(x - x_\ell) v_\ell^m \qquad (4)$$

to the velocity moment of order $m$, where $N_P$ is the number of PIC markers and $N$ represents the total number of physical particles modeled by the PIC algorithm. This is the moment deposition of the standard PIC method.
With these notations, we can now investigate the fundamental mechanism behind the sparse grid PIC method.

## 2.2. Signal reduction at moment deposition in one dimension

The whole idea of the present paper comes from the following simple fact. When a sinusoidal signal with mode number $k$ (integer) is represented by a distribution of markers, if one collects the moments of that distribution on a grid of level $p$ by linear deposition (using shape function (3)), then the sinusoidal signal observed on the collected moment is reduced by the factor [10]

$$f(k, p) = \text{sinc}^2\left(\frac{k}{2^p}\right) = \left[\frac{\sin(\pi k 2^{-p})}{\pi k 2^{-p}}\right]^2, \qquad (5)$$

which should come as no surprise to the reader familiar with signal processing theory. Note that we incorporate the $\pi$ factor in the definition of the cardinal sine. The reason for equation (5) is that the grid density is obtained as the convolution between the shape factor and the marker density. Therefore, in Fourier space, the Fourier coefficients of the marker density are multiplied by the Fourier transform of the shape function [10]. As it turns out, the Fourier transform of the triangular shape function is the square of the cardinal sine function (see also section 7 where other orders of the deposition method are considered).

A less mathematical argument explains qualitatively equation (5) as follows. Assume a physical density given by $\rho(x) = 1 + a \cos(2\pi kx)$, where $k \in \mathbb{N}$ and $a < 1$ (this constraint is necessary only to obtain a non-negative density). Now, assume that this density is represented by a distribution of markers, and that one wishes to use the marker distribution to obtain a numerical estimate of the density $\rho$ on a grid of level $p$. With the deposition algorithm described in the preceding section, it is seen that the markers contributing to the reconstructed value of the density at grid point $x^{(j)}$ are located in the interval $I = [x^{(j)} - h, x^{(j)} + h]$, where $h = 2^{-p}$. Therefore, since they are distributed according to $\rho$, the density of markers in the interval $I$ varies and is not uniformly equal to $\rho(x^{(j)})$. For instance, for the first grid point in $x = 0$, the density of the markers closest to $x = 0$ have a density equal to $1 + a$, but the density of the markers in $x = h$ is lower. In particular, when $k = k_{max} = 2^{p-1}$, the density approaches $1 - a$ when $x$ approaches $h$, which explains why the effect is strongest for fast varying modes.

On a grid of level $p$, the fastest mode is $k_{max} = 2^{p-1}$, for which the reduction factor is equal to $f(2^{p-1}, p) = \text{sinc}^2(1/2) = 0.405$. If $k$ is larger than $k_{max}$, the reduction factor is still given by equation (5), tends to zero as $1/k^2$, and is exactly zero every time $k$ is a multiple of $k_{max}$.

## 2.3. Analytical definition of the filter in two and three dimensions

The signal reduction factor of the one dimensional case can be used in a straightforward manner to compute the signal reduction factor in two and three dimensions. Any periodic signal can be decomposed in products of sine or cosine of each coordinate $x$, $y$ and $z$. Therefore, the reduction factors simply accumulate multiplicatively, so that the reduction factors in two (resp. three) dimensions, for a grid of level $(p_x, p_y)$ (resp. $(p_x, p_y, p_z)$) write

$$f_{2D}(k_x, k_y; p_x, p_y) = f(k_x, p_x)f(k_y, p_y) \qquad (6)$$

$$f_{3D}(k_x, k_y, k_z; p_x, p_y, p_z) = f(k_x, p_x)f(k_y, p_y)f(k_z, p_z) \qquad (7)$$

Note that as a result, the standard PIC method already acts as a low-pass filter, especially in higher dimensions. Indeed, on a grid of level $p$ in a $d$-dimensional code, the fastest modes (which combine fast spatial frequencies in all directions) are reduced by a factor of $(\text{sinc}^2(1/2))^d$. This factor is 0.164 in two dimensions and 0.067 in three dimensions.

Now, let us turn to the sparse grid PIC algorithm, in the case where the finest grid (that on which the evaluation of the moments is desired) is of level $p$. To evaluate a moment, such as, for example, the density,

the standard method is replaced with the following sparse grid combination formula:

$$\rho_{\text{sg},2D}^{(p)} = \sum_{i=1}^{p} \rho_{(i,j)}^{(p)} - \sum_{i=1}^{p-1} \rho_{(i,j)}^{(p)} \qquad (8)$$

$$\rho_{\text{sg},3D}^{(p)} = \sum_{i=1}^{p} \sum_{\substack{j=1 \\ i+j+k=p+2}}^{p+1-i} \rho_{(i,j,k)}^{(p)} - 2 \sum_{i=1}^{p-1} \sum_{\substack{j=1 \\ i+j+k=p+1}}^{p-i} \rho_{(i,j,k)}^{(p)} + \sum_{i=1}^{p-2} \sum_{\substack{j=1 \\ i+j+k=p}}^{p-1-i} \rho_{(i,j,k)}^{(p)} \qquad (9)$$

In these expressions, $\rho_{(i,j)}^{(p)}$ is the density evaluated with a grid of level $(i, j)$, and then interpolated on the grid of level $p$ (and similarly in three dimensions). As explained in ref. [1], this estimate can reduce the PIC noise at a moderate cost in terms of numerical complexity. The increase of numerical complexity is due to the fact that each term of the sums in equations (8)-(9) requires a loop over all the markers of the grid. In two (resp. three) dimensions, the number of terms in the sum is $2p-1$ (resp. $1+3p(p-1)/2$). For instance, a $512 \times 512 \times 512$ ($p=9$) three-dimensional grid requires 109 such loops over the marker distribution, which can represent a large computational overhead.

Here, we show an alternative way of understanding the reason for the reduction of the PIC noise, which has the advantage of explicitly characterizing the filter. Indeed, the fact that the sparse grid acts as a filter is explained in ref. [2], but the filter is not quantified in Fourier space. In addition to providing a more precise understanding of the sparse grid PIC method, this provides an alternative, much more efficient way, of producing the same noise reduction effect. Namely, this can be done by applying the filter directly to the standard PIC estimate of the moments of the distribution function. Although Fourier transforms are required, one avoids looping a large number of times over the marker distribution.

Using the linearity of the Fourier decomposition, we see that a contribution to the total signal having mode numbers $k_x, k_y$ (resp. $k_x, k_y, k_z$ in three dimensions) will be reduced by the factors of equations (6)-(7), when computed on a grid of level $(p_x, p_y)$ (resp. $(p_x, p_y, p_z)$). Also, there will be some signal deformation and spectrum broadening due to the interpolation step. Indeed, take for instance a cosine mode at frequency $2^i$ evaluated on the grid of level $i+1$. This frequency is the fastest compatible with the Nyquist theorem, and it will appear as a triangle signal. Now, when this signal is interpolated (using linear interpolation) on a grid of level $j > i$, it will still appear as a triangle signal. But, on this grid of level $j$, the triangle at frequency $2^i$ Fourier decomposes into a spectrum of modes around the dominant frequency $2^i$. However, in the following discussion, we completely neglect this spectrum broadening phenomenon. We will see with concrete examples that it does not play a dominant role.

Therefore, neglecting the spectrum broadening associated with interpolation, in two (resp. three) dimensions, the mode with mode numbers $k_x, k_y$ (resp. $k_x, k_y, k_z$) is filtered by the sparse grid algorithm according to the following rules:

$$F_{2D}^{\text{sg}}(k_x, k_y; p) = \sum_{i=1}^{p} f_{2D}(k_x, k_y; i, j) - \sum_{i=1}^{p-1} f_{2D}(k_x, k_y; i, j) \qquad (10)$$

$$F_{3D}^{\text{sg}}(k_x, k_y, k_z; p) = \sum_{i+j+k=p+2} f_{3D}(k_x, k_y, k_z; i, j, k) - 2 \sum_{i+j+k=p+1} f_{3D}(k_x, k_y, k_z; i, j, k) + \sum_{i+j+k=p} f_{3D}(k_x, k_y, k_z; i, j, k) \qquad (11)$$
<figure><img src="image_4.png" /><figcaption>Fig. 1. Sparse grid filter in two dimensions for $p=6$ (a) and $p=7$ (b). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

On a given grid of level $p$, the $(2^{p-1} + 1)^d$ different possible values for the mode numbers are reduced by the sparse grid algorithm, according to equations (10)-(11), which can be computed analytically using equations (5)-(7).

## 2.4. Examples in two dimensions

The two dimensional sparse grid filter is plotted in Fig. 1 for grids of level $p=6$ (a) and $p=7$ (b). It is seen clearly that the filter essentially cancels all modes having moderate or high mode numbers in both directions. However, modes that vary rapidly in one direction, but are smooth in the other one, are kept. This property, which is well known to the users of the sparse grid PIC method, means that the algorithm can reduce a large fraction of the noise, while at the same time preserving sharp structures, with the condition that they are more or less aligned with the mesh.

At the same time, the quantitative characterization of the filter allows one to immediately see the downside of the algorithm. One must be cautious not to have physically meaningful modes present in the blue regions of Fig. 1. These modes would disappear from the simulation. This problem is heavily dependent on both the grid size and the geometry, as the following example shows.

Let us assume that a poloidal plane (small section of the torus) of a toroidal plasma is modeled using the sparse grid algorithm. The plasma density can exhibit structures that typically have a well defined angular mode number and a moderate radial extension. Therefore, let us assume that the density reads

$$ \rho(x, y) = 1 + ae^{-(r(x,y)-r_0)^2/(2\sigma^2)} \cos(m\theta(x, y)) \quad (12) $$

where $m=10$, $r_0=0.25$, $\sigma=0.01$, and

$$ x = \frac{1}{2} + r \cos \theta \quad (13) $$

$$ y = \frac{1}{2} + r \sin \theta \quad (14) $$

The spectrum of this density in $(r, \theta)$ variables is sharp. However, it is not the case when expressed in $(x, y)$ variables. Indeed, to produce one single mode with poloidal mode number $m$, a broad spectrum of $(k_x, k_y)$ modes has to be used. Therefore, when the density is collected according to the sparse grid algorithm, a large part of the signal can be significantly deformed, leading to artifacts. Fig. 2 shows these artifacts when the level of the grid is $p=6$. Fig. 2 a) shows equation (12), and Fig. 2 b) the expected result once the sparse grid algorithm is applied to collect the density. To obtain Fig. 2 b), the data of Fig. 2 a) is Fourier

analyzed (with a value of $k_{max} = 2^{p-1} = 32$), then each mode is reduced by the factor given by equation (10), and the signal is reconstructed in real space using this modified spectrum. Fig. 2 c) shows the density collected on the grid of level $p=6$ with a sampling of $N_p = 4 \times 10^7$ markers. Note that in order to draw the markers from this distribution, the Fourier analysis of equation (12) is also used. Cumulative distribution functions can be obtained analytically from the Fourier spectrum in order to draw $x$ randomly from the distribution integrated in $y$, and for each value of $x$, the cumulative distribution function in $y$ is used to draw $y$ randomly. Finally, Fig. 2 d) displays the result of the sparse grid algorithm, equation (8), still with the same grid of level $p=6$. The artifacts predicted by the filter in Fig. 2 b) appear clearly. Note that this confirms our earlier statement that the spectrum broadening associated with interpolation (not taken into account to produce Fig. 2 b)) has little influence on the final results.

In the case of Fig. 2, it is seen that the sparse grid algorithm significantly deforms the signal, instead of reducing its PIC noise. However, let us assume that when the structure of equation (12) appears, the grid actually used in the simulation is of level 10. With $N_p = 4 \times 10^7$, each cell contains roughly 40 markers only, which is not sufficient to obtain a sufficiently good signal to noise ratio with the standard PIC algorithm. However, the sparse grid algorithm is, in this case, very efficient at reducing the noise while at the same time keeping the structure without creating artifacts. This is seen in Fig. 3. The absence of artifact generation in this case is due to the larger grid. In this case, very few mode numbers end up in the region of Fourier space where the modes are reduced to 0.

This illustrates how the sparse grid PIC algorithm can be put to use if one knows the kind of structure to expect in the problem considered. More precisely, one ought to know whether the structures are aligned, or not, with the axes of the grid, and what is their characteristic size. Note that some approaches based on a mathematical error estimate may be able to provide an automatic selection of the best parameters to use for the filter. For example, in ref. [2], which introduced the truncated filter discussed in section 6, the truncation parameter is determined automatically at runtime using the numerical properties of the Fourier transform of the signal at a given instant.

## 3. Properties of the sparse grid filter in two and three dimensions

In this section, we present and discuss the main mathematical properties of the filter in two and three dimensions. The proofs are carried out in the appendices.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Sparse grid deformation of an island structure. Analytical signal (a), sparse grid filtered analytical signal (b), density moment collected with $N_p = 4 \times 10^7$ markers using standard PIC method on a grid of level $p = 6$ (c), and using the sparse grid combination technique with the same level for the finest grid (d).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Sparse grid PIC sampling of the island structure with $N_p = 4 \times 10^7$ markers, on a grid of level $p = 10$ (1024 × 1024).</figcaption></figure>

First of all, $\forall(k_x, k_y, p) \in \mathbb{N}^3$ such that $p \ge 2$ and $k_x, k_y \le 2^{p-1}$, we have $0 \le F_{2D}^{sg}(\mathbf{k}; p) \le 1$, and similarly in three dimensions. This is proven in appendix A.1. The proof is not entirely trivial, and derives from the fact that $\forall x, \text{sinc}^2(x) > \text{sinc}^2(2x)$. As explained in section 7, the property does not hold for all orders of the shape function. The filters are low-pass filters. In particular, we have the following remarkable values in two dimensions

$$F_{2D}^{sg}(2^{p-1}, 2^{p-1}; p) = 0 \qquad (15)$$

$$F_{2D}^{sg}(k_x, 0; p) = \text{sinc}^2(k_x 2^{-p}) \qquad (16)$$

$$\lim_{p \to \infty} F_{2D}^{sg}(\mathbf{k}; p) = 1. \qquad (17)$$

In three dimensions, we have

$$F_{3D}^{sg}(2^{p-1}, 2^{p-1}, 2^{p-1}; p) = 0 \qquad (18)$$

$$F_{3D}^{sg}(k_x, k_y, 0; p) = F_{2D}^{sg}(k_x, k_y; p) \qquad (19)$$

$$\lim_{p \to \infty} F_{3D}^{sg}(\mathbf{k}; p) = 1. \qquad (20)$$

Note that properties (16) and (19) are evidently symmetric with respect to the dimension of the vanishing wave number. Property (15) (resp. property (18)) is true because then, all terms of the sum in equation (10) (resp. equation (11)) vanish. Equation (16) comes from the observation that if $k_y = 0$, then $f_{2D}(k_x, k_y, i, j) = f(k_x, i)$. Equations (17) and (19) are proven in appendices A.2 and B, respectively. Equation (20) follows from the same considerations as the two dimensional case.

At this point, the reader may point out that the formalism is cumbersome, because the filter depends on the level $p$ of the grid. In fact, we show that the filter can be expressed in terms of $d$ generic continuous real functions, defined on $\mathbb{R}_+^{*d}$, where $d$ is the dimensionality of the problem.

Fig. 4 compares the filter for $p = 5$ (a) and $p = 7$ (b). One observes that the filters have a very similar aspect, although the mode resolution is four times coarser in the $p = 5$ case compared to the $p = 7$ case. In fact, when $k_x, k_y$ and $k_z$ are all non-vanishing integers, we have the following interesting property in two and three dimensions, $\forall m \in \mathbb{N}$:

$$F_{2D}^{sg}(2^m \mathbf{k}; p + 2m) = F_{2D}^{sg}(\mathbf{k}; p) \qquad (21)$$

$$F_{3D}^{sg}(2^m \mathbf{k}; p + 3m) = F_{3D}^{sg}(\mathbf{k}; p). \qquad (22)$$

This property is proven in appendix A.2 for the two dimensional case. The three dimensional case follows from identical considerations. Equation (21) explains why Figs. 4 a) and 4 b) look so similar. Using this property, in two (resp. three) dimensions, we can define two (resp. three) functions $G_{2D}^{(0)}$ and $G_{2D}^{(1)}$ (resp. $G_{3D}^{(0)}$, $G_{3D}^{(1)}$ and $G_{3D}^{(2)}$) in order to express the value of the filter for a wave vector $\mathbf{k}$, on a grid of any level $p \ge 2$. The functions are defined for $\mathbf{k} \in \mathbb{R}_+^{*d}$ as follows:
<figure><img src="image_7.png" /><figcaption>Fig. 4. Comparison between the sparse grid filters for $p = 5$ (a) and $p = 7$ (b).</figcaption></figure>

$$G_{2D}^{(i)}(\mathbf{k}) = \lim_{p \to \infty} F_{2D}^{\text{sg}}(2^p \mathbf{k}; 2p + i), \quad i = 0, 1 \qquad (23)$$

$$G_{3D}^{(i)}(\mathbf{k}) = \lim_{p \to \infty} F_{3D}^{\text{sg}}(2^p \mathbf{k}; 3p + i), \quad i = 0, 1, 2. \qquad (24)$$

It is straightforward to see that the limits of equations (23) and (24) exist when the components of $\mathbf{k}$ are all base 2 decimals or dyadic numbers. Indeed, in such a case there exists a critical value $p_c$ such that the components of $2^p \mathbf{k}$ are all integers if $p \ge p_c$. Using property (21) (resp. (22) in three dimensions), we see that the limit is well defined, and $G_{2D}^{(i)}(\mathbf{k}) = F_{2D}^{\text{sg}}(2^p \mathbf{k}; 2p_c + i)$ (resp. $G_{3D}^{(i)}(\mathbf{k}) = F_{3D}^{\text{sg}}(2^p \mathbf{k}; 3p_c + i)$ in three dimensions). The continuity of the $G$ functions on $\mathbb{R}^2$ is proven in appendix A.3 for the two dimensional case. The three dimensional case follows from identical considerations.

The advantage is that once the $G$ functions, which depend only on the dimension but not on the grid size, are known, the value of the reduction factor of any mode number on a grid of arbitrary level can be known. Namely, when $\mathbf{k}, p$ are non-vanishing integers (with the components of $\mathbf{k}$ being less than $2^{p-1}$), then we have

$$F_{2D}^{\text{sg}}(\mathbf{k}; 2p + i) = G_{2D}^{(i)}(2^{-p}\mathbf{k}), \quad i = 0, 1 \qquad (25)$$

$$F_{3D}^{\text{sg}}(\mathbf{k}; 3p + i) = G_{3D}^{(i)}(2^{-p}\mathbf{k}), \quad i = 0, 1, 2. \qquad (26)$$

Therefore, it is sufficient to tabulate a small number of real valued functions in order to find the impact of the sparse grid algorithm on a grid of arbitrary level. This is one of the main results of this paper. In practice, suppose that the maximum level considered in the problem is $p_{max}$, and one wants to be able to have the value of $F_{2D}^{\text{sg}}$ for any even value of $p < p_{max}$. Then one must first compute $G_{2D}^{(0)}$ up to $k_x, k_y = 2^{p_{max}/2-1}$ and with precision $\Delta k = 2^{-p_{max}/2}$. This is done simply, once and for all, by using equation (25), which means the values of $G_{2D}^{(0)}$ at this precision are simply given by $F_{2D}^{\text{sg}}(\mathbf{k}; p_{max}), \forall \mathbf{k} \in [1, 2^{p_{max}-1}]^2$. Then, the values for another choice of $p < p_{max}$, with $p$ even, will be given by picking the right value in the already computed set of values for $G_{2D}^{(0)}$, namely, $F_{2D}^{\text{sg}}(\mathbf{k}; p) = G_{2D}^{(0)}(2^{-p/2}\mathbf{k})$. This is already known since the components of $\mathbf{k}$ are integers so the components of $2^{-p/2}\mathbf{k}$ are multiple of $\Delta k$.

In the next section, we show how a simpler and better performing filter can be constructed, based on the ideas of the sparse grid combination technique. When this new filter is used, the previous discussion about the $G$ functions becomes more of mathematical rather than practical interest.

## 4. The binary sparse grid filter

In this section, we see how it is possible to considerably simplify (and improve) the filters, by studying how the sparse grid combination formulæ, equations (8)-(9), work out in Fourier space directly. A filter, based on the combination formula and taking the values either 0 or 1 in Fourier space (hence the name binary filter), will be defined.

So far, we have approached the problem by looking at the Fourier content of the moments collected using standard PIC techniques. Now, let us consider that it is actually possible to collect directly the Fourier coefficients of a signal represented by PIC markers. Namely, in one dimension for example, one can collect the mode $k$ cosine and sine coefficients of the $m^{th}$ velocity moment of the distribution by the following two sums

$$a_k^c = w_k \frac{N}{N_p} \sum_{\ell=1}^{N_p} v_{\ell}^m \cos(2\pi k x_{\ell}) \qquad (27)$$

$$a_k^s = w_k \frac{N}{N_p} \sum_{\ell=1}^{N_p} v_{\ell}^m \sin(2\pi k x_{\ell}), \qquad (28)$$

where

$$w_k = \begin{cases} 1 & \text{if } k=0 \text{ or } k=2^{p-1} \\ 2 & \text{if } 0 < k < 2^{p-1} \end{cases}. \qquad (29)$$

Incidentally, the corresponding shape function, which can be computed to be

$$S_F(x) = \sum_{k=0}^{2^{p-1}} w_k \cos 2\pi kx, \qquad (30)$$

is delocalized on the entire grid, since all markers contribute to all modes, and a sinusoidal wave is by nature delocalized on the grid. In fact, in order to derive the binary filter, we don't need to consider equations (27)-(30). We only need to consider the Nyquist theorem: if one wants to use the collected Fourier coefficients on a grid of level $p$, one needs to collect the modes up to $k = 2^{p-1}$.

We will restrict ourselves to the two dimensional case (there is no difficulty with the extension to the three dimensional case). The first step is to examine the meaning of $\rho_{(i,j)}^{(p)}$ in the context of direct Fourier decomposition. This quantity is the density evaluated on a grid of level $(i, j)$ and then interpolated on the grid of level $p$. On a grid of level $(i, j)$, the fastest modes available in the $x$ and $y$ directions are respectively $k_x = 2^{i-1}$ and $k_y = 2^{j-1}$. If the interpolation on the grid of level $p$ still uses the Fourier decomposition, rather than bilinear interpolations, $\rho_{(i,j)}^{(p)}$
<figure><img src="image_8.png" /><figcaption>Fig. 5. Selection of retained modes (nodes in contact with colored squares) in the Fourier sparse grid algorithm (p = 4).</figcaption></figure>

is exactly the same as if one had collected all the modes up to $k_x = k_y = 2^{p-1}$, and subsequently set all the modes having $k_x > 2^{i-1}$ or $k_y > 2^{j-1}$ to zero.

From this simple observation, one can easily make sense of equation (8) by rearranging the terms as follows (omitting the (p) exponent to make the notation lighter):

$$ \begin{aligned} \rho_{\text{sg},2D} = & \rho_{(1,p)} - \rho_{(1,p-1)} \\ & + \rho_{(2,p-1)} - \rho_{(2,p-2)} \\ & + \dots \\ & + \rho_{(p-1,2)} - \rho_{(p-1,1)} \\ & + \rho_{(p,1)}. \end{aligned} \quad (31) $$

A term such as $\rho_{(i,j)} - \rho_{(i,j-1)}$ is decomposed as follows. Both terms contain the same modes in the x direction. However, in the y direction, $\rho_{(i,j)}$ contains the modes up to $k_y = 2^{j-1}$, whereas $\rho_{(i,j-1)}$ contains the modes up to $k_y = 2^{j-2}$, so that $\rho_{(i,j)} - \rho_{(i,j-1)}$ contains all the $k_x$ modes up to $k_x = 2^{i-1}$ and all the $k_y$ modes between $k_y = 2^{j-2} + 1$ and $k_y = 2^{j-1}$. The result is seen graphically for $p = 4$ in Fig. 5. In this figure, the nodes represent the modes $k_x \in [0, 2^3]$, $k_y \in [0, 2^3]$. All the nodes in contact with colored squares are kept, all the others are set to zero. The filtering rule can be expressed in terms of a maximum value of $k_y$ for each value of $k_x$. The algorithm to find $k_{y,max}(k_x)$ is as follows:

$$ \begin{aligned} k_{y,max}(k_x = 0) &= 2^{p-1} \\ k_{y,max}(k_x = 1) &= 2^{p-1} \\ k_{y,max}(k_x = 2) &= 2^{p-2} \\ k_{y,max}(k_x = 3-k_x = 4) &= 2^{p-3} \\ k_{y,max}(k_x = 5-k_x = 8) &= 2^{p-4} \\ k_{y,max}(k_x = 9-k_x = 16) &= 2^{p-5} \\ k_{y,max}(k_x = 17-k_x = 32) &= 2^{p-6} \\ \dots \\ k_{y,max}(k_x = (2^{p-2} + 1)-k_x = 2^{p-1}) &= 2^0 \end{aligned} \quad (32) $$

It is interesting to compare the filtering rule (32) with the sparse grid filter derived in the context of standard moment deposition, equation (10). In Fig. 6, the comparison is performed for $p = 6$, where black squares represent the modes that are kept by the binary filter. One sees that the boundary between the modes that are set to zero and the ones that are kept, using rule (32), roughly coincides with the $F_{2D}^{sg} = 0.3$ contour line. A short calculation shows that the fraction of retained degrees

<figure><img src="image_9.png" /><figcaption>Fig. 6. Pcolor of $F_{2D}^{sg}$ together with the retained modes (black squares) in the Fourier formulation for the case $n = 6$. The red contour line is the isocontour $F_{2D}^{sg} = 0.3$.</figcaption></figure>

of freedom, compared to the total number of degrees of freedom, $2^{2p}$, on a two dimensional grid of level $p$, scales as $p^{2-p}$. As expected, most modes are set to zero, and increasingly so when $p$ increases. The gain in terms of PIC noise can become very large, at the expense of not describing a lot of modes. The question is, of course, whether these modes are physical and should actually be treated, or not.

In three dimensions, using the same rule for interpreting terms such as $\rho_{(i,j,k)}^{(p)}$, we obtain similar conclusions. The filter derived from the three dimensional combination formula, equation (9) becomes a binary filter, with coefficients either 0 or 1, that is, some modes are retained, while others are simply discarded. This binary filter, which we shall denote $\bar{F}_{3D}^{sg}$ in the following, acts qualitatively in the same way as the sparse grid filter discussed in sections 2 and 3. Compared to the latter, it has two advantages. The normal sparse grid filter, $F_{3D}^{sg}$, is more difficult to compute as it requires numerous cardinal sine computations. The numerical cost is never quite prohibiting until the grid level exceeds $p = 10$, which is usually sufficient, nonetheless, it can take a few minutes unless the task is parallelized (which is easy). Another, more important, advantage is that in the case of $F_{3D}^{sg}$, the only mode that is not affected has $k_x = k_y = k_z = 0$. Any other mode is reduced by a factor strictly less than one. This is not the case with the binary filter, $\bar{F}_{3D}^{sg}$, where the retained modes do not undergo any deformation.

We now apply the ideas developed in this article to the PIC solver of the magnetohydrodynamic kinetic/fluid hybrid code XTOR-K.

## 5. Application to the XTOR-K hybrid magnetohydrodynamic code

### 5.1. Current moment collection strategy in XTOR-K

The XTOR-K code [24] solves a set of two-fluid visco-resistive magnetohydrodynamics (MHD) equations in interaction with a population of markers evolving in the fluid electromagnetic field. It is a hybrid kinetic/MHD extension of the two-fluid code XTOR-2F [25]. The coupling between the fluid and the marker distribution is carried out through the particle pressure tensor, which is injected in the momentum equation. The fluid grid of XTOR-K is based on the flux coordinates of the underlying equilibrium, computed with the CHEASE code [26]. XTOR-K uses a spectral description in the poloidal ($\theta$) and toroidal ($\phi$) directions, and finite differences in the radial ($s \propto \sqrt{\psi}$, where $\psi$ is the initial equilibrium poloidal flux) direction. The moments of the marker distribution are collected on a specific grid, that is different from the fluid grid. The collection grid is Cartesian in the poloidal plane, which bypasses any issue associated with the singular jacobian at the toroidal grid axis
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Radial structure of the $m=1, n=1$ mode. Filtering before the mapping and with $\bar{F}_{3D}^{sg}$ (a), after the mapping and with $\bar{F}_{3D}^{sg}$ (b), before the mapping and with $F_{3D}^{sg}$ (c), after the mapping and with $F_{3D}^{sg}$ (d).</figcaption></figure>

(which coincides with the magnetic axis at equilibrium). The spectral representation of the fluid code is partially used here, in that the represented toroidal modes are collected one by one on the two dimensional poloidal grid, applying equations (27)-(28). Jointly, the toroidal modes represent the full three dimensional signal. Once the moments (density, velocity and pressure tensor) are known on the collection grid, they are interpolated on the toroidal fluid grid, and finally projected on the set of represented modes. A typical resolution of the fluid grid in the radial, poloidal and toroidal direction is $\ell_{max} \times m_{max} \times n_{max} = 512 \times 64 \times 24$. This typical resolution reflects the fact that the XTOR code (in its purely fluid or hybrid version) is used to study MHD or kinetic MHD instabilities such as the kink [27,28], the tearing [29] or the fishbone mode [24], which have well defined low order angular mode numbers, but have a radial structure exhibiting sharp features around the resonant surface of the instability.

The set of angular modes retained in the simulations is defined by a band of poloidal modes $m$ around each toroidal mode $n > 0$, such that $n - m_{inf} < m < n + m_{sup}$. The $n=0$ case receives a separate treatment, where the poloidal modes are retained up to $m = m_{n0}$. Due to aliasing constraints in the context of nonlinear equations, the largest toroidal mode number (resp. $m_{n0}$) is less than one third of the number of toroidal (resp. poloidal) grid points. For instance when the toroidal resolution is 12 (resp. 24), then the largest toroidal mode is effectively $n=3$ (resp. $n=7$). Therefore, the number of dofs in the angular dimensions is reduced by a factor of roughly three, down from the already moderate value of $m_{max} \times n_{max}$. Since the moments are projected on this set of physical modes, the PIC noise is in part filtered. The problem resides in the radial direction, where the PIC noise is not filtered. The current filtering strategy is to apply a diffusion operator but its main purpose is to clean the noise only close to the axis. The presence of sharp features of the modes at the instability resonant surface (which justifies

the choice of the resolution) means that a strategy of noise reduction based on high frequency filtering (low-pass filter) of the radial profiles only is bound to fail. Indeed, consider what would happen if one would select the low order $(m,n)$ values corresponding to the mode of interest (and its sidebands), and filter the high radial frequencies out. We would also significantly smooth any sharp behavior happening at the resonance, thus strongly modifying the physics. In section 5.3, we will see what strategy we suggest to limit this phenomenon. First, we will examine the issue of energy and momentum conservation in relation with the filtering.

## 5.2. Energy and momentum conservation in XTOR-K and the sparse grid filter

Rigorous conservation of energy and momentum is often less a concern for MHD stability codes, in particular with implicit formulations, than it is for shock problems formulated with explicit schemes [30]. As a matter of fact, it is not easy to write, analytically, a conservative system of even purely fluid extended MHD equations. Naturally, this becomes even more difficult when kinetic physics is added. For example, in general, the noise inherent to PIC codes affects the conserved quantities. A filter such as the one discussed in this paper may also have a detrimental effect on energy and momentum conservation. Consider a specific mode carrying some energy, for instance in the form of internal energy, that is, pressure. From the point of view of the fluid equations, this particle pressure is seen as a moment on a grid, that may be subject to filtering. Following the nonlinear evolution, this energy may be redistributed to different wave numbers, where the value of the filter differs from the original wave number, and some energy will unphysically escape from the system. However, notice that the moment collection step, whatever order it is, already acts as a wave number dependent filter (equation
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Radial structure of the *m* = 1, *n* = 0 mode. Filtering before the mapping and with $\overline{F}_{3D}^{sg}$ (a), after the mapping and with $\overline{F}_{3D}^{sg}$ (b), before the mapping and with $\overline{F}_{3D}^{sg}$ (c), after the mapping and with $\overline{F}_{3D}^{sg}$ (d).</figcaption></figure>

(5)), no matter whether a sparse grid algorithm is implemented or not. Thus, the situation is not necessarily worse off with the filter than it is without it. In XTOR-K, the strategy to avoid that non conservation leads to unphysical situations on the long run is to implement restoring terms. For example, the heat source is implemented as diffusion with respect to an equilibrium, so that any unphysical energy generated by the system always ends up escaping the simulated volume. In this context, it is likely that the effect of an eventual filter on the conservation properties of the model exists, but is not of great concern.

We now explain how the ideas developed in this manuscript suggest filtering strategies to reduce the noise in the radial direction.

### 5.3. Filtering of an XTOR-K distribution

We use an actual marker distribution coming from XTOR-K, in order to evaluate under what conditions the sparse grid filter is able to reduce the PIC noise without deforming the signal. We consider a run where an *m* = 1, *n* = 1 internal kink is simulated in presence of a population of thermal *α* particles of temperature *T*<sub>α</sub> = 2 MeV, while the electron and ion bulk temperature are *T*<sub>e</sub> = 30 keV. The *α* density, *n*<sub>α</sub> = 4 × 10<sup>17</sup> m<sup>-3</sup> in the core is 2% of the bulk density, with similar profile. The *β* of the *α* particles, *β*<sub>α</sub> = 1.2% is close to the fluid *β*, *β*<sub>f</sub> = 1.7% (*β* is the ratio of kinetic to magnetic pressure). In this configuration, the effect of the particles on the kink is significant. The total number of markers is *N*<sub>p</sub> = 172 × 10<sup>6</sup>. The question is whether the radial noise can be reduced on the modes. At *t* = 5000τ<sub>A</sub>, where τ<sub>A</sub> is the Alfvén time, the mode has grown out of the PIC noise and the growth rate of the mode is clearly measurable on the integrated magnetic energies. At this time, we extract the global distribution of markers to compute its moments. We

will be able to study the different poloidal and toroidal mode numbers separately.

The moments are first collected on a cartesian grid in (R, φ, Z) coordinates. As explained in section 5.1, this avoids the issue of the magnetic axis, where the volume of cells tends to zero, increasing the noise. Also, notice that the cartesian grids are slightly larger than the domain of physical interest, so at the boundary of the grids, the signals vanish, which makes periodic boundary conditions suitable. The level of the grid is *p* ∈ {6, 7, 8}. Then, the moments are mapped by linear interpolation to the toroidal coordinates (*s*, *θ*, *φ*) used by the fluid part of the code. The toroidal grid has the same size as the deposition (R, φ, Z) grid (eventually the poloidal and toroidal modes are filtered to retain only the band of simulated modes). We call this step the mapping. The filter can be applied either before or after the mapping. The result will not be the same, because in the first case, the grid is not aligned with the modes present in the distribution, which typically have a well defined periodicity in *θ* and *φ*. Therefore, it would seem more natural to filter after the mapping step. However, the most important modes have low poloidal and toroidal mode numbers *m* and *n*. Therefore, the fast radial frequencies, which contribute to the radial noise, will tend to be well preserved, while they will be cut down for the modes having large *m* and *n*. This is not what we hope to achieve, since we would like the fast frequencies corresponding to noise to be filtered even for low *m* and *n* modes. We see that the order between filtering and mapping matters much, and that it might be preferable to filter before the mapping, where the grid is, in a sense, agnostic regarding the nature of the structures. In particular, we can hope that the sharp features associated with resonance layers will be preserved, at least in part, because they have a broad spectrum.
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Radial structure of the *m* = 2, *n* = 1 mode. Filtering before the mapping and with $\overline{F}_{3D}^{sg}$ (a), after the mapping and with $\overline{F}_{3D}^{sg}$ (b), before the mapping and with $F_{3D}^{sg}$ (c), after the mapping and with $F_{3D}^{sg}$ (d).</figcaption></figure>

For each choice of the order between the filtering and mapping steps, we test both $F_{3D}^{sg}$ and $\overline{F}_{3D}^{sg}$, that is, the filter derived using formula (11) or using the three dimensional counterpart of (32). Finally, the radial structure of any $(m, n)$ mode can be extracted from the filtered signal. For each $(m, n)$ mode, we show how the filtered signal depends on $p$, the level of the grid; every time there are four plots, for each choice of the filter and the order between mapping and filtering. All plots show in black the original signal in XTOR-K. The other curves show the cases $p=6$ (green), $p=7$ (red), $p=8$ (blue). In the top row (subfigures a) and b)) $\overline{F}_{3D}^{sg}$ is used, while in the bottom row, $F_{3D}^{sg}$ is used. In the left column (subfigures a) and c)), the filtering is done before the mapping, while it is done after the mapping in the right column. Figs. 7–10 show the radial structure of the (1, 1), (1, 0), (2, 1) and (2, 2) modes respectively. The other modes are drowned in the noise. Note that our analysis of the noise is qualitative only, because a quantitative study would require to know the underlying true signals, which we do not know.

The bottom rows of these two figures show an interesting behavior. When $F_{3D}^{sg}$ is used, a significant reduction of the signal is observed (Figs. 7, 8, 9 and 10 c)-d)). This reduction of the signal barely shows up for $\overline{F}_{3D}^{sg}$ (Figs. 7, 8, 9 and 10 a)-b)) even at $p=6$ (but we would start to see it for lower values of $p$). The fact that lower values of $p$ lead to stronger deformation of the signal is very natural in view of the properties of the sparse grid filters. The reason for the better robustness of the binary filter $\overline{F}_{3D}^{sg}$ is the following: although it has qualitatively similar properties to $F_{3D}^{sg}$ in terms of noise reduction, in the zone where $F_{3D}^{sg}$ takes up values significantly above zero, $\overline{F}_{3D}^{sg}$ is exactly equal to one (Fig. 6). This means that the retained modes are well preserved for

$\overline{F}_{3D}^{sg}$, whereas they are reduced by some factor in the case of $F_{3D}^{sg}$. This shows that it is largely preferable to use $\overline{F}_{3D}^{sg}$ instead of $F_{3D}^{sg}$.

Now, we compare the left and right columns of the figures (top row), that is, whether the filtering is done before or after the step of mapping to toroidal coordinates. Therefore, we compare the subfigures a) and b) in each figure. In the first case (filtering before the mapping), a good level of filtering is obtained for $p=6$ and $p=7$, even for $p=8$. However, the use of $p=6$ should be prohibited, because we see that the noise reduction comes with a diffusion of the sharp features. Indeed, the MHD modes should be allowed to vary rapidly, close to the resonant surface. The choices of $p=7$ or $p=8$ appear to be a good compromise between noise reduction and good representation of the structures. In other words, it reduces the particle based error without increasing too much the grid based error. In the case when the filter is applied after the mapping, we see that $p=7$ leads to almost no noise reduction, while $p=8$ even seems to increase the noise (actually, this is due to the fact that the additional step of applying a radial diffusion operator, which is used for the XTOR-K data in black, was not used for the sparse grid filtered data). The case of $p=6$, when the filter is applied after the mapping, corresponds roughly to the same level of noise as the case of $p=7$ when the filter is applied before. However, it looks preferable to use the latter, because the noise increases much faster with $p$ in the former case (filter after mapping), so the likelihood of mistuning the algorithm is larger.

In conclusion, the choice of $p=7$ or $p=8$ and of filtering the moments before the mapping appears to be a reasonable choice in XTOR-K. It remains to be fully implemented in the code, in order to perform systematic comparisons. First, the linear growth rates will have to be compared, then, possible differences in the nonlinear physics will have
<figure><img src="image_13.png" /><figcaption>Fig. 10. Radial structure of the m = 2, n = 2 mode. Filtering before the mapping and with $\overline{F}_{3D}^{sg}$ (a), after the mapping and with $\overline{F}_{3D}^{sg}$ (b), before the mapping and with $F_{3D}^{sg}$ (c), after the mapping and with $F_{3D}^{sg}$ (d).</figcaption></figure>

to be tracked before it can be envisaged to routinely use sparse grid filtering.

## 6. Truncated combination schemes

In ref. [2], filtering strategies based on a modification of the combination formulae (8)-(9) are considered. Namely, an integer offset, or truncation, parameter $\tau$ is introduced such that the minimum level of the grids used in the combination becomes $\tau$ instead of 1. Namely, formula (8) is replaced with

$$ \rho_{\text{sg},2D}^{(p,\tau)} = \sum_{i=\tau}^{p} \rho_{(i,j)}^{(p)} - \sum_{i=\tau}^{p-1} \rho_{(i,j)}^{(p)} \quad (33) $$

$$ \begin{aligned} \rho_{\text{sg},3D}^{(p,\tau)} = & \sum_{i=\tau}^{p} \sum_{j=\tau}^{p+1-i} \rho_{(i,j,k)}^{(p)} \\ & - 2 \sum_{i=\tau}^{p-1} \sum_{j=\tau}^{p-i} \rho_{(i,j,k)}^{(p)} \\ & + \sum_{i=\tau}^{p-2} \sum_{j=\tau}^{p-1-i} \rho_{(i,j,k)}^{(p)}. \end{aligned} \quad (34) $$

In ref. [2], $\tau$ is chosen dynamically at run time, depending on an estimation of the value of $\tau$ that minimizes the grid based error. This $\tau$ optimization step involves the Fourier transform of the signal obtained by moment deposition on the finest grid level. Then, it is suggested to deposit the moments on all the sparse subgrids of the algorithm,

and combine them with equation (33). Note that each deposition on a subgrid involves a loop through all the markers of the simulation. As emphasized above, in three dimensions the number of terms in the combination formula scales with $p^2$, where $p$ is the level of the finest grid, with for example 109 terms (109 loops over all markers) for $p = 9$. We suggest that it would be much faster, especially if the Fourier transform of the signal has already been obtained, to simply apply the filter corresponding to equation (33) or its three dimensional counterpart (34). They can be obtained for any level and any value of $\tau$, with the following formula:

$$ \begin{aligned} F_{2D,\text{tr}}^{\text{sg}}(\mathbf{k}; p, \tau) = & \sum_{i=\tau}^{p} \text{sinc}^2(k_x 2^{-i}) \text{sinc}^2(k_y 2^{-j}) \\ & - \sum_{i=\tau}^{p-1} \text{sinc}^2(k_x 2^{-i}) \text{sinc}^2(k_y 2^{-j}) \end{aligned} \quad (35) $$

and similarly in three dimensions. Note that to distinguish this filter from that of equation (10), we have added the subscript tr to indicate that it uses the truncated combination formula. In particular we have $F_{2D,\text{tr}}^{\text{sg}}(\bullet; \bullet, 1) = F_{2D}^{\text{sg}}$. The maps of the filter for order $q = 1$ and grid level $p = 6$, for the cases $\tau = 1$ to $\tau = 4$, are displayed in Fig. 11. In fact, once again the formalism of section 3 can be used to obtain a unified picture of the filters. Let $k_x, k_y, k_z \in \mathbb{N}^*$ and $p \ge 2$, $1 \le \tau \le p-1$. Then, by rearranging the sum (35) (let $i' = i + \tau - 1, j' = j + \tau - 1$) we find

$$ F_{2D,\text{tr}}^{\text{sg}}(\mathbf{k}; p+1-\tau) = F_{2D,\text{tr}}^{\text{sg}}(2^{\tau-1}\mathbf{k}; p, \tau) \quad (36) $$

Therefore, we have

$$ F_{2D,\text{tr}}^{\text{sg}}(2^{\tau-1}\mathbf{k}; p, \tau) = G_{2D}^{(i)}(2^{-r}\mathbf{k}), \quad (37) $$
<figure><img src="image_14.png" /><figcaption>Fig. 11. Filters for order $q = 1$ and grid level $p = 6$, for $\tau = 1$ (a), $\tau = 2$ (b), $\tau = 3$ (c) and $\tau = 4$ (d).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Selection of retained modes (nodes in contact with colored squares) in the Fourier sparse grid algorithm ($p = 4, \tau = 2$).</figcaption></figure>

where $r, i$ are the quotient and remainder of the euclidean division of $p + 1 - \tau$ by 2. Unfortunately, the equality does not hold when the arguments of $F_{2D,\text{tr}}^{\text{sg}}$ are no longer a multiple of $2^{\tau-1}$, that is, in general

$$F_{2D,\text{tr}}^{\text{sg}}(\mathbf{k}; p, \tau) \neq G_{2D}^{(i)}(2^{-(r+\tau-1)}\mathbf{k}). \quad (38)$$

Nonetheless, the preceding equality almost holds in practice, so that $G_{2D}^{(0)}$ and $G_{2D}^{(1)}$ give the general qualitative features of the truncated filters, owing to equation (36), as can be observed in Fig. 11.

Naturally, the truncated combination formula also allows to redefine the binary filter of section 4. For instance, for $p = 4$ and $\tau = 2$, the terms combine in the way represented in Fig. 12 (to be compared with Fig. 5).

## 7. Extension to higher order shape functions

So far, we have restricted the discussion to the case where the shape function is based on the tent function, eq. (2), because this is a very common choice in PIC algorithms. In fact, the tent function is but the order 1 of a series of shape functions built by successive convolutions of the hat function $\tau_0$:

$$\tau_0(x) = \begin{cases} 1 & \text{if } |x| \le 1/2 \\ 0 & \text{if } |x| > 1/2. \end{cases} \quad (39)$$

First, the tent function $\tau_1$ is the convolution product

$$\tau_1(x) = \int_{-\infty}^{\infty} dt \tau_0(x-t) \tau_0(t). \quad (40)$$

More generally, the shape function of order $q$ is given recursively by

$$\tau_q(x) = \int_{-\infty}^{\infty} dt \tau_0(x-t) \tau_{q-1}(t). \quad (41)$$

The Fourier transform of the hat function is the cardinal sine, which explains why that of the tent function is the square of the cardinal sine. More generally, formula (41) immediately leads by induction to the conclusion that the Fourier transform of $\tau_q$ is the cardinal sine to the power $q + 1$. Therefore, the results of this paper can easily be extended to the case of shape functions of any order $q$, if one replaces equation (5) with

$$f(k, p) = \text{sinc}^{q+1}\left(\frac{k}{2^p}\right). \quad (42)$$

This formula has an important consequence. We know that $\forall x$, $\text{sinc}^2(x) > \text{sinc}^2(2x)$, therefore $\forall q \in \mathbb{N}^*$, $\forall x$, $\text{sinc}^{2q}(x) > \text{sinc}^{2q}(2x)$. This property is essential to prove the positivity of the filter when the shape function is built on the tent function, and so this property is
<figure><img src="image_16.png" /><figcaption>Fig. 13. Sparse grid filter for a grid level $p=6$ when the shape function is the hat function (order 0).</figcaption></figure>

preserved for all other odd orders of the shape functions. However, $\mathrm{sinc}(x) - \mathrm{sinc}(2x) = \mathrm{sinc}(x)(1 - \cos \pi x)$, which is of the sign of $\sin \pi x$ for $x > 0$. Therefore, the property that is essential to prove the positivity of the filter is lost for all even powers of the order of the shape function. As a matter of fact, one can check that the filters associated with the shape functions of even order have negative values. For example, the case of order 0 (deposition using the hat function, or, equivalently, deposition on the nearest neighboring grid point) for $p = 6$ ($k_{\max} = 32$) is plotted in Fig. 13. Fortunately, it looks like the minimum (negative) of the filter tends to be reduced in absolute value as $p$ increases. Hence this issue may not be very problematic in practice.

## 8. Summary

In this paper, we have interpreted the action of the sparse grid combination technique on moments of markers distribution, in the context of the PIC Method, as a filter in Fourier space. The essential idea is to remember that the moment collection step is a convolution between the signal represented by the markers and the shape function. The Fourier transform of the shape function is the cardinal sine function elevated to some power, and the filter becomes a combination of sums and products of cardinal sines. We have shown that the filters can be expressed in terms of functions that are independent of the level of the grid. With the quantitative characterization of the filters, it becomes possible to replace the sparse grid combination formula, which involves numerous sums on the markers, with the sparse grid filter, involving only a forward and backward fast Fourier transform. In addition, it is shown how any sparse grid combination formula, including truncated ones, can be translated into a filter with similar qualitative properties, but with values either 0 or 1. The advantage of this binary filter is that it is easier to compute and preserves the large scale structures with more fidelity than the filter derived from the cardinal sine factors. The application of these factors to the hybrid fluid/kinetic magnetohydrodynamic code XTOR-K yields encouraging results.

### CRediT authorship contribution statement

T. Nicolas: Conceptualization, Formal analysis, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing. V. Dubois: Formal analysis, Investigation. Q. Fang: Formal analysis. H. Lütjens: Software, Validation.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgement

The main author would like to thank G. Fubiani and L. Garrigues for interesting discussions about the sparse grid method.

## Appendix A. General properties of the two dimensional filter

$$A.1. \ Proof \ that \ 0 \le F_{2D}^{sg}(k; p) \le 1$$

Let us detail here general results about the two dimensional Filter $F_{2D}^{sg}$ defined in eq. (10). First of all, one can use Eqs. (5), (6) in order to rewrite eq. (10) into

$$F_{2D}^{sg}(k_x, k_y; p) = \sum_{i+j=p+1} \mathrm{sinc}^2(k_x 2^{-i}) \mathrm{sinc}^2(k_y 2^{-j}) \\ - \sum_{i+j=p} \mathrm{sinc}^2(k_x 2^{-i}) \mathrm{sinc}^2(k_y 2^{-j}). \qquad (A.1)$$

Eq. (A.1) is very convenient because it makes it clear that:

* $F_{2D}^{sg}$ can be extended by identification to a function defined on $\mathbb{R}^2 \times \mathbb{N} \setminus \{0,1\} \to \mathbb{R}$. Besides, for all $p \in \mathbb{N} \setminus \{0,1\}$, $F_{2D}^{sg}(\bullet, \bullet; p) \in C^{\infty}(\mathbb{R}^2, \mathbb{R})$. This comes from the fact that the sinc function is infinitely continuous over $\mathbb{R}$.

* $\forall p \ge 2, \forall k_1, k_2 \in \mathbb{R}, F_{2D}^{sg}(k_1, k_2; p) = F_{2D}^{sg}(k_2, k_1; p)$.

* $\forall p \ge 2, \forall k \in \mathbb{R}, F_{2D}^{sg}(0, k; p) = \mathrm{sinc}^2(k^2 p)$.

* $F_{2D}^{sg}(0, 0; p) = 1$.

Now, let us recombine the two sums of eq. (A.1) into one, with one additional term:

$$F_{2D}^{sg}(k_x, k_y; p) = \mathrm{sinc}^2\left(\frac{k_x}{2^p}\right) \mathrm{sinc}^2\left(\frac{k_y}{2}\right) + \sum_{i=1}^{p-1} \left[ \mathrm{sinc}^2\left(\frac{k_x}{2^i}\right) \right. \\ \left. \times \left( \mathrm{sinc}^2\left(\frac{2^i k_y}{2^{p+1}}\right) - \mathrm{sinc}^2\left(\frac{2^i k_y}{2^p}\right) \right) \right].$$

The terms inside the big parentheses are of the form $\mathrm{sinc}^2(x) - \mathrm{sinc}^2(2x) = (\pi x)^2 \mathrm{sinc}^4(x)$. Consequently, the previous expression reduces to

$$F_{2D}^{sg}(k_x, k_y; p) = \sum_{i=1}^{p-1} \mathrm{sinc}^2\left(\frac{k_x}{2^i}\right) \left(\frac{2^i \pi k_y}{2^{p+1}}\right)^2 \mathrm{sinc}^4\left(\frac{2^i k_y}{2^{p+1}}\right) \\ + \mathrm{sinc}^2\left(\frac{k_x}{2^p}\right) \mathrm{sinc}^2\left(\frac{k_y}{2}\right). \qquad (A.2)$$

Eq. (A.2) immediately shows that

$$\forall(k_x, k_y, n) \in \mathbb{R} \times \mathbb{N} \setminus \{0,1\}, F_{2D}^{sg}(k_x, k_y; n) \ge 0. \qquad (A.3)$$

It also implies that $F_{2D}^{sg}(k_x, k_y; p) = 0 \Rightarrow k_x/2^p \in \mathbb{Z}$ or $k_y/2^p \in \mathbb{Z}$.

The same trick used for eq. (A.2) can also be used to compute $H(k_x, k_y; p) = F_{2D}^{sg}(k_x, k_y; p) - F_{2D}^{sg}(2k_x, k_y; p)$:

$$H(k_x, k_y; p) = \sum_{i=1}^{p-1} \left[ \left( \mathrm{sinc}^2\left(\frac{k_x}{2^i}\right) - \mathrm{sinc}^2\left(\frac{2k_x}{2^i}\right) \right) \right. \\ \left. \times \left( \frac{2^i \pi k_y}{2^{p+1}} \right)^2 \mathrm{sinc}^4\left(\frac{2^i k_y}{2^{p+1}}\right) \right] \\ + \left( \mathrm{sinc}^2\left(\frac{k_x}{2^p}\right) - \mathrm{sinc}^2\left(\frac{2k_x}{2^p}\right) \right) \mathrm{sinc}^2\left(\frac{k_y}{2}\right)$$
$$= \sum_{i=1}^{p-1} \left( \frac{\pi^2 k_x k_y}{2^{p+1}} \right)^2 \text{sinc}^4 \left( \frac{k_x}{2^i} \right) \text{sinc}^4 \left( \frac{2^i k_y}{2^{p+1}} \right) \\ + \left( \frac{\pi k_x}{2^p} \right)^2 \text{sinc}^4 \left( \frac{k_x}{2^p} \right) \text{sinc}^2 \left( \frac{k_y}{2} \right).$$

Consequently, $\forall(k_x, k_y, p) \in \mathbb{R}^2 \times \mathbb{N} \setminus \{0, 1\}, H(k_x, k_y; p) \ge 0$. Hence, $F_{2D}^{sg}(k_x/2, k_y; p) \ge F_{2D}^{sg}(k_x, k_y; p)$. Using the symmetry of the function and iterating it, one has the following result: $\forall(k_x, k_y, p, m_1, m_2) \in \mathbb{R} \times \mathbb{N} \setminus \{0, 1\} \times \mathbb{N}^2$,

$$0 \le F_{2D}^{sg}(k_x, k_y; p) \le F_{2D}^{sg} \left( \frac{k_x}{2^{m_1}}, \frac{k_y}{2^{m_2}}; p \right). \quad (A.4)$$

Thus, letting $m_1, m_2 \to +\infty$ and using the fact that $F_{2D}^{sg}$ is continuous, one has $\forall(k_x, k_y, p) \in \mathbb{R} \times \mathbb{N} \setminus \{0, 1\}$,

$$0 \le F_{2D}^{sg}(k_x, k_y; p) \le 1. \quad (A.5)$$

## A.2. Proof of equations (17) and (21)

We start by proving equation (21), which gives a relation between the filters when the grid is refined. Property (17) will come as a consequence. Let us assume that $k_x, k_y \in \mathbb{Z}^*$ and $m, n \in \mathbb{N}$. Then, one can develop

$$F_{2D}^{sg}(2^m k_x, 2^n k_y; p+m+n) \\ = \sum_{i+j=p+m+n+1} \text{sinc}^2(k_x 2^{m-i}) \text{sinc}^2(k_y 2^{n-j}) \\ - \sum_{i+j=p+m+n} \text{sinc}^2(k_x 2^{m-i}) \text{sinc}^2(k_y 2^{n-j}) \\ = \sum_{l=m+1}^{m+p} \text{sinc}^2(k_x 2^{m-i}) \text{sinc}^2(k_y 2^{i-p-m-1}) \\ - \sum_{i=m+1}^{m+p-1} \text{sinc}^2(k_x 2^{m-i}) \text{sinc}^2(k_y 2^{i-p-m}) \\ + \sum_{i \in [1, m] \times [p+m+1, p+m+n]} \text{sinc}^2(k_x 2^{m-i}) \text{sinc}^2(k_y 2^{i-p-m-1}) \\ - \sum_{i \in [1, m] \times [p+m, p+m+n-1]} \text{sinc}^2(k_x 2^{m-i}) \text{sinc}^2(k_y 2^{i-p-m}) \\ = \sum_{i=1}^{p} \text{sinc}^2(k_x 2^{-i}) \text{sinc}^2(k_y 2^{-i-p-1}) - \sum_{i=1}^{p-1} \text{sinc}^2(k_x 2^{-i}) \text{sinc}^2(k_y 2^{-i-p}) \\ + \sum_{l=0}^{m-1} \text{sinc}^2 \left( \underbrace{k_x 2^l}_{\in \mathbb{Z}^*} \right) \left[ \text{sinc}^2(k_y 2^{-(l+p+1)}) - \text{sinc}^2(k_y 2^{-(l+p)}) \right] \\ + \sum_{t=0}^{n-1} \text{sinc}^2 \left( \underbrace{k_y 2^t}_{\in \mathbb{Z}^*} \right) \left[ \text{sinc}^2(k_x 2^{-(l+p+1)}) - \text{sinc}^2(k_x 2^{-(l+p)}) \right], \quad (A.6)$$

with $l=m-i$ and $t=i-p-m-1$ in (A.6). Since the sine function cancels out on $\pi\mathbb{Z}$, sinc is equal to 0 on $\mathbb{Z}^*$. Besides, the first term of (A.6) is precisely $F_{2D}^{sg}(k_x, k_y; p)$. This immediately leads us to the equality in 2D

$$\forall k_x, k_y \in \mathbb{Z}^*, \forall p \ge 2, \forall m, n \ge 0, \\ F_{2D}^{sg}(2^m k_x, 2^n k_y; p+m+n) = F_{2D}^{sg}(k_x, k_y; p). \quad (A.7)$$

Equation (21) is a direct application of eq. (A.7) to the special case $m=n$.

Finally, let us use eq. (A.6) but for $k_x, k_y \in \mathbb{R}$ instead of $\mathbb{N}$. It gives us a weaker, but more general result: $\forall(k_x, k_y, p, m) \in \mathbb{R}^2 \times \mathbb{N} \setminus \{0, 1\} \times \mathbb{N}$,

$$F_{2D}^{sg}(2^m k_x, 2^m k_y; p+2m) \ge F_{2D}^{sg}(k_x, k_y; p). \quad (A.8)$$

Eq. (A.8) can be rewritten into $F_{2D}^{sg}(k_x, k_y; p+2m) \ge F_{2D}^{sg}(2^{-m}k_x, 2^{-m}k_y; p)$. Using again the continuity of $F_{2D}^{sg}$, letting $m \to +\infty$ and considering odd and even values of $p$, one has the result

$$\forall(k_x, k_y) \in \mathbb{R}^2, F_{2D}^{sg}(k_x, k_y; p) \underset{p \to +\infty}{\longrightarrow} 1, \quad (A.9)$$

which proves eq. (17).

## A.3. Continuity of $G_{2D}^{(i)}$

Let us define for all $k_x, k_y \in \mathbb{R}$ and for all $p \ge 1$, $G_p(k_x, k_y) \equiv F_{2D}^{sg}(2^p k_x, 2^p k_y; 2p)$. According to eq. (A.7), $k_x, k_y \in 2^{-m}\mathbb{Z}^* \Rightarrow \forall p \ge 1$, $G_{p+m}(k_x, k_y) = G_m(k_x, k_y)$. Thus, it becomes natural to consider, provided it exists, the limit of the sequence $(G_p)$:

$$G_{2D}^{(0)}(k_x, k_y) = \lim_{p \to \infty} F_{2D}^{sg}(2^p k_x, 2^p k_y; 2p). \quad (A.10)$$

The convergence of $(G_p)$ can be proven as follows. Let us consider the sequence $C_{p,n} \equiv G_{p+n} - G_p$. Using again eq. (A.6), one has

$$C_{p,n}(k_x, k_y) = \sum_{l=0}^{n-1} \left[ \begin{aligned} & \text{sinc}^2(k_x 2^{p+l}) \left( \frac{\pi k_y}{2^{p+l+1}} \right)^2 \text{sinc}^4 \left( \frac{k_y}{2^{p+l+1}} \right) \\ & + \text{sinc}^2(k_y 2^{p+l}) \left( \frac{\pi k_x}{2^{p+l+1}} \right)^2 \text{sinc}^4 \left( \frac{k_x}{2^{p+l+1}} \right) \end{aligned} \right]. \quad (A.11)$$

One can immediately see that this series normally converges on $\mathbb{R}^2$ at fixed $p$ and when $n \to \infty$ since both terms inside the brackets behave like $O(2^{-2l})$. Thus, let us define its limit $C_p = \lim_{n \to \infty} C_{p,n}$. It is continuous thanks to the normal convergence and satisfies the relationship $\forall p \ge 1$, $G_{2D}^{(0)} = C_p + G_p$. As a consequence, $G_{2D}^{(0)}$ is a continuous function that can be expressed as the following series:

$$G_{2D}^{(0)}(k_x, k_y) = \text{sinc}^2(k_x) \text{sinc}^2(k_y) \\ + \sum_{l=0}^{\infty} \left[ \begin{aligned} & \text{sinc}^2(k_x 2^l) \left( \frac{\pi k_y}{2^{l+1}} \right)^2 \text{sinc}^4 \left( \frac{k_y}{2^{l+1}} \right) \\ & + \text{sinc}^2(k_y 2^l) \left( \frac{\pi k_x}{2^{l+1}} \right)^2 \text{sinc}^4 \left( \frac{k_x}{2^{l+1}} \right) \end{aligned} \right], \quad (A.12)$$

where we used $G_1(k_x, k_y) = \text{sinc}^2(k_x) \left(\frac{\pi k_y}{2}\right)^2 \text{sinc}^4\left(\frac{k_y}{2}\right) + \text{sinc}^2\left(\frac{k_x}{2}\right) \text{sinc}^2(k_y)$ and $\text{sinc}^2(x) = \text{sinc}^2(2x) + (\pi x)^2 \text{sinc}^4(x)$.

Finally, $G_{2D}^{(0)}$ satisfies the following properties:

* $G_{2D}^{(0)}$ is solution of eq. (A.10) by construction.

* $G_{2D}^{(0)} \in C(\mathbb{R}^2, \mathbb{R})$ thanks to uniform convergence (it is even infinitely continuous).

* $0 \le G_{2D}^{(0)} \le 1$ thanks to the properties of $F_{2D}^{sg}$.

* $\forall m \in \mathbb{N}, k_x, k_y \in 2^{-m}\mathbb{Z}^* \Rightarrow \forall p \ge 2$, $G_{2D}^{(0)}(k_x, k_y) = F_{2D}^{sg}(2^{p+m}k_x, 2^{p+m}k_y; 2(p+m))$.

* In the particular case $k_x k_y = 0$ (let us take $k_y = 0$), one has

$$C_{p,n}(k_x, 0) = \text{sinc}^2 \left( \frac{k_x}{2^{p+n}} \right) - \text{sinc}^2 \left( \frac{k_x}{2^p} \right) \\ \underset{n \to \infty}{\longrightarrow} 1 - \text{sinc}^2 \left( \frac{k_x}{2^p} \right). \quad (A.13)$$

Thus, $G_{2D}^{(0)}(k_x, 0) = G_{2D}^{(0)}(0, k_y) = 1$.

## Appendix B. General results on the 3D filter

The three dimensional filter $F_{3D}^{sg}$ defined in (11) verifies the same properties as its two dimensional counterpart. The proofs are more cum-
handsome but use the same methods. As in appendix A, let us rewrite eq. (11) with the help of eqs. (5) and (7):

$$
\begin{align}
F_{3D}^{sg}(\mathbf{k}, p) = & \sum_{i+j+l=p+2} \text{sinc}^2\left(\frac{k_x}{2^l}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \text{sinc}^2\left(\frac{k_z}{2^l}\right) \nonumber \\
& - 2 \sum_{i+j+l=p+1} \text{sinc}^2\left(\frac{k_x}{2^l}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \text{sinc}^2\left(\frac{k_z}{2^l}\right) \tag{B.1} \\
& + \sum_{i+j+l=p} \text{sinc}^2\left(\frac{k_x}{2^l}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \text{sinc}^2\left(\frac{k_z}{2^l}\right), \nonumber
\end{align}
$$

with $\mathbf{k} = (k_x, k_y, k_z)$. Eq. (B.1) makes it clear that:

• $F_{3D}^{sg}(\bullet, \bullet, \bullet; p)$ can be identified as an infinitely continuous function defined on $\mathbb{R}^3$ (with $p \ge 3$).

• $F_{3D}^{sg}(\bullet, \bullet, \bullet; p)$ is a symmetric operator.

• $F_{3D}^{sg}(\mathbf{0}; p) = 1$.

In order to retrieve the three dimensional counterpart of the results found in appendix A, one needs to define

$$
P_p(\mathbf{k}) \equiv \sum_{i+j+l=p+1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \text{sinc}^2\left(\frac{k_z}{2^l}\right) \\ - \sum_{i+j+l=p} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \text{sinc}^2\left(\frac{k_z}{2^l}\right). \qquad (B.2)
$$

Now, let us develop eq. (B.2):

$$
\begin{align*}
P_p(\mathbf{k}) &= \sum_{i=1}^{p-1} \sum_{j=1}^{p-i} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \text{sinc}^2\left(\frac{k_z 2^{i+j}}{2^{p+1}}\right) \\
&\quad - \sum_{i=1}^{p-2} \sum_{j=1}^{p-1-i} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \text{sinc}^2\left(\frac{k_z 2^{i+j}}{2^p}\right) \\
&= \sum_{i=1}^{p-2} \sum_{j=1}^{p-1-i} \left[ \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \right. \\
&\qquad \left. \left( \text{sinc}^2\left(\frac{k_z 2^{i+j}}{2^{p+1}}\right) - \text{sinc}^2\left(\frac{k_z 2^{i+j}}{2^p}\right) \right) \right] \\
&\quad + \sum_{i=1}^{p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y 2^i}{2^p}\right) \text{sinc}^2\left(\frac{k_z}{2}\right).
\end{align*}
$$

Using again the $\text{sinc}^2(x) - \text{sinc}^2(2x) = (\pi x)^2 \text{sinc}^4(x)$ identity, one has

$$
P_p(\mathbf{k}) = \sum_{2 \le i+j \le p-1} \left[ \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \right. \\
\left. \left( \frac{\pi k_z 2^{i+j}}{2^{p+1}} \right)^2 \text{sinc}^4\left(\frac{k_z 2^{i+j}}{2^{p+1}}\right) \right] \tag{B.3} \\
+ \sum_{i=1}^{p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y 2^i}{2^p}\right) \text{sinc}^2\left(\frac{k_z}{2}\right).
$$

Now, remarking that $F_{3D}^{sg}(\mathbf{k}; p) = P_{p+1}(\mathbf{k}) - P_p(\mathbf{k})$, we expand $F_{3D}^{sg}(\mathbf{k}; p)$ using eq. (B.3):

$$
\begin{align*}
F_{3D}^{sg}(\mathbf{k}; p) ={}& \sum_{i=1}^{p-1} \sum_{j=1}^{p-i} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \left(\frac{\pi k_z 2^{i+j}}{2^{p+2}}\right)^2 \text{sinc}^4\left(\frac{k_z 2^{i+j}}{2^{p+2}}\right) \\
& + \sum_{i=1}^{p} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{2^i k_y}{2^{p+1}}\right) \text{sinc}^2\left(\frac{k_z}{2}\right) \\
& - \sum_{i=1}^{p-2} \sum_{j=1}^{p-1-i} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2^j}\right) \left(\frac{\pi k_z 2^{i+j}}{2^{p+1}}\right)^2
\end{align*}
$$

$$
\begin{align*}
& \times \text{sinc}^4\left(\frac{k_z 2^{i+j}}{2^{p+1}}\right) \\
& - \sum_{i=1}^{p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{2^i k_y}{2^p}\right) \text{sinc}^2\left(\frac{k_z}{2}\right) \\
&= \sum_{i=1}^{p-2} \sum_{j=1}^{p-1-i} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \left(\text{sinc}^2\left(\frac{k_y}{2^{j+1}}\right) - \text{sinc}^2\left(\frac{k_y}{2^j}\right)\right) \\
& \quad \times \left(\frac{\pi k_z 2^{i+j}}{2^{p+1}}\right)^2 \text{sinc}^4\left(\frac{k_z 2^{i+j}}{2^{p+1}}\right) \\
& \quad + \sum_{i=1}^{p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2}\right) \left(\frac{2^i \pi k_z}{2^{p+1}}\right)^2 \text{sinc}^4\left(\frac{2^i k_z}{2^{p+1}}\right) \\
& \quad + \text{sinc}^2\left(\frac{k_x}{2^p}\right) \text{sinc}^2\left(\frac{k_y}{2}\right) \text{sinc}^2\left(\frac{k_z}{2}\right) \\
& \quad + \sum_{i=1}^{p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \left(\text{sinc}^2\left(\frac{2^i k_y}{2^{p+1}}\right) - \text{sinc}^2\left(\frac{2^i k_y}{2^p}\right)\right) \\
& \quad \times \text{sinc}^2\left(\frac{k_z}{2}\right) \\
&= \sum_{i+j \le p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \left(\frac{2^i \pi^2 k_y k_z}{2^{p+2}}\right)^2 \text{sinc}^4\left(\frac{k_y}{2^{j+1}}\right) \\
& \quad \times \text{sinc}^4\left(\frac{k_z 2^{i+j}}{2^{p+1}}\right) \\
& \quad + \sum_{i=1}^{p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \text{sinc}^2\left(\frac{k_y}{2}\right) \left(\frac{2^i \pi k_z}{2^{p+1}}\right)^2 \text{sinc}^4\left(\frac{2^i k_z}{2^{p+1}}\right) \\
& \quad + \text{sinc}^2\left(\frac{k_x}{2^p}\right) \text{sinc}^2\left(\frac{k_y}{2}\right) \text{sinc}^2\left(\frac{k_z}{2}\right) \\
& \quad + \sum_{i=1}^{p-1} \text{sinc}^2\left(\frac{k_x}{2^i}\right) \left(\frac{2^i \pi k_y}{2^{p+1}}\right)^2 \text{sinc}^4\left(\frac{2^i k_y}{2^{p+1}}\right) \text{sinc}^2\left(\frac{k_z}{2}\right).
\end{align*}
\tag{B.4}
$$

To obtain the first equality in equation (B.4), we have rewritten $\sum_{i=1}^{p-1} \sum_{j=1}^{p-i}$ as $\sum_{i=1}^{p-1} \sum_{l=0}^{p-1-i}$ with $l = j - 1$, and then renamed $l$ with the letter $j$. Eq. (B.4) shows that $F_{3D}^{sg} \ge 0$ and that, as for the two dimensional filter, one has $F_{3D}^{sg}(k_x, k_y, k_z; n) \ge F_{3D}^{sg}(2k_x, k_y, k_z; n)$. Thus, one has the three dimensional counterpart of eqs. (A.5) and (A.9):

$$
\forall \mathbf{k} \in \mathbb{R}^3, \forall n \ge 3, 0 \le F_{3D}^{sg}(\mathbf{k}; n) \le 1 \qquad (B.5)
$$

$$
\forall \mathbf{k} \in \mathbb{R}^3, F_{3D}^{sg}(\mathbf{k}; n) \xrightarrow{n \to +\infty} 1. \qquad (B.6)
$$

In addition, by setting $k_z$ to 0 in eq. (B.4), we find directly expres-
sion (A.2), which shows that

$$
F_{3D}^{sg}(k_x, k_y, 0; p) = F_{2D}^{sg}(k_x, k_y; p). \qquad (B.7)
$$

References

[1] L.F. Ricketson, A.J. Cerfon, Plasma Phys. Control. Fusion 59 (2016) 024002.

[2] S. Muralikrishnan, A.J. Cerfon, M. Frey, L.F. Ricketson, A. Adelmann, J. Comput. Phys. X 11 (2021) 100094.

[3] J.P. Boris, et al., in: Proc. Fourth Conf. Num. Sim Plasmas, 1970, pp. 3–67.

[4] O. Buneman, J. Comput. Phys. 1 (1967) 517.

[5] R.E. Denton, M. Kotschenreuther, J. Comput. Phys. 119 (1995) 283.

[6] S. Jolliet, A. Bottino, P. Angelino, R. Hatzky, T.-M. Tran, B. Mcmillan, O. Sauter, K. Appert, Y. Idomura, L. Villard, Comput. Phys. Commun. 177 (2007) 409.

[7] P. Bratley, B.L. Fox, ACM Trans. Math. Softw. 14 (1988) 88.

[8] R. Sydora, J. Comput. Appl. Math. 109 (1999) 243.

[9] D. Leblond, Simulation des plasmas de tokamak avec XTOR: régimes des dents de scie et évolution vers une modélisation cinétique des ions, Ph.D. thesis, Palaiseau, Ecole Polytechnique, 2011.

[10] C.K. Birdsall, A.B. Langdon, Plasma Physics via Computer Simulation, CRC Press, 2004.
[11] S. Gassama, É. Sonnendrücker, K. Schneider, M. Farge, M.O. Domingues, in: ESAIM: Proceedings, vol. 16, EDP Sciences, 2007, pp. 195–210.

[12] S. Butterworth, et al., Wirel. Eng. 7 (1930) 536.

[13] R.A. Roberts, C.T. Mullis, Digital Signal Processing, Addison-Wesley Longman Publishing Co., Inc., 1987.

[14] M. Griebel, M. Schneider, C. Zenger, 1990.

[15] H.-J. Bungartz, M. Griebel, Acta Numer. 13 (2004) 147.

[16] C. Kowitz, M. Hegland, Proc. Comput. Sci. 18 (2013) 449.

[17] M.M. Ali, P.E. Strazdins, B. Harding, M. Hegland, J.W. Larson, in: 2015 International Conference on High Performance Computing & Simulation (HPCS), IEEE, 2015, pp. 499–507.

[18] K. Kormann, E. Sonnendrücker, in: Sparse Grids and Applications-Stuttgart 2014, Springer, 2016, pp. 163–190.

[19] L. Garrigues, B. Tezenas du Montcel, G. Fubiani, F. Bertomeu, F. Deluzet, J. Narski, J. Appl. Phys. 129 (2021).

[20] L. Garrigues, B. Tezenas du Montcel, G. Fubiani, B. Reman, J. Appl. Phys. 129 (2021).

[21] F. Deluzet, G. Fubiani, L. Garrigues, C. Guillet, J. Narski, ESAIM: Math. Model. Numer. Anal. 56 (2022) 1809.

[22] F. Deluzet, G. Fubiani, L. Garrigues, C. Guillet, J. Narski, J. Comput. Phys. 480 (2023) 112022.

[23] L. Dalcin, M. Mortensen, D.E. Keyes, J. Parallel Distrib. Comput. 128 (2019) 137.

[24] G. Brochard, R. Dumont, H. Lütjens, X. Garbet, Nucl. Fusion 60 (2020) 086002.

[25] H. Lütjens, J.-F. Luciani, J. Comput. Phys. 229 (2010) 8130.

[26] H. Lütjens, A. Bondeson, O. Sauter, Comput. Phys. Commun. 97 (1996) 219.

[27] F.D. Halpern, H. Lütjens, J.-F. Luciani, Phys. Plasmas 18 (2011).

[28] T. Nicolas, H. Lütjens, J.-F. Luciani, X. Garbet, R. Sabot, Phys. Plasmas 21 (2014).

[29] D. Meshcheriakov, P. Maget, H. Lütjens, P. Beyer, X. Garbet, Phys. Plasmas 19 (2012).

[30] L. Chacón, Comput. Phys. Commun. 163 (2004) 143.