Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by the preconditioned nonlinear conjugate gradient method

1. Introduction

• The paper formulates rotating-GPE ground-state computation as constrained energy minimization and reviews approaches based on gradient flow, nonlinear eigenvalue problems, and constrained optimization.Source: §1 · p. 2
• Preconditioned nonlinear conjugate-gradient methods are established tools for orthogonality- and normalization-constrained optimization, including in electronic structure theory.The paper notes that the GPE differs from electronic-structure models in its vortex-resolving requirements and more difficult optimization landscape.Source: §1 · p. 2
• Vortex formation requires discretizations accurate enough to resolve fine details, while the associated shallow energy minima make optimization more difficult.Source: §1 · p. 2
• The combined-preconditioner PCG approach with pseudospectral discretization outperforms the other approaches, especially for very large nonlinearities and in 2D and 3D.The numerical study also examines how the algorithm behaves as rotation speed increases.Source: §1 · p. 3

2. Definitions and notations

• The ground-state problem minimizes the energy E subject to a unit L2 norm, and its first-order condition requires the energy gradient's tangent-space projection to vanish.Source: §2 · p. 3
• The energy gradient and its second derivative are given through the mean-field Hamiltonian and the Hessian, respectively.Source: §2 · p. 3
• The spherical constraint set S contains normalized wavefunctions, and its tangent space consists of directions h satisfying Re ⟨φ, h⟩ = 0.The orthogonal projection onto that tangent space is Mφh = h − Re ⟨φ, h⟩φ.Source: §2 · p. 3
• The Lagrange multiplier associated with the spherical constraint is the chemical potential, making the minimization problem a nonlinear eigenvalue problem.The lowest-eigenvalue property described for linear or specially structured nonlinear problems does not hold here.Source: §2 · p. 3

3. Discretization

• The discretization uses a standard pseudospectral scheme based on FFTs, chosen among other approaches to represent vortex-related fine details accurately.Source: §3 · p. 4
• The wavefunction is truncated to [−L, L]^d with periodic boundary conditions and sampled on a uniform grid with even M points per dimension and spacing h = 2L/M.The description is given in 2D, with extension to other dimensions stated to be straightforward.Source: §3 · p. 4
• Discrete Fourier frequencies are ξp = pπ/L and μq = qπ/L, with indices ranging from −M/2 through M/2 − 1 in each direction.The pseudospectral approximations use Fourier coefficients in the x- and y-directions.Source: §3 · p. 4
• The discrete formulation represents the Laplacian, potential, density, and angular-momentum operator as matrices acting on the discretized wavefunction.The discrete unknown is stored as a complex vector, and a 2D FFT costs O(N log N).Source: §3 · p. 4
• The continuous computations carry over to the discrete setting, with the operators and inner product interpreted as matrices and the standard inner product on C^N.Source: §3 · p. 5

Gradient-flow discretizations

• Projected FE, BE, and CN discretizations decrease energy for sufficiently small ∆t but require a projection step because norm departure is O(∆t^2).Using the unprojected gradient instead gives O(∆t) norm departure; for FE and BE, including λ is equivalent to changing the effective ∆t.Source: §4 · pp. 5–6

Convergence analysis

• The linear iteration converges to the eigenvector associated with A’s largest-modulus eigenvalue, at rate µ_N−1/µ_N, when that eigenvalue is simple.FE is a shifted power method, BE a shift-and-invert approach, and CN uses a generalized Cayley transform.Source: §4 · pp. 6–7
• For elliptic operators, FE’s CFL-like condition forces smaller time steps as the basis grows because λ_N tends to infinity.Implicit BE and CN avoid this explicit-step restriction under the stated assumptions.Source: §4 · p. 7
• BE and CN can achieve arbitrarily small convergence rates by choosing ∆t near the stated eigenvalue-dependent thresholds, whereas FE has a bounded rate.The analysis assumes λ_1 is positive or the stated ∆t conditions hold; otherwise another eigenvalue may have a larger amplification factor.Source: §4 · p. 7

Practical considerations

• Approximating the unknown λ_1 with λ(φ^n) yields Rayleigh Quotient Iteration, which is known to converge cubically and can be viewed as Newton-like.The implicit BE and CN approaches also require globalization for robust nonlinear convergence and linear-system solves; with Fourier pseudospectral discretization, those systems are matrix-free and need iterative solvers.Source: §4 · p. 7
• The paper proposes an explicit method intended to retain BE’s robustness and preconditioning advantages while converging faster and requiring no parameter fine-tuning.Preconditioners can reduce solver iterations, but the BE approach also requires choices such as ∆t and linear-solve precision.Source: §4 · p. 7

5.1. The gradient method

• The projected preconditioned gradient method descends along the negative preconditioned residual, with the direction projected onto the tangent space of the unit-norm constraint.Without preconditioning, this method is identical to the forward Euler method; adaptive stepsize control accepts energy-decreasing steps.Source: §5.1 · pp. 8–9

