Mottness induced superfluid phase fluctuation with increased density

Superfluid phase stiffness is known to suffer from low density, due to the canonical conjugation between carrier density and phase ϕΔnΔϕ ∼ ℏ [1]. Since the density fluctuation is very limited at low density, the phase must then fluctuate significantly, leaving only a small portion in the well-defined phase of the superfluid component. This quantum mechanical limitation on superfluid density provides a natural explanation [1, 2] for the observed reduction of the superconducting transition temperature T_c in the underdoped cuprates (cf δ < 15% of figure 1(a)). As shown schematically in figure 1(a), in the underdoped region (δ < 15%), T_c decreases as the phase stiffness is weakened at lower carrier density (gray dotted line). In that case the superconducting transition temperature T_c would be controlled by the phase stiffness, as opposed to the strength of the pairing.

Zoom In Zoom Out Reset image size

Unexpectedly, recent measurement [3] of penetration depth λ in the overdoped cuprates (δ > 15%) found a surprising reduction of superfluid stiffness upon increasing carrier density. Correspondingly, the low-temperature superfluid phase stiffness (∝λ⁻²) forms a dome shape against doping, qualitatively different from the simple proportionality to carrier density [4] expected in the pairing strength-limited BCS theory. Furthermore, the low-temperature phase stiffness even scales with T_c in the overdoped regime, obeying the same universal Uemura relation [5–7] well-known in the underdoped regime. Since this relation implies a phase coherence-limited superfluidity, this new data apparently reveals that even in the overdoped side, the dominant physics is still the phase fluctuation. Nonetheless, it is not obvious at all why superfluid stiffness should reduce at high carrier density.

In fact, several previous studies already raised similar doubts against the common belief of BCS-like amplitude fluctuation in the overdoped cuprates, and faced the same puzzle how the superfluidity can suffer at higher density. As shown in figure 1(b), angle-resolved photoemission spectroscopy (ARPES) observed [8, 9] a nearly doping independent superconducting gap Δ near momentum (π/2, π/2) over a wide doping range (∼7%–20%), in great contrast to the expected proportionality to T_c in the BCS theory. Similarly, figure 1(d) shows that even earlier, the observed momentum-dependence of superconducting gap-induced weight transfer deviates qualitatively from the expected ${d}_{{x}^{2}-{y}^{2}}$ form of the superconducting order parameter, also from the underdoped regime all the way to the overdoped regime. Both observations suggest more commonality between the underdoped regime and the overdoped one than previously expected. Consistently, VMC calculations [11–13], which ignore phase fluctuation, all found that the pairing still exists beyond 25% doping (black dashed line in figure 1(a)). In this case, the entire superconducting dome would be deep inside the pairing region of the phase diagram, and thus the demise of superconductivity in the overdoped region must be controlled by additional phase fluctuation of unknown origin.

So, the key scientific question is how it is possible to host such a strong superfluid phase fluctuation at such high carrier density? Particularly in the overdoped cuprates, how can the phase quantum-fluctuates even stronger with increased carrier density, as to wipe out superconductivity at roughly δ > 25%? Why is this critical doping level seemingly universal across different families of cuprates and the new nickelate superconductors [14, 15]? What determines this special doping level? What is the origin of the striking similarity between the underdoped and overdoped regime, especially near the quantum critical points (the end of the dome) δ ∼ 5% and 25%? Finally, what essential nature do these puzzles reveal about the unconventional high-temperature superconductivity (HT-SC) in the cuprates?

Here we show that all these important questions can be answered naturally in the scenario of the emergent Bose liquid (EBL) [16–21]. Essentially, the high-energy local Coulomb repulsion that prevents double occupation of carriers in each lattice site at low energy (so-called 'Mottness') leads to a rather large hard core of the superfluid carriers at even lower temperature. Such a large hard core would unavoidably cause a serious jamming in carriers' kinetic motion, and in turn damage the superfluid phase coherence at higher density, a trend opposite to the common lore. Below we demonstrate this general mechanism of stronger superfluid fluctuation at higher density via simple close-pack estimation and Gutzwiller approximation through statistical counting and variational Monte Carlo, all showing a consistent destruction of superfluidity at around 25%–30% carrier density. Specifically for the cuprates and nickelates, this classically intuitive mechanism explains the robust upper limit of doping level in the entire families of cuprate and nickelate superconductors. Combined with the persistent consistency of optical [21], transport [19, 21], quasi-particle (QP) [16, 18, 20], and superconducting [17] properties of the exact same model with corresponding experimental observations, our results reinforce the notion of a new paradigm that the HT-SC in the cuprates (and likely nickelates) is dominated by physics of bosonic superfluidity, as opposed to pairing-strength limited Cooper pairing.

The contents of this paper are organized as follows. Section 2 introduces the EBL model and the basic assumptions. In section 3, we discuss the main results of this study. Section 4 lists important characteristics of EBL model in comparison with other bosonic models. Section 5 gives the connection of our model to the case of the cuprate superconductors.

