Interplay between coherent and dissipative dynamics of bosonic doublons in an optical lattice

We observe the dissipative dynamics of a dense, strongly interacting gas of bosonic atom pairs in an optical lattice, controlling the strength of the two-body interactions over a wide parameter regime. We study how three-body losses contribute to the lattice dynamics, addressing a number of open questions related to the effects of strong dissipation in a many-body system, including the relationship to the continuous quantum Zeno effect. We observe rapid break-up of bound pairs for weak interactions, and for stronger interactions we see doublon decay rates that are asymmetric when changing from attractive and repulsive interactions, and which strongly depend on the interactions and on-site loss rates. By comparing our experimental data with a theoretical analysis of few-body dynamics, we show that these features originate from a non-trivial combination of dissipative dynamics described by a lattice model beyond a standard Bose-Hubbard Hamiltonian, and the modification of three-atom dynamics on a single site, which is generated alongside strong three-body loss. Our results open new possibilities for investigating bosonic atoms with strong three-body loss features, and allow for the better understanding of the parameter regimes that are required to realize strong effective three-body interactions.

We observe the dissipative dynamics of a dense, strongly interacting gas of bosonic atom pairs in an optical lattice, controlling the strength of the two-body interactions over a wide parameter regime. We study how threebody losses contribute to the lattice dynamics, addressing a number of open questions related to the effects of strong dissipation in a many-body system, including the relationship to the continuous quantum Zeno effect. We observe rapid break-up of bound pairs for weak interactions, and for stronger interactions we see doublon decay rates that are asymmetric when changing from attractive and repulsive interactions, and which strongly depend on the interactions and on-site loss rates. By comparing our experimental data with a theoretical analysis of few-body dynamics, we show that these features originate from a non-trivial combination of dissipative dynamics described by a lattice model beyond a standard Bose-Hubbard Hamiltonian, and the modification of three-atom dynamics on a single site, which is generated alongside strong three-body loss. Our results open new possibilities for investigating bosonic atoms with strong three-body loss features, and allow for the better understanding of the parameter regimes that are required to realize strong effective three-body interactions.
Ultracold atomic gases in optical lattices provide a platform for investigating novel many-body coherent and dissipative dynamics [1][2][3][4][5]. In particular, processes such as collisional loss of atoms, which we often seek to avoid, can exhibit features such as the continuous quantum Zeno effect [6][7][8][9], where a strong dissipative process can prevent coherent dynamics from taking place [10,11]. This has resulted in proposals to realize effective three-body interactions through three-body loss [12][13][14][15]. However, there are important open questions related to how these dynamics work when the loss rates become comparable to or larger than the energy band gap, rendering a standard Bose-Hubbard (BH) description insufficient. In this work, we explore the interplay between coherent and dissipative dynamics for a gas of bosonic atom pairs, comparing the dependence of the three-body loss on the interaction strength to theoretical models beyond the standard BH Hamiltonian that combine lattice dynamics with a careful treatment of on-site three-body dynamics that renormalizes the coefficients of coherent interactions and loss. Understanding these effects across different regimes provides a path for future studies of exotic quantum phases induced by strong local three-body loss [16][17][18][19][20][21][22][23].
We make use of the control available in optical lattice systems and study interacting Cesium atoms in the vicinity of a broad Feshbach resonance [24], which allows us to tune the s-wave scattering length a S and therefore the interactions from weak to strong and from repulsive (a S > 0) to attractive (a S < 0). The two-body properties of our system are well understood [25][26][27] and it is straight-forward to prepare an initial state with on-site pairs (often referred to as doublons) [28,29]. Doublons are in general stable in optical lattice systems because of the lack of dissipative phonon modes that can remove energy on short timescales [28,30]. However, by in-creasing the interaction strength one enters a regime where the energy associated with strong three-body dissipative processes [31,32] can easily exceed the band gap in the lattice. This raises questions on how to treat this system with an effective BH model [25,26,[33][34][35], including the question whether off-site loss mechanisms become important in understanding the resulting dynamics, and how effective model parameters would need to be rescaled in that case.
Below, we investigate the effects of strong three-body losses starting with doubly occupied sites and find unexpected phenomena arising from an off-resonant loss process where a single particle tunnels into a virtual triply occupied site. For strong attractive interactions, we observe an interesting nonlinear dependence of the doublon decay rates on a S , and attempt to model this with a BH-type model resulting in values for decay rates that are too large if nearest neighbor losses are included, but too low if they are simply ignored. This indicates that the mixing with higher energy bands is causing a renormalization of both the off-site and the on-site coefficients. When further increasing the attractive interaction strength, we observe a decrease on the decay rate similar to a quantum Zeno type suppression [11]. However, through firstprinciple calculations of the on-site losses induced through mixing with three-body Efimov states [31,32,[36][37][38][39], we find that on-site loss rates are too low to induce this type of suppression and we instead show that the decrease of the decay rates is due to an additional renormalization of the terms responsible for the coherent dynamics.
Our experiments start with the production of an essentially pure Bose-Einstein condensate (BEC) of 1.0 × 10 5 Cs atoms in the lowest hyperfine state in a crossed optical trap [40]. To prepare a comparatively dilute system of doublons, we non-adiabatically load the BEC within 2 ms into a cubic op- Total atom number n tical lattice formed by three retro-reflected laser beams with a wavelength of λ = 1064.5 nm, see Fig. 1(a). Initially, the lattice depths along all three directions are set to V 0 x,y,z = 30 E R , where E R = h 2 /(2mλ 2 ) is the recoil energy with the mass m of the Cs atom. During lattice loading a S is set to 220 a 0 , where a 0 is Bohr's radius. The doublon sample is then purified by Feshbach molecule formation, cleaning and dissociation [41].
Next, we tune a S to a value between −800 a 0 and +800 a 0 by ramping the offset magnetic field B with a rate of ∼ 2.5 G/ms to a value between about 8 G and 37 G, exploiting a broad Feshbach resonance, see Fig. 1(b) [42]. Note that the tuning range includes a zero crossing for a S located at 17.119 G. To induce dynamics of the doublons, we quench the lattice depth along the z-direction within 1 ms to a new value V z , realizing nearly isolated 1D BH systems. After a variable hold time t H we interrupt the dynamics by raising V z in 1 ms back to its original value. Finally, we determine the remaining number of doublons by employing the associationcleaning-dissociation procedure a second time before detecting the atom number after a short time-of-flight.
The doublon number n D shows a strong dependence on a S as can be seen in Fig. 1(c) for four different hold times t H . Clearly, our data exhibits a pronounced asymmetry between strongly attractive and strongly repulsive interactions. We first focus on the region near the zero crossing. A close-up of this region is presented in Fig. 1(d) for four different values of V z and a hold time t H = 100 ms for which the dynamics has reached a steady state. A pronounced dip for n D centered at a S = 0 is evident. In this regime of comparatively weak interactions a standard lowest-band BH description [36,43] should be valid. The dip, which is essentially symmetric as a function of a S = 0, results from rapid dissociation of the doublons when the BH interaction energy |U| and the tunneling energy J become comparable, irrespective of the sign of U. Indeed, when plotting the data as a function of U/J we observe a collapse onto a single curve, see inset to Fig. 1(d). Note that n D remains non-zero at U/J = 0, which is a consequence of the finite size of the 1D systems, resulting in a non-zero probability to find two atoms on the same lattice site.
We stress that in this regime of comparatively weak interactions, the doublons dissociate but atoms are not lost from the sample. This is evident from Fig. 1(e), which depicts the dissociation dynamics for two values of U/J. While the total atom number stays constant, we observe a rapid decay for n D with a rate independent of U/J towards the steady-state doublon number shown in Fig. 1(d). Moreover, the initial decay time τ, obtained from an exponential fit to the data, provides a direct measure of J and obeys the linear relation τ =¯h 4J , see [36]. Comparison with the experimental data as shown in Fig. 1(e) gives good agreement. The small deviations probably arise from finite size effects. Also, the steady-state doublon number is well captured, see the inset to Fig. 1

(d).
Let us now turn to the unexplored regime of strong interactions, where the standard BH description no longer applies. As evident from Fig. 1(c), we observe increased loss of doublons with increasing positive a S , but not so for negative a S . Figure 2(a) shows the typical time evolution of n D for strong repulsive interactions, along with the number of singly occupied sites n S and the total atom number n A . In contrast to the case of weak interactions, the loss of double occupancy is caused by a true loss of atoms from the system. We observe that the number of singly occupied sites increases until all doublons are lost. The ratio of lost doublons to lost atoms is 2.1(2) : 3, the ratio of lost doublons to reappearing atoms at singly occupied sites is 2.1(3) : 1. This represents strong evidence that the dominant loss process is three-body recombination, where two neighboring doublons exhibit a small probability to decay via a triply occupied site, leaving one singly occupied site behind. Note that this mechanism differs from the elastic decay observed for doublons composed of fermions [29].
For a quantitative comparison with a theoretical model, we perform the same experiment in a dense system, i.e. starting from a Mott shell with predominant double occupancy, cleaned from singly occupied sites in the same way as before. This ensures that tunneling dynamics of the doublons themselves do not play any role in the initial dynamics. We estimate the filling fraction to be close to unity for an initial length of about 20 sites for the 1D systems. Figure 2(b) depicts typical decay curves for n D from which we deduce initial doublon decay rates γ D from exponential fits to the data. The decay rates slow down with increased loss as the system gets more dilute. We hence use only the initial parts of the loss curves for the fit.
In this way, we measure γ D for different values of a S and V z . Figure 3(a) shows γ D for a very broad range of interaction strength, from −2000 a 0 to +800 a 0 , revealing the drastic asymmetry between repulsive and attractive interactions that could already be seen in Fig. 1(c). Note that we limit the hold time to 400 ms to avoid tunneling dynamics between neighboring 1D systems. This, however leads to comparatively large uncertainties for small decay rates.
We model this behavior by an extended BH model that accounts for on-site and nearest-neighbor two-body interactions as well as dissipative on-site three-body loss [36]. Note that we neglect terms corresponding to site-to-site three-body loss.
For large values of |a S |, one might expect to obtain large contributions from off-site loss processes. However, this picture only holds for sufficiently small loss rates compared to the band separation energy when the lowest-band Wannier functions provide a good local basis. In the presence of large loss rates, as is the case for large |a S |, the Wannier states mix with states from higher bands and we find that the initial overlaps of nearest-neighboring wavefunctions lead to a very rapid initial decay, although only resulting in a very small loss of atoms. This happens on timescales much faster than the experimental measurements, resulting in renormalized states that are more localized on each lattice site. Atom loss is then driven by on-site decay only, combined with non-resonant tunneling processes, including density-assisted tunneling from nearest neighbor interactions.
To understand the importance of competing on-site and tunneling processes in the dynamics of particle loss, we perform calculations in which we evolve the initial doubly-occupied Mott state using a Lindblad master equation. We incorporate three-body loss, and also single-particle loss terms that account for residual heating induced by the lattice light and/or background gas collisions [26]. We extract the decay rate γ D from an exponential fit to the computed doublon density. In order to identify the roles of different processes, we use two complementary methods to estimate coefficients for our model.
In the first approach, we calculate lattice coefficients using single-particle Wannier functions associated with the lowest Bloch band. On-site loss is modeled as a zero-range threebody recombination with a rate coefficient L 3 , with parameters for Cs fitted to the experimental data from Ref. [24]. As we will see below, the resulting model predicts decay rates at the same order of magnitude as the experimental measurements, capturing qualitative aspects well for positive scattering lengths. However, it fails to predict the finer details for negative scattering lengths due to the presence of a large threebody resonance, for which we clearly must go beyond a standard BH model description.
In our second approach, we determine more accurate onsite parameters by taking the lowest on-site three-body energy level obtained by using the hyperspherical adiabatic representation [31,38]. This allows us to calculate, from first principles, the on-site three-body losses from and energy shifts to the lowest band due to three-body interactions mixing with higher bands [37]. For on-site two-body interactions we also take into account mixing with higher bands using the result presented in Ref. [44] for a single harmonic trap. We find that by incorporating these corrections from strong interactions the main features near the resonance are better approximated. However, we do not modify the nearest-neighbor tunneling and interaction terms, and find that these must also play an important role away from the resonance.
In Fig. 3(a) we compare γ D obtained from the two approaches with the experimental data. We find that for repulsive interactions the decay rates increase with increasing a S , in accordance with the predictions. We also find that the first approach (dotted lines) predicts a minimum in the decay rate around a S = −800a 0 , which coincides with the strong Efimov resonance reported in L 3 [24], and which is not observed in our lattice experiment. Our second theoretical approach (solid lines with symbols), incorporating the effects of the lattice potential on the loss rates, predicts a suppression of the L 3 resonance and we obtain values for γ D that qualitatively match the experiment [36]. This peak suppression is due to the interplay between the Efimov resonance, the mixing with higher bands and the lattice potential, which creates a spread of the particles in momentum space. Indeed, this spread is of the same magnitude as the temperature-induced energy spread that was found to suppress the Efimov resonance in previous measurements [24,45]. In addition, effects of mixing with higher excited states introduce further modifications to the on-site, but also the off-site, loss coefficients producing the double peak feature shown in Fig 3(b).
Furthermore, including corrections from strong interactions (solid lines with symbols in Fig. 3(a)) predicts that, for increasingly large attractive interactions, the decay rates also increase, in sharp contrast to the experimental data. We initially interpreted the trend in the measurements to be due to a quantum-Zeno-type suppression of the losses [11], but the values of the renormalized coefficients in our calculation predict on-site loss rate values that are too low for this to occur [36].
Instead, we believe that this trend can be accounted for through the incorporation of two effects not currently taken into account in our model: mixing between ground and excited on-site three-body states [36], and a further renormalization of the effective tunneling elements [46]. Our experiment enables us to estimate what the values of the rescaled tunneling coefficients,J, would need to be to produce the measured decay rates, accounting for interaction-modified tunneling coefficients [36]. These are shown in Fig. 3(c), and predict a saturation of the magnitudes of the rescaled tunneling rate for negative a S , an effect that is more pronounced for shallow lattices. Additionally, for large positive a S the rescaled tunneling rates are significantly enhanced.
In summary, we have investigated the dynamics of bosonic doublons for a broad range of attractive and repulsive interactions. We have found that for strong two-body interactions (attractive and repulsive) the dynamics are dominated by offresonant decay processes induced by on-site three-body loss, with a large renormalization of the coefficients in both the dissipative and coherent processes, which can be understood by solving the three-body problem on-site. For particularly strong interactions, we enter a further regime, where strong interactions lead to further corrections, which can be attributed to a combination of modified effective tunneling rates, and further mixing of on-site three-body eigenstates.
In this way, our experiments and theoretical analysis both answer existing and (for stronger interactions) open new questions on the most accurate means to model dissipation for atoms in optical lattices. This is important for high-efficiency molecular formation [47], but also opens routes to further explore physics beyond the Bose-Hubbard model, and as well as to identify regimes of true quantum Zeno effect suppression of three-body occupation. This opens routes towards a range of exotic many-body states driven by three-body projections.
We We perform exact diagonalization of small 1D systems of typically 6 atoms on 11 lattice sites with hard wall boundary conditions. To mimic the dilute sample of randomly distributed doublons, we compute the time evolution of each possible initial doublon state and average the time evolution of the doublon probability.

EXTENDED BH MODEL
The 1D extended BH model that we use consists of on-site and nearest-neighbor two-body interactions and on-site threebody loss. In principle there can also be three-body recombination effects between atoms in nearest neighboring Wannier states. However, we obtain much better agreement with the experiment for negative values for a S if we neglect these contributions. This makes physical sense in the limit where the three-body loss rate dominates the dynamics, as any small initial overlap between neighboring states would be very rapidly lost, on timescales much shorter than the experimentally measured dynamics, resulting in more localized on-site wavefunctions. Explicitly, the Hamiltonian is given by (h = 1), The coefficients, calculated through the Wannier functions associated to the lowest Bloch band centered at position z i in the longitudinal direction, w(z − z i ), and r i in the radial direction, w ⊥ ( r − r i ), are given by As before, a s denotes the scattering length, J is the single particle nearest-neighbor tunneling amplitude, and Γ 3 = 2L 3 , where L 3 is the experimentally measured three-body loss parameter [24]. The approximation of restricting the dynamics to 1D is valid as long as the radial trapping frequency is much larger than the longitudinal trapping frequency, ω r ω z . Note that here, we have taken into account three-body Efimov resonances by incorporating the short distance cut-off term, which has to be introduced to regularize the zero-range threebody pseudo-potential [48,49], as a prefactor for a three-body delta-function contact interaction [12,50].
In order to extract the decay rates from this model, we consider four Cs atoms in two lattice sites and we begin with the initial state of two atoms on each site. We then apply a Lindblad master equation, where L = γ 3 /12b 3 i . We have included a single-particle loss rate, H eff = H + H s , where and L s = √S b i . For the loss rate, we use a valueS = 0.006 Hz, which was fine tuned so that this model would reproduce the lowest experimentally measured value. We then calculate the time-dependent expectation value of the total doublon number in the system, and fit an exponential to this function, N (t) /N(0) = exp(−αt), to extract the decay rate α.
Using the master equation we compare two cases. One case where the model coefficients are given by Eq. 2 with the free space L 3 parameter [24], and second case where we include corrections to the on-site lattice parameters due to couplings to higher energy states induced by two-and three-body interactions. For two-body interactions we make use of the analytical results presented in Ref. [44], and we explain how calculations of the three-body states lead us to corrections for the on-site three-body interaction and also give rise to an effective loss rate in the next section. Note that both calculations approximate each lattice site as an isotropic harmonic oscillator: the value of the energy spacing is given by the geometric mean of the three energies of the anisotropic experiment.
If we set the single particle loss rate,S = 0, we can calculate the decay rate through second-order perturbation theory in J U, where E 3B (a s ) and U 2B are the energy shifts due to three and two interacting atoms, respectively. The decay rate is then found through giving us an expression for the decay rate, α = 2Im(∆E (2) ).
We then use this expression to find the rescaled value for , such that the model (with rescaled on-site coefficients) gives rise to the experimentally observed decay rates. Note that we only perform this fit for larger values of |a s | (where the perturbation theory is most valid) because for lower values the calculated decay fits the experimental measurements quite well and is dominated by the single-particle losses.

ON-SITE THREE-BODY ENERGY LEVELS AND DECAY
Our three-body calculations for Cs atoms are performed using the the adiabatic hyperspherical representation [31,38], where the interatomic interactions are given by a two-channel model similar to the one used in Ref. [51], which properly describes the −11-G Feshbach resonance for the F = 3, m F = 3 hyperfine state of Cs [52]. Our methodology, developed in Ref. [53], incorporates the use of Feshbach projectors [54] to represent the fully symmetric three-body spin states as well as the relevant decay channels responsible for atomic losses. In the hyperspherical representation the hyperradius R determines the overall size of the system, while all other degrees of freedom are represented by a set of hyperangles Ω. Within this frame work, the three-body adiabatic potentials U and channel functions Φ are determined from the solutions of the hyperangular adiabatic equation: which contains the hyperangular part of the kinetic energy, expressed through the grand-angular momentum operator Λ 2 and the three-body reduced mass µ = m/ √ 3. In the above equation, our two-channel threshold energies and model potential are given, respectively, by where the background and resonant interactions are V bg (r) = −C 6 /r 6 (1 − λ 6 bg /r 6 ) and V res (r) = −C 6 /r 6 (1 − λ 6 res /r 6 ), with C 6 = 6890.4768 being the van der Waals dispersion coefficient for Cs [55], and interchannel coupling V c (r) = α c exp[−(r − β c ) 2 /2γ 2 c ]. Assuming that ε bg = 0, we choose the parameters λ bg , λ res , α c , β c , γ c and ε res in order to properly reproduce the magnetic-field dependence of the −11-G Feshbach resonance. Our calculations are performed assuming V bg and V res each containing four s-wave molecular states. The values for the energy and corresponding decay are then obtained by solving the hyperradial Schrödinger equation [38] in the presence of a harmonic confinement and leaving the adiabatic channels associated to the diatomic states open, thus allowing the state to decay [37].
The explicit energy and loss values for the lowest and first two excited three-body states are given in Fig. S1. We also include the interaction energy of the initial state in the experiment: two atoms on each lattice site. This allows us to compare the energy shift (see Eq. 5) that the system is required to overcome for each three-body state when the system tunnels from two doubly occupied sites into the triply occupied site, which is the non-resonant process responsible for the atomic loss in the experiment. We can see that for positive and small negative values of a S the three-body state closest to being resonant is the lowest, however, for large and negative values of a S , the most resonant state is the first excited one. This may indicate that this is the most likely state that is populated when the third particle tunnels into the site, which leads us to believe that the remaining discrepancies between our theory predictions and the experimental data presented in the main text could be due to couplings into these excited states.
However, the energy resonance condition is not the whole picture, there is still the question of the rescaling of the effective tunneling rates, which complicates this conclusion. And so it is likely that the correct on-site three-body state to use is a non-trivial superposition between the lowest and the first couple of excited three-body eigenstates. We leave the analysis of the coefficients in this superposition as a future objective.
In Fig. S1(b), we can see that for large negative scattering lengths the renormalized values for the onsite decay rates are two orders of magnitude smaller than those calculated with the free space L 3 parameter (see γ 3 in Equ. 2). This indicates that the onsite losses are too low to induce the quantum Zeno suppression usually expected in this regime.

FUTURE THEORETICAL ANALYSIS
In this section we summarize some of the important approximations that go into deriving our models. The purpose of this is to highlight the additional features that should be incorporated for future theoretical studies into the phenomena presented in the main text. FIG. S1. On-site three-body energies (left) and on-site loss rates (right) calculated using the the adiabatic hyperspherical representation [31,38], for V z = 16E R (blue) and V z = 10E R (red). We include the first, second and third three-body state (solid, dotted, dasheddotted). We also include the (two-body) interaction energy of the initial state of two atoms on each lattice site (left, dashed), calculated using Ref. [44]. And we include the on-site loss rate calculated from the BH model using the free-space L 3 parameter [24] (right, dashed). and the free-space Efimov resonance can allow for additional features on the observed decay rates as displayed in Fig. 3(c).
• Mixing of on-site three-body eigenstates: As shown in Fig. S1, particularly for large and negative values of a S , the energy of the initial state (dashed line) is very close to being resonant with the first excited three-body eigenstate, indicating that in this regime there is a strong possibility that these states will also be important and result in a further rescaling of the on-site loss rates. In the main text, we have attempted to predict the values of the renormalized tunneling terms that couple the initial state into a triply occupied site, but for this calculation we ignored the excited three-body eigenstates. One would then need to calculate the individual matrix elements describing the coupling of the initial state into each of the three-body eigenstates separately in order to better predict the overall decay rates.
• Restricting the dynamics to 1D: The strong coupling with states in higher bands may be enhancing the singleparticle tunneling in the z-direction and potentially increasing tunneling in the other two directions. Such effects could be reducing the validity of our 1D approximation.