Quantum gas microscopy of an attractive Fermi-Hubbard system

TL;DR

This study uses quantum gas microscopy to observe charge-density-wave correlations and superfluidity in an attractive Fermi-Hubbard system, providing insights into pairing phenomena relevant to superconductivity.

Contribution

It presents the first site-resolved measurements of superfluid correlations and charge-density-wave order in an ultracold atomic Fermi gas simulating the attractive Hubbard model.

Findings

Observation of checkerboard charge-density-wave correlations near half-filling.

First evidence of superfluid correlations in a single-band Hubbard system of ultracold fermions.

Charge-density-wave correlations serve as a sensitive thermometer for cooling into superfluid phases.

Abstract

The attractive Fermi-Hubbard model is the simplest theoretical model for studying pairing and superconductivity of fermions on a lattice. Although its s-wave pairing symmetry excludes it as a microscopic model for high-temperature superconductivity, it exhibits much of the relevant phenomenology, including a short-coherence length at intermediate coupling and a pseudogap regime with anomalous properties. Here we study an experimental realization of this model using a two-dimensional (2D) atomic Fermi gas in an optical lattice. Our site-resolved measurements on the normal state reveal checkerboard charge-density-wave correlations close to half-filling. A "hidden" SU(2) pseudo-spin symmetry of the Hubbard model at half-filling guarantees superfluid correlations in our system, the first evidence for such correlations in a single-band Hubbard system of ultracold fermions. Compared to the paired atom fraction, we find the charge-density-wave correlations to be a much more sensitive thermometer, useful for optimizing cooling into superfluid phases in future experiments.

Figures

b,c to $U/t=-5.4$ and $T=0.4t$. Both this value of $U/t$ and the one obtained without corrections to the Hubbard model fall within our experimental uncertainties.

b,c to $U/t=-5.4$ and $T=0.4t$. Both this value of $U/t$ and the one obtained without corrections to the Hubbard model fall within our experimental uncertainties.

Equations

$$\mathcal{H}=-t\sum _{⟨{\text{rr}}^{\prime }⟩,\sigma }\left(c_{\text{r},\sigma }^{†}{c}_{{\text{r}}^{\prime },\sigma }+c_{{\text{r}}^{\prime },\sigma }^{†}{c}_{\text{r},\sigma }\right)+U\sum _{\text{r}}{n}_{\text{r},\uparrow }{n}_{\text{r},\downarrow }-\mu \sum _{\text{r},\sigma }{n}_{\text{r},\sigma }-h\sum _{\text{r}}(n_{\text{r},\uparrow }-n_{\text{r},\downarrow })$$

$$S_{\mathbf{r}}^{-}$$

$$S_{\mathbf{r}}^{+}=(S_{\mathbf{r}}^{-}{)}^{†}$$

$$S_{\mathbf{r}}^z=\frac{1}{2}(c_{\mathbf{r},\uparrow }^{†}{c}_{\mathbf{r},\uparrow }-c_{\mathbf{r},\downarrow }^{†}{c}_{\mathbf{r},\downarrow })$$

$$[S_{\mathbf{r}}^z,S_{\mathbf{r}}^{\pm }]=\pm S_{\mathbf{r}}^{\pm },[S_{\mathbf{r}}^{+},S_{\mathbf{r}}^{-}]=2S_{\mathbf{r}}^z$$

$$[\mathcal{H},S^{\pm }]=[\mathcal{H},S^z]=0$$

$${\eta }_{\text{r}}^{-}$$

$${\eta }_{\text{r}}^{+}$$

$${\eta }_{\text{r}}^z$$

$$[{\eta }_{\text{r}}^z,{\eta }_{\text{r}}^{\pm }]=\pm {\eta }_{\text{r}}^{\pm },[{\eta }_{\text{r}}^{+},{\eta }_{\text{r}}^{-}]=2{\eta }_{\text{r}}^z.$$

$$[\mathcal{H},{\eta }^{\pm }]=\pm (U-2\mu ){\eta }^{\pm },[\mathcal{H},{\eta }^z]=0.$$

$${⟨{\eta }_{\mathbf{r}}^z{\eta }_{\mathbf{r}+\mathbf{a}}^z⟩}_c={⟨{\eta }_{\mathbf{r}}^x{\eta }_{\mathbf{r}+\mathbf{a}}^x⟩}_c$$

$${\Delta }_{\text{r}}^x=c_{\text{r},\uparrow }{c}_{\text{r},\downarrow }+c_{\text{r},\downarrow }^{†}{c}_{\text{r},\uparrow }^{†}$$

$$c_{\mathbf{r},\downarrow }\leftrightarrow (-1{)}^{r_x+r_y}{c}_{\mathbf{r},\downarrow }^{†},c_{\mathbf{r},\uparrow }\leftrightarrow c_{\mathbf{r},\uparrow }.$$

