Dynamical structure factor and a new method to measure the pairing gap in two-dimensional attractive Fermi–Hubbard model

The superfluid state of quantum many-body Fermi gases exhibits a nonzero pairing gap (order parameter) originating from the Cooper pairing. Developing a convenient method to measure this pairing gap is essential for understanding many-body pairing phenomena and dynamical excitations. Currently, the pairing gap is primarily obtained by virtue of various excitation processes, such as the momentum-resolved photo-emission spectroscopy^[1,2] and the radio-frequency spectroscopy.^[3–5] However, these techniques encounter significant challenges in systems with complex band structures.^[6–11] Notably, the dynamical excitations offer a powerful alternative for probing the pairing correlations. These excitations are encoded in the dynamical structure factor, which can be experimentally measured through the two-photon Bragg spectroscopy.^[12–18]

In continuous Fermi gases, the collective modes are typically probed at a small transferred momentum q, whereas the single-particle excitations of the unpaired and paired atoms emerge at a relatively larger q. Specifically, the excitations of the paired atoms correspond to the bosonic molecular excitations. In 2008, Vale et al. investigated the single-particle excitations at large transfer momenta q≫k_F. They observed that the molecular scattering peak gains increasing spectral weights as the interaction is tuned from the Bardeen–Cooper–Schrieffer (BCS) regime to Bose–Einstein–Condensate (BEC) regime, a behavior distinct from the atomic scattering.^[12] Later, they investigated the Goldstone phonon mode and pair-breaking excitation at small q, revealing an interaction-induced suppression of the sound speed.^[13] Subsequent work demonstrated that the bosonic molecular excitations are sensitive to variations in the pairing gap.^[14] Biss et al. experimentally characterized the phonon dispersion across the BCS–BEC crossover using the Bragg spectroscopy.^[15] Theoretical studies have extensively explored the dynamical structure factors of three-dimensional (3D) Fermi gases.^[19–27] For two-dimensional (2D) superfluid Fermi gases, recent experiments using the two-photon Bragg spectroscopy measured the dynamical structure factor across the entire interaction strengths.^[28] Notably, Vitali et al. used the exact quantum Monte Carlo (QMC) method to investigate the spectral weight redistribution between the molecular and atomic excitations.^[29] Several theoretical works on other low-dimensional Fermi gases have investigated the dynamical excitations through the dynamical structure factor.^[30,31]

For the discrete ultracold atomic gases, the optical lattice systems generated by superimposing orthogonal standing waves provide an ideal platform for simulating crystalline physics.^[32,33] These systems are well described by the Bose–Hubbard or Fermi–Hubbard models.^[34–46] The attractive Fermi–Hubbard model, in particular, has drawn significant theoretical interest due to its relevance to the strongly correlated condensed matter systems.^[47–57] Experimentally, this model has been realized in cold atom systems.^[58–63] However, to date, no experimental measurement of the dynamical structure factor has been reported for 2D optical lattices. Vitali et al. numerically investigated the dynamical structure factor of the half-filled attractive Fermi–Hubbard model on a square optical lattice using the QMC simulations.^[64] Their work revealed both the low-energy two collective modes and higher-energy single-particle excitations along the high-symmetry Brillouin zone (BZ) directions. However, the relationship between the roton mode and the pairing gap on an optical lattice remains unclear.

In this paper, we theoretically investigate the doping-dependent dynamical excitations in 2D Fermi superfluids on optical lattices across weak to intermediate coupling regimes. Through detailed analysis of the dynamical structure factor, the key features of both the collective modes and the single-particle excitations are demonstrated, particularly focusing on the roton molecular peak at q = [π,π]. This peak provides a novel strategy for measuring the pairing gap in doped systems with the breaking of particle–hole symmetry, namely, the square of the pairing gap scales linearly with the area under this peak. Notably, this strategy extends to the SOC Fermi gases on an optical lattice, where the complex band structure makes it difficult to directly measure the pairing gap.^[65]

This paper is organized as follows. In Section 2, we use the equations of motion of the Green’s functions to solve the 2D Fermi–Hubbard model in mean field approximation, and self-consistently obtain the chemical potential and pairing gap. In Section 3, we introduce how to calculate the dynamical structure factor with random phase approximation (RPA). We display results of dynamical structure factor at half-filling and compare with the QMC results in Section 4. In Section 5, we present results when the system is away from half-filling, and discuss the hopping dependence of the sound speed and the dynamical excitations at a transferred momentum q = [π,π]. We check the doping dependence of dynamical structure factor in Section 6. Finally we give our conclusions and acknowledgment, and provide some calculation details in Appendix A.

