Fourier Neural Operators and Score-Based Diffusion Models for Scalable Simulation of the Gross-Pitaevskii Equation Beyond Classical Conservation Schemes

Classical time integrators for the Gross-Pitaevskii equation (GPE) face a fundamental tension: methods that conserve both mass and energy, such as CN-FEM and the linear RE-FEM of Besse (2004), deliver reliable approximations in physically complex regimes but at substantial computational cost, while purely mass-conservative linearized schemes (LCN-FEM, Two-Step FEM) are cheaper yet develop energy blow-ups in demanding settings including Bose-Einstein condensates in optical lattices and rotating superfluids. We investigate whether data-driven neural operators and score-based generative models can serve as scalable alternatives that sidestep this tradeoff. We train Fourier Neural Operators (FNOs) as mesh-independent surrogates approximating the solution operator $u_0 \mapsto u(\cdot, t)$ of the GPE $i\partial_t u = -\Delta u + V u + \beta |u|^2 u$, without requiring per-step linear or nonlinear solves. Physics-Informed Neural Networks (PINNs) are trained in parallel, encoding the GPE residual and mass conservation $M[u] = \int_\Omega |u|^2\,dx$ as soft constraints in the loss. Both models are evaluated against all five classical schemes from Henning and Warnegard (2019) on their full benchmark suite: traveling solitons, stationary multi-solitons, BEC dynamics in a 1D optical lattice, vortex patterns in a rotating superfluid, and the Mott insulator to superfluid transition in 2D. A score-based diffusion component learns distributions over physically consistent GPE solution trajectories under varying $\beta$ and $V$, enabling uncertainty quantification over wave function evolution rather than single deterministic predictions. Across smooth test cases, FNO surrogates match the $L^1$ density errors of the linearized RE-FEM while generalizing across interaction parameters without retraining. Score-based models recover the coherent interference pattern in the Mott insulator to superfluid experiment, a regime where the nonlinear methods in Henning and Warnegard (2019) were computationally infeasible at the required spatial resolution of $16 \times 10^6$ degrees of freedom. Evaluation uses $L^1$ density error $\lVert \rho - \rho_{h,\tau} \rVert_{L^1}$ and $H^1$ wave function error $\lVert u - u_{h,\tau} \rVert_{H^1}$ throughout, consistent with the classical benchmark.