$$\mathcal{H}\leftrightarrow -t\sum _{⟨{\text{rr}}^{\prime }⟩,\sigma }\left(c_{\text{r},\sigma }^{†}{c}_{{\text{r}}^{\prime },\sigma }+c_{{\text{r}}^{\prime },\sigma }^{†}{c}_{\text{r},\sigma }\right)-U\sum _{\text{r}}{n}_{\text{r},\uparrow }{n}_{\text{r},\downarrow }-{\mu }^{\prime }\sum _{\text{r},\sigma }{n}_{\text{r},\sigma }-h^{\prime }\sum _{\text{r}}(n_{\text{r},\uparrow }-n_{\text{r},\downarrow })$$

$${\eta }^{\pm }\leftrightarrow S^{\pm },{\eta }^z\leftrightarrow S^z$$

$$t_n=-\frac{4\pi {\hslash }^2{a}_s}{m}\int d^{3}r{w}^{\ast }(\text{r}-\text{d})w^{\ast }(\text{r})w^2(\text{r})$$

Full text

The Fermi-Hubbard model is a fundamental condensed matter model for studying strongly-correlated fermions on a lattice Hubbard (1963); Auerbach . The numerical intractability of the model at low temperatures has spurred much experimental work on quantum simulations using repulsively interacting cold atoms in optical lattices Jördens et al. (2008); Schneider et al. (2008); Greif et al. (2013); Hart et al. (2015); Parsons et al. (2016); Cheuk et al. (2016); Boll et al. (2016); Brown et al. (2016); Cocchi et al. (2016). In contrast, the attractive model has received much less attention from the cold atom community Strohmaier et al. (2007); Hackermüller et al. (2010); Schneider et al. (2012), despite being an ideal playground for exploring short coherence-length superconductivity on a lattice. In the continuum, strongly-interacting attractive fermions have been studied using Feshbach resonances, resulting in the observation of superfluid gases across the crossover from molecular Bose-Einstein condensates (BEC) to Bardeen-Cooper-Schrieffer superfluids (BCS) Ing (2008). Fermionic superfluids have also been prepared in optical lattices close to a Feshbach resonance Chin et al. (2006). However, these systems are not described by a simple Hubbard model due to multi-band couplings and off-site interactions Duan (2005); Carr and Holland (2005); Zhou (2005); Diener and Ho (2006).

In this work, we focus on the 2D attractive Fermi-Hubbard model which has been theoretically studied in detail Hirsch (1985); Scalettar et al. (1989); Micnas et al. (1990); Moreo and Scalapino (1991); Paiva et al. (2004). Our experiments are performed at an interaction energy small compared to the bandgap, where the single-band Hubbard description is applicable. In a grand-canonical ensemble, the Hamiltonian of the system is given by ${\cal{H}}=-t\sum_{\langle\textbf{rr}^{\prime}\rangle,\sigma}\left(c^{{\dagger}}_{\textbf{r},\sigma}c_{\textbf{r}^{\prime},\sigma}+c^{{\dagger}}_{\textbf{r}^{\prime},\sigma}c_{\textbf{r},\sigma}\right)+U\sum_{\textbf{r}}n_{\textbf{r},\uparrow}n_{\textbf{r},\downarrow}-\mu\sum_{\textbf{r},\sigma}n_{\textbf{r},\sigma}$. Here $c^{{\dagger}}_{\textbf{r},\sigma}$ is the creation operator for a fermion with spin $\sigma$ on site r, $n_{\textbf{r},\sigma}=c^{{\dagger}}_{\textbf{r},\sigma}c_{\textbf{r}\sigma}$, $t$ is the tunneling matrix element between nearest-neighbor lattice sites, $U<0$ is the strength of the on-site interaction and $\mu$ is a spin-independent chemical potential. At low temperatures, the fermions undergo a Berezinskii-Kosterlitz-Thouless (BKT) transition to a superfluid phase. As $U/t$ is increased, the superfluid crosses over from a BCS-type superfluid to BEC of hardcore bosons, with the critical temperature reaching a maximum in the intermediate coupling regime. Near this maximum and in the BEC regime, numerical calculations indicate a clear separation between a temperature scale $T^{*}$ at which pairing correlations appear and the BKT transition temperature $T_{c}$ Scalettar et al. (1989).

The phase diagram of the attractive Hubbard model versus filling is shown in Fig. 1a. At half-filling, the model has a “hidden” SU(2) pseudo-spin symmetry, in addition to the ordinary SU(2) spin symmetry Yang and Zhang (1990). Pseudo-spin rotations on site $\textbf{r}=(r_{x},r_{y})$ are generated by the charge-density fluctuation operator $\eta^{z}_{\textbf{r}}=\frac{1}{2}(n_{\textbf{r}}-1)$ and the pairing operators $\eta_{\textbf{r}}^{-}=(-1)^{r_{x}+r_{y}}c_{\textbf{r},\uparrow}c_{\textbf{r},\downarrow}$ and $\eta^{+}_{\textbf{r}}=(\eta_{\textbf{r}}^{-})^{\dagger}$. The local pseudo-spin vector is defined by $(\eta^{x}_{\textbf{r}},\eta^{y}_{\textbf{r}},\eta^{z}_{\textbf{r}})$, where $\eta^{\pm}_{\textbf{r}}=\eta^{x}_{\textbf{r}}\pm i\eta^{y}_{\textbf{r}}$. The Hubbard Hamiltonian is rotationally invariant under global pseudo-spin rotation at half-filling, leading to checkerboard charge-density-wave (CDW) and superfluid correlations of equal strength. The SU(2) symmetry drives the critical temperature to zero. Away from half-filling, the degeneracy is lifted, and superfluid correlations exhibit quasi-long-range order below a non-zero BKT critical temperature while the charge-density-wave correlations remain short-range.

