Resolving Quantum Criticality in the Honeycomb Hubbard Model

TL;DR

The paper resolves the critical exponents of the Gross-Neveu-Heisenberg transition in the honeycomb Hubbard model by pushing projector determinant QMC to unprecedented system sizes using a novel submatrix-T update. It demonstrates that the fermion anomalous dimension $\eta_\psi$ and the correlation-length exponent $\nu$ converge rapidly, while the boson anomalous dimension $\eta_\phi$ requires linear extrapolation to address finite-size effects; the spinless $t$-$V$ model in the Gross-Neveu-Ising class provides an important cross-check that agrees with conformal bootstrap. A sliding-window Bayesian data-collapse workflow is developed to extract $U_c$, $V_c$, $\nu$, $\eta_\phi$, and $\eta_\psi$ with controlled uncertainties, revealing observable-dependent convergence patterns. The results offer state-of-the-art critical exponents and a generalizable methodology for resolving fermionic quantum critical phenomena in strongly correlated systems, with potential implications for graphene-like materials and beyond.

Abstract

The interplay between Dirac fermions and electronic correlations on the honeycomb lattice hosts a fundamental quantum phase transition from a semimetal to a Mott insulator, governed by the Gross-Neveu-Heisenberg (GNH) universality class. Despite its importance, consensus on the precise critical exponents remains elusive due to severe finite-size effects in numerical simulations and the lack of conformal bootstrap benchmarks. Here we try to resolve this long-standing controversy by performing projector determinant quantum Monte Carlo (QMC) simulations on lattices of unprecedented size, reaching 10,368 sites. By developing a novel projected submatrix update algorithm, we achieve a significant algorithmic speedup that enables us to access the thermodynamic limit with high precision. We observe that the fermion anomalous dimension and the correlation length exponent converge rapidly, while the boson anomalous dimension exhibits a systematic size dependence that we resolve via linear extrapolation. To validate our analysis, we perform parallel large-scale simulations of the spinless $t$-$V$ model on the honeycomb lattice, which belongs to the Gross-Neveu-Ising class. Our results for the $t$-$V$ model, including the first QMC determination of the fermion anomalous dimension, show agreement with conformal bootstrap predictions, thereby corroborating the robustness of our methodology. Our work provides state-of-the-art critical exponents for the honeycomb Hubbard model and establishes a systematic finite-size scaling workflow applicable to a broad class of strongly correlated quantum systems, paving the way for resolving other challenging fermionic quantum critical phenomena.

Figures (6)

• Figure 1: Summary of critical exponents $\nu$, $\eta_\phi$, and $\eta_\psi$ reported in the literature and obtained in this work. a, $N=8$ Gross-Neveu-Heisenberg universality class Phys.Rev.D2017ZerfPhys.Rev.B2023LadovrechisPhys.Rev.D2018GraceyPhys.Rev.B2018KnorrArXiv2025Yu-NonequilibriumDynamicsDiracNat.Commun.2025ZengPhys.Rev.B2021OstmeyerPhys.Rev.B2019BuividovichPhys.Rev.X2016OtsukaPhys.Rev.B2015ParisenToldinPhys.Rev.B2025LangPhys.Rev.Lett.2021XuPhys.Rev.B2021LiuPhys.Rev.B2020OtsukaScience2018TangPhys.Rev.B2018Buividovich. b, $N=4$ Gross-Neveu-Ising universality class Phys.Rev.B2018IhrigPhys.Lett.B1992GraceyInt.J.Mod.Phys.A1994GraceyaInt.J.Mod.Phys.A1994GraceyPhys.Rev.B2016KnorrJ.HighEnerg.Phys.2023ErramilliNat.Commun.2025ZengPhys.Rev.B2016HesselmannPhys.Rev.B2016WangNewJ.Phys.2015LiNewJ.Phys.2014WangPhys.Rev.D2020Huffman. The gray horizontal dashed lines divide the plot into three blocks: analytical results (top), QMC results for the same model as in this work (middle), and QMC results for other models in the same universality class (bottom). CB denotes conformal bootstrap J.HighEnerg.Phys.2023Erramilli. See Supplemental Material Tables S5 and S6 for the numerical values and methodological details.
• Figure 2: Phase diagram and crossing-point analysis of the correlation ratio. a, Schematic phase diagram of the interaction-driven semimetal-Mott insulator transition with a Gross-Neveu-Heisenberg critical point. b, Correlation ratio $R_{\mathrm{AFM}}^{(1,0)}$ as a function of $U$ for various system sizes $L$. Horizontal black error bars mark the crossing points $U^\times$ of size pairs $(L,2L)$, determined by quadratic interpolation of the data for each size. Inset: power-law extrapolation of the crossing points.
• Figure 3: Sliding-window data-collapse analysis for the Hubbard model. Columns show three representative fit windows with $L_{\mathrm{max}}=36$, 54, and 72 (see Supplemental Material SM-submatrixLR Section S6 for all windows). The $L$-axis above each panel indicates the sizes included in the fit, with greyed-out sizes excluded. Rows: $R_{\mathrm{AFM}}^{(1,0)}$ collapse yielding $(U_c,\nu)$ (top), $m^2_{\mathrm{AFM}}$ collapse yielding $\eta_\phi$ (middle), and $\langle|G_{AB}|\rangle$ collapse yielding $\eta_\psi$ (bottom).
• Figure 4: Fitting-window dependence of critical parameters for the Hubbard and spinless $t$-$V$ models. a--d, Hubbard model; e--h, $t$-$V$ model. Each point represents the fitted value of $U_c$ (or $V_c$), $\nu$, $\eta_\phi$, or $\eta_\psi$ obtained from a fitting window of six consecutive sizes, with $L_{\mathrm{max}}$ denoting the largest size in each window. Horizontal lines and shaded bands indicate final estimates and confidence intervals; dashed lines show linear extrapolations to $1/L_{\mathrm{max}}\to 0$.
• Figure 5: Schematic illustration of the submatrix-T update algorithm. At a specific time slice, there are $N/k$ proposed local updates in total, which is divided into several delayed blocks of size $n_d$. Within each delayed block, we do not update $T$ (indicated by dashed arrows) but instead perform a recursive update of the submatrix$\Gamma^{-1}$ (indicated by solid arrows in black), which helps in calculating not only the intermediate acceptance ratios $r^{(i)}$ but also the final $T^{(n_a)}$ (indicated by solid arrows in grey), where $n_a$ denotes the total accepted moves in the current delayed block.