Time-Splitting Fourier Pseudospectral Method

• Time-splitting Fourier pseudospectral method is a numerical technique that decomposes complex time-dependent PDEs into simpler subproblems to achieve high spectral accuracy.
• It combines operator splitting with Fourier-based spatial discretization to efficiently handle both linear and nonlinear terms while ensuring stability over long simulations.
• Widely applicable in quantum dynamics, fluid and plasma simulations, and optimal control, it leverages FFTs and advanced boundary treatments to optimize computational performance.

The time-splitting Fourier pseudospectral method is a family of numerical algorithms for the solution of time-dependent partial differential equations (PDEs), characterized by a temporal operator splitting strategy and a spatial discretization based on the Fourier pseudospectral (or spectral collocation) method. These schemes aim to exploit the high accuracy of global spectral methods while efficiently advancing the dynamics of both linear and nonlinear terms by decomposing complex evolution equations into simpler subproblems. The method has found application in a broad range of fields, including quantum mechanics, fluid dynamics, plasma simulation, and general relativity, where it is valued for its spectral or near-spectral accuracy, efficiency, and suitability for advanced parallel architectures.

1. Operator Splitting Formulation and Scope

At the core of time-splitting Fourier pseudospectral schemes is the splitting of a complex evolution equation into a sequence of subproblems, each associated with a distinct operator. This decomposition is typically based on physical or mathematical structure—for example, separating kinetic, convective, potential, or diffusive terms.

Given a general time-dependent PDE,

$\partial_t u = (\mathcal{A} + \mathcal{B})u,$

the method approximates the time propagator for one time step $\tau$ using, for example, Lie-Trotter (first-order) or Strang splitting (second-order) schemes:

• Lie-Trotter: $$u^{n+1} \approx \exp(\tau \mathcal{A}) \exp(\tau \mathcal{B}) u^n$$
• Strang: $$u^{n+1} \approx \exp(\tfrac{\tau}{2} \mathcal{A}) \exp(\tau \mathcal{B}) \exp(\tfrac{\tau}{2} \mathcal{A}) u^n$$

Each substep is solved using a solver appropriate to that component, often exactly or with an efficient explicit scheme. This approach enables the exact or high-fidelity treatment of stiff operators, such as dispersive terms in the Schrödinger equation or convection/diffusion in kinetic models, and supports explicit handling of nonlinearities or source terms.

Applicability

• Quantum dynamics: Time-dependent Schrödinger equation, Dirac equation, Gross-Pitaevskii equation (Serov et al., 2015, Bao et al., 2017, Bao et al., 2023, Lin et al., 14 Jan 2024)
• Fluid and plasma models: Incompressible Navier–Stokes, two-fluid plasma models, turbulence DNS (2002.01392, Morel et al., 2020, Ketcheson et al., 2018)
• Kinetic and open quantum systems: Wigner–Fokker–Planck, Vlasov–Poisson (Yi et al., 14 Sep 2025)

2. Fourier Pseudospectral Spatial Discretization

The spatial discretization relies on the Fourier pseudospectral method, which represents solution fields as periodic trigonometric polynomials: $$u(x, t) \approx \sum_{k=-N/2}^{N/2-1} \hat{u}_k(t) e^{i k x/L}$$ on a uniform grid of $N$ points. Derivatives are computed by spectral differentiation (multiplication by $ik$ or $-(k/L)^2$ in Fourier space), providing "infinite order" accuracy for smooth periodic problems.

• For multidimensional or tensor fields, tensor product bases or double Fourier representations are used (Bruegmann, 2011).
• For non-periodic domains, extensions such as Fourier Continuation techniques are employed to extend non-periodic data to a smooth periodic function on an enlarged domain, restoring spectral accuracy while avoiding the Gibbs phenomenon (2002.01392).
• In higher dimensions or for structural regularity, collocation over rank–1 or rank–$r$ lattices enables efficient spectral transforms with $O(n \log n)$ complexity, and supports high-dimensional problems (Suzuki et al., 2018).

Special attention is necessary for problems with non-smooth data or low-regularity potentials. The extended Fourier pseudospectral (eFP) method uses "windowing" and careful projection/interpolation to maintain optimal convergence rates even when dealing with discontinuous or sharply varying coefficients (Bao et al., 2023, Lin et al., 14 Jan 2024).

3. Splitting Strategies and High-Order Integrators

Sophisticated time-splitting strategies have been developed to accommodate equations with multiple stiff or weakly regular terms:

• Strang splitting: Widely used for second-order accuracy, e.g., in nonlinear dispersive equations (Serov et al., 2015, Einkemmer et al., 2020, Yi et al., 14 Sep 2025).
• Higher-order compact splitting: For example, the $S_{4c}$ fourth-order compact time-splitting for the Dirac equation incorporates a double commutator to achieve fourth-order temporal and spectral spatial accuracy with a compact number of operator exponentials and no negative time steps (Bao et al., 2017).
• Modified Strang Splitting with Transparent Boundary Conditions (TBCs): For dispersive or advection-dominated flows on bounded or unbounded domains, boundary compatibility is ensured by modifying the splitting so inflow conditions are identically satisfied, avoiding order reduction (Einkemmer et al., 2020, Einkemmer et al., 2019).
• Lawson-type and exponential integrators: In cases with low-regularity or time-dependent potentials, as in the Gross-Pitaevskii equation, the time-splitting is combined with an exponential (Lawson-type) integrator. The non-smooth potential is preprocessed using an extended FFT in an enlarged window, and the exponential of the potential is approximated for efficiency and robustness (Lin et al., 14 Jan 2024).

Algorithmic Structure

A typical high-order time-splitting Fourier pseudospectral step consists of:

