Doublon dynamics and polar molecule production in an optical lattice

Polar molecules in an optical lattice provide a versatile platform to study quantum many-body dynamics. Here we use such a system to prepare a density distribution where lattice sites are either empty or occupied by a doublon composed of an interacting Bose-Fermi pair. By letting this out-of-equilibrium system evolve from a well-defined, but disordered, initial condition, we observe clear effects on pairing that arise from inter-species interactions, a higher partial-wave Feshbach resonance and excited Bloch-band population. These observations facilitate a detailed understanding of molecule formation in the lattice. Moreover, the interplay of tunnelling and interaction of fermions and bosons provides a controllable platform to study Bose-Fermi Hubbard dynamics. Additionally, we can probe the distribution of the atomic gases in the lattice by measuring the inelastic loss of doublons. These techniques realize tools that are generically applicable to studying the complex dynamics of atomic mixtures in optical lattices.

P olar molecules with long-ranged dipolar interactions are ideally suited to the exploration of strongly correlated quantum matter and intriguing phenomena such as quantum magnetism, exotic superfluidity and topological phases [1][2][3][4][5][6][7][8] . The recent observation of the dipole-mediated spin-exchange interaction in an optical lattice 9 and the demonstration of the many-body nature of the spin-exchange dynamics 10 mark important steps for the use of polar molecules to study strongly correlated matter. While this initial work was done with a molecular filling fraction of only B5% in a threedimensional (3D) lattice 9 , more recent work has demonstrated a quantum synthesis scheme for molecule production in the lattice that relies on careful preparation of the initial atomic gases 11 . This led to a reduction in the final entropy of polar molecules by a factor of B4, and, correspondingly, a much higher filling fraction of B25% that opens up the possibility for studying non-equilibrium, many-body spin dynamics in a quantum gas of polar molecules where every molecule is connected to others.
The quantum synthesis approach reported in ref. 11 starts by preparing atomic insulator states that depend on atomic interactions, quantum statistics and low temperature 12 . However, realizing the full potential of this approach requires not only control over the atomic distributions, but also a detailed understanding of the molecule creation process.
Here we investigate this important step by leveraging our capability of molecule production in an optical lattice to create a clean system of doublons 13 . This technique allows us to additionally study 3D Bose-Fermi Hubbard dynamics. After creating ground-state molecules, we efficiently remove all unpaired atoms from the lattice and convert the molecules back to free atoms (in their lowest hyperfine states of 9=2; À 9=2 j ifor 40 K and 1; 1 j i for 87 Rb, where F; m F j idenotes the hyperfine state and its projection onto the magnetic field). This realizes a lattice where the sites are either empty or occupied by individual doublons that comprise a pair of bosonic and fermionic atoms. This well-defined initial state allows us to directly address limitations in the molecule creation process by probing the efficiency with which these doublons are converted back to molecules under various experimental conditions that affect atomic tunnelling rates, higher Bloch-band populations and the adiabaticity of a magnetic-field sweep through a higher partial-wave Feshbach resonance. Furthermore, this wellinitialized, non-equilibrium state of a disordered doublon distribution provides an ideal platform to explore the many-body dynamics of a lattice-confined Bose-Fermi mixture in a regime that is beyond the current simulation capabilities.