Here we perform a site-resolved study of the attractive Hubbard model with a fermionic quantum gas microscope Haller et al. (2015); Edge et al. (2015); Omran et al. (2015); Parsons et al. (2015); Cheuk et al. (2015); Yamamoto et al. (2016); Brown et al. (2016) at temperatures above the BKT transition temperature. We measure the in-trap distributions of single and double occupancies and observe charge density wave (CDW) correlations near half-filling. This observation allows us to put a lower bound on superfluid correlations in the system. Furthermore, the CDW correlations enable us to perform accurate thermometry on the attractive lattice gas at the superexchange scale. Thermometry at this scale will be important for future work aimed at observing lattice superfluid phases, including homogeneous superfluids in spin-balanced gases and superfluids with spatially-modulated gaps in spin-imbalanced gases Fulde and Ferrell (1964).

We realize the 2D Fermi-Hubbard model using a degenerate mixture of two hyperfine ground states, $\left\lvert\uparrow\right\rangle\equiv\left\lvert 1\right\rangle$ and $\left\lvert\downarrow\right\rangle\equiv\left\lvert 3\right\rangle$, of 6Li in an optical lattice, where $\left\lvert k\right\rangle$ labels the $k$th lowest hyperfine ground state. A spin-balanced mixture is obtained by optical evaporation in the vicinity of a broad Feshbach resonance centered at 690 G. After the evaporation, the scattering length is set to $-889~{}a_{0}$, where $a_{0}$ is the Bohr radius, obtained by adjusting a bias magnetic field to 305.4(1) G. The mixture is prepared in a single layer of an accordion lattice (for details see Mitra et al. (2016)) and subsequently loaded adiabatically into a 2D square lattice of depth $6.2(2)~{}E_{R}$ , where $t=h\times 1150(50)$ Hz. The lattice depth is chosen experimentally to maximize the observed CDW correlations at half filling (see Supplementary Information). Here, $E_{R}=(\pi\hbar)^{2}/2ma_{latt}^{2}=h\times 14.66$ kHz is the recoil energy, where $a_{latt}=1064~{}\text{nm}/\sqrt{2}$ is the lattice constant. For these parameters, we obtain $U/t=-5.4(3)$ from a bandstructure calculation and a spectroscopic determination of the interaction energy.

We extract the density profiles and correlations in the cloud from site-resolved fluorescence images obtained using quantum gas microscopy techniques (Fig. 1b). After freezing the density distribution in a deep lattice, we shine near-resonant light on the atoms in a Raman-cooling configuration. Light-assisted collisions eliminate atoms on doubly occupied sites, and we measure the singles density $n^{s}=n_{\uparrow}+n_{\downarrow}-2n_{\uparrow}n_{\downarrow}$. To gain the full density information, one needs to measure the doubly occupied sites as well. But they are not directly accessible, so we developed a procedure to selectively image doubly occupied sites (Fig. 1c). After freezing the dynamics, we tag doubly occupied sites using a radiofrequency pulse to transfer atoms in state $\left\lvert 1\right\rangle$ to state $\left\lvert 2\right\rangle$ only on these sites, relying on the interaction energy for spectroscopic addressing. We then push atoms in states $\left\lvert 1\right\rangle$ and $\left\lvert 3\right\rangle$ out of the lattice with a resonant light pulse. We avoid the loss of state $\left\lvert 2\right\rangle$ atoms during the resonant light pulse due to light-assisted collisions by an appropriate choice of a positive scattering length for this step. This procedure gives us access to the doublon density $n^{d}=n_{\uparrow}n_{\downarrow}$. We have measured our doublon detection efficiency to be $0.91(1)$ by analyzing the imaging fidelity of band insulating regions in the cloud (Fig. 1d, see Supplementary Information). The average singles and doubles density profiles obtained from clouds prepared under identical conditions allow us to extract the total density profile $n=n^{s}+2n^{d}$. For most of our measurements, we adjust the atom number to obtain a total density slightly above $n=1$ at the center of the trap to obtain a large region in the cloud near half filling.

