Skip to main content

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

Monte Carlo study of the pseudogap and superconductivity emerging from quantum magnetic fluctuations

Abstract

The origin of the pseudogap behavior, found in many high-Tc superconductors, remains one of the greatest puzzles in condensed matter physics. One possible mechanism is fermionic incoherence, which near a quantum critical point allows pair formation but suppresses superconductivity. Employing quantum Monte Carlo simulations of a model of itinerant fermions coupled to ferromagnetic spin fluctuations, represented by a quantum rotor, we report numerical evidence of pseudogap behavior, emerging from pairing fluctuations in a quantum-critical non-Fermi liquid. Specifically, we observe enhanced pairing fluctuations and a partial gap opening in the fermionic spectrum. However, the system remains non-superconducting until reaching a much lower temperature. In the pseudogap regime the system displays a “gap-filling" rather than “gap-closing" behavior, similar to the one observed in cuprate superconductors. Our results present direct evidence of the pseudogap state, driven by superconducting fluctuations.

Introduction

Even though unconventional and high-Tc superconductivity arises in a diverse set of materials, many of them share similar features in their phase diagram. One prominent feature is a superconducting (SC) dome, which emerges near the termination point of a non-SC phase with either spin or charge order. The second feature is anomalous transport and non Fermi-liquid (nFL) behavior around the putative quantum critical points (QCP). These features have led to the proposal that soft quantum-critical fluctuations of the order parameter serve as the source for the universal behavior and mediate singular interaction that gives rise to superconductivity with nontrivial pairing symmetry, strange metal behavior, and intertwined orders.

In many unconventional superconductors, most notably the cuprates, there is a third feature: the “pseudogap(PG)", a gap-like feature in the fermionic spectrum above the SC phase. Despite decades of investigation, the origin (or origins) of the PG remain intensely debated. One class of proposals names exotic, possibly topological order in the particle-hole channel as the origin1,2,3, while another points to pairing fluctuations in the strong coupling regime4,5,6,7,8,9,10. Substantial numerical efforts have been dedicated to the understanding of PG, see e.g., refs. 1,11,12 and references therein.

The understanding of the coupling between fermionic excitations near the Fermi surface (FS) and bosonic quantum critical fluctuations14,15,16,17,18 is crucial to describe these three features. The development of quantum Monte Carlo (QMC) algorithms for a class of models of this type, pioneered by ref. 19, has created a feasible way to study this physics in an unbiased manner (see the reviews20,21 and references within). In QMC models, FS fermions couple to bosonic fluctuations, representing certain collective modes of low-energy fermions22,23,24,25,26,27,28,29. The bosonic part is bestowed with independent (non-fermionic) dynamics and can be tuned to criticality to mimic the situation in real materials. Crucially, these models are free of the sign-problem plaguing most fermionic QMC, allowing for a realistic test of theory.

In this work, we investigate the PG physics via such a sign-free QMC simulation of fermions near a ferromagnetic QCP. We find robust signatures of a PG above the SC state and are able to study its spectral properties and its interplay with the dynamics of the ferromagnetic degrees of freedom. We also compare the numerical results with several theoretical predictions, and reconcile many key aspects of the two.

Results

Overview

Before going into the details of our work, we present an overview of the essential features of our model and a summary of the main results.

The model we choose to implement is a variant of a quantum critical model, in which the bosons represent critical ferromagnetic(FM) spin fluctuations (a “spin-fermion” model). When looking for a spectral property of the superconductivity, such as a PG, such a model has an advantage over analogous ones, e.g., antiferromagnetic or nematic models (see e.g.,21) because of the simplicity of the momentum structure of the FS (e.g., no hot or cold spots). Furthermore, compared to earlier sign-free QMC studies on ferromagnetism, the coupling strength of our model is stronger in two aspects. First, the spin system is an XY quantum rotor model that is inherently more strongly fluctuating than an Ising model, analyzed earlier30,31. Second, the coupling constant K between the fermionic and bosonic sectors is set to larger values than in previous works. The larger coupling pushes the region of SC fluctuations up to temperatures, where they are discernible in the numerical data. This in contrast with earlier works, where coupling strength was optimized to study normal state properties. As we see below, the larger coupling allows us to reveal the PG behavior.

In the normal state, at low enough temperatures we find in the bosonic sector near the QCP an overdamped dynamics with linear frequency response (z = 2 scaling). This is different from the z = 3 behavior, found in Ising systems, and is a result of a non-conservation of the order parameter in our model.

In the temperature range, where the bosonic susceptibility is linear in frequency, we observe several remarkable features. The uniform susceptibility deviates from Curie-Weiss behavior and actually becomes weaker at smaller T. In the fermionic sector, we find a gap-like feature in the density of states (DOS). Unlike in a BCS superconductor, the size of the gap remains roughly independent on temperature, while the DOS becomes progressively depleted (filled) upon lowering (raising) temperature. Importantly, the scaling behavior of the pairing susceptibility clearly shows that the system is not in a SC state. We thus identify the spectral gap in such a state as a PG.

We note that the “gap-filling" behavior observed in our numerical results has also been observed in tunneling and photoemission experiments on the cuprates32, and has been obtained in a class of γ − models of quantum-critical pairing 6. Our results, obtained from unbiased large-scale QMC simulations, confirm the existence of a PG behavior from pairing fluctuations in a quantum-critical system with itinerant fermions.

The quantum-critical spin dynamics and normal state fermionic properties that we find are consistent with recent theoretical predictions for nFLs at finite temperature, obtained within the modified Eliashberg theory31,33,34. This allows us to benchmark our simulations and extract relevant parameters from the observables (see Supplementary Note 4 for details). The onset temperature for PG behavior, TPG, and SC Tc, extracted from QMC, are consistent with theoretical predictions (see ref. 6 and Methods). Our results therefore provide an attempt to numerically realize the transition from nFL to PG and eventually to superconductivity, lending support to the scenario of pairing fluctuations driven PG phenomena.

