Fourier Split-Step Method

• Fourier split-step method is a family of operator-splitting schemes that decompose evolution equations into linear (Fourier space) and nonlinear (physical space) subflows.
• It employs schemes like Strang factorization to alternate half-steps between representations, achieving second-order accuracy and managing stability through spectral diagonalization.
• Methodological extensions such as filtering, sparse recovery, and GPU implementations enhance its practicality in optical fiber modeling, quantum simulations, and nonlinear PDE applications.

The Fourier split-step method is a family of operator-splitting schemes for evolution equations in which the generator is decomposed into suboperators that are simple in complementary representations, typically a local or nonlinear part in physical space and a dispersive, kinetic, or differential part in Fourier space. In its standard form for nonlinear Schrödinger-type models, one alternates exact or near-exact substeps over a small increment, using FFTs to diagonalize the linear spectral operator and pointwise multiplication to apply the nonlinear phase or potential term. This structure underlies the classical split-step Fourier method (SSFM), the Fourier split operator method, and a large number of later variants developed for sparse recovery, filtering, higher-order splitting, quantum fluctuation propagation, and other contexts (Bauke et al., 2010, Bayindir, 2015).

1. Canonical operator-splitting formulation

In the basic construction, an evolution equation is written as a sum of two operators, and the exact propagator over a small step is replaced by a composition of the corresponding subflows. A standard second-order realization is the Strang factorization

$$\hat U(t+\tau,t) = \hat U_{A_1}\!\left(t+\tau,t,\tfrac12\right)\hat U_{A_2}(t+\tau,t,1)\hat U_{A_1}\!\left(t+\tau,t,\tfrac12\right) + O(\tau^3),$$

which yields an explicit update consisting of a half-step in one representation, a full step in the complementary representation, and a final half-step in the first representation (Bauke et al., 2010).

For nonlinear wave propagation, the formulation is often written as

$$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$

with $\mathcal L$ the linear dispersive operator and $\mathcal N$ the local nonlinear operator (Biancalana, 23 Jun 2026). In the nonlinear Schrödinger equation test problem

$i\eta_t=\eta_{xx}+2|\eta|^2\eta,$

the classical split-step Fourier update is

$$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$

which the paper notes requires two FFTs per time step (Bayindir, 2015).

The same template appears in many model classes. For the cubic nonlinear Schrödinger equation on the torus, the method is written as

$$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$

with exact linear and nonlinear subflows on the Fourier collocation space (Faou et al., 2013). For optical fiber propagation, one likewise splits attenuation and dispersion from Kerr nonlinearity, and in a typical first-order step alternates nonlinear phase rotation in time domain with Fourier-domain dispersion and attenuation multiplication (Li et al., 2020).

2. Spectral realization of the linear and nonlinear substeps

The Fourier component of the method is used because differential operators become algebraic multipliers in spectral space. In the general Fourier split operator formulation, one part of the operator is diagonal in real space and the other is diagonal in Fourier space; after discretization on a lattice of $N$ points, the dominant cost per step is the FFT pair, giving $O(N\log N)$ complexity (Bauke et al., 2010).

For the time-dependent Schrödinger equation in the dipole approximation, the scalar potential is applied as a pointwise real-space phase, while the kinetic term becomes diagonal in momentum space after Fourier transform. For the Dirac equation, the interaction term is local in real space, while the free-particle part is diagonalized in momentum space after a unitary transformation of the free Dirac Hamiltonian (Bauke et al., 2010). This same spectral diagonalization principle is used in fiber optics, where the standard linear frequency-domain transfer function is

$$H_{\text{std}}(f)=\exp\!\left(-\alpha \Delta z + 2j\pi^2\beta_2 f^2 \Delta z\right), \qquad |f|\le \frac{W_s}{2},$$

for the usual SSFM discretization of the optical NLSE (Li et al., 2020).

The spectral substep need not be restricted to the linear operator alone. In the dissipative Kundu–Eckhaus equation, the nonlinear/dissipative update contains the derivative $$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$0, which is evaluated spectrally as

$$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$1

before the subsequent exact Fourier-space dispersive step is applied (Bayindir et al., 2019). In generalized nonlinear propagation equations with derivative nonlinearities, the “Three Operator Story” extends the standard two-operator architecture by introducing an augmented nonlinear operator $$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$2 consisting of a nonlinear coefficient multiplying a derivative operator; that derivative part is then handled by Fourier transformation in the relevant variable so that the exponential propagator remains tractable (Zia, 2016).

Boundary conditions also shape the spectral realization. In the split-step Padé method for guided wave propagation, periodic boundaries use the DFT, homogeneous Dirichlet boundaries use the DST, and homogeneous Neumann boundaries use the DCT; the paper emphasizes that this gives an exact treatment of the vertical operator under homogeneous boundary conditions (Walsken et al., 5 Nov 2025). This suggests that “Fourier split-step” is best understood as a broader pseudospectral paradigm rather than only the periodic FFT implementation.

