A particle-in-Fourier method with semi-discrete energy conservation for non-periodic boundary conditions

We introduce a novel particle-in-Fourier (PIF) scheme that extends its applicability to non-periodic boundary conditions. Our method handles free space boundary conditions by replacing the Fourier Laplacian operator in PIF with a mollified Green's function as first introduced by Vico-Greengard-Ferrando. This modification yields highly accurate free space solutions to the Vlasov-Poisson system, while still maintaining energy conservation up to an error bounded by the time step size. We also explain how to extend our scheme to arbitrary Dirichlet boundary conditions via standard potential theory, which we illustrate in detail for Dirichlet boundary conditions on a circular boundary. We support our approach with proof-of-concept numerical results from two-dimensional plasma test cases to demonstrate the accuracy, efficiency, and conservation properties of the scheme. By avoiding grid heating and finite grid instability we are able to show an order of magnitude speedup compared to the standard PIC scheme for a long time integration cyclotron simulation.


Introduction
The particle-in-cell (PIC) method has become a common approach to numerically solve kinetic equations in plasma physics [4,5,6,7,8].It is physically intuitive, comparatively easy to parallelize, and robust for a wide range of physical scenarios.Traditional PIC schemes, however, are known to be unable to preserve energy [9,10], leading to a non-physical phenomenon known as "grid-heating" and numerical instability for long-time simulations [11].Earlier attempts to improve the conservation properties of PIC fall into several categories: some are based on the Lagrangian formulation in [12], such as [13,14]; others take advantages of implicit integrators, such as [15,16].Structure preserving geometric PIC schemes are also being developed based on a variational formulation [17,18] and discretization of the underlying Hamiltonian structure [19,20].Although many of these schemes have excellent long time stability and conservation properties, they fail to conserve one of either charge or momentum.They also differ significantly from the standard PIC framework and hence may not be as intuitive or easy to transition to from an application and implementation point of view.
Among recent developments for structure-preserving PIC, the particle-in-Fourier (PIF) scheme developed in [1,2] retains the intuitive nature of standard PIC schemes as well as their relative simplicity, and at the same time has excellent conservation properties.By exploiting the non-uniform fast Fourier transform [21,22,23], PIF directly evaluates the Fourier modes of the field efficiently, without any spatial charge spreading and interpolation processes.Both charge and momentum are conserved in this scheme, and energy is conserved in the continuous time limit.The finite-grid instability and mode-coupling effects are also eliminated [24].In [1] the authors show that in finite dimensions PIF is the only scheme which can conserve both energy and momentum simultaneously before time discretization.Even though the exact energy conservation is lost after time discretization with explicit time integrators, the scheme still possesses excellent stability and convergence properties for long time integration simulations as shown in Section 4.2.1.PIF scheme combined with implicit time integration has not been explored yet and could be an interesting research direction to explore in the future.The use of the non-uniform fast Fourier transform reduces the complexity of the algorithm from , where N p is the total number of particles and N m is the number of Fourier modes used in the simulation.Such a Fourier approach has also proved to be successful in other numerical contexts, such as in its application for the immersed boundary method [25].
On the other hand, the choice of a Fourier basis is naturally incompatible with non-periodic boundary conditions, such as free space and Dirichlet boundary conditions.Although particle methods are frequently applied to periodic problems [26,27], most problems of physics and engineering interest cannot be approximated accurately with periodic boundary conditions [26,28,29].One of the strengths of particle methods have been their adaptability to complex geometries and boundary conditions [30,31,28].Limiting PIF to periodic problems would severely limit its range of applicability and its desirability.To handle nonperiodic boundary conditions, one could consider other basis functions, such as Legendre polynomials.Unfortunately, no basis other than the Fourier basis has been found that preserves the spatial translational invariance of the Lagrangian describing the kinetic plasma of interest [1,32].
In this paper, we propose a modified version of PIF scheme enabling the method to adequately handle free space boundary conditions and Dirichlet boundary conditions.We replace the periodic Fourier Laplacian operator 1/k 2 used in PIF with a spatially truncated free space Green's function ĝL (k) proposed in [3], with its truncation radius L determined based on our computational domain and the shape function that we choose.It has been shown in [3,33] that the Fourier kernel ĝL gives us a spectrally accurate solution to the free space Poisson problem for any smooth source ρ(x) on a Cartesian grid.We show that this kernel yields a spectrally accurate solution even if the source comes from Lagrangian particles, which are not aligned with grid points.Therefore, by applying ĝL to our PIF scheme, we are able to obtain free space solutions to the Vlasov-Poisson system with high accuracy and exact charge, momentum and semi-discrete energy conservation.Additionally, by combining the free space PIF algorithm with a Laplace solver based on boundary integral method, we can obtain the solution to any Dirichlet problem.
The remainder of this paper is structured as follows.In Section 2, we revisit the numerical details of the particle-in-cell and particle-in-Fourier algorithms; then we give a brief overview of the free space Poisson solver proposed by Vico, Greengard and Ferrando [3], which plays a key role in our new numerical scheme.In Section 3, we present our strategy of combining the free space Poisson solver and the Laplace solver with the PIF algorithm, followed by an error analysis for the conservation of energy and momentum.In Section 4, we report our numerical results using two-dimensional test cases.In Section 5, we summarize our work and suggest directions for future work.

Prerequisites
Among its various applications, PIC is very often used to solve the kinetic equations that arise in plasma physics [34,9].In this work, we will focus on the kinetic equation describing the evolution of a single species, collisionless, electrostatic plasma immersed in a constant, uniform magnetic field B 0 .The evolution of the distribution function f (x, v, t) of that single-species plasma is described by the following modified Vlasov-Poisson system (here we set the vacuum permittivity ϵ 0 to 1 for simplicity): We note that in equation 2, a background charge density ρ 0 is often subtracted in the context of periodic boundary conditions.In contrast, in this manuscript, we are interested in situations in which the charge density ρ is the physical charge density, which may be the charge density of an isolated non-neutral beam.The only requirement on ρ in this work is that it must have a compact support.

