Emergent relativistic symmetry from interacting fermions on the honeycomb bilayer

We study the phase diagram of interacting spinless fermions on the honeycomb bilayer at charge neutrality using large-scale quantum Monte Carlo simulations. In the noninteracting limit, the low-energy spectrum features quadratically dispersing bands that touch at the corners of the hexagonal Brillouin zone. Weak to intermediate interactions induce a splitting of each of the quadratic band touching points into four Dirac points, located along high-symmetry directions of the reciprocal lattice. Strong interactions lead to the formation of a layer-polarized charge density wave, which spontaneously breaks the $\mathbb Z_2$ layer inversion symmetry and opens an insulating gap in the spectrum. We show that the semimetal-to-insulator quantum phase transition as a function of interaction is continuous and characterized by emergent relativistic symmetry. Our results for the values of the correlation-length exponent $ν$, the order-parameter anomalous dimension $η_ϕ$, and the fermion anomalous dimension $η_ψ$ agree with those of the theoretically predicted 2+1D Gross-Neveu-Ising universality class with eight two-component Dirac fermions within less than 5%\ deviation. We also determine the crossover scale as a function of interaction strength between the nonrelativistic semimetal state at high temperatures, characterized by dynamical critical exponent $z = 2$, and the Dirac semimetal state at intermediates temperatures, characterized by $z=1$. Further reducing the temperature below the crossover scale at a fixed value of the interaction strength above the quantum critical point results in a classical ordering transition in the 2D Ising universality class.

💡 Research Summary

In this work the authors investigate a model of spinless fermions on a Bernal‑stacked honeycomb bilayer at half‑filling using large‑scale determinant quantum Monte Carlo (DQMC) simulations that are free of the sign problem thanks to Majorana reflection positivity. The non‑interacting band structure consists of two quadratically dispersing bands that touch at the inequivalent K and K′ points (quadratic band touching, QBT), yielding a finite density of states at the Fermi level. Weak to intermediate nearest‑neighbor repulsion V is predicted by renormalization‑group analyses to split each QBT point into four Dirac cones located along high‑symmetry directions; strong V is expected to drive a layer‑polarized charge‑density‑wave (CDW) that breaks the Z₂ layer‑inversion symmetry and opens a gap.

Zero‑temperature (projective) QMC simulations on lattices up to 4 × 27² sites are used to locate the quantum critical point (QCP). The authors monitor a charge‑correlation ratio R_c(V) that is RG‑invariant; its size‑dependent crossing yields V_c ≈ 0.900(5). Finite‑size scaling of the crossing points, together with Bayesian extrapolation, provides the critical exponents: correlation‑length exponent ν = 0.938(4) (1/ν = 1.065(48)), bosonic anomalous dimension η_ϕ = 0.856(42) (with dynamical exponent z = 1), and fermionic anomalous dimension η_ψ = 0.0199(46). These values agree with the 2+1‑dimensional Gross‑Neveu‑Ising universality class with eight two‑component Dirac flavors (theoretical estimates: 1/ν ≈ 1.018, η_ϕ ≈ 0.868, η_ψ ≈ 0.0195) within five percent. The agreement demonstrates that, despite the underlying non‑relativistic QBT dispersion, the low‑energy critical theory acquires emergent Lorentz invariance (z = 1).

Finite‑temperature simulations reveal two distinct crossover scales. At low temperature the CDW order melts via a conventional two‑dimensional Ising transition; the authors extract the Ising critical temperature T_Ising by scaling the rescaled charge susceptibility L^{η_Ising‑2} χ with η_Ising = 0.25 and ν_Ising = 1. At higher temperature the fermionic spectrum reconstructs from a Dirac‑like (z = 1) regime to the original QBT (z = 2) regime. This crossover temperature T_cross is identified by fitting the uniform charge susceptibility χ_uni(T) to a linear (Dirac) form at intermediate T and a logarithmic (QBT) form at high T, then locating the temperature where the residuals of the two fits change sign. T_cross grows with increasing V, indicating that stronger interactions push the non‑relativistic regime to higher temperatures. Spectral functions A(k, ω) obtained via maximum‑entropy analytic continuation corroborate the picture: for T < T_cross clear Dirac cones are visible, while for T > T_cross the dispersion becomes quadratic.

Overall, the study provides a comprehensive numerical confirmation that (i) weak interactions split QBT points into Dirac cones, (ii) a continuous semimetal‑to‑insulator transition occurs at V_c, (iii) this transition belongs to the Gross‑Neveu‑Ising universality class with emergent relativistic symmetry, and (iv) the finite‑temperature phase diagram contains both a Dirac‑to‑QBT crossover (z = 1 ↔ z = 2) and a 2D Ising ordering transition. The results are directly relevant for experimental platforms such as gated bilayer graphene, twisted double‑layer systems, or engineered cold‑atom lattices where the interaction strength and temperature can be tuned to explore the predicted quantum critical behavior and the associated emergent symmetries.