A spectral, quasi-cylindrical and dispersion-free Particle-In-Cell algorithm

Remi Lehe; Manuel Kirchen; Igor A. Andriyash; Brendan B. Godfrey; Jean-Luc Vay

Source-linked AI summary

A spectral, quasi-cylindrical and dispersion-free Particle-In-Cell algorithm

Remi Lehe, Manuel Kirchen, Igor A. Andriyash, Brendan B. Godfrey, Jean-Luc Vay

arXiv:1507.04790v2Jul 16, 2015physics.plasm-ph

TL;DR

The paper addresses numerical artifacts and computational cost in PIC simulations, proposing a spectral quasi-cylindrical algorithm based on Fourier and Hankel transforms. Benchmarks show accurate, dispersion-free behavior and avoidance of several artifacts present in finite-difference PIC algorithms.

Problem
Finite-difference PIC algorithms can exhibit numerical dispersion and other artifacts that may affect simulated physics, while full 3D algorithms are computationally costly for close-to-cylindrical problems.
Method
The paper derives and implements a spectral quasi-cylindrical PIC algorithm combining a Fourier transform along z with a Hankel transform in r.
Results
The benchmarks show no spurious numerical dispersion, no zero-order numerical Cherenkov effect in typical lab-frame simulations, and accurate laser force calculation for copropagating electrons.
Takeaways & Limitations
For close-to-cylindrical physical problems, the algorithm provides a lower-cost alternative to full 3D PIC while avoiding several finite-difference numerical artifacts.
Takeaways & Limitations
Spectral algorithms are generally more difficult to parallelize and to implement with specific boundary conditions or a moving window than finite-difference codes.

Abstract

from arXiv · show

We propose a spectral Particle-In-Cell (PIC) algorithm that is based on the combination of a Hankel transform and a Fourier transform. For physical problems that have close-to-cylindrical symmetry, this algorithm can be much faster than full 3D PIC algorithms. In addition, unlike standard finite-difference PIC codes, the proposed algorithm is free of numerical dispersion. This algorithm is benchmarked in several situations that are of interest for laser-plasma interactions. These benchmarks show that it avoids a number of numerical artifacts, that would otherwise affect the physics in a standard PIC algorithm - including the zero-order numerical Cherenkov effect.

1. Representation of the fields and continuous equations

The formalism represents fields through Fourier modes in Cartesian coordinates and a Fourier-Hankel representation in cylindrical coordinates. This transformation decouples spectral modes and enables time integration before transforming fields back to real space.

Cartesian Maxwell equations can be solved by representing fields as sums of Fourier modes, causing different modes to decouple.The Fourier coefficients are then integrated in time and transformed back into real space.
The Fourier-Hankel representation replaces the Cartesian Fourier representation for Maxwell equations written in cylindrical coordinates.It uses Bessel functions and spectral components associated with the cylindrical-field representation.
Cylindrical field components transform differently from the Cartesian z component because their behavior near the axis differs.Scalar fields such as charge density transform like the Cartesian z component.
Substituting the Fourier-Hankel representation into the cylindrical Maxwell equations decouples the different azimuthal modes.The resulting spectral equations have a structure similar to Cartesian spectral equations but differ in signs, factors, and imaginary terms.
The spectral quasi-cylindrical fields can be time-integrated and transformed back to real space using the same general strategy as spectral Cartesian codes.For close-to-cylindrical symmetry, only a few azimuthal modes are typically needed, reducing field manipulation relative to Cartesian 3D arrays.

2. Numerical implementation

The algorithm represents fields on intermediate and spectral grids, transforming between them with Fourier and Hankel transforms. Its discretization uses an evenly spaced radial grid, irregular transverse spectral grids, and matrix-based discrete Hankel transforms within a four-step PIC cycle.

Fields and particles are evolved using a representation that solves Maxwell equations and particle motion in a finite-size simulation box.
Current deposition and field gathering occur locally on an intermediate grid before fields and currents are transformed to spectral space.The spectral representation preserves azimuthal modes computed from particle Cartesian positions.
The transformation to spectral space combines a Fourier transform along z with a Hankel transform along r.The Fourier transform uses an FFT, while the discrete Hankel transform requires a separate discretization choice.
The radial intermediate grid is evenly spaced, whereas the k⊥ spectral grid is irregular because the boundary conditions lead to discrete Bessel modes.Different azimuthal modes use different k⊥ grids because each mode evolves separately.
The discrete Hankel transform is implemented as a matrix operation whose matrices are computed once and reused during the simulation.Its multiplication cost scales as N_r^2, slower than a 1D FFT but faster than the transverse 2D FFT typically used in spectral 3D Cartesian codes.
Charge and current smoothing damps high-frequency spectral components to mitigate increased noise near the cylindrical axis.The noise arises because the cell volume factor decreases as macroparticles approach the axis; smoothing uses a single-pass binomial-filter transfer function.

3. Benchmarks

