W-SLDA Toolkit: A simulation platform for ultracold Fermi gases

TL;DR

The W-SLDA Toolkit tackles the challenge of simulating strongly interacting Fermi superfluids by implementing the Superfluid Local Density Approximation (SLDA) and its time-dependent extension (TDSLDA), enabling self-consistent treatment of normal and anomalous densities across the BCS-BEC crossover, including spin-imbalanced configurations. It provides static (HFB/BdG) and real-time solvers with a flexible energy-density functional framework (BdG, ASLDA, SLDA, SLDAE), regularized to control ultraviolet divergences, and is optimized for hybrid CPU/GPU HPC environments capable of handling fully 3D simulations with up to $10^5$ atoms. Reproducibility is central, with automated reproducibility packs and a template-based workflow that integrates seamlessly with the W-data format and visualization tools, while exploiting translational symmetries to achieve scalable performance. The platform serves as a versatile research infrastructure for ultracold atoms, superconductors, and nuclear-matter physics (via the W-BSk toolkit), bridging microscopic ab initio simulations and experimental realizations through a flexible, HPC-ready, open-source framework.

Abstract

We present the W-SLDA Toolkit, a general-purpose software package for simulating ultracold Fermi gases within the framework of density functional theory and its time-dependent extensions. The toolkit enables fully microscopic studies of interacting superfluid systems across the BCS-BEC crossover, including spin-imbalanced configurations and arbitrary external geometries. It provides both static and time-dependent solvers capable of describing a broad range of phenomena in one-, two-, and three-dimensional settings. In addition, the toolkit incorporates functionality for solving the standard Bogoliubov-de Gennes equations for fermions, extending its applicability to other physical systems such as superconductors. The code is implemented in C with GPU acceleration and is optimized for hybrid CPU/GPU execution on modern high-performance computing platforms. It ensures scalability on leadership-class supercomputers, enabling fully three-dimensional simulations with large atomic numbers, and allows for direct benchmarks of ultracold-atom experimental setups. Its modular architecture facilitates straightforward extensions, user customization, and seamless interoperability with other scientific software frameworks. Furthermore, an extensive collection of practical usage examples is provided through the integrated reproducibility packs functionality, ensuring transparency and reproducibility of computational results.

Figures (5)

• Figure 1: Range of applicability of the functionals implemented within the W-SLDA Toolkit. The horizontal axis represents the dimensionless interaction strength $1/(a_sk_{F})$, while the vertical axis denotes the spin polarization $p = (\rho_a - \rho_b)/(\rho_a + \rho_b)$. The BdG functional is valid in the weak-coupling limit, but it is frequently employed for qualitative studies across the entire BCS--BEC crossover. The ASLDA functional PhysRevLett.101.215301Bulgac2012 is designed to accurately describe spin-imbalanced Fermi gases at unitarity; in the spin-symmetric limit, it reduces to the original SLDA functional Bulgac2007. The SLDAE functional Boulet2022 extends this framework to cover the full range of coupling strengths from the BCS to the unitary regime. By construction, it reproduces the BdG limit in the weak-coupling regime and the SLDA limit at unitarity.
• Figure 2: Types of meshes used in the W-SLDA Toolkit. If the system possesses translational symmetries along two directions ($y$ and $z$), a 1D mesh is employed. If the symmetry is present only along one direction ($z$), functions are represented on a 2D mesh. In the absence of such symmetries, a full 3D mesh is used. Examples are shown in the bottom row, where solutions obtained on the lattice $128\times 32\times 32$ are presented for the 1D harmonic oscillator potential $V_{\sigma}^{\textrm{(ext)}}(x)=\frac{1}{2}\omega_x^2 x^2$, the 2D harmonic oscillator potential $V_{\sigma}^{\textrm{(ext)}}(x,y)=\frac{1}{2}\bigl(\omega_x^2 x^2+\omega_y^2 y^2\bigr)$, and for 3D harmonic potential $V_{\sigma}^{\textrm{(ext)}}(x,y,z)=\frac{1}{2}\bigl(\omega_x^2 x^2+\omega_y^2 y^2 +\omega_z^2 z^2 \bigr)$. The colormaps show the density distribution normalized to the value in the center of the box.
• Figure 3: Schematic illustration of the computation of gradients and the Laplacian of quasiparticle wave functions. The procedure involves one forward Fourier transform and four inverse Fourier transforms. This approach increases memory requirements, as intermediate results must be stored in auxiliary buffers.
• Figure 4: Template-based workflow in W-SLDA: the user copies a template directory to a new location (1), modifies the template files (2), compiles (3), and runs the executable (4).
• Figure 5: Workflow for the family of static codes. Yellow boxes indicate logical blocks of the code, with arrows denoting the order of execution. The process begins with reading input.txt. The section enclosed by the dashed rectangle corresponds to the self-consistent loop, which is repeated until convergence is reached (or the maximum number of iterations is exceeded). The computation terminates at the block labeled Done. Functions shown in purple and blue boxes are located in the files problem-definition.h and logger.h, respectively, and can be customized by the user. The MPI icon marks blocks that involve communication between processes.