An efficient preconditioned conjugate-gradient solver for a two-component dipolar Bose-Einstein condensate

TL;DR

The paper tackles efficient computation of ground states for binary dipolar Bose-Einstein condensates described by the extended Gross-Pitaevskii energy with Lee-Huang-Yang corrections. It introduces a mixture-aware preconditioned nonlinear conjugate-gradient solver formulated on the product-of-spheres normalization manifold, featuring a manifold-preserving analytic line search and two complementary preconditioners—one real-space Hessian-inspired and one Fourier-space kinetic—applied in sequence to balance stiffness across scales. Robustness is enhanced by a Hestenes–Stiefel update with restarts and monotonicity safeguards, enabling monotone energy descent in highly nonconvex landscapes. Benchmark comparisons against imaginary-time evolution show 1–2 orders of magnitude reduction in iterations and typically lower energies, validating time-to-solution advantages and reliable reproduction of textures and phase boundaries in dipolar mixtures. The method is well-suited for large-scale parameter sweeps and phase-diagram mapping, providing a practical bridge between numerical metastable states and experimentally accessible configurations.

Abstract

We develop a preconditioned nonlinear conjugate-gradient solver for ground states of binary dipolar Bose-Einstein condensates within the extended Gross-Pitaevskii equation including Lee-Huang-Yang corrections. The optimization is carried out on the product-of-spheres normalization manifold and combines a manifold-preserving analytic line search, derived from a second-order energy expansion and validated along the exact normalized path, with complementary Fourier-space kinetic and real-space diagonal (Hessian-inspired) preconditioners. The method enforces monotonic energy descent and exhibits robust convergence across droplet, stripe, and supersolid regimes while retaining spectrally accurate discretizations and FFT-based evaluation of the dipolar term. In head-to-head benchmarks against imaginary-time evolution on matched grids and tolerances, the solver reduces iteration counts by one to two orders of magnitude and overall time-to-solution, and it typically attains slightly lower energies, indicating improved resilience to metastability. We reproduce representative textures and droplet-stability windows reported for dipolar mixtures. These results establish a reliable and efficient tool for large-scale parameter scans and phase-boundary mapping, and for quantitatively linking numerically obtained metastable branches to experimentally accessible states.

Figures (3)

• Figure 1: Convergence comparison for a strongly dipolar $^{52}\mathrm{Cr}$ two-component BEC with oppositely polarized components. (a) CG iteration counts $N_{\mathrm{step}}$ with preconditioners $\mathcal{P}_V$ and $\mathcal{P}_\Delta$ over $N_1\in{2,5,8,10}\times10^5$ and $N_1/N_2\in{1,2,5,10}$ (x-axis: $N_1$$(\times10^5)$ on the first line and $N_1/N_2$ on the second; y-axis: $N_{\mathrm{step}}$). $\mathcal{P}_V$ typically converges in fewer steps than $\mathcal{P}_\Delta$. (b) Energy decay of CG with $\mathcal{P}_\Delta$ and $\mathcal{P}_V$ at $N_1=10\times10^5$ and $N_1/N_2=1$; convergence is declared when the maximum residual falls below $10^{-10}$. (c) Iteration counts of CG and ITE for $N_1\in{5,10}\times10^5$ and $N_1/N_2\in{1,2,3,4,5}$ (same axis labeling as in (a)); for fairness, convergence for both methods is defined by an energy variation below $10^{-6}$. (d) Energy decay at $N_1=10\times10^5$ and $N_1/N_2=1$, highlighting the markedly faster convergence of CG compared to ITE.
• Figure 2: Numerical results obtained using the CG method, showing self-organized patterns of a polarized dipolar BEC with $N_1=N_2=3\times10^4$$^{162}$Dy atoms in a $\hat{\bm z}$-aligned magnetic field. The trapping frequencies $\omega_z/2\pi={100,250,300,350,370,425.3}$,Hz generate droplets with $N_d={1,2,3,4,5,6}$, respectively. Parameters: $\mu_1=-\mu_2=10\mu_\text{B}$, $a_{11}=50a_0$, $a_{22}=70a_0$, $a_{12}=150a_0$. Simulations use a $128\times128\times64$ grid with periodic $xy$ boundaries. Results reproduce key features of Ref. Santos2023self.
• Figure 3: Numerical reproduction of Ref. Zhang2024metastable for Dy–Er mixtures. Parameters: $\omega = 2 \pi \times 900$ Hz ($\omega_z = 0.08$ dimensionless); dipole moments $\mu_1 = 10 \mu_B$, $\mu_2 = 7.07 \mu_B$. Balanced case with $N_1 = N_2 = 10^6$ and $\rho_1 = \rho_2 = 625$, using $a_{11} = 118.5 a_0$, $a_{22} = 58.95 a_0$, and $a_{12} = {46.49,57.65,58.5} a_0$. In another configuration ($N_1 = N_2 = 4 \times 10^6$, $a_{11} = 171.6a_0$, $a_{22} = 85.15 a_0$, $a_{12} = 16.7375 a_0$), the parameters expected to produce a ring state instead lead to a stripe structure.