Particle-in-cell method
In order to solve the system of equations ( 1)-( 2), PIC represents the distribution function as particles in a continuous phase space, each particle having a unique position X j and velocity V j .The motion of the particles is described by Newton's law: where E j indicates the self-consistent electric field E evaluated at the particle position X j .
In the PIC method, E is first computed on a grid, before it is evaluated at the particle positions via interpolation in order to advance the particles via Newton's law.Specifically, the charge density carried by the particles is spread onto grid points according to where S is the particle shape function.S is usually chosen as a b-spline function of a certain order, which ensures exact charge conservation and retains positivity of the interpolation at any spatial location.Once the charge density ρ is known at the grid points, one can compute E and φ by solving Poisson's equation on that grid using standard solvers.For example, for periodic boundary conditions and a uniform Cartesian grid, the electric field may be efficiently computed with Fast Fourier Transforms (FFTs): Here, we use the notation k to denote the summation taken over a finite number of modes k.Specifically, if we let k for an arbitrary function f (k), where N m is the number of modes in each dimension.This convention is followed throughout the paper.Once the electric field is known at all grid points, it is evaluated at the particle positions by interpolation with the same shape function S:

Particle-in-Fourier scheme
In order to improve energy conservation in the PIC method, a new scheme called particle-in-Fourier (PIF) has been proposed [1,2].In contrast to the standard PIC approach, the PIF scheme avoids grid heating effects by directly depositing the charge in Fourier space, as well as interpolating the electric field from the Fourier space to the particle locations.Let us assume our problem domain to be Ω = d − 1 2 , 1 2 , where d is the number of dimensions of the problem of interest, and denote N m 1 as the number of Fourier modes in each 1 For simplicity let us assume Nm to be even.
dimension.If we perform a Fourier transform to (4), we get where k is of the form 2π Here, the Fourier coefficients of the shape function Ŝ(k) can be evaluated analytically.Note that ρ(k) in equation (11) is different from the one in (5).This is because the b-spline functions are not band-limited in Fourier space, and in (11) we have to truncate its analytical Fourier spectrum down to N m modes.We then solve for the electric field using equation (7).Once we obtain Ê(k) and φ(k) using the standard spectral method (equations ( 6) and ( 7)), we evaluate the acceleration as Note that both (11) and ( 12) require discrete Fourier transforms of data between equi-spaced points in Fourier space and non-uniformly distributed particle positions.This can be evaluated efficiently using the so-called non-uniform fast Fourier transform (NUFFT) [22,21].Type 1 NUFFT evaluates the summation and its dual, Type 2 NUFFT, evaluates Both Type 1 and 2 NUFFT have a complexity of , which is much faster than a naive direct sum with complexity O(N p N d m ).The NUFFT thus makes the PIF scheme affordable to run long-time particle simulations with high numbers of particles and Fourier modes.

Vico-Greengard-Ferrando free space Poisson solver
Consider Poisson's equation in R d with open (free space) boundary conditions, and where ρ is a source distribution with a compact support Ω. Without loss of generality, let us suppose The solution to Poisson's equation can then be written as where g(r) = 1 4πr for d = 3 and g(r) = −1 2π log(r) for d = 2.The slow decay of the Green's functions g requires the application of fast algorithms for the efficient computation of the convolution in (16).For uniform Cartesian grids, this can be done with the convolution theorem: the Fourier transform of the convolution of the two functions in configuration space is the point-wise product of their Fourier transforms.Equation ( 16) can therefore be evaluated efficiently by computing the Fourier transforms of g and ρ via FFTs, multiplying their Fourier transforms, and computing the inverse Fourier transform of the result, again with an FFT.This is what the Hockney-Eastwood scheme [34,35], commonly used in PIC codes to deal with free space boundary conditions, does, by zero padding the charge distribution in each spatial direction on an interval of the same length as the original domain, and defining a periodic Green's function g on that extended domain.However, the scheme is only secondorder accurate [33], because of the second difficulty arising in (16), namely the fact that the Green's functions are singular at the origin.Hockney and Eastwood propose regularizing the Green's functions according to g(0) = 1.The lack of smoothness of g associated with this choice limits the convergence to second order [36,37].Other elementary regularization procedures have been proposed [38,39], but they do not generally increase the order of convergence because they do not improve the smoothness of the regularized Green's function [36,37].
The Vico-Greengard-Ferrando scheme also relies on a uniform grid, the convolution theorem and FFTs, but significantly improves the accuracy of the computation of the convolution in (16) as compared to Hockney-Eastwood, by addressing the singularity of the Green's functions at the origin in a different manner.
First, note that the Fourier transform of the Poisson Green's functions are in fact known analytically.However, a difficulty arises from the fact that in Fourier space, this transform diverges like 1/|k| 2 at the origin.Vico, Greengard and Ferrando address that difficulty as follows.The key observation is that since we are only interested in the values of φ and E inside Ω, one can subsitute the Green's functions g with truncated Green's functions g L without affecting the results in Ω.Specifically, let where rect(x) is the characteristic function on the unit interval: If we set L > √ d, then the solution φ in Ω is exactly The advantage of truncating the Green's functions in physical space in this manner is that it removes the singularity of their Fourier transforms.Indeed, the Fourier transforms of the truncated Green's functions can also be computed analytically, and are the following smooth functions: Because of the oscillatory nature of the truncated Green's functions ĝL , Vico, Greengard and Ferrando explain that there scheme in principle requires zeropadding to an extended domain Ω with a grid of size (4N g ) d in order to compute the convolution (16) without aliasing errors.That would make their scheme more expensive than the Hockney-Eastwood method.However, they explain that the slightly more unfavourable scaling of the scheme may just be applied to a precomputation step to compute the numerical mollified Green's function T on a grid of size (2N g ) d .Once computed, this numerical kernel T can be used just like the periodic g in the Hockney-Eastwood scheme, to compute φ with the same computational complexity as in the Hockney-Eastwood scheme, via the convolution theorem and FFTs, and with much higher accuracy.Recently a modification taking advantage of the symmetry has been proposed in [40], which improves the memory storage costs of the Vico-Greengard-Ferrando scheme similar to that of the Hockney-Eastwood solver.

Combining PIF with a spectral free space Poisson solver
In this section, we present a novel way to combine PIF with the Poisson solver developed by Vico et al., [3] so that we retain the advantages of both, and can solve problems with free-space boundary conditions.Specifically, we can calculate the electric forces on particles efficiently and accurately for plasmas in free space, while achieving good conservation of charge, momentum, and energy.

