Doublon dynamics of Bose-Fermi mixtures in optical lattices

TL;DR

This paper investigates the out-of-equilibrium relaxation and decay mechanisms of doublons in Bose-Fermi mixtures in optical lattices, providing analytical models, numerical simulations, and experimental relevance.

Contribution

It introduces a cluster expansion approach and analytical expressions for doublon decay, advancing understanding of thermalization in strongly correlated lattice systems.

Findings

Analytical decay expressions for short-time doublon dynamics

Identification of long-time decay mechanisms dependent on quantum statistics

Validation of models through matrix product state simulations

Abstract

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.

Tables

Table 1. Table 1: Typical experimental parameters for the case of fermionic 40 K and bosonic 87 Rb in a λ = 2 ​ π / k = 1064 𝜆 2 𝜋 𝑘 1064 \lambda=2\pi/k=1064\, nm optical lattice. The experimentally adjustable parameters are the lattice depth V B subscript 𝑉 B V_{\text{B}} (1st column, in units of lattice recoils of Rb, E R = ℏ 2 ​ k 2 / 2 ​ m Rb subscript 𝐸 R superscript Planck-constant-over-2-pi 2 superscript 𝑘 2 2 subscript 𝑚 Rb E_{\text{R}}=\hbar^{2}k^{2}/2m_{\text{Rb}} ) and the interspecies scattering length (2nd column, in Bohr).

Equations

$$H_{\mathrm{BF}}$$

$$-U_{\mathrm{BF}}\sum _{\mathbf{i}}{n}_{\mathbf{i},\mathrm{F}}{n}_{\mathbf{i},\mathrm{B}}+U_{\mathrm{BB}}/2\sum _{\mathbf{i}}{n}_{\mathbf{i},\mathrm{B}}(n_{\mathbf{i},\mathrm{B}}-1)$$

$$|\psi (t=0)⟩=\prod _{\mathbf{i}\in {\mathcal{N}}_{\mathrm{occ}}}{a}_{\mathbf{i}}^{†}{c}_{\mathbf{i}}^{†}|\mathrm{vac}⟩.$$

$$P_2(t)=\frac{1}{N}\sum _{\mathbf{i}}⟨\psi (t)|p_{2,\mathbf{i}}|\psi (t)⟩$$

$$P_2(t)\approx P_{\text{sat}}\left[1+\frac{4J_{\mathrm{F}}}{t{U}_{\mathrm{BF}}^2}{J}_1(2J_{\mathrm{F}}t)\cos \left(t\sqrt{4J_{\mathrm{F}}^2+U_{\mathrm{BF}}^2}\right)\right]$$

$$\sum _{c=1}^{n_c}{w}_c{P}_2^{[c]}(t),$$

$$H$$

$$H_0$$

$$H_1$$

$$\Delta E_k^{(2)}$$