In the EBL model [16–21], the charge carriers are assumed to be mostly in a tightly bound two-body state involving the nearest neighboring atoms, below some energy scale, say 150 meV, due to various high-energy physics, for example, bi-polaronic correlation [22], two-dimensional short-range anti-ferromagnetic correlation [23, 24] and/or cancellation of topological spin current [25]. (For example, in cuprates this bosonic two-body state would consist of doped holes of opposite spin residing in neighboring Cu sites, using the half-filled Mott insulator with strong short-range AFM correlation as reference.) Correspondingly, the probability for the fermionic carriers in an unbound one-body state is relatively low below this energy scale. Therefore, most physical properties of the system are dominated by the emerged bond-centered bosonic carriers (cf solid ellipse in figure 2(a)) that reside in the corresponding bond lattice of the original fermionic lattice (black cross in figure 2(a)). Note that this assumption of EBL does not require a complete depletion of the probability for low-energy fermion (or fermionic spectral weight), but only its overall smallness. (More detailed considerations on low-energy fermions are available in appendix A.)

Zoom In Zoom Out Reset image size

In addition, in many strongly correlated materials, intra-atomic electronic repulsion is very strong, well above the energy scale of this effective low-energy model, resulting in an exclusion of double occupation [26–28] (so-called 'Mottness') in the low-energy sector of the fermionic carriers. This gives rise to an 'extended hard-core constraint' of EBL, namely forbidden occupation of the vicinity of a boson by another, as illustrated by the shadowed area in figure 2(a).

The kinetic process of the EBL is dominated by the pivoting motion [16] in the emerged bond-centered lattice. For the specific case of a square lattice of fermions, for example the Zhang–Rice singlets [29] in the cuprates, the bond lattice is a checkerboard lattice in figure 2(a) with two orbitals in the unit cell corresponding to the vertical and horizontal bonds:

$$\begin{equation}H=\sum\limits _{\langle l,{l}^{\prime }\rangle }{\tau }_{l{l}^{\prime }}{b}_{l}^{{\dagger}}{b}_{{l}^{\prime }}\quad (+\;\text{constraint})\end{equation} \tag{ 1 }$$

where b_l denotes the annihilation of a boson located at bond l, and τ_ll' = τ' or τ'' is the fully dressed kinetic process involving the nearest and second-nearest neighboring bonds. In this study, we employ the exact same set of parameters τ' and τ'' previously extracted from ARPES dispersion data [16, 30, 31] (cf appendix C).

It is important to emphasize that the low-energy effective Hamiltonian of EBL is under the above-mentioned extended hard-core constraint, which prevents phase separation related to clustering of the emergent bosons. At low temperature, it also helps to establish the superfluid phase stiffness in the superfluid phase. As demonstrated below, this constraint is the key factor that generates the surprising physical behavior of interest in this study.

Now, since the focus of our study is on the phase stiffness of superfluidity, which involves only the lowest energy excitation, it is advantageous to further integrate out irrelevant (higher-energy) kinetic degree of freedom of equation (1). Figure 2(a) shows that within the bosonic band structure, the lower-energy band is dominated by d-wave symmetry due to the positive τ'. We proceed to decouple the d-wave boson from the s-wave one via a canonical transformation. Practically, this is performed approximately via construction of the d-wave bosonic Wannier orbital corresponding to the lower-energy band in figure 2(a), and then representing the effective Hamiltonian in this decoupled subspace. (See appendix B for more details.) Essentially, this is equivalent to assuming the formation of a well-defined lowest energy QP with a d-wave structure responsible for the superfluidity.

Finally, we arrive at the resulting lowest-energy one-orbital Hamiltonian that controls directly the physics of phase fluctuation,

$$\begin{equation}\tilde{H}=\sum\limits _{\langle i,{i}^{\prime }\rangle }{\tilde{\tau }}_{i{i}^{\prime }}{d}_{i}^{{\dagger}}{d}_{{i}^{\prime }}\quad (+\;\text{interactions}\;\,\&\,\;\text{constraint})\end{equation} \tag{ 2 }$$

corresponds to a larger d-wave Wannier orbital, denoted by ${d}_{i}^{{\dagger}}$ creation operator, now centered at sites i in the original square lattice (cf: figure 2(b)). Note that due to the shift of the Wannier center, the d-wave boson is not strictly local, but instead with a power-law decay similar to the Zhang–Rice singlet [29]. Correspondingly, the d-wave boson can now hop beyond the first and second neighbors with further renormalized ${\tilde{\tau }}_{i{i}^{\prime }}$ given in table C2.

Most importantly, these low-energy superfluid related d-wave bosons now acquire an even more extended hard-core constraint (shadowed area in figure 2(b)) over 16 bonds. As highlighted in figure 2(c), this constraint forbids occupation of 12 surrounding sites by other low-energy d-wave bosons.

