Competitive algorithms for calculating the ground state properties of Bose-Fermi mixtures

TL;DR

The paper addresses computing the ground-state properties of a Bose-Fermi mixture described by a coupled generalized Gross-Pitaevskii equation for bosons and Hartree-Fock equations for fermions, including beyond-mean-field corrections that enable self-bound droplets. It introduces and compares multiple numerical schemes—adiabatic real-time interaction switching, imaginary time propagation with Gram-Schmidt, and iterative eigenvalue approaches with 3D or 1D bases—across four ground-state strategies. The benchmarks on Cs-133 bosons and Li-6 fermions show that the ITP-ITP-GS method yields the lowest ground-state energy with the best efficiency, while RAM limits the feasible basis size for fully 3D implementations. The study demonstrates the formation of self-bound Bose-Fermi droplets under strong attraction and highlights implications for polaron-like physics and boson-mediated fermionic phenomena in ultracold gases.

Abstract

In this work we define, analyze, and compare different numerical schemes that can be used to study the ground state properties of Bose-Fermi systems, such as mixtures of different atomic species under external forces or self-bound quantum droplets. The bosonic atoms are assumed to be condensed and are described by the generalized Gross-Pitaevskii equation. The fermionic atoms, on the other hand, are treated individually, and each atom is associated with a wave function whose evolution follows the Hartree-Fock equation. We solve such a formulated set of equations using a variety of methods, including those based on adiabatic switching of interactions and the imaginary time propagation technique combined with the Gram-Schmidt orthonormalization or the diagonalization of the Hamiltonian matrix. We show how different algorithms compete at the numerical level by studying the mixture in the range of parameters covering the formation of self-bound quantum Bose-Fermi droplets.

Figures (6)

• Figure 1: (Left) Total energy (in units of $\hbar^2/(m_B a_B^2)$) of a trapped Bose-Fermi mixture as a function of imaginary time (or real time for the inset), in units of $m_B a_B^2/\hbar$, for $g_{BF}=-1$ for all studied algorithms. The main frame shows the energy for the ITP-ITP-GS (black solid line), ITP-IEV-3D (green dots), and ITP-IEV-1D (blue circles) methods, while the inset depicts the energy in the case of the A-RTP algorithm. (Right) Simulation acceleration (which is the inverse of the ratio of the wall clock times for the current and single-core computations) as a function of the number of threads. The acceleration is better for the ITP-IEV-1D and ITP-IEV-3D methods.
• Figure 2: Radial bosonic (solid lines) and fermionic (dashed lines) densities for different methods for $g_{BF} = -1$ (left) and $g_{BF} = -3$ (right). The system consists of $N_B = 40$ bosons and $N_F = 4$ fermions and the densities are normalized to one and $N_F = 4$ for bosons and fermions, respectively. When comparing the ITP-IEV-3D method with others for $g_{BF} = -1$, no differences in the density profiles are visually visible. For higher attractive interactions, $g_{BF} = -3$, the situation changes, small differences in the density profiles are visible due to the insufficient number of eigenvectors used in the ITP-IEV-3D method.
• Figure 3: Relative density differences ($(|\psi_{j}^{(B,F)}|^2 - |\psi_{GS}^{(B,F)}|^2) / |\psi_{GS}^{(B,F)}|^2$) for different methods for $g_{BF} = -1$ (left) and $g_{BF} = -3$ (right). The "j" denotes the A-RTP, ITP-IEV-3D, and ITP-IEV-1D algorithms, while "GS" is the reference ITP-ITP-GS method. Solid and dashed lines correspond to fermionic and bosonic densities, respectively.
• Figure 4: (Left) Densities of the system for the ITP-ITP-GS algorithm for different time steps. The solid (dashed) lines denote the fermionic (bosonic) density. Minimal differences can be seen between the red and blue lines. A time step smaller than $\Delta t = 0.01$ does not improve the density profiles. The interaction parameter is equal to $g_{BF} = -1$. (Right) Densities of the system for the A-RTP algorithm for different ramps. The time step is $\Delta t = 0.05$ and the interaction parameter is $g_{BF} = -1$. Red lines show $40\times 10^3$ time unit ramp while blue lines are $20\times 10^3$ time unit ramp. Lower time step or longer ramp does not improve the density profiles.
• Figure 5: Results for the ITP-IEV-3D method for different numbers of basis eigenvectors for $g_{BF} = -3$ (main panels) and $g_{BF} = -1$ (insets). (Left) Densities of the bosons. (Right) Densities of the fermions. Different colors indicate different numbers of basis eigenvectors. The system consists of $N_B = 40$ bosons and $N_F = 4$ fermions.