High-order time splitting methods for the nonlinear Gross-Pitaevskii equation

Technical Summary: High-order Time Splitting Methods for the Nonlinear Gross-Pitaevskii Equation

Problem Statement
The paper addresses the numerical solution of the two-dimensional Gross-Pitaevskii equation (GPE) with a harmonic trapping potential. The GPE models the macroscopic wave function of a Bose-Einstein condensate (BEC). The authors focus on two primary objectives:

1. Computing the ground state of the system, which corresponds to minimizing the Hamiltonian energy under a fixed mass (L²-norm) constraint.
2. Simulating the time evolution of the system, specifically analyzing the long-time preservation of physical invariants (mass and Hamiltonian energy).

The study is restricted to the two-dimensional case ($d=2$) with an isotropic harmonic potential $V(x) = \frac{1}{2}\gamma^2(x_1^2 + x_2^2)$. This choice is motivated by the fact that the linear part of the Hamiltonian admits an explicit spectral decomposition using Hermite functions, allowing for exact computation of the linear flow within a truncated basis. The authors note that while the methods are not restricted to 2D, this setting is particularly suitable for assessing stability and error propagation in high-order splitting schemes due to the complexity of 2D matrix operations compared to the 1D case.

Methodology
The proposed methodology relies on high-order time-splitting schemes with strictly positive coefficients. This distinguishes the approach from classical symplectic methods, which often require negative time steps for orders $q > 2$. Negative steps are incompatible with irreversible evolution problems, such as the gradient flow used for ground state computation.

• Ground State Computation:
  The ground state is found by combining a normalized gradient flow with the time-splitting integrators. The gradient flow is defined as a dissipative evolution (reversing the sign of the time derivative in the GPE).

  • Splitting Strategy: The evolution operator is split into linear and nonlinear subproblems.
    • The linear subproblem (involving the Laplacian and harmonic potential) is solved exactly using a spectral method based on Hermite functions.
    • The nonlinear subproblem is solved by evolving its exact solution.
  • Normalization: Since the gradient flow does not naturally preserve the L²-norm, an explicit L²-normalization step is applied after each iteration of the splitting scheme. This enforces the mass constraint at the discrete level, effectively projecting the flow onto the L²-sphere.

• Time Evolution:
  Once the ground state is computed, the same splitting-based integrators are applied to the original GPE (Hamiltonian dynamics) to simulate time evolution. The linear and nonlinear substeps are solved similarly to the ground state case, but the linear flow involves a complex phase rotation ($e^{-i\mu t}$) rather than decay.

• Order and Complexity:
  The authors utilize symmetric splitting methods of order $q = 2, 4, \dots, 14$. The number of basic steps $S$ required for order $q$ grows quadratically ($S \propto q^2$), contrasting with the exponential growth found in some symplectic alternatives. The methods are noted to be amenable to parallelization.

Key Contributions

1. Adaptation of Positive-Step Splitting: The paper demonstrates the application of high-order splitting methods with strictly positive coefficients (introduced in [15]) to the GPE. This allows for the simultaneous handling of dissipative regimes (ground state search) and Hamiltonian regimes (time evolution) without the instability issues associated with negative time steps.
2. Ground State Algorithm: A robust algorithm is presented that integrates high-order splitting with normalized gradient descent. The inclusion of L²-normalization is identified as essential for the variational consistency of the method.
3. High-Order Benchmarking: The study provides a comprehensive numerical assessment of splitting methods up to order 14 in the presence of a harmonic potential. Previous literature (e.g., [18]) had largely focused on orders up to 6, often in free-space ($V=0$) scenarios.
4. Spectral Implementation: The efficient implementation of the linear flow using a truncated Hermite basis, leveraging the specific properties of the harmonic oscillator, is detailed.

Numerical Results
The authors present three sets of benchmarks:

• Ground State Convergence: The method converges to the ground state under different initializations. Convergence is faster when the initial seed is close to the linear ground state ($h_{0,0}$). The method achieves residuals near machine precision ($10^{-14}$ to $10^{-15}$) and demonstrates that the accuracy is more sensitive to the time step size than the order of the method.
• Invariant Preservation: Long-time simulations show that the methods preserve mass and Hamiltonian energy. While second-order methods ($q=2$) show errors of order $10^{-10}$ to $10^{-12}$, higher-order methods ($q=8$) significantly reduce these errors, reaching the round-off floor ($10^{-14}$) for moderate time steps. The numerical period of the solution matches the theoretical period derived from the chemical potential.
• Cost–Accuracy Trade-off: The study evaluates the computational cost versus accuracy for orders $q=2$ to $14$. Higher-order methods achieve prescribed accuracy with much larger time steps, making them more efficient for high-precision requirements. The observed CPU time ratios relative to order 2 align well with the theoretical quadratic growth prediction ($q(q/2 + 1)$) for lower orders, with deviations at higher orders attributed to computational overhead.

Significance and Claims
The paper claims that the proposed high-order positive time-splitting framework offers a competitive and efficient approach for solving the 2D GPE.

• Robustness: The methods are shown to be stable and accurate for both finding ground states and simulating dynamics.
• Scalability: The ability to reach high orders (up to 14) with manageable computational cost (quadratic growth in steps) suggests these methods are suitable for demanding simulations.
• Accessibility: Notably, all results were obtained on a standard desktop workstation without parallelization or specialized hardware. The authors assert that despite the lack of parallel implementation, the performance metrics (errors and CPU times) are consistent with or competitive to literature results obtained on more powerful setups.
• Limitations: The authors modestly note that the results are provided for reference and that future work could explore adaptive strategies and parallel implementations to further optimize performance. They do not claim the method solves all GPE variants but specifically validates its efficacy for the 2D harmonic case.