The benchmarks test dispersion and wakefield accuracy in vacuum and plasma. PSATD reproduces analytical group velocities without spurious numerical dispersion and closely matches analytical wakefields.

Benchmark design: The benchmarks compare simulated group velocities and wakefields with analytical predictions across vacuum and plasma cases.The tests include vacuum propagation, linear plasma propagation, and a linear laser-wakefield simulation.
3.1. Propagation in vacuum: (c−vg)/c = 1.27×10−4 for the finite-waist vacuum pulse, making the physical velocity difference difficult for PIC codes to resolve.The pulse parameters are w0 = 16 µm and λ = 0.8 µm.
3.1. Propagation in vacuum: PSATD group velocity is practically resolution-independent and agrees closely with the analytical vacuum prediction, unlike PSTD and finite-difference results.PSTD gives vg > c, while the finite-difference algorithm gives vg < c in this example; their relative performance is not generalizable across timestep or cell-aspect-ratio choices.
3.2. Linear propagation in a plasma: In plasma at ne = 10−3 nc, the spectral code remains close to the analytical group velocity at all resolutions, whereas finite-difference results remain resolution-dependent and inaccurate.The spectral code shows a weak residual resolution dependence, likely from current-deposition and field-gathering errors on the finite grid.
3.2. Linear propagation in a plasma: The spectral algorithm obtains the correct group velocity independently of cell aspect ratio, avoiding the costly resolution adjustments used with finite-difference codes.Finite-difference approaches may require a finer longitudinal grid or coarser radial grid, which can be expensive or unsuitable for resolving the target physics.
3.3. Linear laser-wakefield: Analytical and simulated Ez and Ey wakefield profiles overlap very precisely, including Ez sampled near the axis where quasi-cylindrical algorithms are typically noisy.The analytical fields are computed using the simulated longitudinal laser envelope and numerical evaluation of the corresponding integrals.

4. Advantages over finite-difference algorithms

The spectral quasi-cylindrical algorithm avoids numerical artifacts observed in finite-difference simulations while retaining reduced computational requirements for close-to-cylindrical systems.

Spectral algorithms avoid numerical artifacts that can affect finite-difference PIC physics, although they are generally harder to parallelize and configure with boundaries or moving windows.These practical trade-offs accompany their higher accuracy.
The spectral algorithm is free of spurious numerical dispersion, which can otherwise modify the dephasing length of a laser-wakefield accelerator.The effect may be difficult to discern when no analytical formula exists for the dephasing length.
The dispersion-free formulation should suppress zero-order numerical Cherenkov radiation caused by relativistic particles outrunning numerically altered electromagnetic-wave velocities.Such radiation can substantially increase the emittance of an accelerated beam.
766 pC with the finite-difference algorithm and 750 pC with the spectral algorithm were injected, while only the finite-difference bunch showed characteristic high-frequency radiation.The finite-difference spectrum contained the diagnostic double-parabola, and the unphysical field was comparable to the bubble’s focusing fields.
For a copropagating laser and relativistic electron, finite-difference simulations produced resolution-dependent momentum oscillations, whereas spectral simulations produced realistic, nearly resolution-independent oscillations.The spectral result is consistent with accurate Lorentz-force calculation because its E and B fields are not staggered in time.
The algorithm combines reduced computational time and memory from quasi-cylindrical geometry with artifact avoidance from its spectral PSATD formulation.Benchmarks found no spurious numerical dispersion, no zero-order numerical Cherenkov effect in typical lab-frame simulations, and accurate laser force calculation.

Appendix A. Derivation of the spectral quasi-cylindrical representation

The derivation distinguishes Cartesian field components from cylindrical components because cylindrical components are ill-defined on the axis r = 0.

Cartesian components have regular Fourier representations, whereas cylindrical components such as E_r are ill-defined at r = 0 because the azimuthal angle is undefined there.For E = E_0e_x, E_r = E_0 cos(θ), illustrating the axis singularity.

Appendix A.1. Cartesian components

The Cartesian-component derivation expresses fields through Fourier-space representations and cylindrical-coordinate substitutions, yielding the equations used in the spectral quasi-cylindrical formalism.

A Cartesian field component F_u is represented in Fourier space before transforming the transverse variables from (k_x,k_y) to (k_⊥,φ) and (x,y) to (r,θ).Here F is typically E, B, or J, and u is x, y, or z.
The angular expansion introduces Bessel functions and azimuthal harmonics, producing the representations corresponding to Eqs. (6a) and (7a).The expansion reorganizes Cartesian Fourier components into cylindrical-mode contributions.

Appendix A.2. Cylindrical components

Cylindrical field components are derived from Cartesian Fourier-mode coefficients by combining neighboring azimuthal modes into suitable radial-component representations.