3. Order of accuracy, long-time behavior, and instability mechanisms

The method admits both low-order and high-order compositions. For the nonlocal Fowler equation, Lie splitting

$$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$3

is proved first-order accurate in time, with

$$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$4

and the numerics confirm first-order convergence for Lie splitting and second-order convergence for Strang splitting (Bouharguane et al., 2011). In the Wigner(-Poisson)-Fokker-Planck setting, a seven-stage symmetric Strang composition yields second-order accuracy in time and spectral accuracy in phase space (Yi et al., 14 Sep 2025).

Long-time stability is more delicate. For the cubic nonlinear Schrödinger equation on the torus, the split-step Fourier method inherits long-time orbital stability of plane waves under a linear stability condition and a non-resonance condition; in the constant plane wave case, these conditions can be verified under the CFL restriction

$$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$5

(Faou et al., 2013). By contrast, the classical plane-wave von Neumann picture does not capture all regimes relevant in practice. On a plane-wave background for the nonlinear Schrödinger equation, the Fourier split-step threshold recalled in the finite-difference comparison is

$$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$6

but on a soliton background the numerical instability of the Fourier split-step method is highly sensitive to small changes in the step size and computational window, appears only in narrow bands near special frequencies, and strongly violates the frozen-coefficients principle (Lakoba, 2012, Lakoba, 2010).

The contrast with finite-difference split-step methods is explicit. The finite-difference split-step method can become unstable on the background of a soliton even when plane-wave von Neumann analysis would suggest stability, and the resulting unstable modes are supported by the sides of the soliton. The Fourier split-step method, by contrast, is associated with resonance of Fourier harmonics near $$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$7, and the unstable modes are essentially non-localized and involve only a few Fourier harmonics (Lakoba, 2012). A further complication appears in nonlinear Dirac equations: for the massive Gross–Neveu model, the Fourier split-step method exhibits three distinct numerical instability mechanisms, two of which are unconditional, survive $$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$8, and persist in the continuum limit (Lakoba, 2019).

A common misconception is therefore that Fourier split-step stability can be inferred directly from local constant-coefficient reasoning. The supplied analyses indicate the opposite for soliton backgrounds: discrete spectral resonance, periodicity of the numerical domain, and global coupling of Fourier modes are essential parts of the mechanism (Lakoba, 2010).

4. Methodological extensions and modified schemes

A substantial literature modifies the classical algorithm without abandoning the underlying split-step structure. One line of work exploits sparsity. The Compressive Split-Step Fourier Method (CSSFM) keeps only $$\partial_z A=(\mathcal L+\mathcal N)A, \qquad A(z+\Delta z)\approx e^{\frac{\Delta z}{2}\mathcal L}e^{\Delta z\,\mathcal N}e^{\frac{\Delta z}{2}\mathcal L}A(z),$$9 spectral components during time stepping and reconstructs the full $\mathcal L$0-point signal by compressive sensing via $\mathcal L$1-minimization, with the design rule

$\mathcal L$2

for a $\mathcal L$3-sparse signal. In the one-soliton test, SSFM with $\mathcal L$4 and CSSFM with $\mathcal L$5 gave a normalized RMS difference of $\mathcal L$6, and for $\mathcal L$7 time steps the reported timings were $\mathcal L$8 s for SSFM and $\mathcal L$9 s for CSSFM (Bayindir, 2015).

A second modification inserts low-pass filters into the linear steps. The improved SSFM replaces the all-pass linear operator by

$\mathcal N$0

so that spectral components near the Nyquist edges are suppressed before aliasing contaminates the simulation. In the reported optical-fiber experiments, the filtered method reduced NSD by a factor of about $\mathcal N$1 to $\mathcal N$2 at fixed $\mathcal N$3 km, $\mathcal N$4 dBm, and $\mathcal N$5, and the best bandwidths typically lay between $\mathcal N$6 and $\mathcal N$7 of $\mathcal N$8 (Li et al., 2020).

A third direction generalizes the operator algebra itself. The three-operator exponential Fourier method introduces $\mathcal N$9 for constant-coefficient derivative terms, $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$0 for multiplicative nonlinearities, and $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$1 for derivative nonlinearities such as self-steepening. The scheme uses a seven-part symmetrized sequence and is presented as maintaining $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$2rd order error with the additional derivative-nonlinear terms (Zia, 2016). For rotating Bose–Einstein condensates, Magnus expansions are combined with splitting so that the nonautonomous quadratic Hamiltonian is approximated by a product of exponentials computable using only Fourier transforms, effectively requiring six 1D FFTs, or about three 2D FFTs, per step (Bader, 2014).

These modifications do not replace split-step propagation; they alter the representation, the suboperators, or the reconstruction stage. A plausible implication is that the split-step architecture is less a single algorithm than a numerically organized framework in which spectral diagonalization is the stable core.

