Angle-resolved photoemission spectroscopy of a Fermi-Hubbard system

TL;DR

This paper develops ARPES for strongly interacting fermions in optical lattices, enabling the study of pseudogap phenomena and pairing mechanisms in quantum simulators, with benchmarking against quantum Monte Carlo calculations.

Contribution

The authors adapt ARPES to optical lattice systems and demonstrate its effectiveness by measuring the spectral function of an attractive Fermi-Hubbard model, revealing pseudogap behavior.

Findings

Evidence of a pseudogap above the superfluid transition temperature.

Successful benchmarking of ARPES measurements with quantum Monte Carlo simulations.

Potential application to doped repulsive Hubbard models.

Abstract

Angle-resolved photoemission spectroscopy (ARPES) measures the single-particle excitations of a many-body quantum system with both energy and momentum resolution, providing detailed information about strongly interacting materials. ARPES is a direct probe of fermion pairing, and hence a natural technique to study the development of superconductivity in a variety of experimental systems ranging from high temperature superconductors to unitary Fermi gases. In these systems a remnant gap-like feature persists in the normal state, which is referred to as a pseudogap. A quantitative understanding of pseudogap regimes may elucidate details about the pairing mechanisms that lead to superconductivity, but developing this is difficult in real materials partly because the microscopic Hamiltonian is not known. Here we report on the development of ARPES to study strongly interacting fermions in an optical lattice using a quantum gas microscope. We benchmark the technique by measuring the occupied single-particle spectral function of an attractive Fermi-Hubbard system across the BCS-BEC crossover and comparing to quantum Monte Carlo calculations. We find evidence for a pseudogap in our system that opens well above the expected critical temperature for superfluidity. This technique may also be applied to the doped repulsive Hubbard model which is expected to exhibit a pseudogap at temperatures close to those achieved in recent experiments.

Figures

Equations

$$A_{\text{trap}}(k,\omega )=\frac{1}{N}\sum _{i}A(k,\omega +{\mu }_i-\mu ,{\mu }_i),$$

$$G^R(k,t)$$

$$G^R(k,\omega )$$

$$=$$

$$A(k,\omega )$$

$$=$$

$$\rho (\omega )$$

$$⟨n_k⟩$$

$${\int }_{-\infty }^{\infty }d\omega A(k,\omega )$$

$$\mathcal{H}$$

$${\mathcal{H}}^{\prime }$$

$${\mathcal{H}}_o$$

$${\dot{n}}_{kf}$$

$$I(k,h\nu )$$

$$G^R(k,t-t^{\prime })$$

$$I(k,h\nu )$$

$$A(k,\omega )F(\omega )$$

$$I(k,h\nu )$$

$$I(k,{ϵ}_k-\omega -\mu )$$

$$A_{\text{occ, trap}}(k,\omega )$$

$$A_{\text{trap}}(k,\omega )$$

$$A(k,\omega ,n)$$

$$U(x,t|x^{\prime },0)=\sqrt{\frac{1}{2\pi i\sin (\omega t)l^2}}\exp \left[i\frac{1}{2l^2\sin (\omega t)}\left((x^2+x^{\prime 2})\cos (\omega t)-2x{x}^{\prime }\right)\right],$$

Full text

Photoemission spectroscopy measures the occupied single-particle spectral function Damascelli et al. (2003); Chen et al. (2009); Törmä (2016), which describes the allowed energies for a single-particle excitation with given momentum. ARPES has been used to study the presence of a Fermi surface, the lifetime of quasiparticles, superconducting gaps and their symmetries, and surface states in topological materials Vishik et al. (2010); Hasan and Kane (2010). It is illustrative to consider some generic features of the spectral function as we introduce interactions. For a non-interacting system with dispersion $\epsilon_{k}$, the single-particle excitations are eigenstates of the system implying they have a definite energy and infinite lifetime. The spectral function is a delta function, $A(k,\omega)=\delta(\omega-\epsilon_{k}+\mu)$ where $\mu$ is the chemical potential. Turning on weak interactions may create a Fermi liquid state, where the spectral function has a similar form but broadens along the energy axis reflecting the finite lifetime of the quasiparticles. In superconducting systems, more radical changes may occur including the development of a gap separating two disconnected spectral function branches. In weak coupling (BCS) superconductors, the gap vanishes at the critical temperature. However in many strongly-interacting superconducting systems including the high-temperature cuprate superconductors (HTSCs) and the unitary Fermi gas, a depression in the spectral weight at the chemical potential persists in the normal state Mueller (2017); Chin et al. (2004); Schunck et al. (2007); Gaebler et al. (2010); Nascimbène et al. (2011); Feld et al. (2011); Sommer et al. (2012); Murthy et al. (2017). These so called pseudogap states might arise from the same pairing mechanism as the superconducting ground states, and developing a better understanding of their properties may shed light on pairing physics in the HTSCs.