Crucially, we view our finding of a PG as the evidence of a universal mechanism for the formation of a PG from SC fluctuations near a QCP, not limited to the specific model of FM spin fluctuations that we used. We do not claim that we present a model for PG formation in a specific material, but in view of our findings we do expect SC fluctuations to be a contributing factor to PG formation in any system close enough to a QCP, independently of the specific origin of the pairing boson.

The PG, obtained in our work, comes from pairing fluctuations in a situation when the pairing is in turn mediated by a propagator of a FM order parameter. While our model does not directly describe experimental situation in the cuprates, where antiferromagnetic fluctuations are often considered to be a pairing glue, we argue that the mechanism for the PG formation, studied in our work, is a universal phenomenon of the pairing near a quantum-critical point6, and in this sense goes beyond the specific model with FM flucutations. We do expect the SC fluctuations to be a contributing factor to PG formation in any system close enough to a QCP, independently of the specific origin of the pairing glue.

Model

We consider a model of itinerant fermions coupled to SO(2) quantum rotors, as shown in Fig. 1a (rotors are in the middle layer). The model is described by

$$\hat{H}={\hat{H}}_{{{{{{{{\rm{qr}}}}}}}}}+{\hat{H}}_{{{{{{{{\rm{f}}}}}}}}}+{\hat{H}}_{{{{{{{{\rm{qr}}}}}}}}-{{{{{{{\rm{f}}}}}}}}},$$
(1)

where

$$ {\hat{H}}_{{{{{{{{\rm{qr}}}}}}}}} =\frac{U}{2}\mathop{\sum}\limits_{i}{\hat{L}}_{i}^{2}-{t}_{{{{{{{{\rm{b}}}}}}}}}\mathop{\sum}\limits_{\langle i,j\rangle }\cos \left({\hat{\theta }}_{i}-{\hat{\theta }}_{j}\right)\\ {\hat{H}}_{{{{{{{{\rm{f}}}}}}}}} =-{t}_{1}\mathop{\sum}\limits_{\langle i,j\rangle \sigma \lambda }{\hat{c}}_{i\sigma \lambda }^{{{{\dagger}}} }{\hat{c}}_{j\sigma \lambda }-{t}_{2}\mathop{\sum}\limits_{\langle \langle i,j\rangle \rangle \sigma \lambda }{\hat{c}}_{i\sigma \lambda }^{{{{\dagger}}} }{\hat{c}}_{j\sigma \lambda }-\mu \mathop{\sum}\limits_{i\sigma \lambda }{\hat{n}}_{i\sigma \lambda }\\ {\hat{H}}_{{{{{{{{\rm{qr}}}}}}}}-{{{{{{{\rm{f}}}}}}}}} =-\frac{K}{2}\mathop{\sum}\limits_{i\lambda }\left({\hat{c}}_{i\lambda }^{{{{\dagger}}} }{\sigma }^{x}{\hat{c}}_{i\lambda }\cdot \cos {\hat{\theta }}_{i}+{\hat{c}}_{i\lambda }^{{{{\dagger}}} }{\sigma }^{y}{\hat{c}}_{i\lambda }\cdot \sin {\hat{\theta }}_{i}\right).$$
(2)
Fig. 1: Model and Phase diagram.
figure 1

a Sketch of the model in Eq. (1). Deep yellow dots and grids of the top and bottom layers represent fermion degrees of freedom with nearest hopping strength t1 (e.g., \({\hat{c}}_{i}^{{{{\dagger}}} }{\hat{c}}_{j}\)) and next-nearest hopping t2 (e.g., \({\hat{c}}_{i}^{{{{\dagger}}} }{\hat{c}}_{k}\)). Blue arrows and grid in the middle layers denote bosonic parts with an unit vector representing θ [0, 2π) of the rotor on each site. The interaction between two rotors on nearest sites (e.g., θm and θn) is tb. The on-site coupling K between fermions and bosons is shown by the vertical dashed lines. The system size is L × L. b U − T phase diagram of the model obtained from QMC simulation. The inset zooms in to the vicinity of the pseudogap (PG, yellow), ferromagnetic (FM, blue) and superconducting (SC, orange) regions. The blue points on the FM phase boundary are determined by finite size scaling with fixed T or U. Notably, for U = 5.9, as temperature gets lower, the system first enters into the FM phase at T ≈ 0.13, then exits it at T ≈ 0.08. The yellow points of the PG boundary are determined from the onset of a PG in the single-particle spectrum, as shown in Fig. 2. The red points denoting an onset of s-wave superconductivity are determined from the onset of a full gap in the spectrum as well as Kosterlitz-Thouless scaling of the pairing susceptibility. The maximum of SC phase transition temperature Tc is ~0.05. The green points and dashed line, are the phase boundary of the (uncoupled) quantum rotor model13. See the Supplementary Note 3 for additional details as well as a discussion of SC fluctuations above Tc. The errorbars of the points on the FM phase and SC phase boundaries are determined by the data collapse with fixed T or U. For the PG phase boundary, the errorbars come from the uncertainty in identifying the onset of the minimum at ω = 0 for distinct temperatures of DOS.

The first term \({\hat{H}}_{{{{{{{{\rm{qr}}}}}}}}}\) describes a quantum rotor model on a square lattice. Here \({\hat{L}}_{i}\) is the angular momentum of 2D rotor \({\hat{\theta }}_{i}\) at site i. The second term \({\hat{H}}_{{{{{{{{\rm{f}}}}}}}}}\) describes two identical copies of spin-1/2 fermions on a square lattice, with layer index λ = 1 and 2 representing the top and the bottom layers. Fermions in each layer can hop between nearest-neighbor (next-nearest-neighbor) sites with hopping amplitudes t1 (t2), and the chemical potential μ controls the fermion density. The last term \({\hat{H}}_{{{{{{{{\rm{qr}}}}}}}}-{{{{{{{\rm{f}}}}}}}}}\) couples quantum rotors and fermions via an on-site FM interaction that tends to align XY component of a fermion spin with the direction of a rotor on each site.