Radial fields such as F_r, including E_r, B_r, and J_r, are treated generally and derived using the cylindrical-coordinate identities.The derivation applies the same transformation framework used for the Cartesian components.
The radial representation defines combinations such as F̂_−,m = (F̂_x,m−1 + iF̂_y,m−1)/2 from neighboring azimuthal-mode coefficients.Relabeling the summation index organizes Cartesian coefficients into cylindrical components.
The same definitions and method are used to obtain the remaining cylindrical-component relations.

Appendix B. Maxwell equations for the spectral coefficients

The Maxwell equations are transformed into equations for Fourier-Hankel spectral coefficients by separating uncoupled azimuthal and longitudinal modes and handling transverse-wavenumber components through Bessel-function identities and basis properties.

Fourier-Hankel modes with different m and k_z are uncoupled and can be treated separately.
Modes with different k_⊥ may remain coupled through the Bessel functions J_m(k_⊥r) and their derivatives.
Bessel-function relations and rearrangement of the spectral equations yield equations containing one Bessel-function order at a time.
For fixed Bessel order n, the functions J_n(k_⊥r) over different k_⊥ values form a basis, enabling separation of the k_⊥ components.

Appendix C. PSATD scheme, in the Fourier-Hankel representation

The PSATD scheme uses a timestep treatment in which currents remain constant while charge density varies linearly.

Currents are held constant over one timestep, whereas charge density is modeled as linear in time.

Appendix C.1. Expressions for ˆBm

The magnetic spectral coefficients are advanced by combining Maxwell’s equations with the current relation and integrating the resulting differential equations over each timestep using Green-function solutions.

Combining the spectral Maxwell equations with the current equation produces propagation equations for the magnetic coefficients B̂_m.
The resulting second-order differential equations have constant forcing and are solved using Green functions.
During each timestep, the current spectral coefficient Ĵ_m(t) is constant and equal to its timestep value.
The field solution combines the initial field and its initial time derivative through cosine and sine factors.

Appendix C.2. Expressions for ˆEm

The electric spectral coefficients are propagated by analogous equations, with constant currents and linearly varying charge density determining the timestep forcing.

Combining the spectral Maxwell equations with the current relation yields propagation equations for the electric coefficients Ê_m.
Over one timestep, Ĵ_m(t) is constant, so its time derivatives vanish, while ρ̂_m varies linearly in time.
The electric-field equations are integrated with the same cosine–sine solution form used for the underlying differential equations.

Appendix D. Discrete Hankel Transform

The appendix constructs the discrete Hankel transform matrices from constraints based on cavity eigenmodes and special mode cases. It also imposes consistency conditions so forward and inverse transforms recover the initial function.

Transform construction: The discrete Hankel transform uses an evenly spaced real-space grid rather than a grid distributed according to Bessel-function zeros.This choice is intended to avoid inconvenience for current deposition and field gathering.
Transform construction: The transform order n and spectral-grid index m determine the transverse transform cases used in practice.The implementation uses n = m −1, n = m, or n = m + 1.
Matrix determination: Cavity eigenmodes with a perfectly conducting boundary provide N_r different constraints for determining the transform matrix M_n,m.These eigenmodes are required to remain eigenmodes of the PIC cycle.
Special cases: For n = 0 or m = 0, Eq. (D.4) supplies N_r^2 constraints, allowing M_n,m to be completely determined and numerically inverted.The inverse matrix can be extracted directly from the relations before obtaining M_n,m.
Consistency and limitation: The constructed matrices impose that a DHT followed by an IDHT retrieves the initial function exactly, while the method gives satisfying results for m = 1 but not m = 2.Further work is identified to improve the Hankel-transform representation.
Special cases: For n ≠ 0 and m ≠ 0, an empirical additional constraint sets the amplitude associated with J_n(k_0^m r) to zero, completing M_n,m.The corresponding mode has no physical meaning because J_n(k_0^m r) = 0 for every r.

Appendix D.2. Derivation of Eq. (D.3)

The derivation of Eq. (D.3) proceeds by separating cases according to the Hankel-transform order, spectral index, and Bessel-function zero. The relation is established using Bessel-function identities except for a specifically excluded simultaneous case.

Derivation strategy: After the change of variable r = r_max t, Eq. (D.3) becomes an equivalent scaled relation.The transformed relation is then proved case by case.
Case analysis: For n = m and ℓ > 0, the proof uses a Bessel identity together with J_n(α_ℓ^m) = 0.Because n = m, the relevant Bessel-function value vanishes at the specified zero.
Case analysis: For n ≠ m and ℓ > 0, the restricted relation n ∈ {m −1, m, m + 1} reduces the proof to neighboring orders, exemplified by n = m + 1.The n = m −1 case is stated to be similar.
Exceptional cases: The proofs also cover m = 0 and ℓ = 0, while the relation is excluded when n ≠ 0, m ≠ 0, and ℓ = 0 simultaneously.The exceptional case requires separate treatment because the cited relation does not apply there.
Exceptional cases: For n = 0, m ≠ 0, and ℓ = 0, the derivation handles j ≠ 0 and j = 0 separately using Bessel-function relations.This case is treated independently because α_0^m = 0.