GPU Soliton Solver Framework

• Soliton_solver is a GPU-accelerated finite-difference framework that computes topological solitons through energy minimization and arrested Newton flow.
• The framework employs a theory-agnostic, modular CUDA backend to simulate diverse 2D nonlinear field theories, including chiral magnets and Abelian Higgs vortices.
• It integrates real-time, GPU-resident visualization with parallel reduction techniques, supporting both reproducible batch simulations and interactive exploration.

soliton_solver most specifically denotes an open-source, GPU-accelerated finite-difference PDE framework for the simulation and real-time visualization of topological solitons in two-dimensional nonlinear field theories. In the arXiv literature, the same expression also points to a broader class of computational and analytical devices for constructing solitary-wave solutions: finite-difference energy minimizers, Fourier pseudo-spectral time steppers, nonlinear-eigenvalue solvers, determinant-based exact-solution generators, and Hirota or inverse-scattering constructions. The package named soliton_solver is distinguished by a theory-agnostic CUDA backend, modular theory injection, and GPU-resident visualization, so that one numerical framework can be reused across chiral magnets, Abelian Higgs/Ginzburg–Landau vortices, rotating trapped Bose–Einstein condensates, liquid crystals, and related models (Leask, 25 Mar 2026).

1. Definition and conceptual scope

In the software sense, soliton_solver is designed for topological solitons, defined as localized field configurations whose stability is protected by topology: configurations belong to different homotopy classes and cannot be continuously deformed into each other without leaving the finite-energy sector (Leask, 25 Mar 2026). This places it in a different category from solvers aimed at integrable pulse solitons of KdV-, NLS-, or sine-Gordon type, although the computational objective is analogous: resolve localized nonlinear structures with controlled dynamics or controlled stationarity.

In the broader mathematical literature, a solitary wave is described as a solution whose energy travels as a localized packet and preserves this localization in time, while a soliton is a solitary wave with a strong form of stability giving particle-like behavior (Benci et al., 2011). This suggests that “soliton solver” is not a single numerical paradigm but an umbrella term spanning two distinct regimes. One regime emphasizes variationally stable finite-energy defects in nonlinear field theories, which is the domain of soliton_solver proper (Leask, 25 Mar 2026). The other emphasizes exact or asymptotically exact coherent structures in integrable or near-integrable PDEs, often accessed through inverse scattering, Hirota bilinearization, or determinant formulas (Takahashi, 2023, Wang et al., 2023).

2. Numerical core of the GPU package

The package architecture is divided into core/, theories/, visualization/, and examples/. The core/ layer provides finite-difference operators, time-stepping routines, simulation drivers, and GPU memory management. Theories are injected at runtime through a registry that stores canonical theory name, aliases, import path, and version; model-specific modules provide field definitions, parameter packing, initialization routines, CUDA kernels, observables, and optional visualization helpers. The central abstraction is a generic Simulation class, which does not hard-code a particular field theory (Leask, 25 Mar 2026).

The numerical discretization is a rectangular two-dimensional lattice with a fixed halo width and, by default, fourth-order finite-difference stencils with halo width two. The common spatial operators include

$$\partial_x \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 8\phi_{i+1,j} - 8\phi_{i-1,j} + \phi_{i-2,j}}{12\,\Delta x},$$

and

$$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$

with analogous formulas in $y$. Time stepping uses a Courant-like parameter

$C=\frac{\Delta t}{\Delta x},$

with the paper stating that generally $C<1$ for stability and giving the default value $C=0.5$ (Leask, 25 Mar 2026).

The main minimization engine is arrested Newton flow. The discretized energy $E_h[\phi]$ drives the fictitious dynamics

$\ddot{\phi}(t) = - \nabla_\phi E_h[\phi(t)].$

The method starts from rest, integrates the first-order system in $(\phi,\dot\phi)$ with an explicit method such as fourth-order Runge–Kutta, and enforces monotonicity by arresting the flow whenever

$$E_h[\phi^{n+1}] > E_h[\phi^n], \qquad \dot{\phi}^{\,n+1}\leftarrow 0.$$

For constrained theories, projection steps can be inserted during each Runge–Kutta stage, for example to enforce $$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$0 or $$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$1 in magnetization models (Leask, 25 Mar 2026).

3. Supported theories and model class

The package is not tied to one PDE. It targets two-dimensional nonlinear variational field theories with localized finite-energy configurations that are most naturally computed by energy minimization. The built-in theory set spans gauge, scalar, vector, and mixed-field models (Leask, 25 Mar 2026).