In the absence of fermion-rotor coupling, rotors develop quasi-long-range FM order via a Kosterlitz-Thouless(KT) transition13,35. At zero temperature, FM order becomes long range. The KT transition line in (T, U) plane terminates at a QCP at \({(U/{t}_{{{{{{{{\rm{b}}}}}}}}})}_{{{{{{{{\rm{c}}}}}}}}}=4.25(2)\)13,36,37. As we turn on the fermion-rotor coupling, fermion contributions shift the KT phase boundary towards larger U and T. More importantly, the phase transition now involves fermion spins, which at T = 0 also order ferromagnetically below Uc. This allows us to study quantum phenomena near a FM QCP in a metal38. Due to the anti-unitary symmetry and the presence of two copies of fermions, this model can be simulated via QMC techniques without the sign problem (see Supplementary Note 1 for details). This setup then allows us to analyze the universal behavior near a QCP with high numerical accuracy and large system sizes.

We express all quantities in units of tb. In the simulations we set K = 4, t1 = 1, t2 = 0.2 and μ = 0. We varied U and the temperature T and constructed the phase diagram of the model, Fig. 1b, which features a paramagnetic-ferromagnetic transition and several other transitions/phases. The magnetic transition at a finite temperature is of KT type. As U increases, the transition temperature decreases and terminates at a QCP at Uc. The T = 0 transition upon varying U belongs to XY universality class as the coupling to rotors creates an easy plane for fermion spins. Fermion spins order ferromagnetically in the XY plane, breaking a spin-rotational symmetry.

Pseudogap and superconductivity properties

We observe a SC dome around the QCP. Above the dome, we find evidence of PG behavior in the range of T, whose width is comparable to Tc.

First, by measuring correlation functions of Cooper pairs in various pairing channels, we found that the dominant pairing channel is spin-triplet and odd under the interchange between the top and the bottom layers (layer-singlet), i.e., \(\Delta ({{{{{{{\bf{r}}}}}}}})=\frac{1}{\sqrt{2}}({\hat{c}}_{{{{{{{{\bf{r}}}}}}}}1\uparrow }{\hat{c}}_{{{{{{{{\bf{r}}}}}}}}2\downarrow }-{\hat{c}}_{{{{{{{{\bf{r}}}}}}}}2\uparrow }{\hat{c}}_{{{{{{{{\bf{r}}}}}}}}1\downarrow })=\frac{1}{\sqrt{2}}({\hat{c}}_{{{{{{{{\bf{r}}}}}}}}1\uparrow }{\hat{c}}_{{{{{{{{\bf{r}}}}}}}}2\downarrow }+{\hat{c}}_{{{{{{{{\bf{r}}}}}}}}1\downarrow }{\hat{c}}_{{{{{{{{\bf{r}}}}}}}}2\uparrow })\), where 1 and 2 label layers. In the classification of 2D irreducible representations, this is an s-wave gap, as Δ(0) is finite. We verified (see Supplementary Figs. 1 and 2) that the susceptibility in this channel strongly increases when the system approaches a superconducting instability, while the susceptibilities in all other pairing channels remain small and do not increase. This observation is a direct evidence that superconductivity originates from the interaction mediated by soft bosonic fluctuations, associated with the QCP. Indeed, it has long being known that near a FM quantum phase transition, soft dynamical bosonic fluctuations introduce an effective interaction that is attractive in the spin-triplet channel. In the geometry of our model, there are two distinct types of spin-triplet pairing—one is odd under momentum inversion in a layer and even under layer interchange (e.g., p-wave layer-triplet), the other is even within each layer and odd under layer interchange (s-wave layer-singlet). By analogy with previous studies of the pairing mediated by small q fluctuations30, one expects the leading instability to be towards the s-wave layer-singlet, spin-triplet order. The numerical finding of the largest pairing correlations in this channel thus affirms the crucial role of soft FM bosonic fluctuations in the formation of a SC dome.

Second, we obtained the fermionic spectral function, and then integrated it over k − space to obtain the DOS N(ω). For this, we first computed the imaginary-time fermion Green’s function and then converted it to real frequency via stochastic analytic continuation (SAC) method (See Methods for details). We show the results for N(ω) in Fig. 2a. At low T, inside the SC dome, there is clear evidence for an s-wave gap. The data shows that, that as T increases, the magnitude of the gap slightly increases, rather than shrinks, as would be the case in a BCS superconductor. Simultaneously, N(ω) for ω smaller than the gap increases and gradually fills in the states within the gap, ultimately restoring its normal-state value. This phenomenon has been termed gap-filling. It is qualitatively in agreement with experimental observations in many strongly-correlated unconventional superconductors at TTc5,9,10,39, At smaller TTc, the DOS displays gap-closing behavior, like in a conventional BCS superconductor. Guided by the experimental evidence9,10 that gap-filling behavior holds at TTc, we defined the PG region as the one where the DOS gets filled in upon increasing T. We set the lower boundary of this region to where the DOS at the Fermi energy significantly deviates from thermally activated behavior of \({e}^{-\Delta /{k}_{{{{{{{{\rm{B}}}}}}}}}T}\). The upper boundary of the PG region is set at TPG, at which the dip of N(ω) at the Fermi energy becomes invisible. The PG region, obtained this way, is plotted in yellow in Fig. 1b.

Fig. 2: Pseudogap and superconductivity.
figure 2

a Local DOS N(ω) for various temperatures at U = 6 with L = 12. For T = 1/4, far above the PG, the system exhibits a Fermi liquid spectrum. At TPG ≈ 0.1, the SC fluctuations begin to play important role and a noticable gap forms at ω = 0. This gap-forming temperature is consistent with the corresponding intermediate temperature scale in the dynamic bosonic susceptibility χ in Fig. 4b. The gap minimum N(ω = 0) goes down with temperature and eventually reaches zero at T ≈ 0.05, indicating the onset of a SC state as detected by the pairing susceptibility in b. At T < 0.05, fully gapped spectrum corresponds to the behavior of the SC state. b Data collapse of the pairing susceptibility Ps versus temperature at U = 6 for system sizes L = 12, 14, 16, 18 with statistical errors obtained by QMC simulations, consistent with a KT transition. The best fit coefficients are A = 0.75, Tc = 0.048, which is consistent with the temperature of the fully-gapped spectrum in a.