Spatial correlations of the $z$-component of the pseudo-spin near half-filling correspond to CDW correlations, characterized by $\langle\eta^{z}_{\textbf{r}}\eta^{z}_{\textbf{r}+\textbf{a}}\rangle_{c}=\langle\eta^{z}_{\textbf{r}}\eta^{z}_{\textbf{r}+\textbf{a}}\rangle-\langle\eta^{z}_{\textbf{r}}\rangle\langle\eta^{z}_{\textbf{r}+\textbf{a}}\rangle=\frac{1}{4}\langle n_{\textbf{r}}n_{\textbf{r}+\textbf{a}}\rangle_{c}$ where a is the displacement vector between two lattice sites. We extract a closely related quantity, the doublon-doublon correlator $C^{d}(\textbf{a})=4\langle n^{d}_{\textbf{r}}n^{d}_{\textbf{r}+\textbf{a}}\rangle_{c}$, which also exhibits checkerboard order since most of the atoms in the gas are paired. As the attractive gas is compressible for any filling below unit filling, the local density varies across the harmonic trap. Fig. 2a shows the doublon-doublon correlator versus density for the nearest-neighbor and the next-nearest-neighbor, obtained by azimuthally averaging the correlations over the trap. The corresponding trap density profile of the gas is depicted in Fig. 2b. The nearest-neighbor doublon-doublon correlator measured at half-filling is $C^{d}(1,0)=-0.155(6).$ Since $|C^{d}(\textbf{a})|\leq\langle\Delta^{x}_{\textbf{r}}\Delta^{x}_{\textbf{r}+\textbf{a}}\rangle_{c}$ for any filling (see Supplementary Information), our measurements constitute a lower bound for superfluid correlations in the system. Here, the $s$-wave superfluid gap operator on site r is defined by $\Delta^{x}_{\textbf{r}}=c^{\dagger}_{\textbf{r},\downarrow}c^{\dagger}_{\textbf{r},\uparrow}+c_{\textbf{r},\uparrow}c_{\textbf{r},\downarrow}$.

An interesting feature of the next-nearest-neighbor correlator $C^{d}(1,1)$ is that it becomes negative as the average density falls below $\sim 0.4$. This can be understood in the large $U$ limit, where the system can be treated as a gas of hardcore bosons with repulsive nearest-neighbor interactions, leading to negative correlations at distances less than the interparticle spacing. In a recent experiment, a closely related behavior has been observed in the next-nearest-neighbor antiferromagnetic correlations of the $z$-projection of the spin in a spin-imbalanced repulsive Hubbard model Brown et al. (2016). A particle-hole transformation on operators of $\downarrow$ fermions in the Hubbard Hamiltonian, $c_{\textbf{r},\downarrow}\leftrightarrow(-1)^{r_{x}+r_{y}}c^{\dagger}_{\textbf{r},\downarrow}$, changes the sign of the interaction and exchanges the roles of spin-imbalance and doping Ho et al. (2009), leading to this symmetry of the correlators between the two experiments.

We theoretically model our system using determinantal Quantum Monte Carlo (DQMC) Varney et al. (2009) in a local density approximation (LDA) and see very good agreement between theory and experiment for the doublon-doublon correlators and density (Fig. 2). The fits give $U/t=-5.7(2)$ and $T/t=0.45(3)$. The measured temperature is comparable to our recent measurements in a repulsive gas Brown et al. (2016). We have also compared to DQMC in the presence of a spatially varying potential to reproduce the largest experimentally observed density gradients to verify that the LDA holds in our system.

Single-site imaging of doublons allows us to also measure longer range correlations (Fig. 3). We see doublon-doublon correlations up to two sites on a diagonal shown in the correlation matrices $C^{d}(i,j)$. We find good agreement with DQMC calculations corresponding to the experimental fillings, calculated using the same parameters as above. The range of the correlation becomes maximal at half filling as expected from theory.

The previous discussion focused on the upper branch of the Feshbach resonance, used for repulsive Hubbard experiments with lithium Hart et al. (2015); Parsons et al. (2016); Boll et al. (2016); Brown et al. (2016). Since we are studying an attractive model, we also have the choice of using the lower branch of the resonance, previously used in continuum BEC-BCS crossover experiments Ing (2008). It is not clear a priori how the temperature of gases prepared on the two branches would compare. We have explored this question by measuring doublon-doublon correlations in lower-branch gases (Fig. 3c). The least negative scattering length we can access on this branch is $-2800a_{0}$, so we use a lower lattice depth of 4.1(1) $E_{R}$ to reach the same value of $U/t$ as on the upper branch. We obtain nearest-neighbor doublon-doublon correlations of -0.172(5) at half filling, suggesting that we reach similar temperatures for this branch. The correlation matrices for the upper and lower branch do not agree very well, especially for larger site separations, which is expected since the on-site interaction energy is comparable to the vertical lattice confinement, leading to higher-band effects which modify the Hamiltonian. A precise determination of the temperature on the upper branch requires taking these effects into account, an endeavor that is outside the scope of this work.