Theory	Fields	Physical system / solitons
Abelian Higgs / Ginzburg–Landau	$$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$2, $$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$3	cosmic strings, superconducting vortices
Baby Skyrme model	$$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$4	planar analogue of Skyrme solitons
Rotating trapped BEC	$$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$5	superfluid vortices
Chern–Simons–Landau–Ginzburg	$$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$6, $$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$7, $$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$8	anyon superconductors
Magnetization-extended Ginzburg–Landau	$$\partial_x^2 \phi_{i,j} \approx \frac{-\phi_{i+2,j} + 16\phi_{i+1,j} - 30\phi_{i,j} + 16\phi_{i-1,j} - \phi_{i-2,j}}{12\,(\Delta x)^2},$$9, $y$0, $y$1	ferromagnetic superconductors
Micromagnetic chiral magnet	$y$2, $y$3	skyrmions, bimerons
Multicomponent anisotropic Ginzburg–Landau ($y$4)	$y$5, $y$6	fractional vortices, skyrmions
Oseen–Frank chiral liquid crystal	$y$7, $y$8	merons, skyrmions

The common variational formulation is explicit. Representative energies include the Abelian Higgs/Ginzburg–Landau functional

$y$9

the Baby Skyrme energy

$C=\frac{\Delta t}{\Delta x},$0

and the rotating BEC energy

$C=\frac{\Delta t}{\Delta x},$1

with

$C=\frac{\Delta t}{\Delta x},$2

The chiral-magnet and liquid-crystal models additionally couple to scalar potentials through Poisson equations, so the framework covers local and coupled local/nonlocal structure within one solver design (Leask, 25 Mar 2026).

4. GPU residency, rendering, and interactive use

A defining feature of soliton_solver is that the simulation remains GPU resident as much as possible. Device buffers include field variables, fictitious velocities, energy gradients, derivative work arrays, Runge–Kutta stages, scalar energy buffers, and reduction buffers. Observables such as total energy, convergence norms, and extrema are computed by parallel GPU reductions, and only compact partial sums are transferred to the host (Leask, 25 Mar 2026).

The rendering pipeline uses CUDA–OpenGL interoperation. An OpenGL buffer is registered with CUDA, mapped into device address space for each frame, written by a CUDA visualization kernel, and then displayed by OpenGL without staging full image arrays through host memory. The listed software stack includes Python, Numba, CUDA, PyOpenGL, ModernGL, glfw, and cuda-python. The viewer can switch among energy density, order-parameter magnitude, magnetic flux density, magnetic charge density, and number density (Leask, 25 Mar 2026).

The package supports both reproducible batch simulations and interactive exploration. The documented workflow is theory-independent: load a theory from the registry, construct parameters with default_params(...), initialize Simulation(params, theory), initialize the field configuration, and launch the interactive viewer or minimizer. The examples are chosen to demonstrate this reuse. One chiral-magnet application treats a single anti-skyrmion in a Heusler compound with representative parameters $C=\frac{\Delta t}{\Delta x},$3, $C=\frac{\Delta t}{\Delta x},$4, $C=\frac{\Delta t}{\Delta x},$5, $C=\frac{\Delta t}{\Delta x},$6, and $C=\frac{\Delta t}{\Delta x},$7. A rotating trapped BEC example uses $C=\frac{\Delta t}{\Delta x},$8 atoms, $C=\frac{\Delta t}{\Delta x},$9, $C<1$0, $C<1$1, and $C<1$2, and illustrates the transition from Thomas–Fermi ground state to vortex lattice formation (Leask, 25 Mar 2026).

5. Soliton-solver methodologies in the wider literature

The wider arXiv literature shows that soliton solving is methodologically heterogeneous. Exact-solution constructions remain central for integrable or nearly integrable systems. For magnetic excitations in an $C<1$3 spinor Bose–Einstein condensate in a deep one-dimensional optical lattice, exact dark magnetic solitons and elastic two-soliton collisions are derived by reducing the effective anisotropic pseudospin chain to a defocusing NLS-type equation and then applying the Hirota bilinear method; the paper also identifies a critical field $C<1$4 for soliton existence (Li et al., 2010). For two-dimensional supersonic defocusing NLS flow, the known exact single oblique dark soliton can be written in Hirota form with $C<1$5, while the two-soliton ansatz is only approximately consistent, with near-integrability emerging for small angle difference or in the hypersonic limit $C<1$6 (Khamis et al., 2012).

Determinant and inverse-scattering solvers define another branch. The generalized discrete KdV equation admits $C<1$7-soliton Cauchy-type determinant tau functions, and one notable feature is that when $C<1$8, smaller-amplitude solitons may move faster than larger ones (Kanki et al., 2012). The sine-Gordon equation on a stationary elliptic-function background admits reflectionless multi-soliton solutions expressed by determinants of theta functions, with background shift

$C<1$9