Third, to determine the actual SC transition temperature, Tc, we performed scaling analysis of the pairing susceptibility \({P}_{{{{{{{{\rm{s}}}}}}}}}=\frac{1}{{L}^{2}}\int\nolimits_{0}^{\beta }{\sum }_{i}({\Delta }^{{{{\dagger}}} }({{{{{{{{\bf{r}}}}}}}}}_{i},\tau )\Delta ({{{{{{{\bf{0}}}}}}}},0))\), using KT scaling for the pairing susceptibility \({P}_{{{{{{{{\rm{s}}}}}}}}}={L}^{2-\eta }f(L\cdot \exp (-\frac{A}{{(T-{T}_{{{{{{{{\rm{c}}}}}}}}})}^{1/2}}))\) for T > Tc with ηKT = 1/428,40,41. We show the results in Fig. 2b. The data for Ps for various system sizes and temperatures collapse onto a single curve. We fitted the curve by the formula above and extracted Tc = 0.048. This agrees with the lower boundary of the PG region. The upper boundary, TPG, is about twice larger in our simulations, TPG ~ 0.1. We also computed the superfluid density, ρs(T), which has been widely used to estimate Tc in QMC simulations. This is done by detecting the temperature Tρ at which ρs(Tρ) = αTρ, where α is a dimensionless constant41, usually set to 2/π, based on the analysis of the XY model42. This criterion, although qualitatively correct, typically overestimates Tc40. In our case, we found Tc < Tρ ~ TPG. We discuss our analysis of ρs in some length in Supplementary Note 3.

We analyzed the QMC data within the quantum critical theory of itinerant ferromagnets43,44, extended to finite T34 and modified to include two layers of itinerant fermions and superconductivity. We computed fermionic and bosonic self-energies near Uc and found good agreement with the simulations in the normal state (see Supplementary Note 4). We extracted the effective fermion-boson coupling from this comparison, and used it to compute the onset temperature for the pairing within the Eliashberg theory for quantum-critical pairing6. This theory does not differentiate between pair formation and superconductivity, hence the result has to be compared with TPG, extracted from simulations. We obtained theoretical TPG ~ 0.08, quite consistent with TPG ~ 0.1, extracted from QMC data, see Fig. 1b. Further, Eliashberg calculations below TPG show gap-closing behavior at small T and gap-filling behavior at TTPG. The boundary between two regimes has been associated with the actual Tc, based on the analysis of phase fluctuations6. We show this in the phase diagram in Fig. 1b. Based on this comparison, we argue that our unbiased numerical QMC simulations are consistent with the theory and provide strong evidence for PG behavior, originating from preformed pairs above Tc, near a FM QCP in a metal.

We note that previous QMC work (see e.g.,22,23 for an antiferromagnetic model) found an SC dome surrounded by a region of SC fluctuations. These were determined by comparing Tρ determined from the BKT criterion for ρs (see above), and the temperature Tdia at which the system showed diamagnetic behavior, as evidenced by a sign-change of the appropriate current-current correlator. While we too find such fluctuations (see Supplementary Note 3). We stress that the region of gap-filling behavior is predominantly at Tc < T < TPG ~ Tρ, and is therefore distinct from fluctuations on the scale of Tdia. Indeed, from our numerics we observe that TPG < Tdia. We emphasize that in the thermodynamic limit, TPG and Tdia do not correspond to phase transitions, but rather mark crossover regions.

Magnetic dynamics and re-entrance effect

The pairing behavior also has an impact on the magnetic phase transition and the quantum dynamics of the rotors. As shown in Fig. 1b, the phase boundary of the paramagnetic-ferromagnetic transition exhibits a re-entrance behavior at U ~ 5.9, close to the QCP. For example, at U = 5.9, upon reducing the temperature, the system first enters the FM state and then returns to the paramagnetic one, i.e., there is a “back-bending” of the transition line to the FM state. This can also be seen from the Fermi surface behavior. In Fig. 3, we plot the Fermi surface, G(k, τ = β/2) ~ N(k, ω = 0), evolution with temperature. At intermediate temperature T = 0.1, the Fermi surface splits due to the ferromangetic order. However, the split vanishes both either increasing or lowering the temperature. We believe that the re-entrance phenomenon is a consequence of the PG and SC fluctuations, which suppress the fermion DOS and hence the electron-hole contribution to magnetic order45,46. Similar behavior has been seen previously in an antiferromagnetic model22, but no PG was reported there. We emphasize that the paramagnetic-ferromagnetic phase boundary starts to bend to the left roughly at TPG, which is well above the SC dome, indicating that SC fluctuations without phase coherence in the PG region are responsible for the magnetic dynamics.

Fig. 3: Re-entrance.
figure 3

Evolution of the Fermi surface (FS) from non-interacting system with Hf in a, to the nFL FS subjected to strong FM correlation at U = 5.9, T = 0.1 in b, and eventually to the FS in the PG phase at U = 5.9, T = 0.05 in c. The spectral weights are normalized with the same scale. The system size is L = 12 and the twisted boundary condition in the fermion hopping is applied such that the momentum resolution is 4 times larger in both kx and ky directions.

We note, that in the absence of an SC dome, previous works (see e.g.,38,44,47,48) have shown that itinerant FM QCPs are unstable to a first-order transition driven by normal state magnetic fluctuations, which can cut off the FM phase at larger U’s, similar to the behavior seen in our simulations. In our results, we do not observe clear evidence of a first-order magnetic transition, and the correlation of the back-bending with TPG implies that for the parameters that we used the physics is driven chiefly by SC rather than magnetic fluctuations.