5. Principal application domains

The method is central in optical-fiber and waveguide modeling. In coherent optical communication, transmitter and receiver have been optimized end-to-end through an SSFM channel implemented in TensorFlow as a differentiable chain of operations; the study used $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$3 fixed steps per span and reported a threefold increase in achievable distance, from $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$4 to $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$5 km, relative to a standard NFT-based benchmark (Gaiarin et al., 2020). Closely related information-theoretic work studies the finite-dimensional “split-step Fourier channel,” where each segment performs nonlinearity, linearity, and noise addition; for fixed $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$6, the high-power achievable-rate lower bound saturates, while the saturation point grows asymptotically like

$i\eta_t=\eta_{xx}+2|\eta|^2\eta,$7

as the number of segments increases (Keykhosravi et al., 2015).

The same spectral splitting idea extends beyond classical mean-field propagation. The Quantum Split-Step Fourier algorithm propagates the classical field together with Bogoliubov matrices $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$8 and $i\eta_t=\eta_{xx}+2|\eta|^2\eta,$9 for Gaussian quantum fluctuations in nonlinear optical waveguides, while monitoring the commutator-preserving constraints

$$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$0

In the reported resonant-radiation example, the numerical violations $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$1 and $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$2 were $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$3, the radiation-window entropy grew from $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$4 to $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$5 nats, and the effective Williamson-mode number was $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$6 despite occupation of $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$7 FFT bins (Biancalana, 23 Jun 2026).

Outside nonlinear optics, Fourier split-step methods are used for nonlocal conservation laws, dispersive hydrodynamics, quantum kinetic equations, and Bose–Einstein condensates. The Fowler equation study combines an FFT treatment of the linear nonlocal-diffusive subproblem with a finite-difference Burgers solver, proving first-order convergence of Lie splitting and numerically confirming second-order Strang behavior (Bouharguane et al., 2011). For the Kadomtsev–Petviashvili equation, the linear part is advanced with FFTs and the nonlinear Burgers-type part with a semi-Lagrangian discontinuous Galerkin scheme; the paper reports that the approach can outperform numerical methods considered in the literature by up to a factor of five (Einkemmer et al., 2017). For the Wigner(-Poisson)-Fokker-Planck equations, a time-splitting Fourier pseudospectral method is used not only for convergence tests but also to study long-time dynamics and numerical evidence of steady states for non-harmonic external potentials (Yi et al., 14 Sep 2025).

The method is therefore not confined to cubic NLSE propagation. What remains common is the decomposition of the evolution into subflows that are exact or especially simple in spectral variables.

6. High-performance implementations and computational practice

Because the algorithm is dominated by FFTs and pointwise multiplications, it maps naturally onto parallel hardware. GPU implementations using CUFFT and simple real-space and Fourier-space kernels have been shown to accelerate both the Schrödinger and Dirac solvers substantially. For sufficiently large problems, the reported GPU FFT speedup reaches about $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$8, and the total speedup for full split-step propagation reaches about $$\eta(x,t_0+\Delta t)=F^{-1}\!\left[e^{ik^2\Delta t}F\!\left[e^{-2i|\eta(x,t_0)|^2\Delta t}\eta(x,t_0)\right]\right],$$9 relative to a single CPU core on the benchmark platform; the abstract summarizes the overall gain as more than an order of magnitude (Bauke et al., 2010).

Shared-memory parallelism has also been exploited in mixed spectral–implicit variants. For ultra-short laser pulse propagation in air, a parallel $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$0D+$$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$1 split-step Fourier method combines Fourier treatment of the temporal coordinate with a Crank–Nicholson solve in the radial coordinate. On a $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$2-core machine, the method achieved near linear speed-up with efficiency of more than $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$3; for the benchmark case with $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$4, $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$5, and $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$6, the runtime dropped from $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$7 s on one thread to $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$8 s on $$u_K^{n+1}=\Phi_{\mathrm{lin}}^h\circ \Phi_{\mathrm{nonlin}}^h(u_K^n),$$9 threads (Ma et al., 2015).

Implementation details depend on the PDE and on what is diagonal in which basis. In periodic optical-fiber models, each full SSFM step is often nothing more than nonlinear multiplication, FFT, spectral multiplication, and inverse FFT (Li et al., 2020). In guided-wave acoustics, the vertical operator is diagonalized exactly by DST, DFT, or DCT depending on the boundary conditions (Walsken et al., 5 Nov 2025). In more elaborate models, such as ultra-short pulse propagation with self-steepening and stimulated Raman scattering, the nonlinear substep itself is reformulated in Madelung variables and advanced by an upwind method under a CFL condition, while the linear half-steps remain Fourier spectral (Deiterding et al., 2015).

From this computational perspective, the Fourier split-step method is notable for a particular balance: exact or near-exact treatment of structured suboperators, explicit FFT-based data movement between representations, and compatibility with both theoretical analysis and large-scale scientific computing.