so the solver is effectively an exact spectral-data-to-solution map on a periodic lattice background (Takahashi, 2023). For the nonlocal generalized Sasa–Satsuma equation, an improved Riemann–Hilbert method starting from the $C=0.5$0-part of the Lax pair yields compact determinant formulas for simple and high-order $C=0.5$1-soliton solutions under the reflectionless condition (Wang et al., 2023).

A different class consists of direct numerical propagators. The KdV-oriented library sangkuriang uses Fourier pseudo-spectral spatial discretization, adaptive high-order time integration via DOP853, and Numba acceleration on periodic domains, validating isolated propagation, overtaking collisions, and three-body interactions while monitoring $C=0.5$2, $C=0.5$3, and $C=0.5$4 (Herho et al., 17 Jan 2026). The simple-idealized-1d-nlse code solves the focusing 1D cubic NLSE with Fourier pseudo-spectral discretization, DOP853, a $C=0.5$5 dealiasing rule, and an exponential filter, and reproduces single solitons, two-soliton collisions, Akhmediev breathers, and modulation instability (Herho et al., 7 Sep 2025). A KdV Crank–Nicolson solver with predictor–corrector linearization was also compared against an FFT method, with reported $C=0.5$6 for Crank–Nicolson and $C=0.5$7 for FFT in a one-soliton benchmark (Bueno et al., 8 Oct 2025).

Stationary-state solvers form another recognizable category. SWtools treats nonlinear Schrödinger-type equations as nonlinear eigenvalue problems, implementing the spectral renormalization method for prescribed $C=0.5$8 and NSOM for fixed norm $C=0.5$9, with user-specified $E_h[\phi]$0 and linear operator coefficients $E_h[\phi]$1 (Melchert et al., 14 Apr 2025). At a more abstract level, the variational framework of hylomorphic solitons reduces “find a stable soliton” to constrained minimization of $E_h[\phi]$2 at fixed $E_h[\phi]$3, or equivalently to minimization of

$E_h[\phi]$4

so the solver target becomes a $E_h[\phi]$5-compact minimizer with orbital stability (Benci et al., 2011). Taken together, these works suggest that soliton_solver names not one algorithmic doctrine but a family of constructive strategies adapted to distinct notions of soliton.

6. Limitations, distinctions, and recurring misconceptions

The soliton_solver package is explicit about its present scope. It is restricted to two-dimensional theories, targets NVIDIA/CUDA through Numba, and does not yet provide benchmark tables, strong or weak scaling studies, CPU-vs-GPU speedups, convergence-rate measurements, or a systematic validation section. The paper also does not give an exhaustive discussion of boundary conditions or convergence diagnostics, and states that improved diagnostics, additional theories, and further visualization development remain future work (Leask, 25 Mar 2026). As a consequence, its current status is that of a reusable research framework with demonstrated applications rather than a fully benchmarked production code.

A recurrent conceptual distinction concerns what counts as a soliton solver. Exact integrable constructions, pseudo-spectral evolution codes, and finite-energy minimizers solve different mathematical problems. The topological-soliton package soliton_solver computes metastable and stable finite-energy configurations by arrested Newton flow; it does not solve a scattering problem or reconstruct tau functions (Leask, 25 Mar 2026). By contrast, determinant-based or Hirota-based approaches output exact $E_h[\phi]$6-soliton formulas but are usually tied to integrable reductions (Takahashi, 2023, Wang et al., 2023). This suggests that solver choice is inseparable from the ontology of the target object: topological defect, stationary nonlinear bound state, or propagating integrable pulse.

A second distinction concerns the term “soliton” itself. In the generalized fifth-order Hirota–Satsuma coupled KdV study, the reported tanh-polynomial traveling waves are explicitly described as behaving more like kink/antikink-type bounded traveling profiles on nonzero backgrounds than localized $E_h[\phi]$7-pulses decaying to zero (Forozani et al., 2016). Conversely, integrability-sensitive numerical studies increasingly use invariant monitoring, spectral diagnostics, or recurrence analysis to distinguish visually plausible pulses from genuinely soliton-like dynamics. In sangkuriang, for example, conservation errors, amplitude–velocity agreement, spectral entropy, Fisher information, and recurrence measures are all used to support the claim that computed KdV interactions preserve regular integrable phase-space structure (Herho et al., 17 Jan 2026).

The resulting landscape is technically plural. soliton_solver in the narrow sense is a GPU-centered 2D field-theory framework for topological solitons (Leask, 25 Mar 2026). In the broader arXiv sense, it names a family of exact, variational, spectral, and finite-difference methods whose common purpose is the construction, continuation, or dynamical verification of localized nonlinear coherent structures across condensed matter, nonlinear optics, hydrodynamics, and high-energy field theory.