A gridless free space Poisson solver
Since Vico et al. solve Poisson's equation by convoluting the source with a mollified Green's function via the Fourier convolution theorem, we can view their solver as a spectral solver which is inherently compatible with PIF.Indeed, as shown in (11), in PIF we naturally evaluate ρ = F(ρ) required for the Vico-Greengard-Ferrando scheme with the NUFFT.Note that here, the charge density is again just where ( * ) denotes the convolution operation.Since we now impose free space boundary condition for the Poisson solver, it is more natural to use radially symmetric shape functions.Hence, here S(x) is a radially symmetric function, as opposed to tensor product B-spline kernels commonly used in PIC [9].This is particularly convenient for the PIF scheme because one can analytically determine the Fourier coefficients of radially symmetric shape functions.The use of radially symmetric shape functions has two main advantages, namely 1) they preserve rotational invariance, as it should be in the case of free space boundary conditions, and 2) the convolution between the mollified Green's function and the shape function is still radially symmetric, as we will see later in this section.Some sample radially symmetric shape functions that we used in this work are provided in Appendix A.
For the solution to Poisson's equation, we may therefore write: Now, in order to determine the acceleration, we need to perform one more convolution between the fields and the shape function S. We may view the results as potential and electric force fields mollified by the shape function.For the sake of conciseness, let us denote these mollified fields as The accelerations are expressed as For our scheme, we view g L * S and ∇ g L * S * S as our new convolution kernels in replacement of g L , and apply the same Vico-Greengard-Ferrando scheme to compute the acceleration: This can be performed efficiently using Type 1 & 2 NUFFT.
We observe at this point that the explicit computation of the electric and potential fields are not required at any given time step of the PIF simulation, since Eq.( 29) is sufficient to push the particles at any time step.The potential and electric force fields could therefore in principle only be computed as a postprocessing step if desired, when analyzing the results of the simulations.For example, the potential field φ can be computed using the following expression: The simplicity of the formulation described above partially conceals two minor difficulties, in that we need to be careful in choosing the size of the Fourier spectrum for the convolutions, and the radius of the mollified Green's function g L .
Let us start with the second point.Vico et al. choose L = 1.5 > √ d in two dimensions so that the convolution between ρ and g L gives the correct potential inside Ω (here we assume d = 2 but the analysis applies to any dimension).In our solver, instead of g L , the kernel function is ∇(g L * S * S): our goal here is to correctly evaluate ψ and E at any target position inside Ω.Since the radius of the shape function is predetermined and the particle positions are given at each iteration, let us view S * S * j δ(x − X j ) as the source of the field.Its compact support Ω 2R is larger than the original domain Ω by an amount 2R, as shown in figure 1.Following this observation, we choose the size of our mollified Green's function based on the following criterion: where D L (x 0 ) represents the disk of radius L centered on x 0 .In other words, when we translate the center of the mollified Green's function to any possible particle position inside Ω, seen as a target, the distance from that target to any source in Ω 2R should be less than L, so that the mollified Green's function gives the same result as the true Green's function, and ψ and E can be correctly evaluated at particle positions.It is easy to see that this is satisfied if L ≥ √ d+2R, as shown in figure 1.When this condition is set for L, we note that the radius of g L * S is larger than √ d + 3R and the radius of g L * S * S is larger than √ d + 4R.In practice, the radius of the shape function is small compared to √ d, so that the modification to the minimal radius associated with R is not significant.In our two-dimensional examples below, we chose L = 1.5 as for the standard Poisson solver in two dimensions by Vico, Greengard and Ferrando.
Figure 1: Our simulation domain based on the shape function and mollified Green's function chosen.Particles are within the domain Ω (green), and S * S * j δ(x − X j ) has a support Ω 2R (blue).The radius of the Green's function is chosen such that g L is able to cover the entire blue region Ω 2R whenever the center of g L is located inside the green region Ω.To determine the smallest possible radius, we place the center of g L at the upper right corner of Ω.After a simple calculation, we find the smallest radius to be √ 2 + 2R for g L .Finally, the extended domain is Ω (gray).Now, let us address the question of the required Fourier resolution.As discussed in Section 2.3, Vico et al. explain that a zero-padding region extending to Ω = d [−2, 2], which is four times larger than the original domain Ω, is necessary so that the aperiodic convolutions do not lead to aliasing, and so that the oscillatory Fourier transforms of the convolution kernels are properly sampled.Our scheme has the same sampling requirements, unless we implement a precomputation step as suggested in [3] and discussed in the next section.To accurately compute our convolutions in the absence of the precomputation step, we thus upscale the resolution of the Fourier grid by a factor of 4.Although all our operations are done in Fourier space when we do not implement the precomputation step, in the remainder of this article we will sometimes refer to this upsampling in Fourier space as zero-padding of the domain Ω to the extended domain Ω = d [−2, 2], thereby using the same language as Vico et al..

On precomputation
As discussed in Section 2.3, in [3] Vico et al. explain that in the common case in which the simulation domain and its grid remain fixed during the entire simulation, the extended computational grid with (4N g ) d grid points may only be used once, at the start of the simulation.This is to compute the real space representation T on the computational grid of the convolutional kernel for the solution to Poisson's equation, which can be interpreted as a mollified Green's function associated with g L .Once T is computed, the electrostatic potential is calculated for any grid point in the domain Ω with the following convolution (for a two-dimensional setting) where i, j are indices of the grid points.The advantage of precomputing T in this manner is that Eq.( 32) may then be evaluated using standard aperiodic convolution, which only requires zero-padding by a factor of 2. As long as the simulation domain and the discretization remain the same, T can be reused throughout the simulation, and only (2N g ) d grid points are needed for the calculation of the electrostatic potential.Our PIF scheme does not lend itself naturally to the same precomputation method.Indeed, in our algorithm, the locations of the targets at which the potential must be evaluated are not known in advance, and change from one time step to the next.One therefore has two choices.The first choice is to rely on a quadruply extended grid in each direction as discussed above throughout the simulation, guaranteeing the absence of aliasing errors and spectral accuracy for our Poisson solver at each time step.The second choice is to precompute effective convolution operators T 1 and T 2 associated with the two kernels g L * S and ∇g L * S * S via inverse FFTs for a quadruply refined Fourier grid in each direction, as in Vico et al., and then use T 1 and T 2 for standard aperiodic convolutions on a doubly extended grid in each direction for the entire simulation.This second approach is less expensive computationally.However, since T 1 and T 2 only lead to spectrally accurate quadrature rules for target locations matching those of the precomputation, this second approach introduces aliasing errors at the arbitrary target locations of our particles, as with standard FFT-based Poisson solvers [34].Accuracy is therefore reduced in the same way as for these solvers, with second order convergence to be expected.Both choices have merits, depending on the application, and on the number of Fourier modes and of particles for the simulation.In our numerical tests below, we will show results obtained for each of these two options.In most situations, the numerical error of particle based schemes such as ours is dominated by the noise associated with the Monte Carlo nature of the algorithm, and a second-order convergent Poisson solver provides sufficient accuracy.In such cases, relying on T 1 and T 2 and smaller padding regions does not lead to noticeable loss of accuracy, and is the most computationally efficient choice.
We note that the Hockney-Eastwood algorithm [34] also has second order convergence, and for a standard PIC scheme has the advantage that an analytical expression for the Green's function is available in the physical space.One may then think that it could be considered as a simpler alternative to the Vico-Greengard solver with precomputation.However, for our PIF scheme, this advantage is irrelevant, as we have to convolve the Green's function with the shape function in Fourier space.In order to perform the convolution correctly, the usual extension of the Green's function to a 2L domain as in the standard Hockney-Eastwood solver will not be sufficient and the Green's function has to be zero padded adequately depending on the radius of the shape function before the convolution.This partially defeats the simplicity of the Hockney-Eastwood algorithm.On the other hand, the Vico-Greengard solver is a natural option for PIF schemes as the mollified Green's function is directly given as an analytical expression in Fourier space.