Both pairing and CDW correlations may be used as thermometers in an attractive Hubbard gas. To investigate their sensitivity as thermometers, we heated the upper branch system in a controlled fashion by holding it for variable times in the lattice before imaging. Technical noise leads to a linear increase in the temperature of the system. We observe a slight reduction in doublon fraction in a region of half filling for long hold times (Fig. 4a), while the doublon-doublon correlators $C^{d}(1,0)$ and $C^{d}(1,1)$ show a significant change during the same time. These observations illustrate that the doublon fraction is a good thermometer for temperatures on the order of $U$, where CDW correlations are small, while the doublon-doublon correlator is a more sensitive thermometer for temperatures on the order of the exchange $4t^{2}/U$. Fig. 4b,c shows the nearest-neighbor and next-nearest-neighbor doublon-doublon correlators versus the singles fraction at half-filling. Plotting the data this way allows a temperature-independent comparison to DQMC with a single free-parameter ($U/t$), and we see good agreement for $U/t=-5.7$.

In conclusion, we have performed quantum gas microscopy on an attractive atomic Fermi-Hubbard system. We observed CDW correlations and studied their dependence on the lattice filling and temperature. We have shown that these correlations serve as an excellent thermometer in the low temperature regime that will be useful in the quest to reduce the temperature of attractive lattice systems to observe long-range Hubbard superfluids. Above the critical temperature, attractive Hubbard systems will allow experimental exploration of the physics of the pseudogap regime, while at lower temperatures, they will enable studies of the BKT transition in a lattice system. The superfluid correlations inferred in this experiment may be probed more directly using collective excitations of the superfluid order parameter that couple to the conjugate CDW order parameter Demler et al. (1996). Another interesting direction for future work is the study of the spin-imbalanced attractive Hubbard model in 2D, where Fermi surface nesting due to the lattice increases the stability of Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superfluids Moreo and Scalapino (2007); Loh and Trivedi (2010); Kim and Törmä (2012); Gukelberger et al. (2016). Unlike spin-balanced gases investigated in this work, DQMC suffers a fermion “sign” problem in simulating the spin-imbalanced attractive Hubbard model, making theoretical predictions difficult at low temperatures.

I Methods

Preparation of an attractive Fermi gas in an optical lattice The experimental setup is described in detail in the supplement of ref. Brown et al. (2016). We realize the Fermi-Hubbard model using a spin-balanced degenerate mixture of two Zeeman states ($\left\lvert 1\right\rangle=\left\lvert\uparrow\right\rangle$ and $\left\lvert 3\right\rangle=\left\lvert\downarrow\right\rangle$, numbered up from the lowest energy) in the ground state hyperfine manifold of 6Li in an optical lattice.

To create the sample we load a magneto-optical trap (MOT) from a Zeeman slower, then use a compressed MOT stage before loading into a $\approx$ $1\text{\,}\mathrm{mK}$ deep optical trap and evaporating near the $690\text{\,}\mathrm{G}$ Feshbach resonance. For preparing a gas on the upper branch of the resonance, we stop the evaporation before Feshbach molecules form and transfer the atoms to a highly anisotropic ‘light sheet’ trap where the gas undergoes further evaporation to degeneracy at $305.4(1)\text{\,}\mathrm{G}$. At this field, the scattering length is $a_{s}=-889a_{0}$. Next we transfer the gas into the final trapping geometry where a $1070\text{\,}\mathrm{nm}$ beam provides radial confinement and a $532\text{\,}\mathrm{nm}$ accordion lattice with trapping frequency $\omega_{z}=2\pi\times$ $21(1)\text{\,}\mathrm{kHz}$ provides axial confinement (for further details see Mitra et al. (2016)). The spin populations are balanced to within $2.1(9)\%$. We then load the gas into a 2D square lattice with a $25\text{\,}\mathrm{ms}$ long ramp to varying depths from $4-7.5~{}E_{r}$.

For experiments on the lower branch of the resonance, we perform the entire evaporation at the Feshbach resonance. Three-body collisions lead to population of the molecular branch. Before loading into the lattice, the bias field is ramped to 907 G where the scattering length is -2800 $a_{0}$.

Calibration of Hubbard parameters We use a non-separable 2D square lattice with a $752\text{\,}\mathrm{nm}$ spacing formed by four interfering passes of a single vertically polarized beam Brown et al. (2016). We calibrate our lattice depth by measuring the frequencies of the three $d$ bands in a deep lattice using lattice depth modulation, and compare these with a 2D band structure calculation. The inferred depth of the lattice at which our measurements are performed is 6.2(2)\text{\,}~{}$$E_{r}, where $E_{r}=$14.66\text{,}\mathrm{kHz}$$. From that we obtain nearest-neighbor tunneling values $t_{x}=1200$ Hz, $t_{y}=1110$ Hz ($t_{x}/t_{y}=1.08$). The reduction of the lattice depth across the cloud due to the gaussian profile of the lattice beams leads to an increase in the tunneling by $10\text{\,}\%$ at the edge of the cloud compared to the central value. We also calculate a non-zero but negligible diagonal tunneling $t_{11}=$42\text{,}\mathrm{H}\mathrm{z}$=0.04t_{x}$, due to the non-separability of the lattice.

