Modified splitting methods for Gross-Pitaevskii systems modelling Bose-Einstein condensates: Time evolution and ground state computation

The year 2025 marks the 100 and 30 years anniversaries of the discovery of Bose–Einstein condensation and its successful experimental realisation. Inspired by these important research achievements, a conceptually simple approach is proposed to facilitate reliable and efficient numerical simulations. The structure of the underlying systems of coupled Gross–Pitaevskii equations suggests the use of optimised high-order operator splitting methods for dynamical evolution and ground state computation. A second-order barrier, however, prevents the applicability of standard operator splitting methods for both, time evolution as well as imaginary time propagation. An innovative alternative approach accomplishes the design of novel modified operator splitting methods that remain stable under moderate smallness assumptions on the time increments. The core idea is to incorporate commutators of the defining differential and nonlinear multiplication operators, since this permits to fulfill the basic stability requirement of positive method coefficients. Further improvements with respect to convergence at the targeted precision arise from automatic adjustments of the time stepsizes by an inexpensive local error control. The presented numerical experiments confirm the favourable performance of a specific fourth-order modified operator splitting method. Amongst others, it is demonstrated that the excellent mass and energy conservation in long-term evolutions, intrinsic attributes of geometric numerical integrators for Hamiltonian systems, is maintained for a sensible variation of the time stepsizes. Moreover, the benefits of adaptive higher-order approximations in ground state computations are illustrated.

💡 Research Summary

The paper addresses the numerical simulation of multi‑component Bose–Einstein condensates (BECs) modeled by coupled Gross–Pitaevskii (GP) equations. While standard operator splitting methods are attractive for such Hamiltonian systems because they naturally separate the linear kinetic‑potential operator from the nonlinear interaction operator, they encounter a fundamental “second‑order barrier” when applied to non‑reversible flows such as imaginary‑time propagation used for ground‑state computation. Higher‑order splitting schemes inevitably involve negative coefficients, which lead to severe instabilities for the GP system, forcing practitioners to rely on the second‑order Strang splitting as the de‑facto standard.

The authors propose a novel “modified operator splitting” framework that overcomes this barrier. The key observation is that by explicitly incorporating double commutators of the linear differential operator (L) and the nonlinear multiplication operator (N) – i.e., terms of the form (