Dirichlet boundary conditions
Following standard potential theory [41], our free space particle-in-Fourier (FSPIF) scheme can be combined with Laplace solvers in order to simulate systems with Dirichlet boundary conditions.Consider the standard Poisson problem with Dirichlet boundary conditions: where ∂Ω ρ is the boundary of Ω ρ .Let φ P be the free-space solution to the Poisson equation ∆φ P = −ρ discussed previously, and let φ := φ P + φ H , where φ H is a harmonic function satisying the following Laplace problem Then φ solves the desired Poisson problem with Dirichlet boundary conditions, equation (33).
Here, we illustrate this generalization with a simple example, corresponding to two-dimensional plasmas confined within a disk on whose boundary one may apply a potential f , which may depend on the position along the disk boundary.For this particular situation, it is well-known that explicit solutions to equation ( 34) may be constructed via elementary tools of complex analysis.Let us remind the reader of this construction.Suppose that φ satisfies As explained above, we split the solution φ = φ H + φ P , where φ P is the free space solution which we numerically solve according to [3].The other component, φ H , satisfies Laplace's equation on the unit disk: Poisson's formula then gives an explicit expression for φ H [42]: Writing z = [u, v], x = [x, y],the components of the electric field at the particle position corresponding to φ H are The integrals involved in the computation of the electrostatic potential and the electric field can easily be approximated numerically with the trapezoidal quadrature rule, which has spectral accuracy for analytic functions [43].Let us assume that we put N B number of equi-spaced points z j on the boundary, we can then approximate φ H using the following equation: As long as the applied potential f is smooth enough and the charge density ρ is smooth enough for the free-space electrostatic potential φ P to be smooth, and as long as we do not evaluate the electrostatic potential too close to the boundary of the disk [44], we get high order convergence for φ H and E with a simple application of the trapezoidal rule.For situations in which these circumstances are not satisfied, specialized quadrature rules could be implemented to maintain high order accuracy [44].However, in order to determine the acceleration, we need to accurately evaluate ψ H = φ H * S at particle positions.In our PIF framework, this is not trivial, because in order to compute convolutions in Fourier space, we need φH , which is not readily available given the boundary integral form of φ H . Here, we show that this convolution through Fourier space can be avoided if some additional conditions are met.To do so, we make the following assumptions: the shape function S is spherically symmetric, and φ, which is originally defined on the unit disk, is smooth enough to be analytically extended to D 1+R (0) where R is the radius of the shape function.Note that we do not have to evaluate this analytical extension, but we do need to assume its existence for the following derivation.Recall that ψ = φ * S, so that ψ satisfies ∆ψ(x) = −(ρ * S)(x), x ∈ D 1+R (0).
In general, it is a complicated matter to determine the boundary conditions that ψ satisfy, and therefore solve Eq. ( 43).However, if we further assume that ρ has a slightly smaller compact support D 1−R (0), which is often the case in applications [45,46], then we have We can therefore use the mean value property of harmonic functions for φ and the spherical symmetry of the shape function to obtain the following equalities: This identity provides the required boundary condition for ψ(x).We thus conclude that ψ(x) satisfies the following Poisson problem: Let us define ψ P = φ P * S. Then we have with free space boundary conditions.If we call ψ H the solution to the problem then ψ = ψ P + ψ H satisfies Eq.( 49) as desired.Similarly, E = E * S can be found using equation (38) and equation ( 39) by replacing f − φ P with f − ψ P .The potential energy due to φ H can also be simply determined numerically: The simple methods of potential theory we present above do not only apply to the unit disk but to any geometry.For some two-dimensional domains, our approach can be combined with conformal mapping [42,47] for explicit or semi-explicit formulas.For other domains and for three-dimensional problems, standard numerical Laplace solvers may be used.

Algorithm Outline
We combine our gridless free-space Poisson solver with the standard PIF scheme to arrive at the following algorithm, outlined below 1.Before the start of simulation, construct convolutional kernels g L * S and ∇(g L * S * S) in Fourier space, and perform the precomputation step if desired; 2. Initialize the particle positions, velocities and charges; 3. Use Type 1 NUFFT to transform particle positions into an upscaled Fourier grid, in correspondence with the size of the extended domain Ω.Let us relabel k → k/α for the upscaled Fourier grid, then: 4. Find the corresponding Fourier modes of the free space potential φ and acceleration a: 5. Find the free space acceleration of each particle using Type 2 NUFFT: If free space boundary conditions apply to the electric field, then the true acceleration is a j = a P j .Otherwise, if a Dirichlet boundary condition is imposed for the electric field, we need to additionally determine a H j using a Laplace solver, and add a H j to get the true acceleration a j = a P j + a H j ; 6. Push the particles exactly as in standard PIC / PIF using the leapfrog algorithm if no magnetic field is applied, or using the Boris algorithm if there is a constant B 0 .
We make the following observations regarding the outline of the algorithm presented above.First, strictly speaking, the computation of φ in step 4. is not required to push the particles in step 6.The computation of φ and φ may therefore only be included at certain time steps if desired, or not at all, and the electrostatic potential computed as a post-processing step if needed to interpret simulations.Second, the leapfrog and Boris algorithms are suggested as suitable algorithms for particle pushing, but our scheme could be combined with other, possibly higher order, algorithms with suitable momentum and energy conservation properties.

