Pair correlations in the attractive Hubbard model

The mechanism of fermionic pairing is the key to understanding various phenomena such as high-temperature superconductivity and the pseudogap phase in cuprate materials. We study the pair correlations in the attractive Hubbard model using ultracold fermions in a two-dimensional optical lattice. By combining the fluctuation-dissipation theorem and the compressibility equation of state, we extract the interacting pair correlation functions and deduce a characteristic length scale of pairs as a function of interaction and density filling. At sufficiently low filling and weak on-site interaction, we observe that the pair correlations extend over a few lattice sites even at temperatures above the superfluid transition temperature.

The nature of fermionic pairing plays a decisive role in many-body quantum states such as superconductors [1]. For attractive interactions and in the weak-coupling regime, a Bardeen-Cooper-Schrieffer (BCS) description gives rise to large, spatially-overlapping Cooper pairs below the critical temperature T c [2]. Upon increasing attractive interaction, the BCS ground state crosses over to a Bose-Einstein condensate (BEC) of tightly-bound dimers [3]. Above the transition temperature T c , the crossover of a normal state could be characterized by the formation of preformed pairs, which are linked to the emergence of a pairing pseudogap [4][5][6]. In a lattice configuration, the pairing phenomenon can be described in simple terms using the Hubbard model with an attractive s-wave interaction. Understanding the formation of preformed pairs offers a key insight into more complicated pairing mechanisms such as the d-wave pairing in high-temperature superconductors.
Quantum mechanically the pairing behavior between two particles is fully encoded in the second-order pair correlation function g (2) . In position space, g (2) (r) gives the joint probability of finding a pair of particles spaced by distance r. The pair correlation is of particular importance as it is related to various thermodynamical observables including compressibility, pressure and internal energy. In solid materials, g (2) (r) is inferred from the Fourier transform of the density structure factor S(q), which is experimentally accessible by crystallographic methods such as X-ray diffraction and neutron scattering [7,8]. For ultracold atomic gases in optical lattices that emulate the Hubbard model, the occurrence of spatial density and spin correlations with repulsive interactions were detected via high-resolution microscopy or scattering experiments [9][10][11][12][13][14]. More recently, the s-wave pairing pseudogap in the attractive Hubbard model was accessed via photoemission spectroscopy in a quantum gas microscope [15].
In this work, we investigate the formation of pair correlations in the two-dimensional attractive Hubbard model. In particular, we measure the thermodynamic correlation function and make use of the fluctuation-dissipation theorem to obtain an estimate of the pair correlation length.  ↑↑ (r) exhibits fermionic anti-bunching as a correlation hole. The unequal-spin correlation function g (2) ↑↓ (r) depends strongly on the on-site interaction strength U , where attractive (repulsive) interaction raises (lowers) pair correlation at short distance from the classical value of one.
We observe the formation of pairs with increasing attractive interactions above the superfluid critical temperature. The pair correlation length is of great importance as it can distinguish between a pseudogap pairing phase and a Fermi liquid phase above the critical temperature. While past experiments have determined the pair size in a trapped Fermi superfluid across the BCS-BEC crossover [16], it has not been measured for lattice gases.
Due to the Pauli principle, density correlations exist even for ideal fermions. The equal-spin pair correlation function g (2) ↑↑ (r) exhibits a Pauli correlation hole at short arXiv:2005.09121v1 [cond-mat.quant-gas] 18 May 2020 distance, as shown in Fig. 1. This gives rise to an effective repulsion that solely originates from the anti-symmetric nature of the fermionic state. Particles with opposite spin, however, are uncorrelated in an ideal Fermi gas and thus g (2) ↑↓ (r) = 1. For repulsively interacting fermions, the pair correlation between particles with opposite spin is suppressed with respect to the classical value of one [10], which goes hand in hand with the emergence of the Mott insulator. For attractive interactions, the enhancement of the pair correlation signals the formation of pairs.
The static density structure factor S(q) for constant filling n, which incorporates density fluctuations at all length scales, is linked to the Fourier transform of the pair correlation function, and is given by Here, g (2) (r) denotes the full pair correlation function σσ (r) + g (2) σσ (r) [17]. Although theoretically obtaining the full g (2) from Eq. (1) requires access to the density structure factor across all momenta, the limiting case at q = 0 still encapsulates both the Pauli blocking and interacting contributions of g (2) (r). In addition, the fluctuation-dissipation theorem connects the density structure factor at zero momentum to the corresponding thermodynamical susceptibility, in this case, the isothermal compressibility κ = (∂n/∂µ) | T , i.e. S(q = 0) = κT /n [10,18]. Combining this with Eq.
(2) provides us with two important insights. First, the isothermal compressibility entails the competition of Pauli repulsion and on-site attraction, since an external compression induces pressure between particles with equal and unequal spins within the trap. Second, the experimental determination of the compressibility opens access to the pair correlation function. While previous quantum gas experiments have been focusing on measurements of compressibility in relation to transport coefficients such as conductivity [19], its link to pairing has yet to be explored. Our experiment realizes the two-dimensional Hubbard model with attractive interaction on a square lattice, which reads Here c † i,σ (c i,σ ) denotes the fermionic creation (annihilation) operator at site i with spin σ, t is the nearestneighbour tunnelling amplitude, U < 0 represents attractive on-site interaction and µ i is the local chemical potential. We deploy the two lowest hyperfine states of 40 K in the F = 9/2 ground state manifold, serving as the two spin states |↑ = |F = 9/2, m F = −9/2 and |↓ = |F = 9/2, m F = −7/2 . Two-dimensional planes are formed by an optical lattice with a lattice depth of 120 E r , limiting tunnelling, where E r is the recoil energy of the optical light field. The in-plane square lattices are realized at a lattice depth of 6 E r . Here E r = h 2 8ma 2 = h×4.41 kHz is the recoil energy, h is the Plank's constant, a = 532 nm is the lattice spacing and m is the atomic mass. The lattice depth is calibrated by performing lattice modulation spectroscopy, where the laser intensity of the individual lattice beam is periodically modulated. The resonant frequencies indicate the relevant band transitions and therefore offer calibration of the lattice depths. The combination of lattices give rise to an in-plane tunnelling amplitude t = h × 224(6) Hz.
By tuning the magnetic field close to an s-wave Feshbach resonance near 202 G, we realize a wide range of attractive interactions −1 > U/t > −10 [20]. We obtain the on-site interaction energy U based on the analytical solution of two interacting fermions in an axially symmetric harmonic potential with an s-wave δ-pseudopotential [21]. We then employ a lattice-depth dependent correction factor to account for the anharmonicity of the potential around a single lattice site [22] and obtain the values of U/t presented in the manuscript. Using this calculation, we find a very nice agreement for the interaction energy shift between different hyperfine-state pairs measured in radio-frequency (RF) spectroscopy.
We note that we later utilize the local density approximation (LDA) in our analysis in which µ i = µ 0 −V (x, y), where µ 0 is the chemical potential at trap center and V (x, y) is the inhomogeneous confinement due to the optical lattice potential. To obtain this in-plane potential landscape V (x, y), we measure the in-plane trap frequencies ω x = 2π × 19.2(5) Hz and ω y = 2π × 25.9(3) Hz by exciting the dipole oscillation in an non-interacting gas. From the measured trap frequencies, we can infer the relevant lattice beam waists w x = 173(2) µm, w y = 152(1) µm and w z = 105(1) µm. Using these parameters, together with the calibrated lattice depths, we can reconstruct the confinement V (x, y) induced onto the atomic ensembles.
The atomic ensembles are then locked in position by rapidly ramping up the in-plane lattice depths to 60 E r within 1 ms. A single layer of a two-dimensional sample is transferred to another internal state via radio-frequency (RF) tomography in the presence of a vertical magnetic field gradient. We then utilize the interaction shift between singly-and doubly-occupied states to further separate them into two hyperfine states for detection. Next, we deploy high-resolution absorption imaging to obtain the in-situ density profiles of singles S ↑ = n ↑ − n ↑ n ↓ and doubles D = n ↑ n ↓ . We exploit the knowledge of the optical potential V (x, y) and use the local den- (g (2) (r) − 1)d 2 r versus interaction strengths U/t. The black solid line shows the non-interacting expectation at the lowest temperature in the data sets. Data points above the horizontal dashed line signal that the unequal-spin g (2) ↑↓ (r) outweighs its equal-spin counterpart, and vice versa. (c) Interacting contribution (g (2) ↑↓ (r) − 1)d 2 r is obtained by subtracting the equal-spin contribution (g sity approximation to determine the density equation of state n(µ) [23] (also see Appendix A). Thermometry is performed via a chi-squared fitting of n(µ) with determinant quantum Monte-Carlo (DQMC) simulations [24], which allows us to extract the temperature k B T /t (see Appendix B).
Subsequently, we obtain the isothermal compressibility κ = ∂n/∂µ through numerical differentiation of the measured equation of state n(µ). By taking the first derivative of a second-order polynomial fit to a subset of datapoints in n(µ). The fit is performed on the data points over a chemical potential window of h × 700 Hz around the desired µ. For low filling, i.e. n < 0.2, we observe increasing technical noise due to the lower signal-to-noise ratio of the raw data. Thus we deploy an exponential fit in order to avoid fitting negative compressibility. In order to obtain the measurement at the same filling n, as presented in Fig. 2 and 3 across different data sets, we interpolate neighbouring binned data points in µ.
In Fig. 2(a), we depict the density structure factor S(q = 0) derived from the compressibility using the fluctuation-dissipation theorem. We observe that for low filling S(q = 0) exceeds unity, which, quantitatively speaking, indicates particle-bunching. We note that S(q = 0) in Eq. (1) composes of both the equalspin (Pauli) and unequal-spin (interacting) pair correlation functions. At fixed filling n, we observe that the structure factor increases monotonically with interaction strength, suggesting that the interacting pair correlation functions become increasingly dominant.
We analyze the roles of both the equal-and the unequal-spin contributions by plotting the integral of the full g (2) (r) function, computed using Eq. (2), as shown in Fig. 2(b). The horizontal dashed line indicates the point at which the equal and unequal-spin g (2) exactly compensate each other, i.e. g (2) (r) − 1 d 2 r = 0 . The sign of this integral signals the dominant part in the pair correlation function. Thus, it offers a direct indication that the particle bunching observed in the density structure factor S(q = 0) is caused by the fact that the attractive on-site interaction dominates over the Pauli blocking.
In order to quantitatively compare the contributions from the Pauli principle and the attractive interaction, we note that the equal-spin contribution g ↑↑ (r) − 1 d 2 r can be calculated using the tight-binding dispersion relation (k x , k y ) = −2t [cos(k x a) + cos(k y a)] and temperature, see solid line in Fig. 2(b). Strictly speaking, this calculation is exact only for the non-interacting case because interactions in principle modify the dispersion from a simple sinusoidal energy band. However, in the low-filling regime, this estimate remains a faithful approximation since most of the occupied part of the energy band remains harmonic. Subtracting the equal-spin contribution from g (2) (r) − 1 d 2 r, we obtain the interacting  ↑↓ (r) − 1 d 2 r, as shown in Fig. 2(c). We observe that the interacting contribution maintains a similar dependence as S(q = 0). For either decreasing filling or increasing interaction strength, the interacting pair correlation increases. This highlights the parameter space at which the interaction effect on g 2 ↑↓ (r) is most prominent. At low filling, we observe a deviation of the measured pair correlations from the theoretical values, which we attribute to the low signal-to-noise near the low filling regime and the fact that the compressibility is vanishing for n → 0 (vacuum). Both results in insensitivity of density with respect to trapping potential variation and thus leads to increased uncertainties in the compressibility.
To gain further insight into the signature of the pairing, we turn to estimate a length scale up to which the pair correlation extends. We start by noticing that the interacting pair correlation amplitude at r = 0 is g (2) ↑↓ (0) = 4 n ↑ n ↓ n ↑ +n ↓ 2 = D (S ↑ +D) 2 . This implies that the amplitude can be directly obtained from our local density measurement of singly-and doubly-occupied site occupations. We note that in a lattice system, the double occupancy D plays the role of the contact parameter, describing short-range pair correlation [25,26]. Therefore, the non-trivial part of the pairing is reflected in the non-local part of the pair correlation function, which we analyze next.
Although the analytical form of the unequal-spin pair correlation function is not known, an exponential decay e −|r|/ξ is expected to be a good approximation above T c , see inset of inset of Fig. 3(a). Combining the knowledge of the integral g ↑↓ (r) − 1 d 2 r and the amplitude at r = 0, we then infer the characteristic length scale ξ as Eq. (4) renormalizes the measured pair correlation with respect to the on-site contribution. In Fig. 3(a), we plot the estimated pair correlation length ξ as a function of filling n for the same data set in Fig. 2. For low filling n 0.5, we observe a correlation length ξ as large as 0.92(4) a for the lowest interaction strength at U/t = −1.83, where a = 532 nm denotes the in-plane lattice spacing. Although the attractive interaction, therefore the Hubbard model, is purely on-site, we observe that its effect extends beyond the local site, similar to the Pauli blocking leading to beyond-local density suppression.
The trend of ξ as a function of filling can be attributed to two reasons. First, for dilute filling, particles are described by delocalized wave-packets. As temperature decreases, particles tend to explore the bottom of the energy band. This scenario resembles the free-particle case. Therefore, in the superfluid phase, BCS pairs and BEC dimers are expected for weak and strong interaction strengths, respectively. We observe a qualitative agreement to this expectation, despite our temperature being higher than the critical temperature. With increasing filling, the band occupation and thus the influence of the interaction term increases. This results in reduced interparticle spacing, and the latter drives the system away from the weak coupling limit despite the small U . Both contribute to the observed decrease in correlation length.
Second, the localization of particles at high filling means that the continuous integral of Eq. (4) starts to deviate from the discretized sum in a lattice. If the interacting contribution is dominated by the local term g ↑↓ (0), the continuous approximation of Eq. (4) would result in ξ/a ≈ 1/2π ≈ 0.4, which is in agreement with our observation in Fig. 3(a). Upon changing U/t, we observe that the correlation length ξ shrinks as interaction strength increases, as shown in Fig. 3(b). This signals the formation of tightly-bound local pairs and this pairing behavior is most prominent below quarter filling (n 0.5). Above quarter filling, we do not observe a discernible trend of ξ as a function of U/t, as indicated by the lowest data points in Fig. 3(b).
Last but not least, we investigate the temperature dependence of the correlation lengths. For temperatures below the critical temperature and low density, the pairing would be described by a BCS-BEC type behavior. Although the temperatures reached in our experiments remain in the normal phase, we observe a resemblance in the behavior of the pair correlation lengths. As shown in Fig. 3(c), we plot ξ as a function of temperature k B T /t at n = 0.2. The correlation length rises as temperature decreases for the weakly attractive case. In the strongly attractive case, we observe a much less significant trend in temperature dependence. Since the pairing occurs at an energy scale of the interaction U , it is also informative to recast the temperature with respect to |U |. In Fig. 3(d), we compare the temperature with respect to the interaction energy by plotting the correlation lengths as a function of k B T /|U |. Despite a lower achieved k B T /|U | for large interaction strengths, the correlation length remains small due to the energetically favorable dimer state. For weak interaction, the pair correlation length rises at much higher k B T /|U |, signaling the tendency to delocalize and form longer-range pairs.
In conclusion, we investigate the formation of pair correlations across the crossover regime of a normal phase in the attractive Hubbard model. In particular, we observe the competition between Pauli repulsion and on-site attraction. We show that for sufficiently low filling and weak interaction, the pair correlation length extends beyond local distance up to one site. This offers a clear signature for the formation of pairs as a potential precursor of the BCS-BEC crossover in a lattice configuration. Our measurement helps to elucidate the outstanding questions regarding the pairing behavior of a normal state. The approach presented here, in principle, also works beyond simple lattice geometries or on-site interaction, thereby allowing future investigation of long-range correlated systems.
This work has been supported by BCGS, the Using the knowledge of V (x, y), the recorded density profile n(x, y) is mapped to the chemical potential axis under the local density approximation, i.e. µ(x, y) = µ 0 − V (x, y). This allows us to perform iso-potential averaging on in-situ profiles as shown in Fig. 4 and to obtain both quantities as a function of local chemical potential µ. The density equation of state n(µ) can then be computed using n(µ) = 2 [S(µ) + D(µ)], where S and D are the occupation of singly-occupied state ("Singles") and doubly-occupied states ("Doubles"), respectively. Here, we utilize the fact that we are spin-balanced, i.e. n ↑ = n ↓ . The global chemical potential µ 0 and temperature k B T /t can be obtained by a numerical chisquare fit of n(µ) to the result of DQMC simulation.

APPENDIX B: DQMC SIMULATION
For the theoretical prediction shown in the manuscript, we perform determinant quantum Monte Carlo (DQMC) simulation on an 8 × 8 2D Hubbard model in a square lattice, using the Quantum Electron Simulation Toolbox (QUEST) Fortran package [24]. Simulations are performed with a wide range of interactions −10 ≤ U/t ≤ −1, with 1000 warm-up sweeps and 50000 measurement sweeps, and the number of imaginary time slices is set to 25. Typically the average sign of the equal-time Green's function in DQMC exhibits a sign problem for lower filling, especially with strong interaction and low temperature. However, we numerically verify that the sign prob- lem is negligible at our experimental temperature scales. The density structure factor S(q = 0) can be obtained by a summation over all individual density-density correlators.