Doublon dynamics of Bose-Fermi mixtures in optical lattices

We study the out-of-equilibrium dynamics of a dilute, lattice-confined Bose-Fermi mixture initialized in a highly excited state consisting of boson-fermion pairs (doublons) occupying single lattice sites. This system represents a paradigmatic case for studying relaxation dynamics in strongly correlated systems, and provides a versatile platform for studying thermalization and localization phenomena. We provide analytical expressions for the short-time decay of isolated doublons and small doublon clusters due to the competition between tunneling and interparticle interactions. We also discuss a mechanism for long-time decay that crucially depends on the quantum statistics of the particles constituting the doublon, namely, the conversion of pairs of neighboring doublons into an unpaired fermion and a site with a fermion and two bosons. Building on these insights, we develop a cluster expansion method to describe the dynamics in extended systems and compare it to numerically exact matrix product state simulations in one dimension. Finally, we discuss how our predictions can be observed in experiments with ultracold heteronuclear molecules.

Figure 1

(a) Schematic illustration of the cluster expansion method: the initial state of doublons at random positions is partitioned into connected clusters. (b) Relevant states and parameters in the one and two doublon case. Both the bosons and the fermions are in the ground band of the optical lattice. Their tunneling rates are in general different for different masses and polarizabilities. (c) Eigenenergies in the single doublon case for $J_{\text{B}}=0$. The ground state (bound state) is shown as a black (red, blue) solid line for 1D (2D, 3D). Both axes have been scaled by half the width of the band of scattering states (blue shaded area) $2d{J}_{\text{F}}$. In 2D and 3D there is no bound state below a certain critical interaction strength (gray dotted lines). The black dashed line indicates the “bare” bound-state energy $U_{\text{BF}}$.

Figure 2

(a) Time evolution of the doublon fraction for a single-doublon initial state in 1D ($U_{\text{BF}}/J_{\text{F}}=3$, $J_{\text{B}}=0$). The dotted line indicates the saturation value $P_{\text{sat}}$. (b) Dependence of the doublon fraction on the interspecies interaction strength. Dashed lines: analytical prediction for $P_{\text{sat}}$ from perturbation theory. Linear system sizes for the (1D,2D,3D) cases are $L=500,40,13$, respectively. In 1D (for $J_{\text{B}}=0$) the predicted saturation value is exact. The vertical dashed lines indicate the positions of $U_c$ in 2D and 3D, below which no bound state exists [compare Fig. 1]. Note that this point is not defined sharply due to the finite system size. (c) Time evolution of the doublon fraction for initially two doublons on neighboring sites. The cases of two representative parameter sets are shown. Red (upper) and black (lower) solid lines are exact numerical results for these two parameter choices. The dashed lines are results of the analytical formulas in 1D, i.e., perturbative treatment (red, upper) and quasicontinuum model (black, lower). (d) Time evolution of the triplon number. For the analytical results (dashed lines) our assumptions imply that $N_{\text{tri}}\left(t\right)=1-P_2\left(t\right)$. We used a 1D lattice of length $L=14$ with the two doublons initially on the two central sites. Parameters are as follows. Case 1 (black, upper): $J_{\text{B}}/J_{\text{F}}=0.1$, $U_{\text{BF}}/J_{\text{F}}=-10$, $U_{\text{BB}}/J_{\text{F}}=1$; case 2 (red, lower): $J_{\text{B}}/J_{\text{F}}=0.2$, $U_{\text{BF}}/J_{\text{F}}=-30$, $U_{\text{BB}}/J_{\text{F}}=3$.

Figure 3

(a),(b) Long-time evolution of the doublon density in a 1D lattice ($L=1010$) with two initial doublons computed using t-DMRG. Parameters are $J_{\mathrm{B}}=0.1J_{\mathrm{F}}$, $U_{\mathrm{BB}}=J_{\mathrm{F}}$, and $U_{\mathrm{BF}}=10J_{\mathrm{F}}$. In (a) the initial doublons are located on neighboring sites, and there is a fast decay on a time scale of $t\sim J_B^{-1}=10J_{\mathrm{F}}^{-1}$. In (b) the initial doublons are separated by an empty site. Here, the doublons are long-lived and diffuse on a time scale of $t\sim J_{\mathrm{D}}^{-1}=100J_{\mathrm{F}}^{-1}$. In this case the doublons do not decay by forming a triplon and singlon due to the quantum Zeno effect (see text). (c),(d) Time evolution of a system with $L=100$ sites and a random distribution of initial doublons with different filling fractions of $f=5\%$ and $f=20\%$ [panels (c) and (d), respectively]. For low fillings and short times, the exact t-DMRG solution (solid black line, average over 30 initial configurations) is well reproduced by considering dynamics of small connected clusters. The dashed blue (dotted red) line shows weighted averages of evolution of connected clusters with two (three) doublons, respectively.

Figure 4

(a) Energy diagram for states accessible to a two doubloon cluster. For large $U_{\text{BF}}$ the state space can be restricted to the lowest two manifolds as indicated by the red ellipse. (b) These states are connected through different terms in the Hamiltonian. Typically the coupling $J_{\text{F}}$ is large compared to $J_{\text{B}}$. In this case we diagonalize the singlon-triplon subspace (gray box). (c) The resulting problem is a single state (two-doublon state) coupled to a continuum of singlon-triplon states. The gray shading indicates the density of states, which is minimal in the band center in 1D but maximal in the band center in higher dimensions. We have shifted the origin of the energy to be at the center of the band. (d) Relevant states and coupling matrix elements for the three-doublon problem.

Figure 5

Doublon decay in 2D. Fitted decay rate $\mathrm{\Gamma }$ from fitting $P_2\left(t\right)$ with $exp[-\mathrm{\Gamma }t]$ as a function of $U_{\text{BB}}/J_{\text{F}}$ at $J_{\text{B}}/J_{\text{F}}=0.1$ (left) and of $J_{\text{B}}/J_{\text{F}}$ at $U_{\text{BB}}/J_{\text{F}}=1$ (right). Solid lines are a linear fit to the region $1<U_{\text{BB}}/J_{\text{F}}<2d$, and a quadratic fit, respectively. A lattice of size 22 by 23 sites was used with the doublon pair initially sitting on the two central sites.