Concerning the main puzzles of the strong phase fluctuation and diminishing superfluidity density at high carrier density, we will now show that the EBL also provides a simple resolution through its 'extended hard-core constraint'. Such a large extended hard-core implies that the essential kinetic processes will be easily suppressed at moderate density due to blocking and jamming between the low-energy d-wave bosons. Figure 2(d) shows a simple estimation of the special case, in which each of the bosons reach minimum non-zero mobility on average, being able to hop back and forth between two sites. The corresponding density, 2 holes (1 boson)/8 atomic sites = 25%, marks the approximate maximum doping level of superfluidity, above which the d-wave bosons become nearly impossible to move and thus unable to maintain phase coherence. Interestingly, this 25% coincides very well with the experimentally observed end of the superconducting dome in overdoped cuprates in general [32]. Therefore, this mechanism offers a simple and natural explanation of the key puzzle of strong phase fluctuation in the overdoped cuprates.

We employ the well-known Gutzwiller g-factor approximation [33, 34] ${\tilde{\tau }}_{i{i}^{\prime }}\to {g}_{i{i}^{\prime }}{\tilde{\tau }}_{i{i}^{\prime }}$ as the simplest way to capture this renormalization of the kinetic process. The g-factors are calculated via two numerical approaches (see appendix D). First, we count statistically the probability of each hopping under the constraint at various doping level. Figure 3(a) shows that the resulting g-factors are approximately g ∼ 1 − δ/0.3. Second, we calculate the g-factor using VMC based on a noninteracting BEC wavefunction under the constraint. The resulting g-factors in figure 3(b) resemble figure 3(a) very well. As expected, both results decrease rapidly and become rather small beyond 25% doping and diminish around 30%.

Zoom In Zoom Out Reset image size

We now demonstrate the 'dome' shape of T_c with our simple picture. Estimated from the standard BEC [35] of a uniform boson gas, the condensation temperature T_c ∝ n^2/3/m* is a simple function of density n and effective mass m*. Obviously, following $\tilde{\tau }$, 1/m* and T_c are also renormalized by g. Figure 3(c) plots T_c ∝ δ^2/3 (in black) and the renormalized T_c → gT_c (in red). Indeed gT_c is strongly suppressed at higher doping and eventually vanishes around 30%, where the average m* diverges. This δ ∼ 30% upper bound of superfluidity is in excellent agreement with Tl₂Ba₂CuO_6+δ [36] and well consistent with the common $\sim 25\%$ limit of typical cuprates superconductors. If one further incorporates the previously proposed doping dependent bosonic mass that diverges at around 5% due to a level crossing [17], the gT_c in figure 3(d) reproduces the experimentally observed dome shape quite well. In short, the phase fluctuation can indeed grow in the overdoped regime when the suppression of kinetic processes overcomes the increasing density.