An attractive Fermi–Hubbard model in 2D square optical lattices can be described by a Hamiltonian in spatial representation as follows:

where 〈ij〉 means the nearest-neighbor sites on the lattice. Ciσ†(Ciσ) is the creation (annihilation) operator of a particle with spin σ, hopping energy t and chemical potential μ at site i. The Hubbard energy U > 0 is the strength of the on-site two-body attractive interaction. In the following discussions, U is set to be the unit energy, while the lattice length a₀ is used as unit length. Within the mean field theory, the four-operator interaction Hamiltonian can be dealt into a two-operator one with the definition of pairing gap Δ = U〈C_i↓C_i↑〉. The pairing gap Δ can be chosen to be a real number in the ground state. Then the above Hamiltonian is displayed into a mean-field one, whose expression in momentum space reads

where ξ_k = −Ztγ_k–μ and γ_k = (cosk_x + cosk_y)/2. The nearest lattice number satisfies Z = 4 for 2D square lattice.

The above mean-field Hamiltonian can be solved by the equations of motion of the Green’s functions. Here we define the diagonal Green’s function G(k,τ−τ′)=−〈TτCkσ(τ)Ckσ†(τ′)〉 and off-diagonal one Γ†(k,τ−τ′)=−〈TτC−k↑†(τ)Ck↓†(τ′)〉. The diagonal Green’s function is related to the normal particle density and the off-diagonal Green’s function is related to the singlet Cooper pairing information. Their expressions in momentum and energy representation are given by

where Ek=ξk2+|Δ|2 is the quasiparticle spectrum. The chemical potential μ and pairing gap Δ are determined by self-consistently solving the density equation and pairing gap equation

We have set Boltzmann constant k_B = 1, and will consider a typical low temperature T/U = 0.01 (close to zero) in the following analysis.

Generally, the pairing gap Δ exhibits an inverse dependence on the hopping strength t. We plot Δ as a function of t for different particle densities n in Fig. 1. When t is large, Δ asymptotically approaches zero, indicating a phase transition from a superfluid to a normal state. Furthermore, Δ is significantly suppressed as n decreases. It should be noted that Δ under the mean-field theory is overestimated compared with the corresponding QMC results.^[64] For example, at half-filling (n = 1.0) with t/U = 0.25, the mean-field theory yields Δ/U = 0.345 while the QMC predicts Δ/U = 0.1825. Experimentally, the formation and spatial ordering of non-local fermion pairs in an attractive Fermi–Hubbard system are directly observed using two-species degenerate ⁴⁰K atomic gases.^[63]

In this section, we introduce the main idea of the random phase approximation (RPA). The RPA provides an approach to systematically investigate the dynamical excitations by incorporating quantum fluctuations beyond the mean-field theory.

In a superfluid system, there are four different densities, with the normal spin-up/spin-down n^1=〈Ci↑†Ci↑〉, n^2=〈Ci↓†Ci↓〉, and the two anomalous densities describing Cooper pairing n^3=〈Ci↓Ci↑〉, n^4=〈Ci↑†Ci↓†〉. These four densities are coupled with each other through the finite interaction, where any perturbation in one density induces the correlated fluctuations in other densities. Within the linear response theory, the small external perturbation potential V_ext and density fluctuations δn are connected with each other through the full response function χ, namely, δn = χV_ext.

The mean-field theory neglects the contribution from the fluctuation term of interaction Hamiltonian, failing to predict the dynamical excitations in interacting systems. In order to incorporate the quantum fluctuations,^[66–68] the RPA has been proven to be a reliable method for calculating χ beyond the mean-field level. In a 3D BCS–BEC crossover Fermi superfluid, the dynamical excitations obtained from RPA theory even quantitatively agree well with the experimental results.^[15,21,23] Similarly, the 2D RPA incorporating the density fluctuations can qualitatively reproduce the corresponding QMC data.^[30] These results indicate that the RPA approach can provide qualitatively reliable predictions for the dynamical excitations of a 2D lattice Fermi superfluid.

The main idea of the RPA is to treat fluctuation Hamiltonian as part of an effective external potential. This theory establishes the relationship between the full response function χ beyond mean-field level and its mean-field response function χ⁰,

Here G = σ₀ ⊗ σ_x is a direct product of unit matrix σ₀ and Pauli matrix σ_x.

The numerical calculation of mean-field response function χ⁰ is easy to gain, and its expression is given as