1. Transform to Fourier space (using FFT or windowed FFT).
2. Evolve under linear/subproblem (e.g., kinetic term or convolution) – often exact in spectral space.
3. Inverse transform to real space.
4. Evolve under nonlinear/subproblem in real space (explicit, exponential, or Runge–Kutta steps).
5. Repeat transforms as needed for further substeps per the splitting scheme.

4. Filtering, Stability, and Aliasing

In spherical or tensorial geometries, high-frequency instabilities and pole clustering require specialized filtering:

• Spin-weighted spherical harmonic filters (Yn-filters) project tensor fields onto the correct harmonic basis, crucial for stability in double Fourier expansions for tensor evolution on spherical shells. Non-spin-weighted filters are insufficient to prevent exponential instabilities (Bruegmann, 2011).
• Aliasing control: In nonlinear PDEs, especially with polynomial nonlinearities, nonlinear aliasing can be analytically bounded using interpolation and anti-aliasing lemmas.
• Unconditional stability: Many TS-FP schemes are shown to be unconditionally stable (no CFL constraint) under reasonable smoothness and spectral operator choices. Explicit matrix multiplication for differentiation can outperform FFTs for moderate resolution and is especially GPU–friendly (Bruegmann, 2011).

5. Practical Implementation and Parallel Performance

Time-splitting Fourier pseudospectral methods are highly compatible with parallel hardware:

• Matrix–based approaches: For moderate problem sizes (e.g., $N \lesssim 100$), direct matrix multiplication of precomputed differentiation and filtering matrices is preferable to FFTs for speed and simplicity, especially on GPUs (Bruegmann, 2011).
• FFT–based approaches: For large-scale simulation, FFTs provide $O(N \log N)$ scaling. Efficient parallel code, including MPI, OpenMP, CUDA (cuFFT), and custom memory layouts (e.g., slab decomposition and multi–dimensional strided access), enables strong and weak scaling for realistic three-dimensional problems (2002.01392).
• Adaptive time-stepping: Embedded Runge–Kutta pairs (e.g., Bogacki–Shampine BS5(4)) with local error control enable dynamic time step adjustment and significant computational savings for turbulent flows (Ketcheson et al., 2018).
• Periodic and non-periodic boundaries: Fourier Continuation and transparent boundary conditions facilitate the application to both periodic domains (standard pseudospectral) and domains with physical walls, boundaries, or unbounded extent (2002.01392, Einkemmer et al., 2019, Einkemmer et al., 2020).

6. Applications and Numerical Performance

The method has enabled precise and efficient simulation across a spectrum of time-dependent problems:

Equation/Class	Key Features/Results	References
Dirac equation	4th-order compact splitting, norm conservation, uniform parameter accuracy	(Bao et al., 2017)
Nonlinear Klein–Gordon	Reformulation to relativistic NLS, uniform long–time error bounds	(Bao et al., 2020)
Gross–Pitaevskii	eFP/LTSeFP handles low-regularity/oscillatory potentials	(Bao et al., 2023, Lin et al., 14 Jan 2024)
Wigner–Fokker–Planck	7-stage Strang splitting, spectral convergence and steady states	(Yi et al., 14 Sep 2025)
Incompressible fluids	FC pseudo–spectral in walls, efficient explicit integration of turbulence	(2002.01392)
Two-fluid plasma	PSTD–Strang splitting, spectral Maxwell solvers, LWLFn for fluid step	(Morel et al., 2020)
Optimal control ODEs	Fourier integral pseudospectral with edge detection for nonsmooth trajectories	(Elgindy, 2023)

• Numerical experiments repeatedly confirm that the combination of spectral and time-splitting methodology yields high accuracy, strong convergence rates (often spectral in space, second or higher–order in time), and robustness with respect to problem nonlinearity, lack of regularity, and long simulation times.
• New discoveries: Simulation of long-time dynamics in open quantum systems (e.g., Wigner–Fokker–Planck equations) provides numerical evidence for the existence of steady states under general (far-from-harmonic) potentials, extending beyond the reach of current analytical theory (Yi et al., 14 Sep 2025).

7. Limitations and Current Developments

While the time-splitting Fourier pseudospectral framework is powerful, several practical and theoretical challenges shape ongoing research:

• Non-smooth data and coefficients: Standard methods lose accuracy for non-smooth potentials or solutions. Extended windowing, projection/interpolation, and careful quadrature are active areas for improving robustness (Bao et al., 2023, Lin et al., 14 Jan 2024).
• Boundary treatment: High-order accuracy at physical boundaries requires advanced techniques such as Fourier Continuation, dual Petrov–Galerkin spaces, and tailored boundary conditions (TBCs) to avoid spectral pollution and reflection artifacts (Einkemmer et al., 2019, Einkemmer et al., 2020).
• Aliasing in nonlinear evolution: For strongly nonlinear problems, controlling aliasing remains delicate. Analytical anti-aliasing bounds and numerical filtering are required.
• Operator commutator bounds: Higher-order splitting and lattice-based high-dimensional methods necessitate a priori smoothness on coefficients or the solution to ensure observed convergence rates (Suzuki et al., 2018).
• Scalability: While many components are highly parallelizable, the efficiency of matrix-based methods decreases as resolution grows, and data movement can dominate for very large grids.

A plausible implication is that future advances will further integrate adaptive, windowed, and high-order time schemes, potentially combine local and global bases, and expand robust boundary treatments to fully general geometries.

---

In summary, the time-splitting Fourier pseudospectral method provides a rigorous, efficient, and flexible framework for simulating a wide class of time-dependent PDEs, combining operator splitting with spectral spatial discretization, and is supported by extensive theoretical analysis and numerical verification across diverse application domains.