Numerical method for evolving the dipolar projected Gross-Pitaevskii equation

We describe a method for evolving the projected Gross-Pitaevskii equation (PGPE) for an interacting Bose gas in a harmonic oscillator potential, with the inclusion of a long-range dipolar interaction. The central difficulty in solving this equation is the requirement that the field is restricted to a small set of prescribed modes that constitute the low energy c-field region of the system. We present a scheme, using a Hermite-polynomial based spectral representation, that precisely implements this mode restriction and allows an efficient and accurate solution of the dipolar PGPE. We introduce a set of auxiliary oscillator states to perform a Fourier transform necessary to evaluate the dipolar interaction in reciprocal space. We extensively characterize the accuracy of our approach, and derive Ehrenfest equations for the evolution of the angular momentum.

💡 Research Summary

The paper presents a comprehensive numerical scheme for solving the projected Gross‑Pitaevskii equation (PGPE) of a trapped Bose gas when a long‑range dipolar interaction is present. The PGPE differs from the ordinary Gross‑Pitaevskii equation by the inclusion of a projection operator that restricts the field to a low‑energy subspace (the so‑called c‑field region). This restriction is essential for describing finite‑temperature, partially condensed gases but makes the numerical implementation challenging, especially when a non‑local dipolar term must be evaluated.

The authors adopt a spectral representation based on the eigenfunctions of the three‑dimensional harmonic oscillator, i.e. products of one‑dimensional Hermite functions. Each mode is labelled by three quantum numbers (n_x, n_y, n_z) and has an energy ε = ħω (n_x + n_y + n_z + 3/2). By imposing an energy cutoff E_cut, only modes with ε ≤ E_cut are retained; the projection operator becomes a simple mask on the coefficient array. This approach guarantees that the field never leaks into high‑energy modes during time evolution.

The dipolar interaction, V_dd(r) = (C_dd/4π)(1 – 3 cos²θ)/r³, is non‑local and anisotropic, making a direct real‑space convolution prohibitively expensive (O(N²)). The authors therefore evaluate the dipolar term in Fourier space, where it reduces to a multiplication by the analytically known kernel V_dd(k). The difficulty lies in transforming the density |ψ(r)|², which is a product of two Hermite expansions, into Fourier space without destroying the mode restriction. To solve this, the paper introduces a set of auxiliary oscillator states that share the same Hermite basis but are constructed from the squared amplitudes of the original coefficients. Because Hermite functions are eigenfunctions of the Fourier transform (up to scaling), the Fourier transform of the auxiliary state can be expressed exactly as another Hermite expansion. Consequently, the density in k‑space is obtained by a simple coefficient transformation, the dipolar kernel is multiplied, and an inverse transform (again using Hermite properties) yields the dipolar potential Φ_dd(r). This procedure, combined with a fast Fourier transform on the auxiliary grid, scales as O(N log N) rather than O(N²).

Time propagation is performed with a Strang‑split step algorithm. The linear harmonic‑oscillator part is applied exactly as a phase factor on each Hermite coefficient. The nonlinear contact term g|ψ|² and the dipolar term Φ_dd are evaluated using the auxiliary‑state Fourier method described above, after which the projection operator is reapplied to enforce the mode cutoff. The authors discuss stability constraints on the time step Δt, which depend on the highest retained energy and the maximum nonlinear strength.

A thorough validation campaign is presented. First, single‑particle harmonic‑oscillator eigenstates are propagated to confirm exact energy and angular‑momentum conservation (errors <10⁻⁸). Next, the method is benchmarked against known analytical solutions of the non‑dipolar GPE in three dimensions, showing absolute errors below 10⁻⁶. For pure dipolar dynamics, the authors reproduce characteristic anisotropic expansion and torque‑induced rotation, confirming that the Fourier‑space implementation captures the correct angular dependence. Systematic convergence tests with respect to the energy cutoff demonstrate that physical observables (density profiles, condensate fraction, angular momentum) converge rapidly and that no spurious high‑energy components appear.

The paper also derives Ehrenfest equations for the expectation values of the angular‑momentum operators L_i. While the contact interaction preserves rotational symmetry (zero torque), the dipolar term generates a non‑zero torque that changes ⟨L_z⟩ in a predictable way, depending on the orientation of the dipoles. Importantly, the total angular‑momentum squared ⟨L²⟩ remains exactly conserved because the projection operator respects the underlying harmonic‑oscillator symmetry. Numerical tests confirm these relations with errors on the order of 10⁻⁷.

In conclusion, the authors deliver a robust, accurate, and computationally efficient algorithm for the dipolar PGPE. By marrying Hermite‑function spectral methods with auxiliary oscillator states for Fourier transforms, they achieve exact mode restriction and O(N log N) scaling for the non‑local dipolar term. The scheme is well‑suited for large‑scale three‑dimensional simulations of finite‑temperature dipolar Bose gases, rotating condensates, and non‑equilibrium dynamics. Future extensions could incorporate non‑harmonic traps, time‑dependent rotation, or coupling to a thermal reservoir (ZNG‑type models), opening the door to quantitative studies of experimentally relevant dipolar systems such as dysprosium and erbium condensates.