We measure the interaction energy $U$ at the lattice loading field of 305 G using radio frequency spectroscopy. We transfer atoms from state $\left\lvert 1\right\rangle$ to $\left\lvert 2\right\rangle$ and resolve the frequency shift between singly and doubly occupied sites. We determine $U_{13}=\delta U\frac{a_{13}}{a_{13}-a_{23}}$, where $\delta U$ is the measured difference between the singles and doublon peaks and $a_{13}$ ($a_{23}$) is the scattering length at the spectroscopy field for a 1-3 (2-3) mixture. Taking into account a correction due to a significant final state interaction, we obtain $U_{13}=6.6(3)$ kHz. The experimentally measured value agrees with the value determined from band structure calculations of 5.9(1) kHz including higher band corrections to within 10% Idziaszek and Calarco (2005).

Imaging of doublons In the usual scheme of Raman imaging, a singly occupied lattice site produces a fluorescence signal while atoms in a doubly occupied site are lost due to light-assisted collisions. To measure density correlations between doubly-occupied sites, we have developed a new detection scheme. After the gas is loaded adiabatically into the optical lattice, we pin the atoms by increasing the lattice depth to 55(1) $E_{R}$ in 100 $\mu s$ where tunnelling dynamics get frozen. For gases on the upper branch, we adiabatically ramp the field to 594 G where we perform an interaction resolved Landau-Zener sweep (Fig. S1) to selectively transfer 1-3 doublons to 2-3 doublons while not effecting the $\left\lvert 1\right\rangle$ singles. Finally we ramp to a field of 641 G where the 2-3 scattering length is 414 $a_{0}$. We apply resonant pushing pulses of 30 $\mu s$ durations to remove $\left\lvert 1\right\rangle$ and $\left\lvert 3\right\rangle$ atoms, leaving behind with only $\left\lvert 2\right\rangle$ atoms on sites that originally had doublons. The large relative wavefunction of the atoms on a site due to repulsive interactions and a relatively weak vertical confinement significantly reduces the probability of light-assisted collisions during the resonant pulse. On the other hand, in order to get the singles images, we apply no RF or resonant pushing pulses.

For gases on the lower branch, after pinning the atoms, we ramp the field to 725 G where we perform a Landau-Zener sweep on an interaction-resolved transition to the 2-3 upper branch. In order to convert the 2-3 doubles to $\left\lvert 2\right\rangle$-singles with high efficiency, we ramp to 641 G where we apply the resonant pulses to remove $\left\lvert 1\right\rangle$ and $\left\lvert 3\right\rangle$ atoms. This ramp is done quickly (within $\sim$800 $\mu s$) to avoid atom losses due to the crossing of a narrow 2-3 Feshbach resonance around 714 G.

We image the final atom distribution by increasing the lattice depth to 2500 $E_{r}$ within $250~{}\mu$s, ramping up the light sheet to provide axial confinement, and then collecting fluorescence photons during Raman cooling.

To measure the efficiency of imaging doublons, we prepare a band insulator where the filling in the trap center is saturated at two atoms per site. We perform the above-mentioned process to image only the doublons. We measure a combined fidelity (including RF-transfer efficiency and pushing efficiency) of $1-\epsilon_{d}=90.9(5)\%$ of imaging doublons, leading to an underestimation of our doublon-doublon correlator by $(1-\epsilon_{d})^{2}=0.83$ which we correct for. We have also performed the same detection procedure on a Mott insulator, where we expect unit occupancy on the lattice sites. This allows us to extract the probability that a single atom in state $\left\lvert 1\right\rangle$ would get transferred to $\left\lvert 2\right\rangle$ and give a false positive signal of the presence of a doublon. We measure the probability of this process to be only 2.3(3)$\%$.

Raman imaging and reconstruction We perform Raman imaging for $1200\text{\,}\mathrm{ms}$ and collect approximately $1000$ photons per atom. For more details, see Brown et al. (2016). We estimate fidelity errors due to Raman imaging imperfections by taking 40 consecutive images of the same atom cloud and determine the shot-to-shot differences. This leads to a hopping rate during one picture of $0.4(1)\text{\,}\%$ and a loss rate of $1.9(2)\text{\,}\%$. In addition, while holding the atoms in a deep lattice for doublon detection, we lose $2(1)\text{\,}\%$ of the atoms, leading to a net detection efficiency of $\approx$95\text{,}%$$. The densities that we obtain are corrected for the above detection efficiency.

Lower bound on superfluid correlations The full Hubbard Hamiltonian can be written as

[TABLE]

where $c^{{\dagger}}_{\textbf{r},\sigma}$ is the creation operator for a fermion with spin $\sigma$ on site r, $n_{\textbf{r},\sigma}=c^{{\dagger}}_{\textbf{r},\sigma}c_{\textbf{r}\sigma}$, $t$ is the matrix element for tunneling between lattice sites, $U$ is the strength of the on-site interaction ($U<0$ for the attractive model), $\mu$ is a spin-independent chemical potential and $h$ is an effective Zeeman field in the presence of spin-imbalance. For the purpose of this paper, $h=0$, since we only work with a spin-balanced system. Consider first the SU(2) spin symmetry of the problem for $h=0$. The vector spin operator on site $\bf{r}$ is given by $(S^{x}_{\bf{r}},S^{y}_{\bf{r}},S^{z}_{\bf{r}})$. In terms of fermionic creation and annihilation operators, the generators of spin rotations are