Similarly, the low-temperature limit of the phase stiffness λ⁻² would also suffer from the jamming reduced kinetic process. Phase stiffness λ⁻² ∝ n_s /m* is proportional to the superfluid density n_s . Near the above-mentioned mass divergence dictated quantum phase transition points, ${n}_{s}\propto {(1/{m}^{\ast })}^{\beta }\delta$ must vanish with 1/m* with a exponent β > 0, leading to ${\lambda }^{-2}\propto {(1/{m}^{\ast })}^{1+\beta }\delta$. Figure 4(b) shows that this phase stiffness indeed forms a dome shape, vanishing at around δ ∼ 5% and $\sim 30\%$. Particularly, due to the extended hard-core constraint of the EBL, the phase stiffness reduces in the overdoped regime against the growing density, as observed in recent experiments [3, 4]. (Our results' deviation from experiment merely reflects the above-mentioned overestimation of phase coherence at high doping in our simple Gutzwiller treatment.)

Zoom In Zoom Out Reset image size

Furthermore, the effective mass divergence of EBL has an important physical consequence in the relation between λ⁻² and T_c [17]. Estimated from T_c ∝ δ^2/3/m* for a free boson gas, one finds ${\lambda }^{-2}\propto {T}_{\mathrm{c}}^{1+\beta }$. Therefore, distinct from the simple linear relation expected from standard phase fluctuation scenario, λ⁻² vs T_c must show a zero slope as T_c approaches zero (since β > 0). This provides a natural explanation of the super-linear relation in the experimental observations shown in figures 4(a) and (c). In fact, with β = 1 this formula describes very well the observed relation in the entire low-T_c region in both underdoped and overdoped sides. Not surprisingly, the resulting scaling ${\lambda }^{-2}\propto {T}_{\mathrm{c}}^{2}$ agrees with that of the 3 + 1D XY-mode that describes quantum phase fluctuation of superfluidity upon suppressed kinetic energy [38, 39].

Such an apparent 'symmetry' near the quantum phase transition points between the underdoped and overdoped sides is puzzling within the current lore, given that one is supposedly governed by fluctuation in the phase and the other in the amplitude. Especially, the observed physical properties of the normal-state above T_c in these two regions are quite distinct [32, 40, 41]. In our EBL picture, the contrast in non-superfluid properties is easily understood from the rising importance of the jamming effect discussed above. Concerning the higher-energy features, figures 4(d)–(f) show that in the overdoped region the renormalized s-wave (in blue) and d-wave bosons (in red) suffer a significant renormalization, quite distinct from the underdoped and optimally doped region. Yet, the renormalized low-energy coherent d-wave bosonic band near the underdoped (d) and the overdoped transition points becomes heavier (flatter) in a similar fashion (through very different microscopic mechanism though). In other words, the low-energy phenomenon of superfluidity experiences a similar loss of phase stiffness due to the effective mass divergence despite the distinctly different underlying higher-energy physics.

From this perspective, the non-superconducting 'normal state' above 25% doping should have interesting properties. While energetically most favorable locally, the pure d-wave boson would suffer from severe jamming and loss of kinetic energy. It would thus tend to morph into a p-wave boson whose cigar shape allows them to remain mobile at a much larger doping (see appendix E). Therefore, the EBL is expected to behave as an unconventional Bose metal consisting of fluctuating d- and p-wave bosons. Alternatively, since at very high doping, the electronic correlation will eventually become weaker, the emergent bosons might lose their binding strength and start to decompose into weakly bound fermions. Concerning this scenario, we note two very relevant experimental observations: (1) robust paramagnon dispersion at 40% doping [23] indicating that the local magnetic correlation remains strong and (2) linear resistivity at 30% doping [3] indicating that the normal states is still not a fermionic QP based 'normal metal'. Both observations suggest that decomposition of the emergent boson is not immediate after the disappearance of superconductivity.

It is also interesting to realize that in every energy scale our EBL experiences important influence of the 'Mottness' of the underlying doped holes through the suppression of their double occupation at each Cu site due to large local interaction. First, with the help of kinetic energy, it lays the foundation of near neighboring AFM [42] and bi-polaronic [22] correlations that provide a strong tendency to bind doped holes into pairs [26, 43]. Then, it forces the bound pairs to occupy two atomic sites and thus indirectly establishes d-wave form [17]. Finally, it produces the extended hard-core constraint of the EBL that ultimately induces jamming of low-energy d-wave bosons at high enough density. This jamming in turn suppresses the superfluid phase coherence especially with higher density.

Since the dominant physics of jamming inherits from the Mottness at very high energy scale, details of the intermediate- and low-energy scale is therefore not as essential, as long as the generic assumptions of our EBL model remains valid. Interestingly, the recently discovered Nd_1−x Sr_x NiO₂ also demonstrates a superconducting critical doping at around 25% [14, 15], similar to our analysis above. This strongly suggests that the nickelate superconductors can also be described by an EBL, with strong jamming-induced phase fluctuation in the overdoped regime. Future experimental verification of this prediction would be highly valuable.

In short, we show that EBL built from carriers under strong intra-atomic repulsion (Mottness) provides a generic mechanism that reduces superfluid stiffness at higher density due to jamming of its kinetic process at around 25%–30% carrier density. This mechanism naturally explains the recent observation of strong phase fluctuation in the overdoped cuprates [3, 44] and likely to be found also in the nickelates. Assuming that this mechanism is indeed responsible for the puzzling observation in the cuprates, it further suggests that in the entire doping range of the superconducting dome, the unconventional superconductivity can be described by bosonic superfluidity for example of an EBL, without resorting to a crossover into BCS-like amplitude-fluctuating descriptions.

The EBL model assumes a high-energy binding between the nearest neighboring carriers such that the unbinding processes can be integrated out from the low-energy Hilbert space. Therefore, the bosons in EBL model are in a tightly bound two-body state of fermionic carriers of the nearest neighboring position (instead of momentum) with energy spanning far away from the chemical potential (since the binding is assumed strong). Most importantly, with strong short-range magnetic correlation in the background, the two-body state is not necessarily a spin singlet. Due to their emergent nature, they contain the following key characteristics that distinguish them from other boson-like carriers proposed in the literature, such as bosonized Cooper pairs [45–47], resonating valence bond [48–50]:

• (a)  
  They are in a tightly bound two-body state of fermionic carriers (not necessarily a spin singlet, given a strong short-range magnetic correlation).
• (b)  
  They consist of fermionic state with energy spanning far away from the chemical potential, since the binding is assumed strong.
• (c)  
  They are small in size and have a well-defined location, with underlying fermions residing in the nearest neighboring position (instead of a fixed total momentum as in a Cooper pair).
• (d)  
  They have a rod-like shape with two ends at nearest neighboring atoms.
• (e)  
  Correspondingly, they reside in the bonds between atomic sites.
• (f)  
  The bosonic lattice is therefore the bond lattice, different from the underlying fermionic lattice.
• (g)  
  For systems of 2D or higher dimension, the multiple directions of bonds imply multi-orbital nature of EBL.
• (h)  
  The bosons in EBL carry a limited pivoting motion, moving one end at a time. This paves ways to frustration in the kinetic energy and can therefore heavily renormalize the kinetic energy of the first bosonic band (to $\sim 30$ meV for the cuprates).
• (i)  
  They inherit an extended hard-core constraint from the large on-site repulsion between the underlying fermionic carriers.

Discovered almost thirty years ago [51], the exotic phenomenon of HT-SC in the cuprates still remains puzzling to researchers. Conventional superconductivity, a state of the matter that shows no resistance in conducting current, is well described by the standard 'BCS' theory [52] via weakly bound 'Cooper pairs' that fluctuate in amplitude. On the other hand, the HT-SC in cuprates shows qualitatively different behavior. For example, in the weakly hole doped ('underdoped') region, the isotope effect increases dramatically and yet the corresponding superconducting transition temperature T_c decreases [53] instead. The observed superconducting gap Δ₀ in the cuprates is typically significantly larger than the canonical value of twice the transition temperature $\sim 2{k}_{\mathrm{B}}{T}_{\mathrm{c}}$ in BCS theory [54, 55]. During the phase transition, the measured specific heat shows two clear kinks 10 K apart [56], qualitatively different from a standard second-order phase transition from the BCS theory. Furthermore, in the underdoped region, the low-temperature specific heat shows no T² contribution expected from the observed d-wave quasiparticles, but only a dominant T³ instead [56]. In addition, the observed upper critical field H_c2 does not saturate at low temperature and sometimes even exceeds the Pauli limit [57]. These qualitative distinct features indicate clearly that the HT-SC in the cuprates is of a different nature.

The assumptions of the EBL model are actually quite applicable given the strong high-energy spin and lattice correlations observed in the cuprates [22, 24, 42]. Particularly, the EBL model is in good consistency with the recent observations of bosonic behaviors in the cuprates [58–60], and is obviously the most intuitive one to account for the observed dominant role of phase fluctuation in the superconducting phase transition of the cuprates [1, 3].

Below we revisit several recent studies that demonstrated persistent agreement of optical [21], transport [19, 21], QP [16, 18, 20], and superconducting [17] properties of EBL with corresponding experimental observations. The fact that these successful applications of EBL are all via the same Hamiltonian with the same set of parameters strongly supports the scenario in which charge-related physics in the cuprates can be understood intuitively by a simple EBL.

5.1. Doping dependent superconducting gap

Previously this EBL model was shown to offer a simple explanation for the disappearance of superfluidity in cuprates at δ ∼ 5%, as the effective mass of the emergent boson diverges [17]. (The same study also suggested a scenario for the lack of superfluidity below 5% doping.) In a related study, a second kind of 'superconducting gap' was found in the QP spectrum, resulting from coherent kinetic scattering against the Bose–Einstein condensation (BEC) of the EBL [16]:

$$\begin{equation}{{\Delta}}_{\mathbf{k}}(T)\sim {f}_{d}(k)\cdot ({\varepsilon }_{k}-{\varepsilon }_{0})\sqrt{{n}_{0}(T)}\end{equation} \tag{ 3 }$$

where ${f}_{d}(k)=\frac{1}{2}[\mathrm{cos}({k}_{x})-\mathrm{cos}({k}_{y})]$ gives the momentum dependent d-wave form factor of the order parameter and n₀(T) the temperature-dependent condensation density. ɛ_k − ɛ₀ denotes the QP energy measured from the low-energy band center, and signifies the kinetic origin of the gap (as opposed to the pairing potential). Beyond the 'Fermi arc', where ɛ_k moves from the chemical potential toward the band center, it provided additional momentum dependence [16] that reproduced nicely the ARPES measurements [10] at various doping levels (cf figure 1(c)).

We first observe that this model actual 'predicted' naturally the above-mentioned nearly doping independent superconducting gap in ARPES measurements [8, 9]. Since the QPs on the Fermi arc are at the chemical potential μ, ɛ_k − ɛ₀ is given by μ − ɛ₀ ∝ 1 − 2δ upon hole doping into the quasi-2D singlet band of the cuprates. Away from the ends of the dome, in the absence of strong quantum fluctuation, the condensation density n₀(T = 0) is roughly proportional to the superfluid density, n_s (T = 0) ∝ T_c. Together, this gives a weakly doping dependent near-node gap scale ${V}_{{\Delta}}\propto (1-2\delta )\sqrt{{T}_{\mathrm{c}}}$, very similar to the experiments shown in figure 5, but in great contrast to the stronger dome shape of T_c. Particularly, notice that the 1 − 2δ factor shifts the maximum of V_Δ to a significantly lower doping δ ∼ 0.12, compared to T_c. Similar to the above-mentioned deviation from d-wave form of the gap, this again reflects the kinetic nature of the anomalous scattering gap absent in typical weak-coupling pairing scenario.

Zoom In Zoom Out Reset image size

5.2. Transport, quasiparticle and superconducting properties

In summary, within the scenario EBL with extended hard-core constraint, we propose a generic mechanism to explain the puzzling reduction of superfluid phase stiffness and transition temperature at high carrier density. Essentially, the strong intra-atomic repulsion (the Mottness) leads to a large extended hard-core constraint for the emergent bosons that can be easily jammed in their motion at high density and thereby host weakened superfluid stiffness. Our results are well consistent with the experimental observations in some superconductors, such as the cuprates and the nickelates. Furthermore, together the persistent consistency of the EBL with many key spectral and transport observations of the cuprates, this new finding adds credibility to the notion that EBL model can provide a good description for the charge channel of the low-energy electronic structure of the cuprates.

Note added: after initial submission of our manuscript, the shot noise experiment proposed in our original manuscript was performed [59], whose data indicated indeed a 2e-charge quanta, subject to further independent experimental confirmation. In the final version of this manuscript, we therefore remove this suggestion and instead incorporate this new development in the reference.

This work is supported by National Natural Science Foundation of China (NSFC) under Grant Nos. 12274287 and 12042507, and Innovation Program for Quantum Science and Technology No. 2021ZD0301900. FY acknowledges support from NSFC No. 12074031 and 11674025.

All data that support the findings of this study are included within the article (and any supplementary files).

It is a common question whether one is allowed to employ a bosonic picture like the EBL model to systems with experimentally observed low-energy fermionic QPs, such as the cuprates [40, 62]. As shown in figure A1(c), if one would consider only the BEC limit corresponding to an overwhelming binding energy E_B, the low-energy fermions indeed should be completely depleted with a clean gap of the scale of E_B − W around the chemical potential in the fermionic one-body spectral function. (Here W denotes the bandwidth of unbound fermions.) However, in general, with intermediate binding strength (E_B slightly larger than W/2) illustrated in figure A1(b) [63] one often finds a small probability of the unbound fermions (in blue) being occupied and forming the Fermi surface. Even though, fermions in the system are actually predominately part of the well-defined bound bosonic state, as reflected by the larger spectral weight in red below the chemical potential. In this scenario, the rare gapless fermionic QPs' contributions to low-energy physical properties are generally overwhelmed by those of the more likely bosonic carriers. Thus, a bosonic description of the low-energy physics, such as the EBL model used in this study, would be the most convenient even in the presence of rare gapless unbound fermions.

Zoom In Zoom Out Reset image size

Importantly, in contrast to the weak influence of the rare unbound fermions on the bound fermion, the latter is however expected to have significant impact on the former. Particularly at low temperature, contrary to the depletion of typical bosonic excitations such as phonon and magnon, the number (or more precisely the probability) of these bound bosons is basically fixed. Therefore, scattering against the bosonic carrier can lead to many unusual non-Fermi-liquid behaviors of the rare unbound fermions that cannot result from phonon or magnon, such as kink structure [18], Fermi arc and pseudogap [20] (cf section 5 of the main text). In essence, the unbound fermions can sensitively reflect the influence of bosonic carriers and serve as good 'probes' to the more essential bosonic carriers.

Since the focus of this study is on the superfluid stiffness, our main interest is the lowest-energy fluctuation around the superfluid. We therefore aim to construct a low-energy QP whose local structure absorbs the influence of the high-energy orbitals. Conceptually, this is similar to integrating out τ' to construct a local Wannier orbital with a procedure similar to that of the Zhang–Rice singlet, where the nearest neighbor hopping processes between two oxygen are integrated [29]. As long as the energy of resulting local object is lower enough than other state, this approximation is controlled.

In practice, it is easily achieved numerically by simply constructing the d-wave Wannier function from the lower-energy band in figure 2(a). The resulting energy of d-wave boson is ɛ_d = τ'' − 2τ' = −41.2 meV which is lower than the other states, such as p-wave boson with energy ɛ_p = −τ'' = −25.8 meV (with the parameter τ' = 33.5 meV, τ'' = 25.8 meV, see table C1). With the Wannier orbital, the lowest-energy effective one-band Hamiltonian equation (2) can be obtained by representing equation (1) within the subspace of d-wave boson, similar in spirit to the construction of the t–J model from the Zhang–Rice singlet [29].