$$\displaystyle=\left\{\begin{array}[]{ll}\frac{8J_{\rm F}^{2}}{L}\sum_{l\text{ odd}}\frac{-\sin^{2}(\pi l/L)}{U_{\rm BF}+2J_{\text{F}}\cos(\pi l/L)}&\text{if }k=0\\ \frac{8J_{\rm F}^{2}}{L}\frac{\sin^{2}(\pi k/L)}{2J_{\rm F}\cos(\pi k/L)+U_{\rm BF}}&\text{if }k\neq 0,k\text{ odd}\\ 0&\text{else.}\end{array}\right.$$

$$|\Delta {\varphi }_k^{(1)}⟩=\sum _l|{\varphi }_l⟩\frac{⟨{\varphi }_l|H_1|{\varphi }_k⟩}{E_k^{(0)}-E_l^{(0)}}$$

$$\displaystyle=\left\{\begin{array}[]{ll}2J_{\text{F}}\sqrt{2/L}\sum_{l\text{ odd}}\frac{-\sin(\pi l/L)}{U_{\text{BF}}+2J_{\text{F}}\cos(\pi k/L)}\ket{l}&\text{if }k=0\\ 2J_{\text{F}}\sqrt{2/L}\frac{\sin(\pi k/L)}{2J_{\text{F}}\cos(\pi k/L)+U_{\text{BF}}}\ket{0}&\text{if }k\text{ odd}\\ 0&\text{else.}\end{array}\right.$$

$$P_2(t)=|⟨\psi (0)|\psi (t)⟩{|}^2={\left|\sum _k{e}^{-i{E}_k^{(2)}t}{\left|⟨{\varphi }_k^{(1)}|0⟩\right|}^2\right|}^2$$

$$={\left|\frac{e^{-i{E}_0^{(2)}t}}{N_0^2}+\sum _{k\text{ odd}}\frac{2e^{-i{E}_k^{(2)}t}}{L{N}_k^2}{\left|\frac{2J_{\text{F}}\sin (\pi k/L)}{2J_{\text{F}}\cos (\pi k/L)+U_{\text{BF}}}\right|}^2\right|}^2$$

$$P_2(t)\approx$$

$$\approx$$

$$H_{\mathrm{free}}=-\sum _{\alpha =B,F}\sum _{x_{\alpha }=1}^L{J}_{\alpha }\left(|x_{\alpha }⟩⟨x_{\alpha }+1|+|x_{\alpha }⟩⟨x_{\alpha }-1|\right)$$

$$⟨R,r|$$

$$=$$

$$-J_F{e}^{iKR}\left(e^{-iK/2}{\psi }_K(r+1)+e^{iK/2}{\psi }_K(r-1)\right)$$

$$H_{\mathrm{free}}{\psi }_K(r)={ϵ}_K(k){\psi }_K(r)$$

$$H$$

$$H_0$$

$$H_1$$

$$\Delta E^{(2)}_{k}=\left\{\begin{array}[]{ll}\frac{8J_{\rm B}^{2}}{L}\sum_{l}\frac{\sin^{2}(\pi l/L)}{-U_{\rm BB}-2J_{F}\cos(\pi l/L)}&\text{if }k=0\\ \frac{8J_{\rm B}^{2}}{L}\frac{\sin^{2}(\pi k/L)}{2J_{K}\cos(\pi k/L)+U_{\rm BB}}&\text{if }k\neq 0\end{array}\right.$$

$$|\Delta {\varphi }_k^{(1)}⟩=$$

$$\displaystyle\left\{\begin{array}[]{ll}2J_{\rm B}\sqrt{2/L}\sum_{l}\frac{\sin(\pi l/L)}{-U_{\rm BB}-2J_{\text{F}}\cos(\pi k/L)}\ket{l}&\text{if }k=0\\ 2J_{\rm B}\sqrt{2/L}\frac{\sin(\pi k/L)}{2J_{\text{F}}\cos(\pi k/L)+U_{\rm BB}}\ket{0}&\text{if }k\neq 0\end{array}\right.$$

$$P_2(t)=$$

$$\approx$$

$$|{\varphi }_k^{\pm }⟩$$

$$E_k^{\pm }$$

$$H_1=2J_B\frac{1}{\sqrt{2}}\left(|1^{+}⟩+|1^{-}⟩\right)⟨0|+2J_B|{\varphi }_L⟩⟨0|+h.c.$$

$$|{\varphi }_0^{(1)}⟩$$

$$\displaystyle+\frac{\sin\left(\frac{\pi l}{L}\right)}{U_{\rm BB}+E_{l}^{-}}|\phi_{l}^{-}\rangle\biggr{)}-2J_{\rm B}\frac{1}{U_{\rm BB}}\ket{\phi_{\rm L}}$$

$$|{\varphi }_{l,\pm }^{(1)}⟩$$

$$|{\varphi }_{\mathrm{L}}^{(1)}⟩$$

$$E_0^{(2)}$$

Full text

I Introduction

Understanding the relaxation dynamics of quantum systems out of equilibrium is vital to many outstanding problems in physics. In particular in the case of strongly interacting many-body systems and in the absence of a quasiparticle context, the study of such relaxation has profound consequences. Topics of interest range e.g. from eigenstate thermalization in quantum statistical physics Deutsch (2018); D’Alessio et al. (2016); Nandkishore and Huse (2015) to fundamental questions in nuclear physics Berges et al. (2018) or cosmology Kofman et al. (1997). Due to the many-body nature of the problem, the theoretical analysis is very challenging and efficient numerical approaches are often restricted to one-dimensional (1D) systems Orús (2014); Schollwöck (2011); Verstraete et al. (2008). During the last decades it has become possible to realize strongly interacting closed quantum systems under highly controlled conditions using cold atoms trapped in optical lattices Bloch (2018); Bloch et al. (2012); Lewenstein et al. (2007). In particular, experiments realizing tunable Fermi Jördens et al. (2008); Esslinger (2010) and Bose-Hubbard models Greiner et al. (2002); Bakr et al. (2009) recently revealed a large range of interesting non-equilibrium phenomena such as many-body localization, strongly correlated multi-fermion scattering processes, or the direct observation of repulsively bound pairs (see Gross and Bloch (2017) for a recent review).

Hubbard models describing lattice-confined mixtures of bosons and fermions are less well studied despite the fact that they also exhibit a rich phenomenology of relaxation dynamics. The early experimental efforts in trapping Bose-Fermi mixtures Günter et al. (2006); Ospelkaus et al. (2006a, b) fueled theoretical investigations of equilibrium properties of these mixtures, as well as those of the resulting dipolar bosonic or fermionic molecules into which the atom pairs are assembled Lang et al. (2008); Danzl et al. (2010). For example, it was shown that Bose-Fermi mixtures can be used to stabilize a supersolid phase and charge density wave order, as well as quantum phases of composite fermions Titvinidze et al. (2008); Albus et al. (2003); Büchler and Blatter (2003); Lewenstein et al. (2004). In addition, Bose-Fermi mixtures have been of interest in the context of cold molecule creation, where a number of groups has successfully created ultracold heteronuclear molecules Ospelkaus et al. (2006c); de Miranda et al. (2011); Reichsöllner et al. (2017); Seeßelberg et al. (2018); Bohn et al. (2017).

Here we study the dynamics in a Bose-Fermi mixture prepared in an out-of-equilibrium state in which bosons and fermions are paired into doublons on different sites. The dynamics of a single doublon after a quench have been studied extensively theoretically for Bose-Bose Winkler et al. (2006); Denschlag and Daley (2006); Petrosyan et al. (2007); Wang et al. (2008); Valiente and Petrosyan (2009); Javanainen et al. (2010); Deuchert et al. (2012), Fermi-Fermi Valiente et al. (2010); Nguenang and Flach (2009); Ohashi (2008); Hofmann and Potthoff (2012) and Bose-Fermi doublons Piil et al. (2008). However, many aspects of the few and many-body physics of the doublons remain un-explored. We provide a description of the decay dynamics of isolated and small clusters of doublons which will be essential for benchmarking the dynamics observed in quantum simulators.

We find that in the case of Bose-Fermi doublons a decay channel involving the formation of triplons, i.e. sites occupied by a fermion and two bosons, facilitates the decay of doublons on neighboring sites – a process which is not present for fermionic doublons. We study these processes and their importance on various time-scales using both analytical and numerical methods. We use our analytical results for small doublon clusters to efficiently calculate the short-time dynamics. We show that for the experimentally relevant case of strongly imbalanced tunneling rates the system can be understood in terms of the dynamics of highly mobile fermions in a lattice of defects formed by the heavier bosons. In the case where the inter-species interaction strength $U_{\text{BF}}$ is the dominant energy scale, the coupling between different energy bands separated by multiples of $U_{\text{BF}}$ can be treated perturbatively as illustrated in Fig. 1(b). In this limit, the decay of adjacent doublons into singlon-triplon pairs can also be understood in terms of a perturbative treatment. While the dynamics of Bose-Fermi mixtures in high dimensional extended geometries is in general numerically intractable, we show that in certain parameter regimes and on short time scales, cluster expansion techniques [see Fig. 1(a)] provide a good approximation to the full solution. We benchmark our approximate methods with exact t-DMRG simulations White and Feiguin (2004); Daley et al. (2004); Schollwöck (2011) in 1d and using various analytical and numerical techniques.

The paper is structured as follows: In section II we introduce our model of Bose-Fermi mixtures in optical lattices. In section III we present our results on single and few doublon dynamics and the cluster expansion model. Section IV provides conclusions and an outlook on future directions, including a discussion of realistic experimental parameters for a possible implementation in a mixture of fermionic potassium and bosonic rubidium atoms. Details of our analytical calculations are given in the appendix.

II Model

In the following we study the non-equilibrium dynamics of a Bose-Fermi mixture in a $d$-dimensional cubic lattice of linear size $L$. Atomic motion is restricted to the ground band of the lattice potential. The system is thus described by the single-band Hamiltonian

[TABLE]

where $a_{\mathbf{i}}(a_{\mathbf{i}}^{\dagger})$ and $c_{\mathbf{i}}(c_{\mathbf{i}}^{\dagger})$ are bosonic and fermionic annhilation(creation) operators for a particle on the $d$-dimensional lattice site $\mathbf{i}$, respectively. The notation ${\langle\mathbf{i\,j}\rangle}$ indicates a summation over all neighboring sites and $J_{\rm F}$ and $J_{\rm B}$ are the tunneling rates for the fermionic and bosonic species. The occupation number operator for a particle at site $\mathbf{i}$ is given by $n_{\mathbf{i},B}=a_{\mathbf{i}}^{\dagger}a_{\mathbf{i}}$ ($n_{\mathbf{i},F}=c_{\mathbf{i}}^{\dagger}c_{\mathbf{i}}$). $U_{\rm BB}$ denotes the on-site interaction strength between the bosonic particles, while $U_{\rm BF}$ is the inter-species on-site interaction. The individual terms of Hamiltonian (II) are shown schematically in Fig. 1(b). Inspired by recent experiments with fermionic 40K and bosonic 87Rb atoms Covey et al. (2016), we consider repulsive intra-species interactions and attractive inter-species interactions. We have defined the inter-species interaction term in the Hamiltonian with a negative sign such that both $U_{\rm BF}$ and $U_{\rm BB}$ take positive values.

The system is initially prepared in a product state where each site is either empty or occupied by a doublon, i.e. a Bose-Fermi pair:

[TABLE]

Here $\mathcal{N}_{\rm occ}$ denotes the initially occupied sites, and $\ket{\rm vac}$ is the state of the empty lattice. Our main observable is the normalized fraction of remaining doublons as function of time

[TABLE]

where $p_{2,{\bf i}}=a_{\bf i}^{\dagger}c_{\bf i}^{\dagger}\ket{\rm vac}\bra{\rm vac}c_{\bf i}a_{\bf i}$ are projection operators onto local single doublon states. The sum runs over all $L^{d}$ sites of the lattice, and $\ket{\psi(t)}=\exp(-{\rm i}Ht)\ket{\psi(0)}$ is the time-evolved state ($\hbar\equiv 1$). This choice of initial state and observable is motivated by recent experiments which have shown that both the preparation of an initial state of doublons and the detection of the doublon fraction is possible using magneto-association and dissociation of the doublons into Feshbach molecules Covey et al. (2016), see also Sec. IV.

In general, dynamics under $H_{\text{BF}}$ leads to a decay of $P_{2}(t)$. In the following we identify the main processes contributing to this decay by breaking down the full dynamics into dynamics of small clusters with single, two and three doublons [Fig. 1(a)].

III Results

III.1 Single doublon

The decay dynamics of a heteronuclear doublon are governed by the relative magnitudes of $J_{\rm B}$, $J_{\rm F}$ and $U_{\rm BF}$, as well as the dimensionality of the system. Qualitatively, for large $U_{\rm BF}$ ($U_{\text{BF}}\gg J_{\rm B},J_{\rm F}$) doublons are bound and their decay is suppressed. For smaller $U_{\rm BF}$ the dimensionality of the system plays an important role. In contrast to 1d systems which support a bound-state for all $U_{\rm BF}\neq 0$, 2d and 3d systems support a bound state only if $U_{\rm BF}>U_{c}$, with a critical interaction strength $U_{c}$, as shown in Fig 1(c). Hence for higher dimensional systems the doublon fraction decays to zero for $U_{\rm BF}<U_{c}$ and to a constant for larger attractive interactions. In 1d the problem of a single doublon can be solved analytically using the Greens function formalism described in  Denschlag and Daley (2006); Piil et al. (2008). In the limit where $J_{\rm B}\ll J_{\rm F}$, i.e. when the boson tunneling is negligible, the problem further simplifies and one can use a perturbative approach to find a compact expression for the doublon fraction (see appendix A for details),

[TABLE]

where $J_{1}(x)$ is the first Bessel function. It can be seen from the above expression that after initial oscillations at frequency $\sqrt{4J_{\rm F}^{2}+U_{\rm BF}^{2}}$, the doublon fraction saturates to $P_{\text{sat}}=U_{\rm BF}^{2}/(U_{\rm BF}^{2}+4J_{\rm F}^{2})$ after a time $t_{\text{sat}}\sim 1/(2J_{\rm F})$, see Fig. 2(a). This means that for $U_{\text{BF}}\gg J_{\text{F}}$, the doublon fraction saturates to a non-zero value while for $U_{\text{BF}}\lesssim 2J_{\text{F}}$ a significant doublon decay is observed as shown in Fig. 2(b). This is expected from a spectral analysis of the relevant state space shown in Fig. 1(c). Dividing the relevant states into bound states with $E_{\rm B}<-2dJ_{\rm F}$, connected to the initial state in the limit $J_{\text{F}}\rightarrow 0$, and scattering states, we expect that these two manifolds become strongly mixed as soon as the bare bound state energy $U_{\text{BF}}$ enters the band of scattering states, which happens at $U_{\text{BF}}=2J_{\text{F}}$ in 1d. The expression (4) can be generalized to the case of finite $J_{\text{B}}$ by separating the problem into center-of-mass and relative motion before applying the perturbative treatment to the relative motion part, see appendix A. In higher dimensions it is not possible to obtain an analytical expression for $P_{2}(t)$. However, in the limit $J_{\rm F}\ll U_{\rm BF}$ one can obtain a good approximation for $P_{\text{sat}}$ by replacing $J_{\text{F}}$ by $\sqrt{d}J_{\text{F}}$ in the above expression. The statement that strong doublon decay occurs as soon as bound and scattering states overlap energetically ($U_{\text{BF}}\lesssim 2dJ_{\text{F}}$) still holds.

In the regime where $J_{\rm F}\gg J_{\rm B}$ the short-time behavior of systems of many doublons will be dominated by single doublon decay, as we will see when discussing the cluster expansion approach. This is consistent with the experimental observations made in Covey et al. (2016). The following analysis goes beyond this regime and takes into account additional decay channels that open through interactions between doublons.

III.2 Two and three doublons

The single doublon decay described in the previous section occurs at a rate $J_{\rm F}$ and determines the short-time dynamics of the system. In the case of two doublons two additional processes have to be considered: i) If two doublons initially occupy neighboring sites, the tunneling of the boson results in the production of a triplon (two bosons and one fermion) and a singlon (lone fermion), cf. Fig. 1(b). The singlon can subsequently move through the lattice freely. This process occurs at a rate $\sim J_{\text{B}}$. ii) Another important indirect process is doublon tunneling to neighboring sites, which occurs at an effective tunneling rate $J_{\rm D}=2J_{\rm F}J_{\rm B}/U_{\rm BF}$.

In the case that $U_{\text{BF}}$ is the largest scale in the problem, the doublon tunneling is slow compared to all other processes mentioned so far (see also the discussion below). Moreover, in the regime $J_{\text{F}}\ll U_{\text{BF}}$, the single doublon process of the fermion tunneling off the doublon site is strongly suppressed. Thus the process relevant to the decay of two neighboring doublons is the formation of the triplon-singlon complex and the subsequent tunneling of the singlon.

Under these approximations, the dynamics is restricted to the initial two-doublon state and the singlon-triplon states. In order to find analytical expressions for the for the doublon fraction, we employ the perturbative approach used in the single-doublon case, supplemented by a quasi-continuum treatment of the scattering states (see appendix A for details). Here, the two-doublon configuration takes the role of the bound state while the singlon plus triplon sector represents the scattering states. The energy splitting and coupling between the two sectors are given by $U_{\text{BB}}$ and $J_{\text{B}}$, respectively [see Fig. 1(b)]. The width of the band of scattering states is $4dJ_{\text{F}}$. Thus, in analogy with the single-doublon case, doublons are stable if the two sectors are energetically well separated ($U_{\text{BB}}\gg 2dJ_{\text{F}}$) and decay if the bound and scattering states overlap ($U_{\text{BB}}<2dJ_{\text{F}}$). In the former regime a perturbative treatment yields saturation value of $P_{\text{sat}}=U_{\rm BB}^{2}/(U_{\rm BB}^{2}+8J_{\rm B}^{2})$ in 1d. In the latter regime, a quasi-continuum treatment of the scattering states shows that the doublon fraction decays exponentially with rate $\Gamma=8J_{\text{B}}^{2}/J_{\text{F}}\sqrt{1-[U_{\text{BB}}/(2J_{\text{F}})]^{2}}$ (see appendix A). In higher dimensions the behavior is qualitatively the same: If $U_{\text{BB}}\gg 2dJ_{\text{F}}$ two doublons on neighboring sites are essentially stable and $P_{2}(t)$ quickly saturates at a value near unity, while for $U_{\text{BB}}<2dJ_{\text{F}}$, the doublon pair decays exponentially with a rate $\Gamma\propto J_{\text{B}}^{2}/J_{\text{F}}$. For $J_{\text{B}}\ll J_{\text{F}}$ the transition between both regimes is sharp.

Figures 2(c) and (d) show the time evolution of the doublon and triplon fraction for two sets of parameters in 1d. This confirms that for $U_{\text{BB}}<2dJ_{\text{F}}$ the doublon fraction decays quickly while for $U_{\text{BB}}>2dJ_{\text{F}}$ it saturates. The triplon number increases in the same way the doublon fraction decreases, showing that the doublon decay is indeed caused by the formation of triplons. It should be noted that in our perturbative treatment the sum of the doublon fraction and the triplon fraction is always 1, as a direct result of the Hilbert space truncation.

In the case of three doublons on neighboring sites we observe that two among them can form a singlon plus triplon configuration, which leads to a decay of initially three doublons to one doublon in the regime $U_{\text{BB}}<2dJ_{\text{F}}$ (see appendix A for details).

Having discussed the decay dynamics of a cluster of doublons initially sitting on neighboring sites we now discuss the case of two doublons separated by an empty site and $J_{\rm B}\ll J_{\rm F}\ll U_{\rm BF}$. This discussion will justify the use of only connected clusters in our cluster expansion approach. In the limit $J_{\rm F}\ll U_{\rm BF}$ a single doublon is stable ($P_{\rm sat}\sim 1$), while two neighboring doublons may decay by forming a singlon-triplon configuration at rate $\Gamma$. Since doublon tunneling is strongly suppressed ($J_{\text{D}}\ll U_{\rm BF}$), the separated doublon pair is governed by the single doublon decay and thus stable. The doublons simply propagate on a time-scale given by $J_{\rm D}^{-1}$. This striking dependence on the initial doublon positions is confirmed in exact calculations shown in Figs. 3(a) and (b).

While one might expect to observe a decay process on a time scales $\sim J_{\rm D}^{-1}$ caused by decay after two doublons hopped on neighboring sites, this is not what we observe in the numerical computations of Fig. 3(b). The doublons remain stable in a much longer time scale than $J_{\rm D}^{-1}$. This can be understood by noting that $J_{\rm D}\ll\Gamma$, which places the system in a “quantum Zeno regime”. In this regime, the coherent pair tunneling happens at a much slower rate compared to the decay process at rate $\Gamma$ which connects the two-doublon state to a continuum of scattering states consisting of a fixed triplon and a mobile singlon. As a result the decay process for two doublons initially separated by an empty site is suppressed. A perturbative expansion in $J_{\text{D}}/\Gamma$ shows that the effective doublon decay rate is given by $\Gamma_{\rm eff}\sim J_{\rm D}^{2}/\Gamma$. For the parameters considered in Fig. 3(b) this corresponds to a decay process occurring on a time-scale of $\sim 200J_{\rm F}^{-1}$, consistent with our numerical observations. The discussed effects make it plausible why one may accurately describe the dynamics of a dilute system by considering solely the dynamics of connected few-doublon clusters.

III.3 Extended systems: Cluster expansion

Next, we consider extended systems of randomly placed doublons at a given filling fraction. An upper limit on the propagation speed of particles is given by the maximal group velocity of the most mobile particles, typically the fermions, $v_{g}=\partial E/\partial k\leq 2J_{\text{F}}$ (1d) limiting the speed at which correlations can spread in the system. Thus, in order to calculate the value of a local observable at site $i$ after time $t$ it is always sufficient to simulate a volume of radius $r=2J_{\text{F}}t$ around site $i$. If the system shows some kind of localization length $\xi<2J_{\text{F}}t$ a smaller simulation volume might be sufficient. For a translationally invariant system this is equivalent to simulating a box of linear size $L>2J_{\text{F}}t$ with periodic boundary conditions and evaluating the observable at all sites of the box. Thus short-time dynamics can always be simulated exactly.

As described above, we know that in the case of large $U_{\text{BF}}$ the decay rate of two initially separated doublons is small. Thus, in this case a viable ansatz is to consider only the dynamics of connected clusters up to a certain size. This cutoff size can be small for low densities. For example, in 3d, at a filling fraction of $f=10\%$ the fraction of doublons that are part of a connected cluster of size larger than $3$ is about $1.5\%$. In this case the process dominating the doublon decay, is the triplon formation from nearest neighbor doublons. This cluster expansion is schematically illustrated in Fig 1(a).

In Fig. 3(c) and (d), the validity of this phenomenology is illustrated. There we compare an exact t-DMRG calculation in a 1d box with $L=100$ sites and initial small filling fractions $f=5\%$ and $f=20\%$. To model the exact evolution we compare these results to weighted averages

[TABLE]

where $P_{2}^{[c]}(t)$ denotes the evolution of a cluster with $c$ connected doublons and $w_{c}$ the probability of the occurrence of the cluster in the initial state. Clearly this simple approach reproduces the exact calculation in the limit of short times and small filling fractions.

IV Conclusions & Outlook

In conclusion we have studied the non-equilibrium dynamics of a lattice confined Bose-Fermi mixture initialized in a state consisting of doublons randomly placed on the lattice sites at a given dilute filling fraction. We found that the short-time relaxation dynamics can be understood in terms of the decay of connected doublon clusters. We have derived analytical expressions for the decay of clusters of one, two, and three doublons valid when the inter-species interaction is the dominating energy scale. We found that the relaxation of the doublon population is governed by multiple different decay time scales resulting from different processes such as single-doublon decay, doublon-pair decay via triplon formation, and decay after doublon tunneling at even longer times. We verified our cluster expansion model by numerically exact t-DMRG simulations in 1d.

Our choice of initial state and observable is motivated by the fact that both can be realized straight forwardly in experiment using the ability to associate and dissociate Feshbach molecules. In Ref. Covey et al. (2016) a lattice gas of ground state KRb molecules was prepared where any remaining unpaired atoms can be removed from the lattice. After exciting the molecules to a weakly bound state (Feshbach molecules) they can subsequently be dissociated by an adiabatic ramp of the magnetic field, resulting in a configuration with all sites either empty of occupied by a doublon. The interspecies interaction strength can be adjusted by varying the endpoint of the magnetic field ramp. The remaining doublon fraction after a free evolution time can be determined by reversing the preparation process. Molecules only form on sites that are still occupied by a doublon. Unbound atoms can again be discarded and the molecules detected. Experiments that have taken advantage of the efficient conversions of pairs of the molecular constituents to molecules to create a low entropy lattice gas of polar molecules have recently been reported by a number of groups Moses et al. (2015); Takekoshi et al. (2012, 2014); Park et al. (2012, 2015); Gröbner et al. (2016); Molony et al. (2014), opening up a new direction for experimental investigation of non-equilibrium dynamical properties of Bose-Fermi mixtures.

In order to show that the regimes studied above are indeed experimentally relevant, we provide typical parameters for the case of 40K and 87Rb atoms in an optical lattice of wavelength $\lambda=1064\,$nm as used in Ref. Covey et al. (2016). The different masses and polarizabilities lead to different tunneling rates for the two species. We express all parameters in units of $J_{\text{F}}=J_{\text{K}}$. $U_{\text{BF}}$ can be adjusted independently of the other parameters by varying the magnetic field near an inter-species Feshbach resonance. Increasing the lattice depth leads to slower tunneling and enhanced on-site interactions. Table 1 shows that all the different regimes are experimentally accessible, and also that $U_{\text{BF}}$ is indeed always the largest parameter and $J_{\text{B}}\ll J_{\text{F}}$ in this case. Thus a time resolved measurement of the doublon fraction should reveal the effects found here, and in addition provide a means to explore the long-time dynamics of large 3d ensembles, which are inaccessible to numerical simulations.

The rich relaxation dynamics observed already for the restricted parameter regime considered here, and the prospect of experimental realization, call for a more systematic study of this system, addressing more parameter regimes and longer-time dynamics. Such a study would be quite challenging and require novel numerical and analytical methods for going beyond 1d and short times. Nevertheless, our results already showed a hierarchy of energy scales at which relaxation processes occur. Considering that such hierarchical relaxation is also found in spin glasses and soft matter physics, we take this as a hint for a rich phenomenolgy of relaxation dynamics in lattice Bose-Fermi mixtures yet to be uncovered. Adding an additional disorder potential, this system might be suited for studying localization phenomena in a new class of Hubbard models accessible to state-of-the-art quantum simulation experiments.

Appendix A Details of analytical calculations of few-doublon dynamics

In this appendix we present details of the derivation of our analytical results for single, two and three-doublon clusters. In each scenario we calculate the doublon fraction, $P_{2}(t)$ from Eq. (3). For two and three doublons, the initial state is a connected cluster where the doublons are initialized on neighboring sites.

A.1 Single doublon physics

In this section we describe two approaches which yield analytic expressions for single doublon dynamics in 1d. These approaches will then be used to also obtain analytical results for more than one doublon within certain approximations.

We first apply standard perturbation theory to the single doublon problem. The prescriptive nature of this method allows for straight forward extensions to larger number of doublons and higher dimensions. We start by considering the limit where one constituent of the doublon is stationary and can be approximated as a fixed defect. We then consider the case where both constituent atoms are mobile. Here we are able to build on the expressions derived in the prior approach and simply replace certain parameters.

A.1.1 Perturbation Method: Stationary Defect

Consider a system described by Hamiltonian (II) with $J_{\text{B}}=0$ on a ring of $L$ sites with periodic boundary conditions, where site $L$ is identified with site 0. We use basis states $\{|n\rangle\}$, with $n=0,\dots,L-1$ corresponding to a boson on site $i=0$ and the fermion on site $n$. In this basis the system is described by the Hamiltonian

[TABLE]

where we treat $H_{1}$ as a perturbation.

The eigenstates of $H_{0}$ are $\ket{\phi_{0}}=\ket{0}$, forming the bound state, and $\ket{\phi_{k}}=\sqrt{2/L}\sum_{n}\sin(\pi kn/L)\ket{n}$ which are the scattering states. The corresponding unperturbed eigenenergies are given by $E^{(0)}_{0}=-U_{\rm BF}$ and $E^{(0)}_{k}=2J_{\text{F}}\cos(\pi k/L)$.

The effect of $H_{1}$ on the eigenenergies is non-vanishing only at second order. The second-order perturbative corrections to the energies are given by

[TABLE]

The first order corrections to the eigenstates are given by

[TABLE]

With this we find that the remaining doublon fraction is given by

[TABLE]

where $N_{0}$ and $N_{k}$ are normalization factors. In the limit $U_{\rm BF}\gg 2J_{\rm F}$, we can approximate $2J_{\text{F}}\cos(\pi k/L)+U_{\text{BF}}\approx U_{\text{BF}}$ in the denominator. Neglecting terms of order $L^{-1}$ this leads to $N_{0}=1+2J_{\text{F}}^{2}/U_{\text{BF}}^{2}$, $N_{k}=1$, $E_{0}^{(2)}=-U_{\rm BF}-2J_{\rm F}^{2}/U_{\rm BF}$, and $E_{k}^{(2)}\approx E_{k}^{(0)}=2J_{\rm F}\cos(\pi k/L)$. With this we find

[TABLE]

where $\Delta E_{k}=E_{0}^{(2)}-E_{k}^{(0)}=-U_{\text{BF}}-2J_{\text{F}}^{2}/U_{\text{BF}}-2J_{\text{F}}\cos(\pi k/L)$ and $P_{\text{sat}}=1/(1+4J_{\text{F}}^{2}/U_{\text{BF}}^{2})$. This is the perturbative solution shown as a dashed line in Fig. 2(a) and (b).

A.1.2 Perturbation Method: Relative Motion

We now include the motion of the boson ($J_{\text{B}}>0$). We first consider the non-interacting Hamiltonian describing the motion of the boson and fermion in the lattice. We denote the position of the boson in the lattice with $x_{B}$ and the position of the fermion with $x_{F}$. The non-interacting part of the Hamiltonian thus reads

[TABLE]

We now introduce relative, $r=x_{B}-x_{F}$, and center-of-mass coordinates, $R=\left(x_{B}+x_{F}\right)/2$, a well as total and relative quasimomenta denoted by $K$ and $k$, respectively. Then, using the ansatz $\langle R,r|\Psi\rangle=\sqrt{1/2\pi}e^{iKR}\psi_{K}(r)$ for the relative wave function at a given momentum $K$ of the center-of-mass motion (with the obvious labeling $\ket{R,r}$ of the position basis states) we find

[TABLE]

The action of $H_{\rm free}$ on the relative coordinate wavefunction $\psi_{K}(r)$ can be written more compactly using the quantities $J_{\pm}=\frac{J_{F}\pm J_{B}}{2}$:

[TABLE]

with $\epsilon_{K}(k)=E_{K}\cos\left(\phi_{K}+k\right)$. Here we have introduced $E_{K}=-4\sqrt{J_{+}^{2}\cos^{2}\left(K/2\right)+J_{-}^{2}\sin^{2}\left(K/2\right)}$ and $\phi_{K}=\arctan\left(\frac{J_{-}}{J_{+}}\tan\left(K/2\right)\right)$.

We can now rewrite the Hamiltonian considered above for the stationary defect, interpreting the basis states $\{|n\rangle\}$ referring to the relative distance between the fermion and the boson, with $J_{F}\to E_{K}/2$. Applying a perturabtive treatment valid for $U_{\text{BF}}\gg J_{\text{B}},J_{\text{F}}$ results in a refined expression $P_{\rm sat}\to U_{\rm BF}^{2}/\left(U_{\rm BF}^{2}+8(J_{+}^{2}+J_{-}^{2})\right)$.

A.2 Clusters of two and three doublons

In the following we continue to assume that $U_{\text{BF}}$ is the largest energy scale such that we can restrict the Hilbert space to the manifold containing no lone bosons, indicated by the red ellipse in Fig. 4(a). Consider a chain of $2L-2$ sites, where the sites are indexed by integers in the intervals $\left[-L+1,-1\right]$ and $\left[1,L-1\right]$ (omitting [math] for symmetry). We denote the initial state in which the two doublons are placed on the nearest neighboring sites -1 and 1 by $\ket{0}$. The rest of the basis states are denoted by $\ket{n}$, $n=1,\dots,L-1$. These states are symmetrized states where the triplon and the fermion are separated by $n$ sites, cf. Fig. 4(b). Using this basis and in the limit of vanishing $J_{\rm D}$ the Hamiltonian takes the form

[TABLE]

We observe that by replacing $U_{\text{BB}}\rightarrow U_{\text{BF}}$ in $H_{0}$ and $J_{\text{B}}\rightarrow J_{\text{F}}/2$ in $H_{1}$ one obtains the single doublon Hamiltonian discussed above (except for the periodic boundaries). Thus the problem is reduced to a single state $\ket{\phi_{0}}=\ket{0}$ at energy $-U_{\text{BB}}$ coupling to a band of states $\ket{\phi_{k}}=\sqrt{2/L}\sum_{n}\sin(\pi kn/L)\ket{n}$ of width $2J_{\text{F}}$ [see Fig. 4(c)].

A.2.1 Perturbative treatment

If the the bound state and the band of scattering states are energetically separated ($U_{\text{BB}}>2J_{\text{F}}$), we can again employ the perturbative approach. The perturbative corrections to the eigenenergies and eigenfunctions of $H_{0}$ are given by

[TABLE]

and

[TABLE]

Using these expressions we find that the doublon fraction is given by

[TABLE]

with $\Delta E_{k}=-U_{\rm BB}-4J_{\rm B}^{2}/U_{\rm BB}-2J_{\text{F}}\cos(\pi k/L)$ and $P_{\text{sat,2}}=1/(1+8J_{\rm B}^{2}/U_{\rm BB}^{2})$, where we have replaced the denominator in the energy corrections by $U_{\rm BB}$.

In the following we will extend the perturbative method to three doublons. We consider a chain of $2L-3$ sites. The initial state consists of three doublons in three neighboring sites, denoted by $\ket{0}$. This state connects to two disconnected manifolds of states. The first set of states are labeled as $\ket{n}$ with $n=1,\cdots,L-1$, where $n$ labels the distance between the fermion and the triplon, and are shown in Fig. 4(d) panel $(i)$. These states are connected to the states $\ket{\eta}$ (see panel $(ii)$ in Fig. 4(d)) if the triplon and the doublon on its neighboring site exchange places. For these states $\eta=1,\cdots,L-1$ labels the distance between the fermion and the doublon. We use $\ket{n}$ and $\ket{\eta}$ to denote the symmetric combination of these two families of states. The second energetic manifold is given by the state $\ket{L_{\rm s}}$, the symmetric combination of the states $(iii)$ and $(iv)$ in Fig. 4(d). In this manifold the fermionic atom cannot move and thus these states, along with $\ket{0}$ form bound states.

In analogy with the previous perturbative treatments, we begin by determining the eigenenergies and eigenfunctions of $H_{0}$. The unperturbed localized states are $|\phi_{0}\rangle=\ket{0}$ and $\ket{\phi_{\rm L}}=\ket{L_{\rm s}}$ with $E_{0}=-U_{\rm BB}$ and $E_{\rm L}=0$. respectively. The unperturbed scattering states and their corresponding energies are given by

[TABLE]

where $|n^{\pm}\rangle=\frac{1}{\sqrt{2}}\left(|n\rangle\pm|\eta\rangle\right)$, $n=1,\dots,L-1$. We now proceed to find the perturbative corrections to the wavefunctions and the energies with the perturbation Hamiltonian

[TABLE]

To first order the eigenfunctions are given by

[TABLE]

and the second order energy corrections are

[TABLE]

We can now proceed as usual to find the time-dependence of the doublon fraction:

[TABLE]

where $\Delta_{l}^{\pm}=U_{\rm BB}+\frac{12J_{\rm B}^{2}}{U_{\rm BB}}-E_{l}^{\pm}$ and $\Delta_{L}=U_{\rm BB}+\frac{4J_{\rm B}^{2}}{U_{\rm BB}^{2}}$. The above expression saturates to $P_{\rm sat}=1/\left(1+16J_{B}^{2}/U_{\rm BB}^{2}\right)$ where we have used the approximation $U_{\rm BB}+E_{l}^{\pm}\approx U_{\rm BB}$.

A.2.2 Quasi-continuum model

We now examine the decay of two neighboring doublons in the regime $U_{\text{BB}}<2J_{\text{F}}$, i.e., when the two-doublon state lies inside the continuum of singlon-triplon states, see Fig. 4(c). In this case perturbation theory breaks down since the spacings between the levels of the singlon-triplon band are small compared to the perturbation $H_{1}\sim J_{\text{B}}$. We expect an exponential decay of the two-doublon state population and our objective here is to find its decay constant. We will use a treatment known as Bixon-Jortner quasi-continuum in quantum optics Bixon and Jortner (1968). We use the notation defined in the previous section.

The Hamiltonian in the basis of bound and scattering states defined above reads

[TABLE]

Note that the couplings $W_{k}$ and detunings $\Delta_{k}$ are real-valued. The time-dependent Schrödinger equation becomes

[TABLE]

where $b(t)$ and $c_{k}(t)$ are the time-dependent amplitudes of states $\ket{0}$ and $\ket{\phi_{k}}$, respectively. Initially the system is in state $\ket{0}$, thus $b(0)=1$ and $c_{k}(0)=0$. Applying a Laplace transform [$\bar{f}(s)=\mathcal{L}[f(t)]=\int_{0}^{\infty}e^{-st}f(t)dt$] this set of equations becomes

[TABLE]

Solving the second equation for $\bar{c}_{k}(s)$ and substituting into the first, we obtain

[TABLE]

Now we have to invert the Laplace transform to obtain $b(t)$.

At this point we make the following approximation: Since $J_{\text{B}}\ll J_{\text{F}}$, the state $\ket{0}$ only significantly couples to a small window within the band of states $\ket{\phi_{k}}$ [see Fig.4(c)]. If this window is small compared to the band width $2J_{\text{F}}$ we can assume that the energy differences between the states $\ket{\phi_{k}}$ are constant (linearize the spectrum, $\Delta_{k}\propto k_{0}-k$) and that the couplings $W_{k}$ are constant over the relevant range of $k$’s. Moreover, for a dense (quasi-continuous) spectrum, we can assume that there is always a resonant state $k_{0}$ defined by $\Delta_{k_{0}}=2J_{\text{F}}\cos(\pi k_{0}/L)+U_{\text{BB}}=0$. In quantum optics this is known as the Bixon-Jortner quasi-continuum Bixon and Jortner (1968). (Originally it was reported in the context of radiationless decay in molecules.) It is a good approximation in the center of the band, i.e., for $U_{\text{BB}}\ll 2J_{\text{F}}$. Near the band edges it can be problematic. We thus have

[TABLE]

where

[TABLE]

With these approximations the sum in $\bar{b}(s)$ becomes

[TABLE]

In the limit of large system size $L\rightarrow\infty$, i.e., $\Delta\rightarrow 0$, we can use $\lim_{x\rightarrow\infty}\coth(x)=1$, to obtain $\bar{b}(s)=(s+\Gamma/2)^{-1}$, with $\Gamma=2\pi W^{2}/\Delta$, which is obviously the Laplace transform of $b(t)=\exp(-\Gamma t/2)$. Note that $W^{2}/\Delta$ is independent of $L$.

To summarize the final result, we have

[TABLE]

with

[TABLE]

This model was used to produce the black dashed lines in Figs. 2c and d.

A.2.3 Higher dimensions

Consider the problem of 2 doublons on neighboring sites in $d$ dimensions ($d>1$). We again neglect all states that are detuned by at least $U_{\text{BF}}$ from the initial state. Analogous to the 1d case we can split the problem into the two-doublon state plus a continuum of singlon-triplon states. As in the case of a single doublon in higher dimensions, the problem is that the eigenstates of the singlon-triplon block cannot be found analytically. They are the eigenstates of a free particle on a $d$-dimensional lattice with a defect, i.e., a single site that is not accessible (the site that is occupied by the triplon). Unlike in the single doublon problem, the coupling of the bound state to the scattering states is $J_{\text{B}}$, which is different from the tunneling rate $J_{\text{F}}$ that determines the width $4dJ_{\text{F}}$ of the band of scattering states. Note that the triplon is in principle also mobile for $d>1$, since the fermion (singlon) can propagate to any of the $2d$ neighboring sites of the triplon and a boson can subsequently tunnel to that site to form two doublons on sites different from the initial doublon sites, which can then again decay into a triplon on the neighboring site of the original triplon site plus a singlon. This process leads to a center-of-mass motion of the cluster.

However, for the understanding of the short-time dynamics, it is sufficient to consider a single two-doublon state plus a continuum of scattering states. The global properties of the band are not changed significantly with respect to a perfect lattice by the addition of a defect, in particular the density of states and the width of the band stays the same. The crucial parameter is still the position of the two-doublon state with respect to the band: If it is outside the band ($U_{\text{BB}}>2dJ_{F}$), $P_{2}(t)$ will quickly saturate and $1-P_{\text{sat}}\propto J_{\text{B}}^{2}$, the square of the coupling matrix element. $P_{\text{sat}}$ will depend on $U_{\text{BB}}/J_{F}$ in a complicated way, and for $U_{\text{BB}}/J_{F}\gg 2d$ we have $1-P_{\text{sat}}\propto J_{\text{B}}^{2}/U_{\text{B}}^{2}$, however, in this limit $1-P_{\text{sat}}$ will be extremely small and the doublon pairs are essentially stable. Thus, for $J_{\text{B}}\ll J_{\text{F}}$, the doublon pairs become stable very quickly as soon as $U_{\text{BB}}>2dJ_{F}$.

In the regime $U_{\text{BB}}<2dJ_{\text{F}}$, where the bound state lies inside the band of scattering states, we can again employ the Bixon-Jortner quasi-continuum approximation, which means that we assume that the bound state couples to a continuum of equally spaced states (spacing $\Delta$) with a constant coupling strength $W$. If this is a valid description of the continuum, then in the limit of large $L$, the bound state population will decay like $\exp[-\Gamma t]$, with $\Gamma=2\pi W^{2}/\Delta$. In a higher-dimensional situation, for a given energy of the bound state, there are in general several resonant states ($\Delta_{k}=0$) and each of those states couples to the bound state with a different $W_{k}$. What we can do is to define $W^{2}$ as the average squared coupling of the scattering states at a given energy and $\Delta$ as the mean splitting between neighboring states in an energy interval, i.e., as the inverse of the density of states $\rho(E=-U_{\text{BB}})=\int d\mathbf{k}\delta(E(\mathbf{k})+U_{\text{BB}})$. Calculating the density of states involves an elliptical integral and cannot be solved analytically. However, it is clear that the density of states will be proportional to $1/J_{\text{F}}$ and thus $\Delta\propto J_{\text{F}}$. The squared coupling will be proportional to $J_{\text{B}}^{2}$. Thus $\Gamma=(J_{\text{B}}^{2}/J_{\text{F}})f(u)$, where $u=U_{\text{BB}}/J_{\text{F}}$. The function $f(u)$ depends on where in the band the bound state is located ($0<u<d$) and will decrease to zero at the band edge.

In Fig. 5 we show the fitted decay rates $\Gamma$, obtained by fitting the function $\exp[-\Gamma t]$ to the numerically obtained $P_{2}(t)$ as a function of the rescaled model parameters. We find that an exponential decay is indeed a good fit except close to the edges of the band $|U_{\text{BB}}|\approx 2dJ_{\text{F}}$. We observe that $\Gamma$ decreases to zero monotonically in the interval $0<U_{\text{BB}}/J_{\text{F}}<2d$ and is approximately linear for $U_{\text{BB}}/J_{\text{F}}>1$. In addition, $\Gamma\propto J_{\text{B}}^{2}/J_{\text{F}}$ is very well satisfied for $J_{\text{B}}/J_{\text{F}}\lesssim 0.5$. These statements also apply for $d=3$.

Acknowledgements.

We thank Jun Ye, Jacob Covey, Steven Moses, and Michael Wall for helpful discussions. This work is supported by the AFOSR grant FA9550-18-1-0319 and its MURI Initiative, the ARO single investigator award W911NF-19-1-0210, the NSF PHY1820885, NSF JILA-PFC PHY-1734006 grants, and by NIST. M.G. acknowledges support from the DFG Collaborative Research Center SFB1225 (ISOQUANT). J.S. is supported by the French National Research Agency (ANR) through the Programme d’Investissement d’Avenir under contract ANR-11-LABX-0058_NIE within the Investissement d’Avenir program ANR-10-IDEX-0002-02.