Energy and momentum conservation analysis
In this section, we provide energy and momentum conservation analyses for our free space particle-in-Fourier scheme, before and after time discretization.

Energy conservation
In [2] the authors provided a proof for the energy conservation of the PIF scheme under the continuous time limit with periodic boundary conditions.Here we generalize this proof to free space boundary conditions.Without loss of generality, we let q = m = 1.We define the potential energy as U E , and the kinetic energy U K , at time t in terms of the Fourier modes of the relevant physical quantities as follows: Taking time derivative of the first quantity U E : where G is the Green's function we choose depending on the boundary condition: G = 1 k 2 for periodic boundary conditions, and G = ĝL for free space boundary conditions.Notice that since ρ is a real function, ρ(k) = −ρ(−k); Also, since the Green's function g is symmetric by construction, G(k) = G(−k).Therefore, we have the following derivation: where the last line comes from the fact that Ŝ(k) is real.Finally, we take the time derivative of the kinetic energy: One can easily check the identity from equations ( 70) and (75).Furthermore, we show that with the standard leapfrog time integrator and in the sole presence of electric forces, the energy conservation error has second order convergence for both periodic and free space boundary conditions.We begin by noticing that Eq.( 21) and Eq.( 22) take the following form in Fourier space: where the superscript n denotes the quantities at time step n.We therefore immediately find the potential energy to be: Since the leapfrog scheme uses the particle velocities at half time steps, we write: where the expression for a n s was introduced in equation ( 29).If we use the midpoint rule to define V n , then we have the following relationship: é .
(81) We can then compute particle positions at time n + 1 in terms of positions and velocities at time n: We can plug this expression into the equation (78) for the potential energy: Now, for ∆t small, Substituting this equality in equation (85) gives The velocities at the next time step, V n+1 s , are calculated as: Note that we again used the midpoint rule for the last equality.This equality gives us the total kinetic energy at time n + 1 to second order in ∆t: Adding equations ( 99) and (89) we conclude that we have second order convergence for energy conservation with respect to ∆t: This result, showing second order convergence for energy conservation independently of the number of particles and modes, is not shared by classical particle-in-cell scheme.In a standard explicit particle-in-cell scheme, the grid heating effect and finite grid instability will cause energy to increase without bound over time if the spatial grid does not resolve the Debye length.We note that for our numerical tests, we will consider physical situations that involve other fields, such as an external magnetic field B 0 or a homogeneous potential field φ H derived from Dirichlet boundary conditions.By choosing suitable numerical algorithms, we can preserve the second order convergence of energy conservation in these situations.In the presence of a magnetic field B 0 , one may choose the Boris algorithm to update the particle positions and velocities, which achieves exact energy conservation if E = 0 and at least second order otherwise [48].In the presence of a homogeneous potential field φ H to satisfy Dirichlet boundary conditions, one can use the boundary integral method mentioned in Section 3.3 to compute the full electric field.The mathematical proof of energy conservation for these situations is left as future work; the numerical convergence of energy conservation in these cases will however be confirmed in the next section.

Momentum conservation
In this section, we analyze the momentum conservation of the free space particle-in-Fourier scheme.Without loss of generality we assume the mass of a particle is m = 1.We define the total momentum of the system at time t n−1/2 as By equation 80, we are able to determine the momentum at the next time step: éé .
(102) Therefore, the difference between the momentums at two adjacent time steps are Finally, we make note that are all even functions of k, whereas k = −(−k) is an odd function of k.We therefore conclude that the summation in equation 104 equals zero, and thus the total momentum of the system is conserved over time.