The key physics of interest in this study, namely enhanced superfluid phase fluctuation at higher density, is dominated by the extended hard-core constraint inheriting from the high-energy physics of Mottness. Therefore, it is completely independent of the parameters of equation (1). Nonetheless, this section provides the parameters for interested readers.

Table C1 lists the parameters of equation (1) used in this study. These parameters were initially extracted [16] from ARPES experiments [30, 31] on La_2−δ Sr_δ CuO₄, and kept fixed in all the studies [16–21, 61] using the same EBL model. Other families of cuprates might have slightly different parameters, but we do not expect any significant deviation.

Table C2 gives examples of ${\tilde{\tau }}_{i{i}^{\prime }}$ of the d-wave boson at δ = 15%, corresponding to the band in figure 2(b). They are derived from integrating out the s-wave boson corresponding to the higher-energy band in figure 2(a).

Table C2. Examples of kinetic parameters ${\tilde{\tau }}_{i{i}^{\prime }}$ of d-wave boson in equation (2) for δ = 15%.

i − i'	(1,0)	(1,1)	(2,0)	(2,1)	(2,2)
${\tilde{\tau }}_{i{i}^{\prime }}$ (meV)	−5.19	−6.02	3.61	1.20	−0.24

D.1. Statistical counting of Gutzwiller g-factor