In addition, we measured the inverse dynamical bosonic susceptibility of the rotors across different regions of the phase diagram. Our results are summarized in Fig. 4, showing data for three representative U at various temperatures. To study the dynamics, we subtract the static part of the inverse susceptibility and focus on the spin polarization χ−1(q, ω) − χ−1(q = 0, ω = 0) (for details of what follows see Supplementary Note 2 and 3). Deep in the FM phase, Fig. 4a, we find an ω2 dependence (dynamical exponent z = 1). This is similar to that of the bare rotor model, and indicates that the fermionic contribution to the dynamics is negligible because of the spin gap. Similarly, deep in the Fermi liquid phase, Fig. 4c, we find an ω2 dependence, except at the lowest frequencies, which furthermore extrapolates to a nonzero value. The saturation is readily understood as resulting from the non-analyticity of the Lindhard function which implies χ−1(q = 0, ω → 0) ≠ χ−1(q → 0, ω = 0) at weak coupling. In the quantum-critical regime, Fig. 4b, we find a qualitatively different behavior indicating strong fermionic correlations.

Fig. 4: Magnetic dynamics.
figure 4

Inverse bosonic dynamic susceptibility χ versus ωn in three different regions, at U = 3, 6, 8, corresponding a in the FM phase, b in PG and SC phases, and c disordered phase. a log-log plot for various system size L = 6, 8, 10, 12, 14, each of which includes various β = 12, 16, 20, 24. Red line is a quadratic line of χ−1 ~ ω2 for low frequency part ωn < 1. b log-log plot for various β = 10, 16, 20, 24, 30 with L = 12. At temperature T = 0.1 (β = 10), the fermions are in the quantum critical regime, the bosonic susceptibility is the linear function of ω, as indicated by the orange line. When temperature gets lower, the fermion goes into PG phase, prompting the bosonic scaling behavior to deriviate from linear function. And upon entering the SC phase, the χ−1 ~ ω2 (the blue line, a guide to the eye) as in the FM phase. c Bosonic susceptibility in the disordered phase at U = 8 plotted for system size L = 6, 8, 10, 12 with various β = 12, 16, 20, 24. At high frequency all data points successfully merge together, as indicated by the red line quadratic in ω for ω > 1. At low frequency, the χ−1 ~ ω as indicated by the blue line, which is a guide to the eye, due to the non-conserved rotor order parameter.

First, at higher frequencies we find a linear frequency dependence (z = 2), which does not saturate to a finite value. This is surprising at first glance, since Landau damping for a ferromagnet has an ω/q form (z = 3) rather than linear ω. We note that, in purely electronic models χ−1(q, ω) is required to be non-analytic at any coupling strength due to spin conservation, and non-analytic behavior was seen previously in simulations of Ising-ferromagnets. However, in our simulations of the XY model, the order parameter is not conserved, leading to linear frequency dependence, in direct contrast to Ising model studied in refs. 30,31.

Second, at lower Matsubara frequencies accessible at lower temperatures, the ω2 behavior is again restored even in the quantum-critical region. As discussed above, this is again a direct result of the formation of a gap—this time the pseudogap, which depletes the low-energy fermion density of states and reduces the fermionic feedback on the bosons.

From the analysis above, we see that the spin dynamics is consistent with the quantum-critical behavior and PG physics.

Discussion

In this work, we performed a large-scale quantum Monte Carlo simulation of a FM spin-fermion model. We reported direct spectral and thermodynamic evidence of the formation of a PG prior to the SC transition. Within such a PG phase, the temperature evolution of the fermion spectral gap exhibits a gap-filling behavior, in sharp contrast with that of a conventional superconductor. Moreover, we found that the dynamics of the spin fluctuations display a different behavior than the well-known Landau damping behavior with z = 3.

Remarkably, we were able to reconcile all these features with theoretical predictions of Eliashberg theory and its generalization to the γ-model. Experimentally, PG phases have been observed in various unconventional superconductors49,50, most notably the cuprates51. Our results imply that a PG arising from strong dynamical fluctuations should be ubiquitous in quantum-critical metals, and we expect this to be a fruitful direction for future research.

Methods

QMC simulations and data analysis

We employ the determinant quantum Monte Carlo (DQMC) method20,30 to simulate the Hamiltonian in Eq. (1). The quantum rotor model plays the role of the auxiliary field in the conventional DQMC and the quantum rotor model can be efficiently simulated with non-local update scheme developed in our previous work13. For each realization of the rotor in space-time, the fermion determinant is evaluated with the kinetic part and the coupling part of the Hamiltonian included as the configurational weight and the Markov chain of the Monte Carlo process is carried out according the weight. Detailed measurements of the physical observables are given in the Supplementary Note 2.

In order to obtain real-frequency spectral functions, the SAC scheme is employed to obtain the spectral function N(ω) from the imaginary-time correlation function G(τ),

$$G(\tau )=\int\nolimits_{-\infty }^{\infty }d\omega \frac{{e}^{-\omega (\tau -\beta /2)}}{2\cosh (\beta \omega /2)}N(\omega )$$
(3)

It is known that the problem of inverting the Laplace transform is equivalent to find the most probable spectra N(ω) out of its exponentially many suggestions to match the QMC correlation function G(τ) with respect to its stochastic errors, and such transformation has been converted to a Monte Carlo sampling process52,53,54. This QMC-SAC approach has been successfully applied to quantum magnets and interacting fermion systems ranging from the simple square lattice Heisenberg antiferromagnet55 to deconfined quantum critical point and quantum spin liquids with their fractionalized excitations56,57,58 and to the continuum model of twisted bilayer graphene and benchmarked with the exact solution at the chiral limit59,60.

Theoretical analysis

We analyzed the QMC data for fermionic and bosonic response using the modified Eliashberg theory, which is a low energy effective dynamical theory for itinerant fermions near a QCP at finite temperatures. The theory accepts as parameters the static properties of a coupled fermion-boson system near a QCP, e.q. fermion bandstructure, bosonic susceptibility, etc., and computes the dynamical response of the system in terms of the fermionic self energy Σ(k, ωn) and bosonic self energy (polarization) Π(q, Ωn), taking into account the low energy excitations near the FS. It accounts for deviations from the canonical T → 0 quantum critical behavior, e.g., deviations from the \(\Sigma \sim {\omega }_{n}^{2/3}\) nFL self energy, and from the Landau damping Π ~ Ωn/(vFq) as discussed in the main text. For details on the method see refs. 31,34.