Results
Preparing the doublons. The experiment proceeds in steps as depicted schematically in Fig. 1. To prepare the doublons, we create a sample of molecules in their ro-vibrational ground state in the lattice as described in ref. 11 and then remove unpaired atoms with resonant light, so that all lattice sites are either empty or contain a single molecule. We then transfer the ground-state molecules back to a weakly bound Feshbach molecule state, followed by a magnetic-field (B) sweep to above the resonance to create a clean system of doublons. The solid black line in the upper panel of Fig. 1 shows schematically B relative to the s-wave Feshbach resonance (dashed line) that is used to manipulate the atomic inter-species interactions and to create molecules. After this preparation, the doublons are left to evolve in the lattice for a variable time t. Our measurement then consists of sweeping B to below the resonance to associate atoms into Feshbach molecules and determining the fraction of K atoms that form molecules. Specifically, we measure the molecule number using the following protocol. We first apply radio frequency (rf) to spin-flip the unpaired K atoms to another hyperfine state, which renders the unpaired K atoms invisible for subsequent molecular detection. We then sweep B back above the resonance to dissociate the molecules, and measure the number of resulting K atoms by spinselective resonant absorption imaging. The conversion efficiency is determined by dividing this molecule number by the total number of K atoms measured when we do not apply the rf. d-wave Feshbach resonance. We begin by investigating a narrow d-wave Feshbach resonance 14-17 that is located less than 0.1 mT above the 0.3-mT-wide s-wave resonance that is used for making molecules (Fig. 2a). With a pair of atoms confined in the same lattice site, the on-site density is orders of magnitude higher than that in ordinary optical traps, and thus this narrow Feshbach resonance can adversely affect the magneto-association process, where B is swept down from above the s-wave resonance to create molecules. Crossing the d-wave resonance too slowly will produce d-wave molecules, which will not be coupled to the ground state by the subsequent STIRAP laser pulses, as the process is weak and off resonance. If B is swept sufficiently fast to be diabatic for this narrow resonance (but still slow enough to be adiabatic for the broad s-wave resonance), crossing the d-wave resonance has no impact; however, the high effective densities at each site in an optical lattice can make it challenging to sweep fast enough. Although we study here specific resonances for the K-Rb system, the possibility of having to cross other Feshbach resonances and the issue of sweep speeds are general to magneto-association of atoms in an optical lattice.  Starting with a mixture of K, Rb and doublons (the smaller blue ball, the larger red ball and the pair grouped with grey background, respectively) in a 3D lattice, we sweep the magnetic field from above the s-wave Feshbach resonance (at B res s ¼ 54:66 mT) to below the resonance to create Feshbach molecules. These molecules are then transferred to their ro-vibrational ground state via STIRAP (stimulated Raman adiabatic passage). After unpaired atoms are removed with resonant light, the STIRAP process is reversed to transfer the ground-state molecules back to Feshbach molecules. The field is then swept above B res s to dissociate the molecules and create doublons. After holding for a time, t, at B hold , we measure the conversion efficiency when sweeping the field below B res s to re-form Feshbach molecules. To detect molecules, we use a rf pulse to spin flip the unpaired K atoms to a dark state (ball with black dashed edge) before dissociating the Feshbach molecules and imaging K atoms. The bottom panel illustrates possible dynamics of the doublons during B hold . As shown schematically, lattice sites populated with a K and a Rb atom have an interaction energy shift U 0 KRb . The K tunnelling energies in the lowest and first excited bands are denoted by J 0 K and J 1 K , respectively. Rb tunnelling happens at a slower rate since it experiences a deeper lattice.
In the experiment illustrated in Fig. 1, we investigate the d-wave resonance by varying the rate, _ B, and the final value, B hold , of the sweep that creates doublons. We then measure the subsequent molecule conversion efficiency after t ¼ 1 ms using a fast 1.68-mT ms À 1 magneto-association sweep. Figure 2b illustrates relevant states, above and below the resonance, for two atoms in a lattice site 18,19 : at low fields these are the s-wave molecule ( 1 ), d-wave molecule ( 5 ) and unbound atoms ( 4 ), and at high fields these are unbound atoms in the ground band of the lattice ( 2 ) and atoms with a band excitation in their relative motion ( 3 ). For simplicity, we illustrate states for a harmonic potential whose trap frequency o is the same for both atoms, with eigenstates of relative motion denoted by u ¼ 0; 1; 2. The dashed arrows show the diabatic ( 1 -2 ) and adiabatic ( 1 -3 ) trajectories for the dissociation of s-wave Feshbach molecules when crossing the d-wave resonance, while the solid arrows show the diabatic trajectories ( 2 -1 and 3 -4 ) for the subsequent, fast magneto-association sweep. Figure 2c shows the measured molecule conversion efficiency as a function of _ B when sweeping across the d-wave resonance from 54.56 to 56.24 mT (there are no other resonances in this field range), while Fig. 2d shows the effect of the final field B for a relatively slow, 0.018-mT ms À 1 , sweep. The data are taken for lattice depths of V latt ¼ 35E R (circles) and 30E R (diamonds), where E R ¼ : 2 k 2 /(2m Rb ) is the recoil energy for Rb, m Rb is the Rb atom mass, k ¼ 2p=l and l ¼ 1064 nm. For our highest sweep rates, or when B hold is below the d-wave resonance, the measured molecule conversion efficiency is near unity. This high conversion of doublons [20][21][22] is crucial for the quantum synthesis approach to producing molecules with a high filling fraction in the lattice. The near-unity conversion also provides an excellent starting point for diagnosing potential limitations to molecule production, and the data in Fig. 2c,d clearly show the negative effect that the d-wave resonance can have on magneto-association in the lattice.
The lines in Fig. 2c,d show fits used to extract the width (D d ) and position of the d-wave resonance. We use a Landau-Zener formalism 23 where the probability to cross the d-wave resonance diabatically, and therefore create s-wave Feshbach molecules in the subsequent magneto-association step, is where A depends on the on-site densities and the Feshbach resonance parameters. By approximating the sites in the deep optical lattice as harmonic oscillator potentials, we extract D d where o HO is the harmonic trap frequency for relative motion of the two atoms (see Methods) and L HO ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ' =ðmo HO Þ p is the harmonic oscillator length with the doublon reduced mass m (ref. 24) (we note that the right-hand side of equation (26) in this reference is missing a factor of p). From an exponential fit (line in Fig. 2c), A ¼ 0.110(7) mT ms À 1 , and using a background scattering length of a bg ¼ À 187(5)a 0 (ref. 25), where a 0 is the Bohr radius, we extract a width of D d j j ¼ 9:3ð7ÞÂ10 À 4 mT. By fitting an error function (line) to the data in Fig. 2d, we determine the location of the resonance to be 54.747(1) mT, which is consistent with previous experiments where atom loss was observed 14,16 . ARTICLE The precise determination of the width of the d-wave resonance allows us to gauge its significance in molecule creation. Our typical sweep rate of 0.34 mT ms À 1 for magneto-association, which has remained the same since the first creation of KRb molecules in an optical lattice 22 , gives B70% probability of being diabatic when crossing the d-wave resonance. This suggests that we create a substantial fraction of d-wave molecules that are dark to our detection ( 5 in Fig. 2b). These d-wave molecules may have played a role in limiting the lattice filling fraction for polar molecules achieved in ref. 11.
Short-time tunnelling dynamics. Tunnelling dynamics of doublons in the lattice 26 can also affect molecule production. In the quantum synthesis approach, achieving a high lattice filling for molecules requires not only the preparation of a large fraction of lattice sites that have doublons, but also that these doublons are not lost due to tunnelling and/or collisions prior to conversion to molecules. In our system, K feels a lattice depth that, in units of recoil energy, is 2.6 times weaker than for Rb due to differences in atomic mass and polarizability. Consequently, K tunnels faster than Rb. While a sufficiently deep lattice can prevent tunnelling of both K and Rb, practically this may not be possible in all cases, especially for polar molecule production using two atomic species that have large differences in mass and polarizability. Figure 3 illustrates doublon dynamics due to the interplay between tunnelling and interactions, which we control by varying the lattice depth, interspecies scattering length a K-Rb and band population. The fraction of doublons that remain after t is essentially equal to the measured molecule conversion efficiency described above. We note that for a K-Rb 4 À 850a 0 , the B sweep crosses the d-wave Feshbach resonance with a _ B that varies from 0.5 to 1.9 mT ms À 1 . Using our measured width of the d-wave resonance, the data presented in Fig. 3 have been multiplied by a factor that increases the doublon fraction to account for the finite _ B when crossing the d-wave resonance. Figure 3a shows the effect of the lattice depth for t ¼ 1 ms at three different values of B hold , corresponding to different values of a K-Rb . This timescale is relevant for both molecule production and K tunnelling dynamics. We observe that the remaining doublon fraction is highly sensitive to the lattice depth for weak interspecies interactions, for example, a K-Rb ¼ À 220a 0 , with a lower doublon fraction for shallower lattices that exhibit higher tunnelling rates. For stronger interactions, the dependence on lattice depth becomes less significant and almost disappears in the strongly interacting regime, for example, a K-Rb ¼ À 1,900a 0 . Similar behaviour is observed if we fix the lattice depth but vary the interspecies interactions, as shown in Fig. 3b.
Modelling. The data in Fig. 3 clearly show evidence of decay of doublons due to tunnelling that is affected by both the lattice depth and interspecies interactions. We can model these doublon dynamics with the following Hamiltonian: The doublon fraction is plotted for three lattice depths as a function of the scattering length. (c) The doublon fraction for 1.68 mT ms À 1 sweeps, t ¼ 1 ms and a K-Rb ¼ À 220a 0 is shown as a function of the lattice depth for the case of higher excited-band fraction (squares) and lower excited-band fraction (circles). (d) Band-mapping images of the initial K gas are shown for the two different initial temperatures, where image i corresponds to the red circle data points and ii corresponds to the green square data points in (c). Each image is the average of three measurements. The colour bar indicates the optical depth (OD). Below the images, we display the OD for a horizontal trace through the image, with averaging from À :k to þ :k in the vertical direction. All error bars represent 1 À s standard error.
where Z ¼ 0 and 1 denote, respectively, the ground and the first excited lattice bands. The first and second terms are the kinetic energy of the K and Rb atoms, respectively. Here, a i a y i is the bosonic annihilation (creation) operator for a Rb atom at lattice site i in the lowest band, and c i;Z c y i;Z is the fermionic annihilation (creation) operator for a K atom at lattice site i and band Z. We use i; j h i to indicate nearest-neighbour hopping between sites i and j with matrix element J Z a with a ¼ K or Rb. The third term describes the inter-species on-site interactions with matrix element U Z KRb . The last term is the on-site intra-species interaction between ground-band Rb atoms with strength U 0 RbRb , with n 0 Rb;i as the occupation of site i. The tunnelling rates and interaction energies are calculated for a particular V latt and a K-Rb (ref. 27). For example, for V latt ¼ 10E R , J 0 K =h ¼ 386 Hz, J 0 Rb =h ¼ 38:9 Hz. The solid curves in Fig. 3a,b show the calculations based on the Hamiltonian given in equation (1), where we have neglected Rb tunnelling by setting J 0 Rb ¼ 0. We start with a single doublon, evolve the K for a hold time t, and then extract the doublon fraction from the probability that the K atom remains on the same site as the Rb atom. In this treatment, we ignore the role of the magnetic-field sweeps. Calculations for a single doublon (solid lines), where the initial decay scales as 1 À 12ðJ 0 K =U 0 KRb Þ 2 , agree well with the data, except at doublon fractions below B30%, where the disagreement arises from the finite probability in the experiment that a K atom finds a different Rb partner. Simulating a Gaussian distribution of doublons with 10% peak filling accounts for this effect (dashed lines) (see Methods). The good agreement of these calculations with the data shows that tunnelling of K, which is suppressed for deeper lattices, is the dominant mechanism for the reduction of the doublon fraction at short (B1 ms) times. The on-site interaction with Rb suppresses the K tunnelling when the interaction energy becomes larger than the width of the K Bloch band 28 .
When studying doublon dynamics measured for two different initial atom conditions, we find indirect evidence for excited-band molecules. Here, we compare results for our usual molecule preparation using atomic insulators to a case where we start with a hotter initial atom gas mixture at a temperature above that for the Rb Bose-Einstein condensation transition. Using a band-mapping technique, we measure the initial population of K in the ground and first excited band, as shown in Fig. 3di and ii (see Methods). We find that 11(2)% of the K atoms occupy the first excited band for the colder initial atom gas (these conditions are similar to those in ref. 11 and are used in all the measurements described in this work, except for the green squares in Fig. 3c). When starting with the hotter atom gas, we measure a significantly higher K excited-band population of 31(6)%. When looking at doublon dynamics for these two cases (Fig. 3c), we observe a lower doublon fraction for the hotter initial gas for V latt r25E R . These data are taken for 1.68 mT ms À 1 sweeps, t ¼ 1 ms and a K-Rb ¼ À 220a 0 .
The lower doublon fraction can be explained by excited-band K atoms, which have a high tunnelling rate (J 0 K =h and J 1 K =h are 89.3 and 1110 Hz, respectively, for V latt ¼ 25E R ). The presence of excited-band K atoms suggests that the B sweeps for magnetoassociation (and dissociation) couple excited-band K atoms (plus a ground-band Rb atom) to excited-band Feshbach molecules. Moreover, the data suggest that the conversion efficiency for the excited-band Feshbach molecules is still high for V latt r25E R since the observed difference in the initial excited-band K atoms is similar to the observed difference (roughly 20%) in the doublon fraction (Fig. 3c). Since, in our preparation scheme, the doublons are directly formed from the dissociation of ro-vibrational ground-state molecules, these results further indicate that a polar molecule sample prepared from a finite-temperature atom gas can contain a small fraction of molecules in an excited motional state in the lattice. We also observed a Rb excited-band population of 31(5)% after loading the thermal gas in the lattice; however, even for the excited band, the off-resonant Rb tunnelling is slow compared to the 1-ms time scale of the measurements presented in Fig. 3.
The green dashed curve in Fig. 3c shows the theoretical results for a K excited-band fraction of 24%. For comparison, the red solid curve, which is the same as the red curve in Fig. 3a, includes no excited-band population. The estimated excited band fraction ignores the effects of harmonic confinement on tunnelling, which are more significant for the hotter initial atom gas, where the resulting molecular cloud is also larger. For the hotter initial atom gas, the green dashed curve overlaps the data at the shallower lattice depths, but deviates from the measured doublon fraction at larger lattice depths (the excited band fraction of the initial K gas is independent of lattice depth). This may be expected since in the limit of a very deep lattice and a fully adiabatic magneto-association sweep, one expects that only the heavier atom (Rb) in excited bands (plus a ground-band K atom) will couple to excited-band Feshbach molecules. Future studies of the magneto-association process in a lattice for systems such as K-Rb where centre-of-mass and relative motion are coupled 29 would be interesting and relevant to polar molecule preparation.
Long-time tunnelling dynamics. In Fig. 4, we present data taken for t up to 40 ms, in order to look for the effects of Rb tunnelling. Measurements of the remaining doublon fraction are shown for two lattice depths (10E R and 15E R ) and two values of a K-Rb ( À 910a 0 and À 1900a 0 ). In Fig. 4, the doublon fraction has been normalized by the measured value for t ¼ 1 ms in order to remove the effect of the shorter-time dynamics that are presented in Fig. 3a,b. Similar to the shorter-time dynamics, at the longer hold times we observe a reduction in the doublon fraction that is suppressed for a deeper lattice and for strong inter-species interactions. Modeling these dynamics is theoretically challenging, and the lines in Fig. 4 are exponential fits that are intended only as guides to the eye. Compared to doublons composed of identical bosons 13 or fermions in two-spin states 30 , the heteronuclear system has the additional complexities of two and 15E R (blue diamonds and triangles) for either a K-Rb ¼ À 1,900a 0 (circles, diamonds) or a K-Rb ¼ À 1,900a 0 (squares, triangles). The lines are fits to an exponential decay and are intended only as guides to the eye. (Inset) The doublon decay can involve tunnelling of doublons through empty sites (i) prior to loss by the Rb tunnelling process illustrated in (ii). All error bars represent 1 À s standard error.
particle masses, two tunnelling rates and two relevant interaction energies. For example, for large a K-Rb , the interspecies interactions will strongly suppress Rb tunnelling from a doublon to a neighbouring empty site. Similarly, tunnelling of a doublon to an empty lattice site is a slow second-order process at the rate J pair ¼ 2J 0 Rb J 0 K =U 0 KRb due to the energy gap of U 0 KRb (Fig. 4  inset i). However, Rb tunnelling between two neighbouring doublons, which creates a triplon (Rb-Rb-K) on one site and a lone K atom on the other site (Fig. 4 inset ii), may occur on a faster time scale due to a much smaller energy gap of U 0 RbRb , which is smaller than the K tunnelling bandwidth. While the theoretical description is complicated, we observe that the time scale of the doublon decay roughly matches 1= 2pJ pair À Á .
Measuring atomic distributions with doublon detection. The studies discussed thus far demonstrate that the Feshbach molecule conversion that we use to detect doublons could potentially underestimate the doublon fraction. For example, the conversion efficiency of doublons containing excited-band atoms is complicated to calculate and is likely to be less than 1. In addition, the efficiency of converting doublons to Feshbach molecules depends on the magnetic-field sweep rate, and, as shown in Fig. 2c, a very slow sweep does not always yield a unity conversion efficiency. Finally, Feshbach molecules can suffer losses from inelastic collisions with other Feshbach molecules or unpaired K atoms 31 , which could reduce the measured number. Given these factors and the importance of measuring the doublon fraction as a powerful diagnostic for optimizing molecule production from ultracold atoms in a lattice, we have implemented a second, complementary approach for measuring the doublon fraction using inelastic collisional loss in the initial atomic mixture, without the molecular purification step. In our system, inelastic collisions are initiated by transferring the Rb atoms from the 1; 1 j i to the 2; 2 j i state. Collisions of the 2; 2 j i Rb atoms with K can result in spin relaxation back to the Rb F ¼ 1 manifold. At B ¼ 55 mT, the 2; 2 j i state is higher in energy by h Â 8.1 GHz; this inelastic collision releases a large amount of energy compared to the trap depth and therefore results in atom loss from the trap. At a collision energy corresponding to 1 mK, the calculated inelastic collision rate using the coupled channels model of ref. 15 is b ¼ 6Â10 À 12 cm 3 s À 1 , and using the on-site densities in a V latt ¼ 25E R lattice, the resulting doublon lifetime is B2 ms.
In Fig. 5a,b, we show example data for the number of Rb atoms as a function of time after a 2.1-ms rf sweep that transfers Rb atoms to the 2; 2 j i state. We observe a fast loss on the time scale of a few ms, followed by slower loss. We attribute the fast loss to inelastic collisions of Rb atoms in lattice sites shared with K, and the slow loss to tunnelling of atoms followed by inelastic collisions. The dashed lines in Fig. 5a,b show a fit to the sum of two exponential decays with different time constants. We can extract the fraction of Rb that is lost on the short timescale from the fits. We compared this technique with Feshbach molecule formation, and found that the two measurements generally agree.
As a further demonstration of the inelastic collision technique, we use this to probe the initial atomic distribution in the lattice before molecule formation, providing quantitative information on the Rb Mott insulator. Figure 5c shows the fraction of Rb atoms that are lost quickly from a V latt ¼ 25E R lattice after the Rb atoms are transferred to the 2; 2 j i state. For these data, we vary the initial number of Rb atoms that form a Mott insulator in the optical lattice prior to molecule creation. In Fig. 5c, the blue diamonds correspond to the data shown in Fig. 5a,b. For the data shown in circles, the fraction lost is determined by comparing the Rb number measured before to that measured 8 ms after the rf transfer. The solid curve shows a calculation of the expected loss for a Mott insulator with a temperature T=J 0 Rb ¼ 15 and a total radial harmonic confinement of 33 Hz for Rb. At low Rb number, where we expect only one Rb atom per site in the Mott insulator, the fraction lost is just the K filling fraction, assuming no double occupancy for K. For higher Rb number, double (and eventually triple and higher) occupancy in Mott shells causes a reduction in the fractional loss under the assumption of one Rb and one K lost per inelastic collision. The shaded area indicates a 10% uncertainty in the harmonic trapping frequency and 30% uncertainty in T. From the fit, we extract a K filling fraction of 0.77 (2), which is in excellent agreement with the measured peak K filling reported in ref. 11. We note the previously measured fraction of Rb converted to Feshbach molecules at low Rb atom number was significantly less, at about 0.  panels (a,b). At low Rb number, where the Mott insulator has one Rb atom per site, the fraction lost should be equal to the fraction of sites that have a K atom. As the Rb filling increases and the second Mott shell becomes populated, the fraction lost decreases. This technique yields both the filling of K and a measure of Rb atom number that corresponds to the onset of double occupancy of the Rb Mott insulator. All error bars represent 1 À s standard error. disagreement may now be attributed to the large number of K atoms present in the lattice after molecule formation, which can induce losses through inelastic collisions, and to the effect of the d-wave resonance when making molecules, as discussed above.