Photoemission spectroscopy with cold atoms has been used to study the pseudogap regime in strongly interacting Fermi gases without a lattice. The idea was first proposed using Raman techniques Dao et al. (2007), but was realized using momentum-resolved radiofrequency (rf) spectroscopy Stewart et al. (2008). Previous experiments have explored the BCS-BEC crossover in 3D Stewart et al. (2008); Gaebler et al. (2010); Sagi et al. (2015) and 2D Feld et al. (2011); Fröhlich et al. (2012) continuum systems. These experiments inferred the presence of a pseudogap from spectral features that exhibit “back-bending” near the Fermi surface at temperatures above $T_{c}$, in analogy to the back-bending observed in the BCS dispersion. This interpretation has been criticized because back-bending can also arise due to other effects, including universal short range physics expected in continuum systems Schneider and Randeria (2010). Exact theory techniques were not available for these systems, therefore experiments compared with various $T$-matrix approximation schemes Loktev et al. (2001); Perali et al. (2002, 2011) which become exact for weak interactions.

Here we extend photoemission spectroscopy to strongly-correlated lattice systems. We take advantage of quantum gas microscopy techniques to obtain high signal-to-noise data using a single layer 2D system, avoiding the complication of integration over a third axis of momentum space. Combining photoemission spectroscopy with quantum gas microscopy poses a technical challenge because photoemission spectroscopy relies on a large field of view for time-of-flight measurements to reach momentum space, whereas quantum gas microscopy is limited to a small field of view over which the atoms can be pinned during imaging. We solve this problem by developing a new measurement protocol which consists of four parts, illustrated in fig. 1a. First, we excite atoms from the interacting system to a non-interacting final state using a radiofrequency pulse. Next, we perform band mapping to adiabatically connect quasimomentum states in the lattice to real momentum states. Then, we allow the atoms to expand in a harmonic trap for a quarter period which maps the atomic momentum distribution to the real space distribution Murthy et al. (2014). Finally, we turn on a deep optical lattice to freeze the atoms in place in preparation for site resolved imaging. The momentum detection portion of this protocol has been discussed in a recent proposal for measuring the spectral function in 1D systems using lattice modulation techniques Bohrdt et al. (2018).

We benchmark the ARPES technique by studying the pseudogap regime in the attractive Hubbard model. There has been comparatively little experimental work on this system Strohmaier et al. (2007); Hackermuller et al. (2010); Schneider et al. (2012); Mitra et al. (2017), but it has a number of advantages as a testing ground for our new protocol. Unlike the repulsive model, the attractive model does not have a sign problem, meaning it is possible to make high-precision theory calculations using DQMC at low temperatures. Furthermore, the 2D attractive model exhibits phenomenology similar to other strongly interacting superconducting systems. In particular, it exhibits a BCS-BEC crossover with increasing on-site interaction strength and a Berezinskii-Kosterlitz-Thouless (BKT) transition to an $s$-wave superconducting state at low temperatures Paiva et al. (2004). At weak interactions the ground state is BCS-like, and pairing is destroyed at the critical temperature for superconductivity, $T_{c}$. As the interaction strength is increased, pairing persists above the critical temperature and is destroyed at some higher temperature $T^{*}>T_{c}$. Intermediate temperatures, $T_{c}<T<T^{*}$, exhibit a pseudogap which appears as a depression in either the density of states or spectral function at the Fermi energy Singer et al. (1996, 1998, 1999). As the interaction strength increases, the pairing changes character from momentum-space to real-space and the pseudogap temperature increases monotonically. In the BEC limit, the pseudogap represents a regime of real-space bosonic pairs which develop phase coherence below $T_{c}$.

In this experiment, we study an attractive Hubbard system with interactions ranging from intermediate coupling, $U/t=-4$, to strong coupling, $U/t=-8$. We observe the emergence of a pseudogap as we increase the interaction at fixed temperature, $T/t\sim 0.5$. This is well above the theoretically calculated maximum $T_{c}$ for this model $T/t\sim 0.15$ Paiva et al. (2004). At strong coupling and hotter temperatures, we observe the presence of two branches in the occupied spectral function due to the coexistence of single and paired atoms.