[TABLE]

These operators obey the usual commutation relations

[TABLE]

In the absence of an effective Zeeman field, the operators $S^{z}=\sum_{r}S^{z}_{\mathbf{r}}$ and $S^{\pm}=\sum_{r}S^{\pm}_{\mathbf{r}}$ satisfy

[TABLE]

implying the SU(2) spin-symmetry of the Hubbard model. In other words, the Hamiltonian is invariant under a global rotation of the spin degree of freedom. To demonstrate the other “hidden” symmetry of the Hubbard model, we define a new set of generators for rotations of the “pseudo-spin” on a site given by

[TABLE]

It can be easily verified Zhang (1990) that they obey the following commutation relations

[TABLE]

The operators $\eta^{z}=\sum_{r}\eta^{z}_{\mathbf{r}}$ and $\eta^{\pm}=\sum_{r}\eta^{\pm}_{\mathbf{r}}$ satisfy

[TABLE]

regardless of the sign of the interaction $U$. Exactly at half-filling, $\mu=U/2$ and all the pseudospin generators in Eq. S7 commute with $\cal{H}$ similar to the spin generators in Eq. S4, meaning that the Hamiltonian is invariant under global rotations of the pseudo-spin.

The pseudo-spin on a site can be visualized on the Bloch sphere like a regular spin. If it points up (down) the site contains a doublon (hole), while if it lies in the equatorial plane, the site contains an equal superposition of a doublon and hole with a complex relative phase determined by the azimuthal angle. In the limit of large attractive interactions, the spin-balanced Hubbard Hamiltonian can be approximated as a Heisenberg Hamiltonian with antiferromagnetic interactions between the pseudo-spins, leading to charge-density-wave and superfluid correlations corresponding to $z$ and $x,y$ antiferromagnetic pseudo-spin correlations respectively. At half-filling, the pseudo-spin rotational symmetry implies

[TABLE]

where $\eta^{x}_{\textbf{r}}=(\eta^{+}_{\textbf{r}}+\eta^{-}_{\textbf{r}})/2$. Defining the s-wave pairing operator in the $x$ direction as

[TABLE]

one concludes that the density correlations $C^{n}(\textbf{a})=\left\langle n_{\textbf{r}}n_{\textbf{r+a}}\right\rangle_{c}=4\left\langle\eta^{z}_{\textbf{r}}\eta^{z}_{\textbf{r+a}}\right\rangle_{c}$ and pairing correlations $C^{\Delta}(\textbf{a})=\left\langle\Delta^{x}_{\textbf{r}}\Delta^{x}_{\textbf{r+a}}\right\rangle_{c}=4\left\langle\eta^{x}_{\textbf{r}}\eta^{x}_{\textbf{r+a}}\right\rangle_{c}$ are equal in magnitude at half-filling.

Deviation from half-filling introduces an effective Zeeman field along $z$ that couples to the pseudo-spin. This leads to canted antiferromagnetic pseudo-spin correlations, with stronger correlations in the direction orthogonal to the field, i.e. the superfluid correlations become stronger than the density correlations. Therefore, measurement of charge-density-wave density correlations in the attractive Hubbard model provides a lower bound on the superfluid correlations at any filling. In our experiment we measure a more accessible quantity, the doublon correlator $C^{d}(\textbf{a})=4\left\langle n^{d}_{\textbf{r}}n^{d}_{\textbf{r+a}}\right\rangle_{c}$. The doublon and density correlators become equal in the limit of low temperatures and large interactions. However, the doublon correlator $C^{d}$ still provides a bound on the superfluid correlator $C^{\Delta}$ as we have verified numerically for the entire range of interactions, temperatures and fillings studied in the experiment (Fig. S2).

Relationship to earlier work on repulsive spin-imbalanced Hubbard systems On a bipartite lattice, one can define a particle-hole transformation given by following mapping of the annihilation operators

[TABLE]

Under this transformation, $n_{\uparrow}\leftrightarrow n_{\uparrow}$ and $n_{\downarrow}\leftrightarrow 1-n_{\downarrow}$ and the Hubbard Hamiltonian (S1) transforms to (up to a constant energy offset)

[TABLE]

where, $\mu^{\prime}=h-U/2$ and $h^{\prime}=\mu-U/2$. In addition, the pseudo-spin rotation generators map on to the regular spin rotation generators under this transformation

[TABLE]

In other words, this transformation provides a mapping between the attractive and the repulsive Hubbard models, with the roles of doping and spin-imbalance interchanged Ho et al. (2009). The doped attractive Hubbard system studied in this work maps onto a spin-imbalanced repulsive Hubbard system. Therefore, the $z$-spin correlations measured in our previous work Brown et al. (2016) are closely related to the CDW correlations measured in this work.