Discussion
Our investigation of heteronuclear doublons and their conversion to molecules by magneto-association reveals the important roles played by the lattice depth for both atomic species, the interspecies interactions, the population in excited motional states of the lattice and the magnetic-field sweep rate. The doublon dynamics uncovered in this study provides insights into the universal mechanism of their decay and atomic mixture dynamics in a 3D optical lattice, and allows preparation of optimal conditions for producing polar molecules. The highly nonequilibrium state of doublons that we use for these studies also provides an intriguing system for exploring the Hubbard dynamics of a Bose-Fermi mixture, where the behaviour of the many-body system can depend on two different tunnelling rates and two different interaction strengths 28,32 . This system sets the stage for performing benchmarking experiments in 1D for theory, and investigating the thermalization of an isolated many-body quantum system, including novel phases such as quasicrystallization and many-body localization in higher dimensions [33][34][35][36] .

Methods
Optical trapping potentials. The preparation of the atomic gas in a 3D lattice, with a wavelength of 1,064 nm, as well as the creation of ground-state polar molecules, follows the procedures described in ref. 11. The lattice is superimposed on a crossed-beam optical dipole trap that is cylindrically symmetric. The dipole trap alone has an axial trap frequency of o z ¼ 2pÂ180 Hz in the vertical direction and a radial trap frequency of o r ¼ 2pÂ25 Hz for Rb. The measured optical trap frequencies for K are 2pÂ260 Hz and 2pÂ30 Hz.
Width of the d-wave resonance. The scattering lengths reported in Fig. 3 have been calculated using a KÀRb ðBÞ ¼ a bg ½1 À D s =ðB À B res s Þ with a bg ¼ À 187 (5) We note that ref. 17 predicts D 0 d ¼ À 6:3Â10 À 4 mT, which is larger than our determination of D 0 d ¼ À 2:0ð2ÞÂ10 À 4 mT. In this determination, we ignore the coupling between the centre of mass and relative motion that arises from the fact that K and Rb experience different trapping potentials in the optical lattice. We use an effective trap frequency o HO ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi that governs the dynamics in the relative coordinate. Here, m Rb and m K are the masses of the Rb and K atom, respectively. The trap frequency for Rb is given by o Rb ¼ 2ðE R =' Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi V latt =E R p , and for the 1,064 nm lattice the trap frequency for K is o K % 1:4o Rb . For V latt ¼ 35E R , o HO ¼ 30:4 kHz.
Density distribution of doublons. The dashed lines in Fig. 3a,b have been obtained by random sampling of initial doublon positions according to a Gaussian probability distribution of the filling fraction. A peak filling of 10% and widths s x ¼ s y ¼ 6.5s z ¼ 21 sites have been used, corresponding to NE2,000 sites occupied with a doublon. The experimentally determined cloud sizes are slightly larger (s x ¼ 25-42 sites), but we confirmed that the resulting doublon fraction is converged with respect to the cloud size. In-situ absorption images of the cloud are consistent with a Gaussian distribution of 5-10% peak filling. Tunnelling of Rb is neglected in the model, where initially each K atom is localized on a site containing a Rb atom and the doublon fraction is defined as the probability to find the K atom on a site with Rb after the evolution time t.
Band mapping. To measure the excited-band fraction of the initial K atoms, we use a band-mapping technique (Fig. 3d). Starting with the K atoms in the 3D lattice plus optical dipole trap potential, we turn off the lattice in 1 ms and allow the the K gas to expand in the optical dipole trap for a quarter trap period 38 . We image the cloud with a probe beam that propagates along the vertical direction.