The dimension of χ⁰ demonstrates the coupling situation among the four density channels. These 16 matrix elements are obtained through the corresponding density–density correlation functions derived from the previously defined Green’s functions. Due to all possible symmetries of the system, only 6 of these matrix elements are independent, i.e., χ110=χ220, χ120=χ210=−χ330=−χ440, χ310=χ320=χ140=χ240, χ410=χ420=χ130=χ230. The symmetry of the matrix is closely related to the symmetry of the Green’s functions, such as Γ^†(k,ω) = Γ(k,ω) = Γ(k,−ω), which leads to χ120=χ210, where χ120=−〈TτC↑†(r,τ)C↑(r,τ)C↓†(r′,τ′)C↓(r′,τ′)〉. Based on Wick’s theorem, χ120=−Γ†(r−r′,τ−τ′)Γ(r′−r,τ′−τ). Their expressions are listed in Appendix A.

The total density response function χ_n is defined by χ_n ≡ χ₁₁ + χ₁₂ + χ₂₁ + χ₂₂, and its expression is given as

where

According to the fluctuation–dissipation theory, the density dynamical structure factor S(q,ω) is connected to the imaginary part of the density response function χ_n by

where q and ω are respectively the transferred momentum and energy. δ is a small positive number in the numerical calculation (usually we set δ = 0.003).

We first study the dynamical structure factor of 2D attractive Fermi–Hubbard model at half-filling (n = 1). By analyzing the density dynamical structure factor S(q,ω) under different transferred momenta q, both the collective excitations and single-particle excitations of Fermi atomic gases can be obtained. We calculate the energy and momentum dependencies of S(q,ω) and present its spectral weight distribution along the high-symmetry directions in the BZ for varying hopping strengths t/U, as shown in Fig. 2. In the low-energy region, S(q,ω) displays sharp peaks, signaling the presence of two gapless collective modes. The first mode, emerging from q = [0,0], increases nearly linearly with momentum along both [0,0] → [π,0] and [0,0] → [π,π]. This behavior is characteristic of a phonon mode originating from the spontaneously U(1) symmetry breaking of pairing gap (or order parameter). The slope of this mode at q = [0,0] defines the sound speed. The second mode is the roton mode appearing near q = [π,π]. Its emergence can be attributed to a global pseudospin SU(2) symmetry breaking.^[68,71,72] Owing to the particle–hole symmetry of the Hamiltonian (1) at n = 1, the fluctuations of superfluid and charge density wave (CDW) become degenerate, leading to the appearance of the gapless roton mode. As the momentum increases, the phonon mode gradually merges into the single-particle excitations, acquiring a finite expansion width because of the scattering processes with the single-particle excitations. The existence of these two collective modes has been confirmed by the QMC simulations.^[64]

In the high-energy regime, the dynamical excitations are dominated by the single-particle continuum band arising from the pair-breaking effect. The energy required to break Cooper pairs at a certain transferred momentum q, given by E_k+q + E_k, is determined by the quasiparticle spectrum E_k. This pair-breaking process generates a continuous excitation band, with its minimum energy threshold min[E_k+q + E_k] marked by a green dotted line in Fig. 2(b). This horizontal green line corresponds to the transferred energy ω = 2Δ, representing the minimum pair-breaking threshold and providing an experimental strategy for measuring the pairing gap. However, this measurement strategy becomes unreliable in systems with complex band structures (e.g., under SOC), where the threshold varies with q. As evident in Fig. 2, increasing t shifts the threshold toward the low energies due to the suppression of the pairing gap. Furthermore, compared with the RPA results, the QMC simulations reveal that the upper single-particle excitation branch has significantly lower energy, highlighting the RPA’s limitations in strong interacting systems. To address this, we introduce the quasiparticle coherent weight Z_F^[74–77] and a normal-state full Green’s function incorporating many-body interaction: g(k,ω) = 1/(ω – ξ_k – Σ(k,ω)), where the self-energy Σ(k,ω) is decoupled as Σ(k,ω) = Σ_e(k,ω) + ωΣ_o(k,ω). The Z_F is defined as ZF−1(k,ω)=1−Σo(k,ω). Under the static limit approximation, Z_F = Z_F(k = k_F,ω = 0) by taking the Fermi momentum k = k_F, simplifies the full Green’s function to g(k,ω)=ZF/(ω−ZFξk)=ZF/(ω−ξ¯k), where ξ¯k=ZFξk is the renormalized energy spectrum. For the non-interacting systems, Z_F = 1. Increasing the interaction strength reduces Z_F, lowering the upper branch energy. In superfluids, Z_F also suppresses the pairing gap, explaining the smaller pairing gap observed in QMC compared with the RPA results. In this paper, we do not discuss the effect of Z_F.