We applied the theory to our QMC data, both to verify our assumptions on the normal state of the system and to extract the effective fermion-boson coupling. In the bare theory, the coupling \(\bar{g} \sim {K}^{2}\), but it is renormalized by fermions with energies of order of the bandwidth, so it should be extracted by fitting from the QMC data. We present results for U = 6 which is almost above the QCP in Supplementary Fig. 13, showing good agreement between theory and data. For details of the fitting procedure and a discussion of the quality of the fits are presented in Supplementary Note 4. We found \(\bar{g}=6.3\pm 0.2\), representing about a 20% renormalization of the bare vertex K, which is consistent with earlier works34.

Finally, we used the obtained \(\bar{g}\) to predict TPG within Eliashberg theory (the γ − model). Our model corresponds to γ = 1/3. The analytical prediction for TPG can be found in ref. 6, and details of the conversion from our \(\bar{g}\) to the γ − model parameters are in the Supplementary Note 4. We found TPG ≈ 0.08, in good agreement with the QMC TPG ~ 0.1.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability

All numerical codes in this paper are available upon reasonable request to the authors.

References

  1. Scheurer, M. S. et al. Topological order in the pseudogap metal. Proc. Natl Acad. Sci. 115, E3665 (2018).

    MathSciNet  CAS  PubMed  PubMed Central  MATH  Article  Google Scholar 

  2. Sachdev, S. Topological order, emergent gauge fields, and fermi surface reconstruction. Rep. Prog. Phys. 82, 014001 (2018).

    ADS  MathSciNet  PubMed  Article  CAS  Google Scholar 

  3. Wu, W. et al. Pseudogap and fermi-surface topology in the two-dimensional hubbard model. Phys. Rev. X 8, 021048 (2018).

    CAS  Google Scholar 

  4. Emery, V. J. & Kivelson, S. A. Importance of phase fluctuations in superconductors with small superfluid density. Nature 374, 434 (1995).

    ADS  CAS  Article  Google Scholar 

  5. Keimer, B., Kivelson, S. A., Norman, M. R., Uchida, S. & Zaanen, J. From quantum matter to high-temperature superconductivity in copper oxides. Nature 518, 179 (2015).

    ADS  CAS  PubMed  Article  Google Scholar 

  6. Wu, Y.-M., Abanov, A., Wang, Y. & Chubukov, A. V. Interplay between superconductivity and non-fermi liquid at a quantum critical point in a metal. ii. Phys. Rev. B 102, 024525 (2020).

    ADS  CAS  Article  Google Scholar 

  7. Wang, H., Chudnovskiy, A. L., Gorsky, A. & Kamenev, A. Sachdev-ye-kitaev superconductivity: Quantum kuramoto and generalized richardson models. Phys. Rev. Res. 2, 033025 (2020).

    CAS  Article  Google Scholar 

  8. Dahm, T. et al. Strength of the spin-fluctuation-mediated pairing interaction in a high-temperature superconductor. Nat. Phys. 5, 217 (2009).

    CAS  Article  Google Scholar 

  9. Reber, T. J. et al. Prepairing and the “filling” gap in the cuprates from the tomographic density of states. Phys. Rev. B 87, 060506 (2013).

    ADS  Article  CAS  Google Scholar 

  10. Kanigel, A. et al. Evolution of the pseudogap from fermi arcs to the nodal liquid. Nat. Phys. 2, 447 (2006).

    CAS  Article  Google Scholar 

  11. Gull, E., Parcollet, O. & Millis, A. J. Superconductivity and the pseudogap in the two-dimensional hubbard model. Phys. Rev. Lett. 110, 216405 (2013).

    ADS  PubMed  Article  CAS  Google Scholar 

  12. LeBlanc, J. P. F. et al. Solutions of the two-dimensional hubbard model: Benchmarks and results from a wide range of numerical algorithms. Phys. Rev. X 5, 041041 (2015).

    Google Scholar 

  13. Jiang, W., Pan, G., Liu, Y. & Meng, Z. Y. Solving quantum rotor model with different Monte Carlo techniques. Chin. Phys. B 31, 040504 (2022).

  14. Abanov, A., Chubukov, A. V. & Schmalian, J. Quantum-critical theory of the spin-fermion model and its application to cuprates: Normal state analysis. Adv. Phys. 52, 119 (2003).

    ADS  CAS  Article  Google Scholar 

  15. Metzner, W., Rohe, D. & Andergassen, S. Soft fermi surfaces and breakdown of fermi-liquid behavior. Phys. Rev. Lett. 91, 066402 (2003).

    ADS  CAS  PubMed  Article  Google Scholar 

  16. Metlitski, M. A. & Sachdev, S. Quantum phase transitions of metals in two spatial dimensions. i. ising-nematic order. Phys. Rev. B 82, 075127 (2010a).

    ADS  Article  CAS  Google Scholar 

  17. Metlitski, M. A. & Sachdev, S. Instabilities near the onset of spin density wave order in metals. N. J. Phys. 12, 105007 (2010b).

    Article  CAS  Google Scholar 

  18. Lee, S.-S. Recent developments in non-fermi liquid theory. Annu. Rev. Condens. Matter Phys. 9, 227 (2018).

    ADS  CAS  Article  Google Scholar 

  19. Erez Berg, S. S. & Metlitski, M. A. Sign-problem-free quantum monte carlo of the onset of antiferromagnetic in metals. Science 338, 1606 (2012).

    ADS  MathSciNet  PubMed  MATH  Article  CAS  Google Scholar 

  20. Xu, X. Y. et al. Revealing fermionic quantum criticality from new monte carlo techniques. J. Phys. Condens. Matter 31, 463001 (2019a).

    CAS  PubMed  Article  Google Scholar 

  21. Berg, E., Lederer, S., Schattner, Y. & Trebst, S. Monte carlo studies of quantum critical metals. Annu. Rev. Condens. Matter Phys. 10, 63 (2019).

    ADS  CAS  Article  Google Scholar 

  22. Schattner, Y., Gerlach, M. H., Trebst, S. & Berg, E. Competing orders in a nearly antiferromagnetic metal. Phys. Rev. Lett. 117, 097002 (2016).

    ADS  PubMed  Article  CAS  Google Scholar 

  23. Gerlach, M. H., Schattner, Y., Berg, E. & Trebst, S. Quantum critical properties of a metallic spin-density-wave transition. Phys. Rev. B 95, 035124 (2017).

    ADS  Article  Google Scholar 

  24. Bauer, C., Schattner, Y., Trebst, S. & Berg, E. Hierarchy of energy scales in an o(3) symmetric antiferromagnetic quantum critical metal: A monte carlo study. Phys. Rev. Res. 2, 023008 (2020).

    CAS  Article  Google Scholar 

  25. Xu, X. Y. et al. Monte carlo study of lattice compact quantum electrodynamics with fermionic matter: The parent state of quantum phases. Phys. Rev. X 9, 021022 (2019b).

    CAS  Google Scholar 

  26. Liu, Z. H., Pan, G., Xu, X. Y., Sun, K. & Meng, Z. Y. Itinerant quantum critical point with fermion pockets and hotspots. Proc. Natl Acad. Sci. 116, 16760 (2019).

    ADS  MathSciNet  CAS  PubMed  PubMed Central  MATH  Article  Google Scholar 

  27. Chen, C., Xu, X. Y., Qi, Y. & Meng, Z. Y. Metal to orthogonal metal transition. Chin. Phys. Lett. 37, 047103 (2020).

    ADS  CAS  Article  Google Scholar 

  28. Chen, C., Yuan, T., Qi, Y. & Meng, Z. Y. Fermi arcs and pseudogap in a lattice model of a doped orthogonal metal. Phys. Rev. B 103, 165131 (2021).

    ADS  CAS  Article  Google Scholar 

  29. Gazit, S., Assaad, F. F. & Sachdev, S. Fermi surface reconstruction without symmetry breaking. Phys. Rev. X 10, 041057 (2020).

    CAS  Google Scholar 

  30. Xu, X. Y., Sun, K., Schattner, Y., Berg, E. & Meng, Z. Y. Non-fermi liquid at (2+1)d ferromagnetic quantum critical point. Phys. Rev. X 7, 031058 (2017).

    Google Scholar 

  31. Xu, X. Y., Klein, A., Sun, K., Chubukov, A. V. & Meng, Z. Y. Identification of non-fermi liquid fermionic self-energy from quantum monte carlo data. npj Quantum Mater. 5, 65 (2020).

    ADS  Article  Google Scholar 

  32. Vishik, I. M. Photoemission perspective on pseudogap, superconducting fluctuations, and charge order in cuprates: a review of recent progress. Rep. Prog. Phys. 81, 062501 (2018).

    ADS  MathSciNet  CAS  PubMed  Article  Google Scholar 

  33. Chubukov, A. V., Betouras, J. J. & Efremov, D. V. Non-landau damping of magnetic excitations in systems with localized and itinerant electrons. Phys. Rev. Lett. 112, 037202 (2014).

    ADS  PubMed  Article  CAS  Google Scholar 

  34. Klein, A., Chubukov, A. V., Schattner, Y. & Berg, E. Normal state properties of quantum critical metals at finite temperature. Phys. Rev. X 10, 031053 (2020).

    CAS  Google Scholar 

  35. José, J. V., Kadanoff, L. P., Kirkpatrick, S. & Nelson, D. R. Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model. Phys. Rev. B 16, 1217 (1977).

    ADS  Article  Google Scholar 

  36. Hasenbusch, M. & Török, T. High-precision monte carlo study of the 3d xy-universality class. J. Phys. A: Math. Gen. 32, 6361 (1999).

    ADS  MathSciNet  MATH  Article  Google Scholar 

  37. Meng, Z. Y. & Wessel, S. Phases and magnetization process of an anisotropic shastry-sutherland model. Phys. Rev. B 78, 224416 (2008).

    ADS  Article  CAS  Google Scholar 

  38. Brando, M., Belitz, D., Grosche, F. M. & Kirkpatrick, T. R. Metallic quantum ferromagnets. Rev. Mod. Phys. 88, 025006 (2016).

    ADS  Article  Google Scholar 

  39. Reber, T. J. et al. The origin and non-quasiparticle nature of fermi arcs in Bi2Sr2CaCu2O8+δ. Nat. Phys. 8, 606 (2012).

    CAS  Article  Google Scholar 

  40. Paiva, T., dos Santos, R. R., Scalettar, R. T. & Denteneer, P. J. H. Critical temperature for the two-dimensional attractive hubbard model. Phys. Rev. B 69, 184501 (2004).

    ADS  Article  CAS  Google Scholar 

  41. Costa, N. C., Blommel, T., Chiu, W.-T., Batrouni, G. & Scalettar, R. T. Phonon dispersion and the competition between pairing and charge order. Phys. Rev. Lett. 120, 187003 (2018).

    ADS  CAS  PubMed  Article  Google Scholar 

  42. Pokrovsky, V. Properties of ordered, continuously degenerate systems. Adv. Phys. 28, 595 (1979).

    ADS  Article  Google Scholar 

  43. Rech, J., Pépin, C. & Chubukov, A. V. Quantum critical behavior in itinerant electron systems: Eliashberg theory and instability of a ferromagnetic quantum critical point. Phys. Rev. B 74, 195126 (2006).

    ADS  Article  CAS  Google Scholar 

  44. Maslov, D. L. & Chubukov, A. V. Nonanalytic paramagnetic response of itinerant fermions away and near a ferromagnetic quantum phase transition. Phys. Rev. B 79, 075112 (2009).

    ADS  Article  CAS  Google Scholar 

  45. Moon, E. G. & Sachdev, S. Competition between spin density wave order and superconductivity in the underdoped cuprates. Phys. Rev. B 80, 035117 (2009).

    ADS  Article  CAS  Google Scholar 

  46. Moon, E. G. & Sachdev, S. Quantum critical point shifts under superconductivity: Pnictides and cuprates. Phys. Rev. B 82, 104516 (2010).

    ADS  Article  CAS  Google Scholar 

  47. Belitz, D., Kirkpatrick, T. R. & Vojta, T. Nonanalytic behavior of the spin susceptibility in clean Fermi systems. Phys. Rev. B 55, 9452 (1997).

    ADS  CAS  Article  Google Scholar 

  48. Green, A. G., Conduit, G. & Krüger, F. Quantum order-by-disorder in strongly correlated metals. Annu. Rev. Condens. Matter Phys. 9, 59 (2018).

    ADS  CAS  Article  Google Scholar 

  49. Kasahara, S. et al. Giant superconducting fluctuations in the compensated semimetal fese at the bcs-bec crossover. Nat. Commun. 7, 12843 (2016).

    ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

  50. Oh, M. et al. Evidence for unconventional superconductivity in twisted bilayer graphene. Nature 600, 240 (2021).

    ADS  CAS  PubMed  Article  Google Scholar 

  51. Lee, P. A., Nagaosa, N. & Wen, X.-G. Doping a mott insulator: Physics of high-temperature superconductivity. Rev. Mod. Phys. 78, 17 (2006).

    ADS  CAS  Article  Google Scholar 

  52. Sandvik, A. W. Stochastic method for analytic continuation of quantum monte carlo data. Phys. Rev. B 57, 10287 (1998).

    ADS  CAS  Article  Google Scholar 

  53. Beach, K. S. D., Identifying the maximum entropy method as a special limit of stochastic analytic continuation. Preprint at https://arxiv.org/abs/cond-mat/0403055.

  54. Sandvik, A. W. Constrained sampling method for analytic continuation. Phys. Rev. E 94, 063308 (2016).

    ADS  PubMed  Article  Google Scholar 

  55. Shao, H. et al. Nearly deconfined spinon excitations in the square-lattice spin-1/2 heisenberg antiferromagnet. Phys. Rev. X 7, 041072 (2017).

    Google Scholar 

  56. Sun, G.-Y. et al. Dynamical signature of symmetry fractionalization in frustrated magnets. Phys. Rev. Lett. 121, 077201 (2018).

    ADS  CAS  PubMed  Article  Google Scholar 

  57. Ma, N. et al. Dynamical signature of fractionalization at a deconfined quantum critical point. Phys. Rev. B 98, 174421 (2018).

    ADS  CAS  Article  Google Scholar 

  58. Zhou, C. et al. Amplitude mode in quantum magnets via dimensional crossover. Phys. Rev. Lett. 126, 227201 (2021).

    ADS  CAS  PubMed  Article  Google Scholar 

  59. Zhang, X., Pan, G., Zhang, Y., Kang, J. & Meng, Z. Y. Momentum space quantum monte carlo on twisted bilayer graphene. Chin. Phys. Lett. 38, 077305 (2021).

    ADS  CAS  Article  Google Scholar 

  60. Pan, G., Zhang, X., Li, H., Sun, K. & Meng, Z. Y. Dynamic properties of collective excitations in twisted bilayer Graphene. Phys. Rev. B 105, L121110 (2022).