D.2. Variational Monte Carlo (VMC) calculation of Gutzwiller g-factor

Due to the extremely large many-body Hilbert space of bosonic systems, we can only afford to evaluate the Gutzwiller g-factor using an approximate wave function. Luckily, as argued above, the results should be quite insensitive to the choice of wave function.

We calculate the ratio between the expectation values of the hopping term $\langle {d}_{i}^{{\dagger}}{d}_{i\prime }\rangle$ in the 'Gutzwiller-projected' BEC state and that in the mean-field BEC state as the Gutzwiller g-factor, g_ii',

$$\begin{equation}{g}_{ii\prime }=\frac{\langle \mathrm{B}\mathrm{E}\mathrm{C}\left\vert {P}_{\mathrm{G}}{d}_{i}^{{\dagger}}{d}_{i\prime }{P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\rangle }{\langle \mathrm{B}\mathrm{E}\mathrm{C}\left\vert {d}_{i}^{{\dagger}}{d}_{i\prime }\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\rangle }.\end{equation} \tag{ D.1 }$$

Here, the P_G represents the extended hard-core constraint imposed on the bosonic hole-pairs, and the BEC mean-field state |BEC⟩ is expressed by the following normalized formula,

$$\begin{equation}\vert \mathrm{B}\mathrm{E}\mathrm{C}\rangle \equiv \frac{1}{\sqrt{N!}}{({d}_{\mathbf{k}=0}^{{\dagger}})}^{N}\vert 0\rangle =\frac{1}{\sqrt{N!}}{\left(\frac{1}{\sqrt{M}}\sum\limits _{i=1}^{M}{d}_{i}^{{\dagger}}\right)}^{N}\vert 0\rangle ,\end{equation} \tag{ D.2 }$$

where M and N represent the total site number and boson numbers respectively, with N = Mδ/2.

The denominator of equation (D.1) can be easily obtained as,

$$\begin{equation}\begin{aligned}\hfill \langle \mathrm{B}\mathrm{E}\mathrm{C}\left\vert {d}_{i}^{{\dagger}}{d}_{i\prime }\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\rangle & =\frac{1}{M}\sum\limits _{\mathbf{k}}{\text{e}}^{-\text{i}\mathbf{k}\cdot ({\mathbf{R}}_{i}-{\mathbf{R}}_{i}\prime )}\langle \mathrm{B}\mathrm{E}\mathrm{C}\left\vert {d}_{\mathbf{k}}^{{\dagger}}{d}_{\mathbf{k}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\rangle \hfill \\ \hfill & =\frac{1}{M}\langle \mathrm{B}\mathrm{E}\mathrm{C}\left\vert {d}_{0}^{{\dagger}}{d}_{0}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\rangle =\frac{N}{M}=\delta /2.\hfill \end{aligned}\end{equation} \tag{ D.3 }$$

The numerator of equation (D.1) can be evaluated as,

$$\begin{equation}\langle \mathrm{B}\mathrm{E}\mathrm{C}\left\vert {P}_{\mathrm{G}}{d}_{i}^{{\dagger}}{d}_{i\prime }{P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\rangle =\sum\limits _{\alpha }{\left\vert \left\langle \alpha \left\vert {P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\right\rangle \right\vert }^{2}\sum\limits _{\beta }\left\langle \alpha \left\vert {d}_{i}^{{\dagger}}{d}_{i\prime }\right\vert \beta \right\rangle \cdot \frac{\left\langle \beta \left\vert {P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\right\rangle }{\left\langle \alpha \left\vert {P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\right\rangle }\equiv \sum\limits _{\alpha }{P}_{\alpha }{B}_{\alpha },\end{equation} \tag{ D.4 }$$

where ${P}_{\alpha }\equiv {\left\vert \left\langle \alpha \left\vert {P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\right\rangle \right\vert }^{2}$ represents the weight of each configuration $\left\vert \alpha \right\rangle$ in the Gutzwiller-projected wave function P_G|BEC⟩, and ${B}_{\alpha }={\sum }_{\beta }\left\langle \alpha \left\vert {d}_{i}^{{\dagger}}{d}_{i\prime }\right\vert \beta \right\rangle \cdot \frac{\left\langle \beta \left\vert {P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\right\rangle }{\left\langle \alpha \left\vert {P}_{\mathrm{G}}\right\vert \mathrm{B}\mathrm{E}\mathrm{C}\right\rangle }$ represents the measurement for the configuration $\left\vert \alpha \right\rangle$. The weight P_α and the measurement B_α for any configuration |α⟩ are easily obtained. From equation (D.2), one easily finds that the weights P_α for all the configurations |α⟩ which are permitted by the extended hard-core constraint are equal, while those for the constraint-prohibited configurations are zero. The value B_α of any constraint-permitted configuration |α⟩ is chosen as follows: defining $\vert \beta \rangle ={d}_{i\prime }^{{\dagger}}{d}_{i}\vert \alpha \rangle$, if configuration |β⟩ is permitted by the extended hard-core constraint then B_α = 1; otherwise B_α = 0. These formulae provide an appropriate start-point for the following Monte-Carlo calculations.

In the Monte-Carlo calculation, we start from an arbitrary configuration |α₁⟩. We then randomly select a particle in configuration that let the particle hop to any other hard-core-constraint-permitted position on the lattice to obtain a second configuration |α₂⟩. This completes an update. Continuing the updates, we obtain a configuration series |α₁⟩, |α₂⟩, |α₃⟩, .... After about 10⁴ Monte-Carlo steps thermalization can be realized. In successive Monte-Carlo steps, one begins to extract the measurement B_α for the configuration |α⟩, and to average these B_α to obtain the numerator defined by equation (D.4). To avoid auto-correlation, we take one measurement after each N Monte-Carlo steps. About 10⁵ measurements are performed to attain convergence. The final result for the Gutzwiller g-factor is given by equations (D.1), (D.3) and (D.4).

Since we are interested in the superfluid phase stiffness in this study, our main study focuses only on the jamming of the d-wave emergent bosons. In this section we briefly discuss the possible phase beyond the overdoped superfluid regime, where the jamming of the d-wave emergent bosons costs them too much kinetic energy.

Let us first recall that at low density, the freely propagating d-wave bosons gain the most kinetic energy, 2ɛ_d = 2τ'' − 4τ', compared to those of other local symmetries. This energy can be approximately split into two contributions: locally forming a d-wave form (a four-site Wannier function in figure 2(b)) and its propagation in the system. The former contribution can be easily estimated by solving a local four-site problem and turns out to be ɛ_d , exactly half of the total kinetic energy. Similarly, a p-wave emergent boson would have a total kinetic energy 2ɛ_p = −2τ'', half of which is also gained by forming local p-wave symmetry. Therefore, when τ' > τ'', 2ɛ_d < 2ɛ_p and d-wave form is dominant.

As discussed in the main text, as the density increases upon higher doping, due to the large hard core the d-wave boson starts to suffer from a significant jamming effect. Represented by the Gutzwiller factor g_d , a crude estimation of the jamming-induced renormalized kinetic energy can be made: (1 + g_d )ɛ_d . Similarly, for the p-wave boson, the energy reduces to (1 + g_p )ɛ_p .

Interestingly, as shown in figure E1(a) due to the narrower shape of the p-wave bosons, they naturally suffer a lot less from the jamming effect than the d-wave bosons. Consequently, the p-wave boson can survive a much larger doping, approximately δ = 50% based on the estimation in figure E1(b). Indeed, figure E1(c) shows a much slower decay in g_p against δ from a simple statistical counting. Therefore, one can expect that at high-enough doping, the d-wave symmetry should no longer be dominant, since the bosons would be forced to adapt to p-wave whenever they need to pass each other in short distance. Figure E1(d) shows that for a fixed τ' and τ'' (from optimal doping), (1 + g_p )ɛ_p will eventually become lower than (1 + g_d )ɛ_d . Therefore, assuming that the EBL remains intact beyond the superconducting dome, the system should become a Bose metal consisting of a mixture of p- and d-wave bosons. Notice that such a Bose metal should experience a different structure of main phase fluctuation compared with the recently proven [19] homogeneous Bose metal in the underdoped region, and is an interesting candidate for another homogeneous Bose metal. Such a strong phase fluctuating EBL should present non-Fermi liquid behaviors as well, as found experimentally in references [3, 23]. The emergence of p-wave boson is also consistent with the observation of a similar direction of the charge order found in beyond underdoped and beyond overdoped systems [64].

Zoom In Zoom Out Reset image size