Comparison between PIC and PIF schemes
The computational cost of PIF schemes scales as O((|logε|+1 , where ε is the tolerance used for the NUFFT.In contrast, the computational cost of a standard PIC scheme with cloud-in-cell shape function and an FFT-based field solver scales as O(2 d N p + N d m logN d m ).Thus the standard PIC scheme has a computational complexity similar to the coarsest NUFFT PIF scheme with the tolerance of ε = 10 −1 .In other words, for the same total number of particles and number of modes, PIF schemes are costlier than the standard PIC scheme.However, the PIF scheme has advantages in the following respects: 1.As mentioned in [2], we can use pre-computed shape functions of arbitrary order in PIF schemes without any increase in their computational cost, whereas for PIC schemes the computational cost increases with the order of the shape function.2. We can choose the tolerance ε of the NUFFT in a PIF scheme depending on the desired levels of conservation in energy and charge for the problem of interest.In practice, even with the same number of particles and grid points/modes, PIF schemes can only be slightly costlier than the PIC schemes as shown in Section 4.2.1.3. The mesh size in explicit PIC schemes have to be on the order of the Debye length in order to avoid grid heating.This typically requires a much higher number of grid points than the number of modes considered in PIF schemes.Moreover, in order to have the same number of particles per cell, the total number of particles have to be increased as well.As a result of these practical constraints, one ends up having a much higher number of particles and grid points in an explicit PIC scheme compared to the PIF scheme, which in turn leads to increased time to solution.As an example we shown in Section 4.2.1 that for comparable accuracies in the quantities of interest, the PIF scheme leads to more than an order of magnitude speedup compared to the explicit PIC scheme we implemented for comparison.4. PIF schemes are spectral particle methods, which belong to the general class of geometric, structure preserving PIC schemes [1].There are other geometric PIC schemes [17,18,19,20] which also offer excellent long time conservation and stability properties.Compared to them, PIF schemes retain the intuitive and ease of implementation nature of standard explicit PIC schemes.A more detailed comparison between PIF schemes and other geomteric PIC schemes in terms of accuracy and performance is beyond the scope of current work and will be carried out elsewhere in the future.
In terms of parallelization strategies, usually domain decomposition is used for PIC schemes where both the particles and grid points are divided between the cores.PIF schemes are relatively more global in nature than PIC schemes.Hence, as mentioned in [2] they are more suited for particle decomposition parallelization strategy where only the particles are divided between cores and each core carries all the Fourier modes.However, this strategy has a bottleneck in both memory and computing costs when the simulation requires a lot of Fourier modes (as they are duplicated and not decomposed).In order to alleviate this issue recently a space-time parallelization strategy has been proposed in [49] which shows promise for large scale 3D-3V PIF simulations on extreme-scale computing architectures.

Numerical Tests
In this section, we first test the numerical accuracy of our gridless free space Poisson solver with the method of manufactured solutions.We then verify our free space PIF implementation in a two-dimensional beam dynamics problem which has been extensively studied analytically [50,51].Finally, the accuracy of our PIF scheme for Dirichlet boundary conditions based on FSPIF and the boundary integral method is also examined.For NUFFTs needed for these examples, we rely on fiNUFFT, an efficient library presented in [52].The source codes for all our numerical results can be found at https://github.com/Nigel-Shen/Free-Space-Particle-In-Fourier.

Free space Poisson problem
In this test case, we place 30 random particles in the unit box, and solve the free space Poisson problem using our gridless solver.We choose the shape function to be a truncated Gaussian with standard deviation σ = 1/100 and truncation radius R = 1/8 (see Appendix A), so that ρ is almost C ∞ within the domain Ω.We choose the radius of the mollified Green's function to be L = 1.75, and let the extended box, which includes the zero padded regions, be  Notice that without performing the precomputation step, the solver achieves spectral accuracy, whereas with the precomputation step, it has second order convergence.
Let us denote the particle positions as X j .Both ρ and φ can be determined analytically: where Ei is the exponential integral function defined as We vary the number of modes N m and evaluate the numerical solution φ(x) on a refined grid of resolution 128×128 using Fourier interpolation (i.e., padding zeros in the spectral space before inverse FFT), and compare it with the analytical solution.
We begin by plotting the analytical solution φ and compare it with the numerical solutions obtained when N m = 16, 24, 32 and 40 in Figure 2. We then plot the absolute error ||ε|| 2 = ||φ − φ num || 2 as a function of the number of modes N m on a log-log scale, shown in Figure 3. Spectral accuracy can be clearly observed as the L 2 norm decays faster than any polynomial order in the case without precomputation step.If we do a precomputation step before solving the system, the solver now is second order if evaluated at arbitrary positions inside the domain, for the reason explained in Section 3.2.

Infinitely long beam confined by a guiding magnetic field, in free space
In this example, we consider an infinitely long charged particle beam in the z direction, in free space, and confined by a strong magnetic field aligned with the z axis.This situation has been considered as an insightful model for understanding the space charge dynamics in high intensity cyclotrons, and in particular the process of beam axisymmetrization [53,50,51].We let our domain be as the computational box including the zero padded regions.We normalize the spatial length unit by the Debye length λ D = ε 0 T /qρ, and the temporal unit by the inverse plasma frequency ω −1 p = ε 0 m/qρ, where the particle mass m and charge q are each normalized to one.In order to confine our charged particles, we apply an external constant magnetic field B 0 = 300z.The magnitude of B 0 has been chosen in such a way that the force due to the magnetic field is approximately 50 to 100 times larger in magnitude compared to the force due to the electric field.The initial particle distribution is an elongated normal distribution in configuration space, and Maxwellian in velocity space: where We fix the number of particles to be 40,000, and use a second-order radially symmetric b-spline function for the shape function (see Appendix A).The radius of the shape function is set to be R = 1/N m where N m is the Fourier resolution.We first perform the simulation with N m = 128 in FSPIF and compare it with the result obtained using the standard PIC method.For FSPIC, the charges are spread onto the grid points using the second order tensor-product b-spline function, the standard Vico solver [3] is applied for spectrally accurate field solutions at the grid points, and the force is interpolated back to the particles using the same shape function as the one for the spreading.As shown in Figure 4, beam spiraling and axisymmetrization are observed in both our FSPIC and FSPIF numerical simulations, as expected, in agreement with the results in [50,51].We use negatively charged particles in all our simulations and hence Figure 4 shows the negative of charge density.The same is true for all the rest of the figures where we show charge density from the test cases.In order to demonstrate that the number of modes in Figure 4 is sufficient to resolve the spatial scales of this problem, we present in Figure 5 the absolute values of the Fourier modes for N m = 128, and it is clear that this Fourier resolution is sufficient to capture the part of the Fourier spectrum that has a significant impact on the evolution of the beam.We measure the total energy E over the full duration of the simulation for different choices of ∆t. Figure 6 shows the energy fluctuation over time for ∆t = 0.001, 0.0005, 0.00025.Notice that the fluctuation is smooth and has a period equal to the cyclotron period, which suggests that the error is mainly caused by the strong magnetic field.Since the Boris algorithm is energy-conserving if E = 0 and at least second order convergent in the presence of an electric field, the global energy error has a second order of convergence, i.e., ε(E) ∼ ∆t 2 , as can be seen in Figure 7.Moreover, from Figure 7 we observe that even with the precomputation step, the energy still converges quadratically, and in fact the error is smaller than in the case without precomputation by a constant factor.
We also compare FSPIF and FSPIC with the same number of modes for FSPIF as the number of grid points for FSPIC in each dimension (N m = N g = 32), the same number of particles N p = 40, 000, and the same time step size ∆t = 0.0005, and show the numerical results in Figure 8.Clearly, the use of a grid in FSPIC causes the energy to increase over time, whereas it is very stable with FSPIF.This will be more pronounced for certain long time integration simulations as we show in the next section.

Grid heating in a long time integration simulation
In this section, we consider a long time integration simulation, and compare FSPIC and FSPIF schemes in terms of the accuracy of the physical quantities of interest, as well as time to solution.We consider a setting slightly different from the one in Section 4.2, with the domain Ω = [−1, 1] × [−1, 1], external magnetic field B 0 = 40z, total charge Q = −20 and the standard deviations in equation 114 scaled corresponding to the domain as σ x = 1/15 and σ y = 1/5.We compare an FSPIF simulation with N m = 64 modes with two FSPIC simulations, one with N m = 64 grid points in each direction, and one with N m = 512 grid points in each direction.We take ∆t = 0.01, total number of particles 40, 000 and integrate until final time T = 1000 for FSPIC and FSPIF with N m = 64, whereas for FSPIC with N m = 512 grid points in each direction, we only integrate until T = 160 to reduce computing time requirements.FSPIF with precomputing is used for the simulations.
We compare the energy conservation between FSPIC and FSPIF schemes in Figure 9.While the energy error is very stable for FSPIF with N m = 64 in Figure 9(b), it is two orders of magnitude larger and increases continuously for FSPIC scheme with N m = 64 in Figure 9(a) due to aliasing and grid heating.For FSPIC with N m = 512 the slope of the error increase is much smaller in Figure  The grid heating and lack of energy conservation in FSPIC scheme are reflected physically in the charge density and phase-space plots in Figures 10  and 11.We can clearly see that the beam radius for the FSPIC scheme with N m = 64 is larger in Figure 10(a) as compared to the beam radius in Figure 10(b) for the FSPIF scheme with N m = 64.The beam radius for FSPIC with N m = 512 in Figure 10(c) is comparable to the FSPIF scheme, however, since the total number of particles is not increased, it has much more statistical noise.We observe similar results in the x−phase space plot in Figure 11(c) where the x−emittance is comparable for FSPIF with N m = 64 and FSPIC with N m = 512, whereas for FSPIC with N m = 64 the semi-minor axis is twice larger.The effects are much more pronounced at T = 1000 (this corresponds to 10 5 time steps, which is typical for many production-level cyclotron simulations), where the semi-minor axis is now approximately three times larger in Figure 11(d) and the corresponding beam radius is also larger as shown in Figure 11(a).
Even though in the earlier tests we used a strict tolerance of 10 −14 close to machine precision for the NUFFTs in the FSPIF scheme, it can be made lenient depending on the time step size in the explicit time integrator so that the energy conservation is not affected.For the chosen time step size of ∆t = 0.01 we found from numerical experiments that the NUFFT tolerance of 10 −4 produces indistinguishable results from that of 10 −14 .We note here that while the momentum conservation is unaffected by the tolerance of NUFFT, charge conservation is proportional to it and for this test a charge conservation error of O(10 −4 ) did not change the results in any significant way.Now the times to solution for the FSPIC and FSPIF schemes are compared in Table 1.The computations are performed in one node of CPU partition of JURECA DC supercomputer at Jülich Supercomputing Centre.Each node of JURECA DC has 128 cores of AMD EPYC 7742, 2.25GHz processors with 512GB RAM.The computations are neither parallelized nor optimized explicitly and the default settings of fiNUFFT as well as other libraries such as numpy and scipy are used.While the FSPIF scheme with N m = 64 and NUFFT tolerance of 10 −4 is slightly more expensive than the FSPIC scheme with N m = 64 (approximately 1.1 times), it is approximately 12 times faster than the FSPIC scheme with N m = 512 which is needed for comparable accuracies in physical quantities of interest such as beam radius and emittance.Thus even in 2D-2V, FSPIF gives significant speedup compared to the FSPIC scheme with the fine resolutions necessary to avoid grid heating.We expect much larger benefits in 3D-3V simulations due to the increased number of grid points and particles.

Expanding plasma in free space
We consider an expanding plasma in free space in this section, and verify the momentum conservation property of the PIF scheme numerically.The    We consider a larger domain and a smaller total charge than the previous section so that the expanding plasma does not leave the domain.Similar to the previous section, we take ∆t = 0.01, total number of particles 40, 000, N m = 64 and integrate until T = 1.0.The tolerance for NUFFT is taken as 10 −4 and precomputation is used for the FSPIF scheme.
In Figure 12 we show the initial and final charge densities as well as the relative error in total momentum over time.We get momentum conservation close to double precision accuracy as seen from Figure 12(c) which confirms our theoretical result in Section 3.5.2.

Infinitely long beam confined by a guiding magnetic field and a perfectly conducting cylinder
In this section, we consider the same infinitely long beam confined by a guiding magnetic field aligned with it, except that the beam is now also surrounded by a perfectly conducting cylinder.The examples we study serve to verify our implementation of the Dirichlet boundary conditions for our PIF solver, which follows the scheme discussed in Section 3.3.Physically, they are motivated by experiments by Noah Hurst and collaborators [46].
We begin by first checking the accuracy of our boundary integral solver using the method of manufactured solutions.We choose our analytic solution to be a Let the coordinate of a boundary point be described by a single variable θ where sin(θ) = x and cos(θ) = y, then the boundary data for our Laplace equation is: We vary the number of boundary points and use the trapezoidal rule to evaluate equation on a 256 × 256 equi-spaced grid that covers the entire domain D 1 (0).Due to the singular integrand, the numerical value u num does not converge in the vicinity of the boundary.This behavior is expected, and advanced quadrature schemes have been developed to address this issue [44].However, the implementation of such schemes is not required if the beam remains sufficiently far from the conductor during its evolution, and is beyond the scope of this work.We therefore only measure the error ε(x) = |u num −u|(x) within a slightly shrunk domain D 0.9 (0), and plot them in Figure 13.Overall, the solver provides spectral accuracy inside D 0.9 (0) away from the boundary, which matches the spectral convergence of the trapezoidal rule for periodic analytic functions [54,43].Now we combine our boundary integral solver with the free space PIF algorithm as described in Section 3.3, and test its numerical properties.We consider the same initial charge distribution as in Section 4.2.In order to keep our computation in the unit box, in this test we use D 0.5 (0) as our domain, instead of the unit disk as in Section 3.3.When solving for the fields and accelerations, we add the harmonic function from the boundary integral method in order to satisfy Dirichlet boundary conditions at the boundary, C 0.5 (0).We perform tests for two Dirichlet problems, described as   .We choose 128 boundary points to evaluate φ H and E H using the method described in Section 3.3.We choose the number of modes to be Nm = 128 and the number of particles to be 40, 000.The beam drifts upwards in the y direction over time, as expected.
problem f = f 2 , we expect the beam to receive an extra force along the y axis, so the non-neutral beam will drift upward in the y direction.
Our simulations, shown in Figures 14 and 15, confirm our intuition.The global relative energy error measured with the L ∞ norm is plotted as a function of the time step size ∆t in Figure 16.As in the free space case, the energy error converges in second order with respect to ∆t.

Conclusion
We present a novel Particle-In-Fourier (PIF) scheme which can be applied to simulations with free space boundary conditions while maintaining semi-discrete energy conservation.We show that, by using a mollified Green's function kernel, we can obtain highly accurate solutions to the free space Poisson equation, even for non-equispaced data.Combining this Poisson solver with the PIF scheme allows us to simulate the evolution of plasmas described by the Vlasov equation with free space boundary conditions.We prove the energy and momentum conservation of the scheme in semi-discrete and fully discrete settings and validate it with numerical results.
Our approach can be extended to general Dirichlet boundary conditions on general domains, via standard potential theory.As an illustration of this, we presented an algorithm combining the free space PIF scheme with the boundary integral method on a unit disk for general Dirichlet boundary conditions.We show similar accuracy and energy conservation properties as those of our free space PIF scheme.Numerical implementation and validation of three-dimensional examples in arbitrary geometries are left as future work.
The high accuracy of the field solver and conservation of energy observed in our numerical tests show potential and advantages for long-time three-dimensional plasma simulations based on particle codes.It is however important to remember that in most cases, the error in a particle code is dominated by the noise of particle sampling rather than numerical error due to the field solver.Therefore, in order to further increase the accuracy and performance of our scheme, noise reduction strategies are of particular interest.Noise reduction is typically achieved by running simulations with a larger number of computational particles, taking advantage of scalable massively parallel computing.The design of efficient parallelization strategies for our PIF schemes therefore represents a natural direction for future work.In a complementary way, sparse grids have been used to enhance the performance of PIC by drastically reducing particle sampling error for a given number of computational particles [55,56].The development of counterparts to sparse grid PIC schemes in the PIF setting appears to be another promising research direction.

Figure 2 :
Figure 2: Potential φ from the analytic solution given by equation (111) with particles represented as blue circles (top) and numerical error of our free space Poisson solver using different numbers of modes (bottom).We choose the radius of the Green's function to be L = 1.75, and the solver uses a (4Nm) 2 box.

Figure 3 :
Figure 3: Global absolute errors in L 2 norms versus the number of modes in each dimension Nm, calculated in an upscaled physical grid of size 512 × 512 using Fourier interpolation.Notice that without performing the precomputation step, the solver achieves spectral accuracy, whereas with the precomputation step, it has second order convergence.

Figure 4 :
Figure4: Charge density ρ at different times of an infinitely long non-neutral beam confined by a strong constant and uniform magnetic field, simulated using 40,000 particles and Nm = Ng = 128 using free space PIC and PIF methods.For the FSPIF scheme, the precomputation step is not performed.

Figure 5 :
Figure 5: Absolute values of Fourier modes when Nm = 128 for the same nonneutral beam simulation as in Figure 4.These figures indicate that most high-frequency modes are not very significant, and there is no obvious mode-coupling or aliasing effects.

Figure 6 :
Figure 6: Total energy versus time for the full duration of the simulation and for different time step sizes ∆t = 0.001, 0.0005, and 0.00025.The periodic variation of the energy has a period equal to the cyclotron period, indicating that the energy fluctuation is mostly caused by the strong magnetic field.

Figure 7 :
Figure7: Convergence of the global relative energy error (measured by the L∞ norm) for different time step sizes ∆t for our free space particle-in-Fourier simulation using the Boris algorithm.The Fourier resolution is set to Nm = 32.Second order convergence is observed whether or not we perform a precomputation step.Interestingly, the precomputation step reduces the energy error by a constant compared to the errors without precomputation.

Figure 8 :
Figure 8: Comparison of the total energy fluctuation as a function of between our PIC and PIF implementations.Our PIC scheme uses the Vico solver to compute the free space Poisson electrostatic potential on the collocated grid.The time step size is ∆t = 0.0005, and the number of modes in PIF and grid points in PIC is set to be Ng = Nm = 32.

9
(c) compared to Figure9(a) but the values are still an order of magnitude higher than for the FSPIF scheme in Figure9(b).

Figure 9 :Figure 10 :
Figure 9: Comparison of the total energy error as a function of time between FSPIC and FSPIF schemes for the simulations presented in Section 4.2.1.

Figure 11 :
Figure 11: Phase space and charge density comparisons at T = 160 and T = 1000 for FSPIC and FSPIF schemes for the simulations presented in Section 4.2.1.

Figure 12 :
Figure 12: Charge density at initial time (left), T = 1.0 (center) and relative error in total momentum conservation (right) for an expanding plasma in free space with the FSPIF scheme.

Figure 13 :
Figure 13: Convergence of the boundary integral method as a function of the number boundary points N B in the domain D 1 (0).The method is spectrally accurate inside the domain except for the regions in the vicinity of the boundary.

f 2 (
x) = y ∀x ∈ C 0.5 (0).(118)The evolution of the system under these two different Dirichlet boundary conditions over time is plotted in Figures14 and 15with Fourier resolution N m = 128, where the blue scatter dots are the 128 boundary points where we evaluate the boundary condition.For the first Dirichlet BC problem f = f 1 , we expect the solution to be similar to the free-space solution; for the second Dirichlet BC

Figure 14 :
Figure14: Charge density ρ at different times an infinitely long non-neutral beam confined by a strong constant and uniform magnetic field and a grounded cylinder.We choose 128 boundary points to evaluate φ H and E H using the method described in Section 3.3.The number of modes is Nm = 128, and the number of particles is 40,000.The solution is similar to the vacuum solution in Figure4, as expected.

Figure 15 :
Figure15: Charge density ρ at different times of an infinitely long non-neutral beam confined by a strong constant and uniform magnetic field and a cylinder with biasing potential f 2 (x) = y .We choose 128 boundary points to evaluate φ H and E H using the method described in Section 3.3.We choose the number of modes to be Nm = 128 and the number of particles to be 40, 000.The beam drifts upwards in the y direction over time, as expected.

Figure 16 :
Figure16: Convergence of the global relative energy error for time step sizes ∆t for our particle-in-Fourier simulation with Dirichlet boundary conditions f 1 (x), corresponding to a grounded, perfectly conducting cylinder.Here, the time integrator is the Boris algorithm.The Fourier resolution is set to Nm = 32.As for the free space version of our scheme, second order convergence is observed.

Table 1 :
Scheme No. of grid points/modes CPU time (s) Time to solution for FSPIC and FSPIF schemes with 40, 000 particles, ∆t = 0.01 integrated up to T = 160.FSPIF scheme with 64 2 modes and FSPIC scheme with 512 2 grid points have comparable accuracies in quantities of interest however FSPIF is approximately 12 times faster than FSPIC.