The band width of the single-particle excitations W increases with the hopping term t, specifically, W = 8t. Along [0,0] → [π,π] direction, S(q,ω) exhibits a distinct upper boundary, marked by the blue arrow in panel (b), described by

The lower contour, marked by the black arrow in panel (b), is obtained as

The red-dashed line represents the dispersion of E_d. Physically, E_d does not correspond to a collective mode but is analogous to the 1D lower boundary of E_d = 4tsin(q) at U = 0.^[69,70] In the Heisenberg model, E_d corresponds to the des Cloizeaux–Pearson (dCP) dispersion.^[70] Here, due to the pairing gap, the prefactor is reduced to 1.83, yielding E_d = 1.83sin(q). Notably, E_d(q) exhibits double periodicity compared with E_u. Both E_u and E_d arise from the single-particle excitations, which can be understood by ℏω_kq = ξ_k+q – ξ_k. Along [0,0] → [π,π], this simplifies to ℏω_kq = 8tsin(k + q/2)sin(q/2). When sin(k + q/2) = 1, the maximum excitation occurs, yielding ℏωkqmax=8tsin(q/2)=Eu. At n = 1, the Fermi momentum k_F = π/2. When k = k_F, the excitation energy reduces as ℏω_kq = 4tsin(q) = E_d.

Doping modifies the particle density and shifts the Fermi energy, thereby significantly impacting the dynamical excitations. The particle density is given by n = 1 – δ, where δ denotes the doping concentration. We analyze the energy and momentum dependencies of S(q,ω) at n = 0.8 (away from the half-filling). In Fig. 3, we present the contour plots of S(q,ω) along the high-symmetry directions in the BZ for varying hopping strengths. Three key differences emerge when comparing n = 0.8 with the half-filling, namely, the molecular excitations at q = [π, π], the split of E_d, and the doping-dependent variation of the sound speed. These differences are elaborated in the following three subsections.

Compared with Fig. 2, the roton mode minimum at [π, π] moves above the zero energy when n = 0.8. Away from the half-filling, the particle–hole symmetry of the system is broken, resulting in opening a gap in this roton mode. In optical lattices, this gapped roton mode remains well separated from the single-particle excitation band and is interpreted as the molecular Cooper-pair excitations here.^[71] To show this clearly, Fig. 4 displays S(q,ω) as a function of ω for varying t, with insets highlighting the atomic excitation regime. Obviously, S(q,ω) exhibits two distinct features. First, a sharp low-energy peak corresponds to the excitation of the bosonic molecules from a molecular condensate. Second, a broad high-energy single-particle excitation band arises from the atomic (particle–hole) excitations. As t increases (U decreases), the weight of the molecular peak decreases while the atomic excitation band increases quickly.

The spectral weight of the molecular peak can be quantified by integrating its intensity over frequency, yielding the peak area A_peak. In Fig. 5, we plot A_peak (pink dotted line) as a function of the hopping strength t, compared with the square of the pairing gap Δ² (green solid line). Our results reveal that A_peak has almost the same t dependence as Δ². Both quantities exhibit particularly large values at small t and decrease gradually as t increases from the intermediate to the weak coupling regime. This correlation indicates that Δ² governs the molecular excitation at q = [π, π], suggesting that Δ can be experimentally measured by detecting S(q = [π, π],ω). It is worth noticing that this measurement method is universal. For 3D cubic optical lattices, our calculations confirm that A_peak at the roton mode q = [π, π, π] remains proportional to Δ², demonstrating dimensional independence. Furthermore, this method extends to the SOC Fermi systems, where conventional methods for measuring Δ are hindered by the complex band structures.^[65]

Doping also significantly modifies the atomic excitations. Compared with the half-filling case in Fig. 2, the boundary E_d splits into two distinct branches upon doping. To elucidate this effect, we calculate S(q = [π/2,π/2],ω) as a function of ω for varying doping concentrations at t/U = 0.4 in Fig. 6.

In the low-energy region, a sharp peak emerges as a signature of collective mode excitations, while a broader higher-energy band corresponds to the atomic excitations. At n = 1.0 (red solid line), a characteristic peak at ω/U = 1.83 (marked by a black arrow) corresponds to E_d. However, for n = 0.8 (red dashed line) and 0.6 (blue dotted line), E_d splits into two branches: U_d and D_d, which also exist in the normal state. This splitting is not induced by the interaction term Hint=−∑k(Δ*Ck↓C−k↑+H.c.) in Eq. (1). Moreover, as doping increases, the energy gap between U_d and D_d increases. Along [0,0] → [π,π], the Fermi momentum is k_F = π/2 at n = 1.0 and decreases with doping, thereby altering the single-particle excitation spectrum.