Dependence of CDW correlations on the interactions We studied the dependence of the CDW correlations on $U/t$. The system was prepared on the upper branch at a scattering length of $-889~{}a_{0}$ and $U/t$ was varied by changing the lattice depth. Fig. S3 shows the nearest-neighbor and next-nearest-neighbor correlator vs. density for different lattice depths. The data is compared to DQMC results with $U/t$ and $T/t$ used as fit parameters. A slight increase in the temperature is observed for larger lattice depths. The measured correlations initially increase as the lattice depth is increased, but do not show significant variation for lattice depths between $6.2~{}E_{R}$ and $7.4~{}E_{R}$. At larger depths, it is experimentally difficult to obtain a half-filling condition near the center of the trap because of the increasing radial confinement from the lattice beams.

Corrections to the Hubbard model We investigated the importance of corrections to the Hubbard model in our analysis of the experimental data. For this purpose, two effects were taken into account:

• •

Effects of next-nearest-neighbor tunneling From a tight-binding calculation, we obtain for a lattice depth of $6.2~{}E_{R}$ (used for upper branch experiments) tunnel couplings $t_{11}/t\sim 4\%$ and $t_{20}/t\sim-10\%$, where $t_{ij}$ is the tunneling matrix element to the (1,1) and (2,0) neighboring site, respectively. We studied the effect of including these terms in the Hamiltonian on the CDW correlations (Fig. S4) and found it to be negligible within our experimental errors.

• •

Density-dependent tunneling The simple Hubbard model we use in this paper assumes that the tunneling matrix element for an atom is independent of the occupancy of the site it is tunneling to. This approximation starts to break down as the scattering length $a_{s}$ becomes large, and it may be expected that this effect is important in our system since most of the atoms are in pairs. The density-dependent hopping can be expressed in a tight-binding model as Jürgensen et al. (2014)

[TABLE]

where $w(\textbf{r})$ is the Wannier function of an atom on a site and d is the displacement vector to a neighboring site. For the upper branch experiments, $t_{n}/t\sim-8\%$ and for the lower branch experiments $t_{n}/t\sim-16\%$. To study the effect of the density-dependent tunneling on the measured CDW correlations and singles fraction, we performed Numerical Linked Cluster Expansion (NLCE) simulations where the tunneling was replaced by $t_{\sigma,\textbf{r},\textbf{r}^{\prime}}=t+(n_{-\sigma,\textbf{r}}+n_{-\sigma,\textbf{r}^{\prime}})t_{n}$, where if $\sigma$=$\uparrow$, $-\sigma$=$\downarrow$ and vice versa. In other words, the tunneling of a spin $\sigma$ from site r to $\textbf{r}^{\prime}$ depends on the total density of the opposite spin at the two sites between which the tunneling process is taking place. The strongest effect was about $10\%$ on the singles fraction and that would change the fit value of the data in Fig. 4b,c to $U/t=-5.4$ and $T=0.4t$. Both this value of $U/t$ and the one obtained without corrections to the Hubbard model fall within our experimental uncertainties.

Determinantal Quantum Monte Carlo calculations The DQMC simulations shown in the paper were all performed using a Fortran 90/95 package called QUantum Electron Simulation Toolbox (QUEST) developed and maintained by R. T. Scalettar et. al. Varney et al. (2009). For a spin-balanced system with attractive interactions ($U/t<0$) the calculations do not suffer from a fermion sign problem.

Simulations for Fig. 2 and Fig. 3 of the main text were performed on a square lattice of 8$\times$8 sites with $U/t=-5.7$ and a chemical potential ($\mu/t$) varying from -3 to 1.5 with $\mu=0$ representing half filling. The inverse temperature $\beta=Ld\tau$ was split into $L=40$ imaginary time slices, with an interval $d\tau=0.0556$. To obtain higher statistics, the simulations were run over 100,000 passes. For Fig. 4 of the main text, $U/t$ was fixed at -5.7 and $\mu$ to zero. The temperature was scanned by varying $d\tau$ from 0.0046 to 0.083 for a fixed $L=40$ and each point was averaged over 100,000 passes.

The DQMC simulations in the paper are performed using homogeneous systems. We rely on the local density approximation (LDA) for comparison to the experiment where the density varies slowly due to the harmonic trapping potential. For data shown in Fig. S3, the density gradient at $n=1$ varied from 0.07 atoms/(site)2 for a depth of $4.1~{}E_{R}$ to 0.14 atoms/(site)2 for a depth of $7.4~{}E_{R}$. To verify that the LDA holds for our system, we have performed DQMC simulations in the presence of a linearly varying chemical potential along one direction of the 2D Hubbard lattice to reproduce the maximal density gradients observed in the experiment. The results of this calculation showed that the LDA holds to excellent approximation for our experimental parameters. However, we note that violation of the LDA had been predicted for such a system at much lower temperatures and higher gradients Assmann et al. (2012).