The Quantum Split-Step Fourier Algorithm for Nonlinear Optical Waveguides
A new numerical framework from Biancalana (arXiv:2606.24643) propagates a classical mean field and the Bogoliubov matrices U, V of its quantum fluctuations through the same split-step operator used in nonlinear SchrΓΆdinger simulations, then reconstructs the full reduced Gaussian state of any chosen spectral window β symplectic eigenvalues, Williamson modes, von Neumann entropy, and purity all follow from one Bogoliubov pair. Applied to soliton-driven resonant radiation, QSSF shows the selected band steadily thermalises toward a mixed Gaussian state whose purity is governed by only a handful of dominant Williamson modes.
Paper: arXiv:2606.24643
arXiv:2606.24643 β The Quantum Split-Step Fourier Algorithm for Nonlinear Optical Waveguides. Biancalana, Fabio.
Classical and quantum, propagated together
The classical split-step Fourier (SSF) method has been the workhorse of ultrafast optics for three decades: alternate a half-step under the nonlinear operator in real space with a full step under the dispersive operator in Fourier space, and the symmetric Strang splitting makes the integrator second-order accurate at near-machine-precision cost. The QSSF paper keeps that integrator intact but enlarges the state vector. Each frequency bin now carries a Bogoliubov pair of operators $(\delta\hat a_k, \delta\hat a_k^\dagger)$, evolved by a $2\times 2$ block of the matrices $U(z)$, $V(z)$ that commute under the nonlinear Hamiltonian the same way the classical mean field does.
What the extra width buys is symplecticity for free. Because the nonlinear step is symmetric under the Bogoliubov transformation (the commutator $[U,V]$ contracts against the sympletic form $\Sigma$ in exactly the way the nonlinear SchrΓΆdinger equation is symmetric under particle-number reversal), the discretised evolution preserves the canonical commutation relations to numerical precision β no artificial squeezing, no spurious heating, no drift in the vacuum eigenvalues. That is what allows QSSF to read the second moments of the field as a physical covariance matrix rather than as a numerically corrupted one.
From $U$, $V$ to a reduced Gaussian state
Once $U(z)$ and $V(z)$ are known for every propagation distance, the second-moment matrix $\sigma$ of any chosen spectral window is reconstructed by tracing the Bogoliubov-evolved operators against the vacuum initial state. The window can be one resonant-radiation sideband, the soliton core, or any linear combination of bins β QSSF evaluates the window matrix on the fly, so the diagnostics do not require a separate simulation pass.
Diagonalising the window's $\sigma$ in the basis that symplectic-rotates it to a normal form yields the Williamson modes: the canonical modes whose excitations $\nu_k$ measure, exactly, how far each mode has been driven away from the vacuum. $\nu_k=1$ on every mode is vacuum; $\nu_k>1$ means the mode is squeezed; $\nu_k<1$ is anti-squeezed. The reduced state itself factorises as a tensor product of one thermal-mode state per Williamson mode, so the full Gaussian picture is in hand once the $\nu_k$ spectrum is known.
Entropy and purity as propagation diagnostics
Both the von Neumann entropy of the reduced state and its purity $\mathrm{Tr}(\rho^2)$ are explicit functionals of the symplectic eigenvalues. So once QSSF has run, the entropy-versus-distance curve, the purity floor of a chosen band, and the rate at which a sideband thermalises all come out without an extra post-processing step. For soliton-driven resonant radiation, the diagnostic is sharp: the selected band acquires a steadily increasing $S_{\mathrm{vN}}$ with propagation distance, while the rest of the spectrum stays close to vacuum. The loss of purity of the band is direct, quantitative evidence of entanglement with the modes it scattered into.
Because QSSF is lossless by construction, the entropy growth is not numerical heating β it is the physical entanglement entropy between the selected window and its complement. That is the diagnostic the abstract leans on: a Gaussian window is mixed because the photons it counts are quantum-correlated with photons it does not count, and QSSF makes that accounting exact.
Why a few Williamson modes carry the radiation
The most surprising numerical finding of the paper is structural. The resonant-radiation band occupies many tens β sometimes hundreds β of Fourier bins, and one might expect a comparable number of squeezed Williamson modes. The opposite is true. Across the parameter sweep, the bulk of the entropy and the bulk of the impurity are concentrated in just a handful of Williamson modes; the rest of the $\nu_k$ spectrum sits at the vacuum floor within numerical precision.
QSSF therefore turns what is normally a high-dimensional correlation problem into a low-dimensional one. Engineers modelling quantum correlations in nonlinear frequency conversion, supercontinuum generation, or multimode squeezed-light formation in chip-scale waveguide platforms can read off the effective number of dominant modes from a single spectrum plot and treat the rest as vacuum. That is what the paper means by an "information-theoretic diagnostic": not just measuring entropy, but locating where it lives.
Where QSSF slots in next to existing methods
Truncated Wigner and the positive-$P$ representation are the standard tools for quantum propagation in nonlinear optics, and QSSF is not a replacement for either β it is a different observable. Wigner and $P$ give you trajectory-level observables (single-shot homodyne traces, g2 functions), at the cost of stochastic sampling. QSSF gives you the exact Gaussian observables β symplectic eigenvalues, Williamson modes, entropy, purity β at the cost of staying in the Gaussian sector. For diagnostics that are linear or quadratic in the second moments, which is most of the information-theoretic vocabulary of modern quantum optics, QSSF is the cheaper, more transparent tool.
The integration with the existing classical split-step Fourier code is the practical hook. Anyone already running SSF on a soliton or supercontinuum simulation can graft QSSF on as a post-processing layer: re-run the saved $U$, $V$ matrices through QSSF and ask how entangled the output is with the modes you did not measure. The classical infrastructure β split-step propagators, absorbing boundaries, adaptive step sizing β carries over unchanged.
Takeaway
QSSF is the right tool when the question is Gaussian β symplectic eigenvalues, Williamson modes, von Neumann entropy, purity of a spectral window β and the system is lossless. It propagates the Bogoliubov matrices $U$, $V$ alongside the classical mean field inside the same split-step Fourier loop, so the cost over a classical SSF run is one extra $2\times 2$ block per bin per step. It is exact within the Gaussian sector, preserves the canonical commutation relations by construction, and exposes a low-dimensional structure in resonant-radiation entanglement that was hidden inside high-dimensional covariance matrices. For ultrafast waveguide platforms β where the classical SSF integrator is already the production tool β QSSF is the cheapest way to turn pulse propagation into a quantum-information diagnostic.
- **Bogoliubov-preserving integrator.** The same Strang splitting that integrates the nonlinear SchrΓΆdinger equation also evolves $U$, $V$ with the same order and the same symplectic invariant.
- **Exact Gaussian diagnostics.** Symplectic eigenvalues, Williamson modes, entropy, and purity are direct functionals of the propagated covariance matrix; no sampling, no truncation.
- **Low effective dimension.** Resonant-radiation bands with many Fourier bins are dominated by a handful of Williamson modes β QSSF quantifies exactly how many.
- **Lossless Gaussian regime.** Where Wigner and $P$-rep would introduce sampling noise or non-Gaussian truncation, QSSF stays exact as long as the Gaussian-sector approximation holds.