5.2. The conjugate gradient method

• The conjugate-gradient method adds a Polak–Ribière update and restarts as a gradient step when its coefficient is negative, reverting to a gradient step if its direction is not descending.Numerical tests report that descent-direction checks and stepsize control are useful early; near a minimum, the direction is always descending and the prescribed stepsize decreases energy.Source: §5.2 · pp. 9–10

5.3. Stopping criteria

• Minima may form a continuum because of symmetries, making a criterion based on distance to a particular minimum too restrictive; large-rotation residual tolerances can also prolong computation without changing energy much.The paper uses an energy-based stopping criterion in its 2D and 3D cases.Source: §5.3 · p. 10

5.4. Convergence analysis

• The gradient method has convergence guarantees with a bounded-below energy and optimal stepsizes, while conjugate-gradient guarantees require suitable restarts; no asymptotically optimal convergence bound is known for the constrained, nonquadratic, preconditioned problem.The expected iteration dependence on the condition number is proportional to its square root for conjugate gradient and to the condition number for gradient descent.Source: §5.4 · pp. 10–11
• The preconditioned Hessian guides expected convergence rates, but unpreconditioned Hessian eigenvalues grow with grid refinement and confining-domain size, motivating adequate preconditioning.The largest eigenvalues grow as O(h^-2) with mesh refinement and O(L^p) with domain size for potentials growing like |x|^p.Source: §5.4 · p. 11

5.5. Preconditioners

• The kinetic preconditioner is efficient in Fourier space and reduces resolution dependence on a fixed domain, but can be inefficient for growing domains or strong nonlinear interactions; the potential preconditioner has the complementary stability pattern.The kinetic preconditioner is an inverse shifted Laplacian, while the potential preconditioner is diagonal in real space; the latter deteriorates as resolution increases.Source: §5.5 · pp. 11–12
• Combined preconditioners have small eigenvalues of size O(L^-p + h^2), so conditioning deteriorates as both domain size and spatial resolution increase.Their asymptotic condition numbers are described as the better of the kinetic and potential preconditioners.Source: §5.5 · p. 12
• Combined preconditioners perform similarly, with the symmetrized version potentially more stable in some circumstances, but may add Fourier-transform cost and lack a theoretical explanation for their observed behavior.They act asymptotically like the better of the kinetic and potential preconditioners, but may be less efficient before the asymptotic regime.Source: §5.5 · pp. 12–13

6.1. Numerical results in 1D

• The conjugate-gradient method is always more efficient than the other methods by factors varying between one and two orders of magnitude; its iteration count grows only like h^-1, owing to Krylov-like properties and optimal stepsizes.The gradient/Forward Euler method’s iteration count grows like h^-2, while Backward Euler grows like h^-1 and becomes more efficient as h decreases.Source: §6.1 · p. 14
• Across domain and mesh sizes, the combined preconditioner PC performs best, and PCGC is usually most efficient, followed by PGC and BEC.Observed near-independence from h and L is attributed to the specific initial guess, which may not excite slowly convergent modes enough to reveal the expected dependence.Source: §6.1 · p. 15, Figure 4 · p. 16
• For fixed η and h, PCGC converges much faster than PGC; both are almost independent of h, but PGC is more sensitive to η.The comparison uses L = 128 and varies η from 0 to 10^5.Source: §6.1 · p. 15, Figure 5 · p. 16

6.2. Numerical results in 2D

• A relatively large stopping tolerance for the error criterion can require long computation even when energy changes little, particularly at large ω.The experiments therefore fix the error tolerance for subsequent examples.Source: §6.2 · p. 15
• In the 2D rotating tests, PCGC outperforms PGC; for PCGC, iteration counts increase with ω at fixed η, and dependence on η strengthens as ω increases.At ω = 0, iteration counts show only slight dependence on η.Source: §6.2 · pp. 15, 17, Figure 7 · p. 17
• For the tested 2D cases, PCGC usually reaches the lowest-energy stationary state from initial data of type (d) or (d̄), and multigrid is more robust than fixed-grid initialization.The reported robustness concerns CPU time and the possibility of obtaining a stationary state with lower energy.Source: §6.2 · pp. 18–19, Table 2 · p. 19, Table 3 · p. 19

6.3. Numerical results in 3D

• For four challenging 3D cases, PCGC converges in 1,403–21,971 seconds using Thomas-Fermi initial data and multigrid.The cases span η = 100–10,000 and ω = 1.4–3.Source: §6.3 · p. 19

• The method computes stationary states of the Gross-Pitaevskii equation with fast rotation and large nonlinearities.Source: §7 · p. 19
• Numerical examples in 1D, 2D, and 3D show the method is robust and accurate, and more efficient than standard approaches.Source: §7 · p. 20
• A simple multigrid approach can further accelerate the method and improve robustness through its initial data.Source: §7 · p. 20
• The authors describe direct extension to more general GPE systems, including highly nonlinear 3D cases with very large rotations.Source: §7 · p. 20