Download references

Acknowledgements

We thank R.M. Fernandes, M.H. Christensen, Y. Schattner, E. Berg, and X. Wang for valuable discussions. W.L.J. thanks Z. Liu for the support of the code, and G. Pan for the helpful suggestions. W.L.J., Y.Z.L. and Z.Y.M. acknowledge support from the RGC of Hong Kong SAR of China (Grant Nos. 17303019, 17301420, 17301721 and AoE/P-701/20), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB33000000), the K. C. Wong Education Foundation (Grant No. GJTD-2020-01) and the Seed Funding “Quantum-Inspired explainable-AI" at the HKU-TCL Joint Research Centre for Artificial Intelligence. We thank the Center for Quantum Simulation Sciences in the Institute of Physics, Chinese Academy of Sciences, the Computational Initiative at the Faculty of Science and the Information Technology Services at the University of Hong Kong and the Tianhe platforms at the National Supercomputer Centers in Tianjin and Guangzhou for their technical support and generous allocation of CPU time. The authors also acknowledge Beijng PARATERA Tech CO.,Ltd.(https://www.paratera.com/) for providing HPC resources that have contributed to the research results reported within this paper. The work by A.V.C. was supported by the Office of Basic Energy Sciences, U.S. Department of Energy, under award DE-SC0014402. A.K. and A.V.C. acknowledge the hospitality of KITP at UCSB, where part of the work has been conducted. The research at KITP is supported by the National Science Foundation under Grant No. NSF PHY-1748958. Y.W. is supported by startup funds at the University of Florida and by NSF under award number DMR-2045871.

Author information

Authors and Affiliations

Authors

Contributions

A.K., Y.W., K.S., A.V.C. and Z.Y.M. initiated the work. W.L.J. and Y.Z.L. developed the program and performed the QMC calculations, W.L.J., Y.Z.L. and Z.Y.M. carried out the numerical data analysis. A.K., Y.W., K.S. and A.V.C. performed the theory analysis. All authors wrote the manuscript together.

Corresponding author

Correspondence to Zi Yang Meng.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Nature Communications thanks the anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Jiang, W., Liu, Y., Klein, A. et al. Monte Carlo study of the pseudogap and superconductivity emerging from quantum magnetic fluctuations. Nat Commun 13, 2655 (2022). https://doi.org/10.1038/s41467-022-30302-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1038/s41467-022-30302-x

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Search

Quick links

Nature Briefing

Sign up for the Nature Briefing newsletter — what matters in science, free to your inbox daily.

Get the most important science stories of the day, free in your inbox. Sign up for Nature Briefing