The sound speed of the Goldstone phonon mode is defined by its low-energy dispersion in the limit ω → 0 and q → 0, c_s = ω/|q|. The sound speed depends sensitively on the interaction strength, parameterized by t/U. Figure 7 plots c_s as a function of t, revealing that c_s decreases with decreasing t from the intermediate to the weak coupling regimes. This trend is qualitatively consistent with experimental observation in ⁶Li Fermi gases.^[13,28] The qualitative behavior of c_s can be understood through its proportionality to the root-mean-square Fermi velocity 〈vF2〉, which will be introduced in Fig. 9.

Here we discuss the interplay between dynamical excitations and the doping concentration in the Hubbard model. Figure 8 displays S(q,ω) along the high-symmetry directions of the BZ for (a) n = 0.6 and (b) n = 0.4 with t/U = 0.4. Our results reveal two key doping-dependent features: the molecular peak at q = [π,π] moves to the higher energy, enlarging the roton gap; the sound speed exhibits a pronounced doping dependence. To clarify this, we plot c_s as a function of n in Fig. 9(a). Our results show that the sound speed exhibits a non-monotonic doping dependence. It initially increases and then decreases as n increases. This behavior can be qualitatively understood within a weak interaction theory,^[73] where c_s is closely related to the root-mean-square Fermi velocity, cs=〈vF2〉[1−UN(0)]/2, with N(0) denoting the density of states (DOS) at the Fermi energy. The DOS is computed as N(0)=(2π)−2∫d2kδ(ξk), and is proportional to the Fermi surface length or n. The Fermi velocity v_F = ∂ξ_k/∂k|_k=[kFx,kFy] is evaluated at the Fermi wave vector determined by ξ_k=[kFx,kFy] = 0, with k_F, μ and Δ obtained self-consistently for a given n (see Eq. (4)). We compute the averaged squared Fermi velocity 〈vF2〉=1N0∑kFvF2, where the sum runs over N₀ = 120 points along the Fermi surface. We plot 〈vF2〉 as a function of n in Fig. 9(b). It is shown that 〈vF2〉 has qualitative agreement with the doping dependence of c_s. At half-filling, v_F exhibits strong anisotropy across the BZ. As n decreases, the Fermi surface contracts, reducing the v_F anisotropy and aligning the system’s physical properties with those of continuum Fermi gases. In the low-momentum region, the cosine function can be expanded as cosk ≈ 1 – k²/2, yielding that the energy spectrum of the optical lattice can be approximated as ξk≈Zt(kx2+ky2)−tZ−μ, which matches the continuum form ξk=(kx2+ky2)/(2m)−μ. Moreover, as n → 0, N(0) → 0, and the sound speed simplifies to cs=〈vF2〉/2, consistent with ideal continuum Fermi gases.

In conclusion, the doping and hopping dependencies of the dynamical structure factor in the 2D attractive Fermi–Hubbard model were investigated using RPA theory. Two collective modes emerge: the phonon mode at small transferred momenta and the roton mode at transferred momentum q = [π,π] regime. The roton mode corresponds to Cooper pair molecular excitations. First, the area of the molecular excitation peak at q = [π,π] scales with the square of the pairing gap at fixed doping. This observation motivates a universal strategy for measuring the pairing gap, being applicable to both 2D/3D and SOC optical lattices. While this result is obtained from the weak coupling to the intermediate coupling regime within the RPA, we conjecture that it remains valid even for strong interactions. Because this strategy does not depend on dimension and SOC interaction. Furthermore, the relation between the roton mode and the pairing gap is confirmed by both Zhang’s theory and Dyke et al.’s Bragg experiments. It is also essential to verify this strategy beyond the RPA. Second, the atomic excitation peak splits into two branches when the system deviates from half-filling, with the splitting magnitude increasing with doping. Third, at a given doping, the sound speed is suppressed by interaction strength and governed by the root-mean-square Fermi velocity.

The mean-field response function χ⁰ of 2D interacting Fermi atoms in a square optical lattice is numerically calculated, and all 6 independent matrices elements of χ⁰ are

The corresponding functions in above equations Fk,q(1), Fk,q(2), Fk,q(3), Fk,q(4) are defined as

where

f(E_k) and f(E_k+q) are Fermi distributions.