To begin, we prepare an equal mixture of hyperfine states $\left\lvert\uparrow\right\rangle=\left\lvert 3\right\rangle$ and $\left\lvert\downarrow\right\rangle=\left\lvert 1\right\rangle$ in a 2D trap at a field of $\sim$120\text{,}\mathrm{mG}\text{/} where the scattering length $a_{13}\sim$-800\text{\,}\mathrm{a_{o}}\text{/} O’Hara et al. (2000), where $\mathrm{a_{o}}\text{/}$ is the Bohr radius and $\left\lvert i\right\rangle$ denotes the $i$th lowest hyperfine state of lithium. We adiabatically load a 2D optical lattice with lattice spacing $a=$752\text{,}\mathrm{nm}\text{/}to depths between $5$ and $10\text{\,}\mathrm{E_{r}}\text{/}$, where\mathrm{E_{r}}\text{/}$/h=$14.7\text{,}\mathrm{kHz}\text{/} is the lattice recoil energy. We then apply a $1\text{\,}\mathrm{ms}\text{/}$ radiofrequency (rf) pulse to transfer $\sim 15\%$ of the $\left\lvert\uparrow\right\rangle$ atoms to a third state $\left\lvert f\right\rangle=\left\lvert 2\right\rangle$ (fig. [1](#S0.F1)a, left panel). Final state interactions are negligible because the scattering lengths $a_{12}\sim$0\text{\,}\mathrm{a_{o}}\text{/} and $a_{23}\sim$25\text{,}\mathrm{a_{o}}\text{/} are small. We determine the initial quasimomentum distribution of the probed atoms by performing band mapping and subsequently expanding for a quarter period in a harmonic trap. To achieve this, we first turn on a harmonic trapping potential with trapping frequency $\omega=(2\pi)$360\text{\,}\mathrm{Hz}\text{/} in a few microseconds. Then, we ramp the lattice potential to zero linearly in $200\text{\,}\mathrm{\SIUnitSymbolMicro s}\text{/}$ before expanding for $T/4=$700\text{,}\mathrm{\SIUnitSymbolMicro s}\text{/}$$. The real space density distribution after the quarter period expansion reflects the initial quasimomentum distribution of the $\left\lvert f\right\rangle$ atoms according to $\tilde{n}_{f}(k)=n_{f}(r)l^{4}$, where $l=\sqrt{\hbar/m\omega}$ is the harmonic oscillator length and $k=r/l^{2}$.

We obtain the energy dependence of the spectral function by measuring the density profile at a range of rf frequencies. For each frequency, $\nu$, we obtain a two-dimensional image, $n_{f}(r,\nu)$ as shown in fig. 1b. Each pixel in the image is associated with a different quasimomentum. To improve our signal-to-noise, we spatially average using the symmetry of the square. The final density distribution of state $\left\lvert f\right\rangle$ atoms is proportional to the occupied single-particle spectral function, $n_{f}(r,\nu)\propto A(\mathbf{k},\epsilon_{\mathbf{k}}-h\nu-\mu)F(\epsilon_{\mathbf{k}}-h\nu-\mu)$ sup , where $F$ is the Fermi function. The rf photon carries negligible momentum, therefore the initial momentum is the same as the final momentum. The energy of the final state is given by the tight-binding dispersion $\epsilon_{\mathbf{k}}/t=4-2\cos(k_{x}a)-2\cos(k_{y}a)$, where $t$ is the tunneling energy, and the rf photon energy $h\nu$ leading to an initial energy $\epsilon_{\mathbf{k}}-h\nu$. This transition is illustrated schematically in momentum space (fig.1a, left panel).

The density profiles change shape as the rf frequency is scanned, a characteristic of an interacting system. At large red detuning relative to the non-interacting transition, the transferred density is more heavily weighted towards the edges of the the Brillouin zone (fig. 1b, leftmost). As the detuning is decreased, the density becomes peaked at the center and the transfer increases (center). Decreasing the detuning further leads to less transfer but does not significantly change the shape (rightmost). To further quantify the energy dependence of the spectral function we focus on individual quasimomenta and plot the corresponding energy distribution curves (EDCs) in fig. 1d. We compare our experimental results with determinantal quantum Monte Carlo (DQMC) calculations at $U/t=-7.5$ and $T/t=0.55$, which are shown in fig. 1c and d and find good agreement. The full DQMC spectral function (fig. 1d) exhibits a depression near the central chemical potential, indicating a pseudogap at this temperature and interaction.

Having demonstrated the viability of this photoemission spectroscopy technique, we focus on elucidating the properties of the attractive Hubbard system by exploring the evolution of the spectral function with interaction strength and temperature. First, we consider the low temperature regime $T/t\sim 0.5$ for a range of interactions $U/t=-4$, $-6$, and $-7.5$. For each interaction strength we determine the temperature, interaction, and global chemical potential by fitting in-situ densities and correlators to DQMC results sup . We determine the occupied spectral function from the density according to $A_{\text{occ}}(k,\omega)=n_{f}(k,\epsilon_{k}-\omega-\mu)$ and obtain the profiles shown in fig. 2. For the weakest interaction we consider, $U/t=-4$, the spectral weight is shifted to negative energy relative to the non-interacting dispersion and there is significant spectral weight near the chemical potential. As we increase the interaction, the spectral weight shifts further from the non-interacting dispersion and there is less spectral weight near the chemical potential.

We quantify the distribution of spectral weight by fitting each EDC to a Lorentzian profile and extracting the center position and half width half maximum (HWHM) values. We emphasize that the resulting graph cannot be interpreted as a dispersion for two reasons. First, there are trap averaging effects and second even in a homogeneous system the Fermi function shifts the position of the peak in $A_{\text{occ}}(k,\omega)$ relative to $A(k,\omega)$, causing significant deviation from the true dispersion near and above the Fermi surface sup . For $U/t=-4$, the center of the spectral weight reaches the chemical potential at a point along the $X-M$ line in momentum space (fig. 2a). For stronger interactions, $U/t=-6,-7.5$ the center of spectral weight no longer approaches the chemical potential. This pulling back of the spectral weight at the chemical potential is a signature of pseudogap behavior. The spectral width increases weakly with the interaction. At $U/t=-4$, the HWHM is $\sim 1.2t$ with little dependence on the quasimomentum. As we increase the interaction the width increases to $\sim 2t$ and $\sim 2.5t$ at $U/t=-6$ and $U/t=-7.5$ respectively. The measured width is larger than the Fourier broadening which is $1\text{\,}\mathrm{kHz}\text{/}$, corresponding to $\sim 0.7t$, $\sim 0.8t$, and $\sim 1t$ for these interactions. The width is not intrinsic but due to the spread of chemical potentials across the harmonic trap sup .

We compare our measurements with DQMC calculations which account for both the presence of the trap and Fourier broadening of the rf pulse. These results, shown in the lower panels of fig. 2, display similar trends to the experiment. We determine the proportionality constant between the experimental data and the theoretical spectral function using a least squares fit. After scaling the experimental data, the theory agrees well at low quasimomentum sup . At larger quasimomentum near the $X$ point, the experimental spectral function is systematically smaller than the theory. This may be due to imperfections in the band mapping process which are more significant near the edge of the band. A pseudogap appears as a depression in the spectral function at the chemical potential in the theory results for $U/t=-6$ and $-8$ sup . The pseudogap is not present at $U/t=-4$, consistent with the calculations of Singer et al. (1998). We determine the peak locations using a peak finding algorithm, which identifies multiple peaks more efficiently in the high quality numerical data.

The spectral function is broadened by the presence of the harmonic trapping potential which causes the system to sample a range of chemical potentials. We account for this effect in our theory results using the local density approximation. The trap spectral function is given in terms of the homogeneous system spectral functions by

[TABLE]

where $N$ is the number of sites in our trap, $\mu$ is the global chemical potential, and $\mu_{i}$ is the local chemical potential at site $i$. The trap causes broadening of spectral features and tends to obscure signatures of gaps or pseudogaps. We illustrate this effect for an experimental density profile in fig. 3. The theoretical homogeneous system spectral functions shown in fig. 3b exhibit a pseudogap for densities $n>0.2$. When these functions are combined using eq. 1, the pseudogap feature is broadened over the spread of chemical potentials where the uniform system exhibits a pseudogap (fig. 3c). This effect is also visible in fig. 1c, where the pseudogap feature appears shifted away from the central chemical potential. Before comparing these results with experiment, we account for Fourier broadening of the rf pulse by convolving the trap averaged spectral function with a Gaussian with standard deviation $h\times$1\text{,}\mathrm{kHz}\text{/}$$.

We also consider the evolution of the spectral function with increasing temperature at the strongest interaction, $U/t=-7.5$. At the lowest temperature we find little weight at the Fermi level, fig. 4a. As the temperature is increased, a second peak emerges which approximately follows the non-interacting dispersion, fig. 4b,c,d. At the same time, the lower peak shifts to more negative energy and the spectral function is gapped at all momenta. At the hottest temperatures (fig. 4c,d), the Fermi weighting approaches a constant value and occupied spectral function is proportional to the full spectral function, $A_{\text{occ}}(k,\omega)\sim A(k,\omega)/2$. The upper feature is due to unpaired atoms and the lower feature due to pairs.

We fit the EDCs to either one Lorentzian for the coldest temperatures or two Lorentzians for the hotter temperatures. We find that temperature has a weak effect on the peak widths. The HWHMs of the upper and lower peak are $\sim 1.5t$ and $\sim 2.5t$ respectively. The theory results show that the dispersion obtained from the occupied spectral function nearly coincides with that obtained for the true spectral function $T/t\sim 0.5-3$. They also show that a pseudogap is present at $U/t=-8$ up to $T/t=5$ sup . The large value of $T^{*}$ is characteristic of being in the strong coupling regime.

In conclusion, we have developed a technique to measure the single-particle spectral function in Hubbard systems using a quantum gas microscope. In combination with recent transport experiments Brown et al. (2019); Nichols et al. (2019), this extends the toolkit for measuring dynamical properties of these systems, which is so far rather limited compared with the techniques available to study HTSCs and other real materials. Our experiment paves the way for future studies of the single-particle spectral function in the doped repulsive Hubbard model, a quantity which is difficult to calculate using exact theory techniques because of the fermionic sign problem. Such an experiment might elucidate properties of the Hubbard system at high temperature where dynamical mean field theory studies have suggested quasiparticles survive into the linear resistivity regime Xu et al. (2013); Deng et al. (2013) or test recent parton theories which may describe the doped antiferromagnetic regime Grusdt et al. (2018). At lower temperatures currently beyond the reach of exact theory methods, this technique could be used to search for the symmetry of the pseudogap and superconducting states.

I Acknowledgments

Acknowledgements.

This work was supported by the NSF (grant no. DMR-1607277), the David and Lucile Packard Foundation (grant no. 2016-65128), and the AFOSR Young Investigator Research Program (grant no. FA9550-16-1-0269). TPD and EWH acknowledge support from the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under Contract No. DE-AC02-76SF00515. Computational work was performed on the Sherlock cluster at Stanford University.

II Supplementary Material

III The single-particle spectral function

The single-particle spectral function comes from the single-particle retarded Green’s function, which is defined by

[TABLE]

This function gives information about the single-particle excitations of a many body system. The first term in the anticommutator, $\left<c_{k}(t)c^{\dagger}_{k}(0)\right>$, describes what happens if we create a particle excitation with momentum $k$, watch this excitation propagate until time $t>0$, and then measure what portion of the excitation survives. The second term, $\left<c^{\dagger}_{k}(0)c_{k}(t)\right>$, can be interpreted similarly but for holes instead of particles.

Information about the allowed energies for an excitation with momentum $k$ is contained in the spectral (Lehmann) representation of the Green’s function, which we arrive at by Fourier transforming and expanding using a basis of eigenstates $\left\lvert n\right\rangle$ with energies $\epsilon_{n}$ and partition function $Z$

[TABLE]

Here $\eta\in\mathbb{R}^{+}$ and we always consider the limit $\eta\to 0^{+}$. The terms $G^{+}$ and $G^{-}$ come from the anticommutator in $G^{R}$ and can again be associated with creating a particle and hole excitation respectively.

We define the single-particle spectral function as

[TABLE]

$A_{\text{unocc}}$ represents particle injection and $A_{\text{occ}}$, the occupied spectral function, represents photoemission.

The spectral function counts single-particle excitations at a given momentum and frequency. In this sense it may be thought of as a momentum resolved density of states. More precisely, we have the following relationship between the density of states, $\rho(\omega)$ and the spectral function

[TABLE]

The spectral function is also closely related to the momentum distribution function

[TABLE]

Finally, the spectral function is normalized to unity

[TABLE]

IV Linear Response

The Hamiltonian of our system is given by

[TABLE]

where $\mathcal{H}^{\prime}$ describes the radiofrequency probe with Rabi frequency $\Omega$ and $\mathcal{H}_{o}$ is the Fermi-Hubbard Hamiltonian on a square lattice with $\epsilon_{k}=-2t\left[\cos(k_{x})+\cos(k_{y})-2\right]$. We choose the band minimum of the non-interacting dispersion as zero energy. This implies that half-filling occurs at $\mu/t=4+U/2$.

In the experiment, we measure the total density of atoms transferred by the perturbation portion of the Hamiltonian acting for some time $\delta t$ on our system. Assuming that the time is short enough and the perturbation is weak enough, we can apply linear response theory and the atom density is proportional to the atom current. The atom current operator is given by

[TABLE]

The expectation value of the momentum resolved current for rf photon energy $h\nu$ is given by the Kubo formula as

[TABLE]

Supposing that the final state $\left\lvert f\right\rangle$ is unoccupied and does not interacting with the initial state $\left\lvert\uparrow\right\rangle$, this expression can be simplified Chen et al. (2009); Törmä (2016) to arrive at

[TABLE]

where $\mu$ is the chemical potential of the interacting system, and $A$ is the spectral function of the interacting system described by $\mathcal{H}_{o}$. This expression assumes that the energy of the final state $\left\lvert f\right\rangle$ is larger than the initial state. If the energy of the final state is smaller, then we should replace $h\nu$ with $-h\nu$ in the above expression.

In this expression $\epsilon_{k}$ is the energy of the outcoupled atom, implying that $\epsilon_{k}-h\nu$ is the initial energy. We must additionally subtract $\mu$ from the argument because the spectral function is defined such that $\omega=0$ represents the system at the chemical potential.

This expression implies that the momentum resolved current is a measure of the occupied spectral function, however the presence of $\epsilon_{k}$ in the energy argument means that a measurement at a single frequency $\nu$ is a complicated surface in $k-\omega$ space.

We can rewrite eq. S17 to give us the occupied spectral function at equal $\omega$,

[TABLE]

V Trap averaging

Due to the harmonic trapping potential, the density of our Hubbard system develops a spatial dependence. Working in the semiclassical local density approximation (LDA), we model this as an effective spatial variation of the chemical potential according to $\mu_{i}=\mu-V_{i}$, where $V_{i}$ is the potential at site $i$. The current is now given by

[TABLE]

where we explicitly include the dependence of $A$ on $\mu$.

Determining the spectral function for the trapped system is more complicated than in the uniform system (eq. S18) because the convention is to reference $\omega$ to the chemical potential. This choice of reference is of little consequence when working at fixed density, but when working with multiple densities we must use a consistent reference point. In this case the physically meaningful quantity is $\omega+\mu_{i}$. In the trapped system our energy reference is the chemical potential (or the central chemical potential in LDA language), meaning that $\omega+\mu_{i}-\mu$ becomes the argument of our spectral functions

[TABLE]

As we expect eq. S17 to hold for the spectral function for the trap $A_{\text{trap}}$, we identify

[TABLE]

where normalizing by the total number of sites in the trap ensures that the eqs. S8, S9, and S10 continue to hold.

In the experiment we have direct access to the density on each site, $n_{i}$. Therefore, using DQMC we effectively compute $\mu(n)$, the chemical potential as a function of the atomic density.

We compute the spectral function at the densities $n=0.2-1.8$ in steps of $0.2$ and for each density in the trap we use the closest computed density to determine the chemical potential.

VI EDCs experiment to theory comparison

It is instructive to compare the energy distribution curves for the experiment and theory data shown in the main text. Presenting the data in this way allows a direct comparison of theory and experiment. Here we present the EDCs corresponding to fig. 2 in the main text in fig. S1, and the EDCs corresponding to fig. 4 in the main text in fig. S2.

VII Full spectral functions

Determining the presence of a pseudogap in the experimental system is complicated by the fact we measure the occupied spectral function. However, in the theory we can directly see the pseudogap open in the full spectral function. Here we present the full trap averaged spectral functions corresponding to the occupied spectral functions shown in figs. 2 and 4. Here, unlike in the main text, we do not use any Gaussian broadening.

In fig. S3, a pseudogap is visible for $U/t=-6$ and $-8$ at lower temperatures $T/t\sim 0.5$. No pseudogap is visible for $U/t=-4$. This figure corresponds to fig. 2 in the main text.

In fig. S4 a pseudogap is visible for $T/t<3$ and appears to have disappeared by $T/t=5$. This figure corresponds to fig. 4 in the main text.

VIII Dispersion from $A(k,\omega)$ versus $A_{\text{occ}}(k,\omega)$

The dispersion inferred by determining the peak in spectral weight at each momentum is significantly affected by the presence of the Fermi function. Near the chemical potential this causes the dispersion of the occupied spectral function to differ significantly from the real dispersion. This effect is most pronounced outside of the pseudogap regime because in that case there is significant spectral weight at the chemical potential. To illustrate this effect, we perform the peak-finding procedure described in the text using both the full and occupied spectral functions. We show only peaks with heights greater than $0.02$. The results are shown in fig. S5.

IX DQMC calculation of the single-particle spectral function

We perform DQMC simulations on the attractive Hubbard model Blankenbecler et al. (1981); White et al. (1989) with interaction strengths and temperatures $U/t=-3.7,T/t=0.48$; $U/t=-6,T/t=0.5$; $U/t=-7.5;T/t=0.55$; $U/t=-8,T/t=0.44$; $U/t=-8,T/t=1.04$; $U/t=-8,T/t=1.2$; $U/t=-8,T/t=2.83$; and $U/t=-8,T/t=5$. For each parameter set, we consider filling levels of $\left<n\right>=0.2,0.4,0.8,$ and $1$. The chemical potential is tuned to achieve the desired doping level to within an accuracy of $O(10^{-4})$. The imaginary time interval $[0,\beta]$ is discretized into steps of around $0.05$ or less, resulting in negligible Trotter errors for our simulations. All simulations were on square clusters of size $16\times 16$. We run $160$ independently seeded Markov chains with $1000$ warmup and $5000$ measurement sweeps each, giving a total of $8\times 10^{5}$ sweeps for each doping and temperature. Unequal time measurements are performed on every other sweep.

To ensure numerical stability in computing the equal-time Green’s functions, we use the prepivoting stratification algorithm as described in Tomas et al. (2012), allowing up to 10 matrix multiplications before performing a QR decomposition. The unequal time Green’s functions are constructed using the Fast Selected Inversion algorithm described in Jiang et al. (2016), with blocks corresponding to the product of matrices from 10 time steps.

We perform maximum entropy analytic continuation (MaxEnt) Jarrell and Gubernatis (1996) to extract the single particle spectral function from the imaginary time Green’s functions measured in DQMC. The MaxEnt optimization utilizes Bryan’s algorithm as described in Jarrell and Gubernatis (1996). For the choice of model function, we use a non-informative flat function, and select the parameter $\alpha$ using the prescription given in Bergeron and Tremblay (2016). We have checked that the results are nearly identical when using the Classic formulation of MaxEnt with an annealing method for model functions Jarrell and Gubernatis (1996).

To determine the spectral function for densities above half-filling, we take advantage of the symmetry

[TABLE]

where $Q=(\pi,\pi)$.

X Fitting the temperature and interaction

To determine the interaction, temperature, and chemical potential of our gas we perform a simultaneous fit of the quantities $C_{\uparrow}(1)$, $C_{\downarrow}(1)$, $n_{s}$, $C_{s}(1)$, $d$, $C_{d}(1)$ versus total density to DQMC data with $U/t$ and $T/t$ as free parameters. Here $C_{\uparrow}(1)=\left<n_{\uparrow}(i_{x},i_{y})n_{\uparrow}(i_{x}+1,i_{y})\right>-\left<n_{\uparrow}(i_{x},i_{y})\right>\left<n_{\uparrow}(i_{x}+1,i_{y})\right>$ and similarly for the other correlators. By fitting these various correlators versus density we remove the need to precisely character the trapping potential by taking advantage of the known equation of state. The procedure is similar to that described in Brown et al. (2019), but using additional quantities such as the doublon-doublon correlator and the single-single correlator allows simultaneous determination of both $U/t$ and $T/t$ and leads to higher precision fitting. The result of an exemplary fit is shown in fig. S6. In this figure, the DQMC data has been corrected using imaging fidelities and the experimental correlators are the values measured in the experiment. We used a doublon transfer efficiency of $0.9$ and an imaging fidelity of $0.96$. Note that the doublon imaging efficiency is the product of these two efficiencies. Densities are corrected by a single factor of the relevant efficiency, and correlators are corrected by the efficiency squared. The fit in fig. S6 yields a temperature of $T/t=0.38$ and interaction of $U/t=-4.2$.

XI Band mapping simulation

As our band mapping time, $200\text{\,}\mathrm{\SIUnitSymbolMicro s}\text{/}$, is not small compared with the quarter period expansion time of $700\text{\,}\mathrm{\SIUnitSymbolMicro s}\text{/}$, we perform a 1D simulation of the procedure following McKay et al. (2009) to check what effect atomic motion during the band mapping has on the final density profiles.

We first find the lowest eigenstates of a single-particle Hamiltonian with potential given by the sum of a lattice and harmonic trap in real space. Then we change the lattice depth linearly from the initial depth to zero while time evolving the initial eigenstates using the Crank-Nicholson method. We simulate the quarter period expansion in the harmonic trap evolution using the propagator

[TABLE]

where $l=\sqrt{\frac{\hbar}{m\omega}}$ is the harmonic oscillator length.

The results of this procedure are shown in fig. S7 for an initial temperature of $T/$\mathrm{E_{r}}\text{/}$=1/30$, a lattice depth of $6.5\text{\,}\mathrm{E_{r}}\text{/}$, a harmonic trapping frequency of $\omega_{i}=(2\pi)$200\text{,}\mathrm{Hz}\text{/} and a band mapping time of $200\text{\,}\mathrm{\SIUnitSymbolMicro s}\text{/}$. The expansion trap frequency is $\omega_{f}=(2\pi)$360\text{\,}\mathrm{Hz}\text{/} and one quarter period is $\sim$700\text{,}\mathrm{\SIUnitSymbolMicro s}\text{/}$$.

We compare the results of the band mapping procedure to the initial quasimomentum distribution, which we determine by projecting the eigenstates of the initial trap on to the Bloch eigenstates of the lattice. We determine the scaling factor between the two using a least-squares fit as we do for the experimental data. This fit typically requires less than a $10\%$ rescaling of the density. We find that the largest discrepancy occurs at the edges of the Brillouin zone where the gaps are smallest. The average fractional discrepancy between the quasimomentum distribution and the scaled band map density distribution is less than $10\%$.

In the experiment the trap parameters are $\omega=(2\pi)$160\text{,}\mathrm{Hz}\text{/}, $t/h=$1400\text{\,}\mathrm{Hz}\text{/}, and $5\text{\,}\mathrm{E_{r}}\text{/}$ at $U/t=-4$; $\omega=(2\pi)$185\text{,}\mathrm{Hz}\text{/}, $t/h=$1230\text{\,}\mathrm{Hz}\text{/}, and $5.5\text{\,}\mathrm{E_{r}}\text{/}$ at $U/t=-6$; and $\omega=(2\pi)$200\text{,}\mathrm{Hz}\text{/}, $t/h=$1000\text{\,}\mathrm{Hz}\text{/}, and $6.5\text{\,}\mathrm{E_{r}}\text{/}$ at $U/t=-8$.

XII Magnetic field stabilization

To reach $120\text{\,}\mathrm{mG}\text{/}$ we null Earth’s field in the $x$ and $y$ directions using stable current sources. Due to the high magnetic field sensitivity of the $\left\lvert 2\right\rangle-\left\lvert 3\right\rangle$ rf transition at low fields, we must actively stabilize the field in the $z$ direction. We measure the magnetic field using a Stefan Mayer FLC100 fluxgate sensor Dedman et al. (2007), and control it using a Helmholtz coil of radius $\sim$0.5\text{,}\mathrm{m}\text{/}$$ and a home built servo driving a power MOSFET.

We test the stability of our field by measuring the width of the $\left\lvert 2\right\rangle-\left\lvert 3\right\rangle$ transition using a polarized gas of $\left\lvert 3\right\rangle$ atoms in a deep lattice, ($\sim$60\text{,}\mathrm{E_{r}}\text{/}), which freezes out tunneling. We fit the spectra to a Lorentzian line shape and find a half width half max of $\sim$300\text{\,}\mathrm{Hz}\text{/}. The full width corresponds to a field fluctuation of $0.3\text{\,}\mathrm{mG}\text{/}$.

We assess the long term stability of our field by measuring the resonance frequency of the polarized gas before and after each experimental data set. We see frequency drifts of $\sim$1\text{,}\mathrm{kHz}\text{/}$$ ($0.5\text{\,}\mathrm{mG}\text{/}$), which is of the same order as our Fourier broadening.