Absence of a Spin Liquid Phase in the Hubbard Model on the Honeycomb Lattice

Sandro Sorella; Yuichi Otsuka; Seiji Yunoki

Source-linked AI summary

Absence of a Spin Liquid Phase in the Hubbard Model on the Honeycomb Lattice

Sandro Sorella, Yuichi Otsuka, Seiji Yunoki

arXiv:1207.1783v2Jul 7, 2012cond-mat.str-el

TL;DR

The paper tests whether the honeycomb-lattice Hubbard model hosts a stable spin liquid between a semi-metal and antiferromagnetic insulator. Using auxiliary-field Monte Carlo on clusters up to 2592 sites, it instead finds strong support for a direct and continuous transition between these phases.

Problem
Whether a stable spin liquid exists between the semi-metal and antiferromagnetic-insulator phases of the honeycomb-lattice Hubbard model remains debated.
Method
The study uses auxiliary-field Monte Carlo simulations on symmetry-preserving periodic clusters containing up to 2592 sites.
Results
The possible spin-liquid region shrinks to 3.8t ⪅ U/t ⪅ 3.9t, while charge-gap opening coincides with antiferromagnetic order, supporting a direct and continuous SM-to-AFMI transition.
Takeaways & Limitations
The calculations favor the conventional scenario in which antiferromagnetic order appears concomitantly with insulating behavior rather than across a broad spin-liquid phase.
Takeaways & Limitations
The broader claim about SU(2)-invariant spin-1/2 models is not established universally, and further numerical study is required to identify ingredients that stabilize spin liquids in realistic electronic models.

Abstract

from arXiv · show

A spin liquid is a novel quantum state of matter with no conventional order parameter where a finite charge gap exists even though the band theory would predict metallic behavior. Finding a stable spin liquid in two or higher spatial dimensions is one of the most challenging and debated issues in condensed matter physics. Very recently, it has been reported that a model of graphene, i.e., the Hubbard model on the honeycomb lattice, can show a spin liquid ground state in a wide region of the phase diagram, between a semi-metal (SM) and an antiferromagnetic insulator (AFMI). Here, by performing numerically exact quantum Monte Carlo simulations, we extend the previous study to much larger clusters (containing up to 2592 sites), and find, if any, a very weak evidence of this spin liquid region. Instead, our calculations strongly indicate a direct and continuous quantum phase transition between SM and AFMI.

Introduction

The paper examines whether the half-filled honeycomb-lattice Hubbard model hosts a stable spin liquid between a semimetal and an antiferromagnetic insulator. Larger-cluster simulations challenge the previously reported broad spin-liquid interval and instead support antiferromagnetic order emerging with insulating behavior.

Introduction: A spin liquid is a Mott insulator that is not adiabatically connected to a band insulator and preserves all symmetries at zero temperature.
Introduction: The central issue is whether spin degrees of freedom can avoid Bose-Einstein condensation and crystallization, both necessary for stable spin-liquid behavior without long-range order.
Introduction: The half-filled honeycomb Hubbard model describes electrons at one per site and is used here to study its insulating ground state.
Introduction: The unfrustrated model is known to evolve from a semimetal to a classically Néel-ordered antiferromagnetic insulator as U/t increases.
Introduction: Meng et al. reported a possible spin-liquid phase for 3.4 ⪅ U/t ⪅ 4.3 between the semimetal and antiferromagnetic insulator using clusters up to 648 sites.
Introduction: Using clusters up to 2592 sites, this study finds the possible spin-liquid region reduced to 3.8t ⪅ U/t ⪅ 3.9t, if it exists, and supports a conventional SM-to-AFMI transition.

Results

The simulations test magnetic order, spin gaps, and charge correlations across interaction strengths and cluster sizes. They find finite antiferromagnetic order, a critical interaction near Uc/t = 3.869, and charge-gap onset coincident with the magnetic transition, supporting a direct continuous SM–AFMI transition.

Results: Clusters with N = 2L^2 sites and periodic boundaries preserve lattice symmetries, with L reaching 36, while auxiliary-field Monte Carlo evaluates ground-state observables.
Results: At U/t = 4, finite-size scaling of long-distance spin correlations yields a finite antiferromagnetic order parameter in the thermodynamic limit.
Results: The extrapolated spin structure factor and maximum-distance spin correlations remain consistent with antiferromagnetic long-range order, contrasting with a spin-disordered ground state.
Results: β ≃ 0.8 describes the approximately linear scaling of the staggered magnetization, close to the Hartree-Fock prediction β = 1.
Results: Uc/t = 3.869 ± 0.013 is the best estimate for the interaction where antiferromagnetic order melts, below the previously reported value ≈4.3.
Results: The spin gap is zero within statistical errors across the studied U/t values, with an error as small as 0.004t at U/t = 4.
Results: Charge correlations change from power-law to exponential behavior around Uc, with uncertainty below 0.1t, indicating charge-gap opening at the antiferromagnetic transition.

Discussion

The authors argue that trivial ground-state phases in this unfrustrated model favor antiferromagnetic order once a charge gap opens, while acknowledging limits to generalizing this criterion.

Discussion: The Marshall sign rule makes the ground-state phases trivial in the no-double-occupancy sector of this unfrustrated model.The authors contrast this with highly nontrivial, entangled phases in established spin-liquid models.
Discussion: Finite-size extrapolations of SAF and Cs(Lmax) give consistent antiferromagnetic order estimates, with the extrapolated Cs(Lmax) providing an accurate lower bound.The fits remain stable when the largest or smallest cluster is removed, and the two estimates agree within two standard deviations.
Discussion: Magnetic long-range order is expected to occur once the charge gap becomes finite because Bose-Einstein condensation is difficult to avoid.
Discussion: The best-fit critical exponent is β = 0.80 ± 0.04, while the estimated transition lies at Uc/t = 3.869 ± 0.013.The critical Uc remains between 3.8 and 3.9 for the compared exponent choices.
Discussion: The authors caution that highly nontrivial and entangled ground-state phases are expected for a true two-dimensional spin liquid, but this observation is not established generically beyond the stated scope.They specifically note the importance of restricting the discussion to SU(2)-invariant models and call for further numerical study.

Methods

The simulations use auxiliary-field Monte Carlo with optimized trial states and direct singlet–triplet energy calculations, while resampling quantifies extrapolation errors.

Methods: The calculation uses different left and right trial wave functions, including a Gutzwiller projection to improve efficiency.The projection contains a variational parameter g that reduces statistical errors for suitable values.
Methods: The spin gap is evaluated directly as Δs = E(S = 1) − E(S = 0) by separately simulating singlet and triplet sectors.Improved estimators reduce the statistical errors without encountering a negative sign problem.
Methods: Different symmetry choices for the left and right trial states accelerate convergence while addressing degeneracies at momenta K and K′.
Methods: Resampling generates normally distributed fictitious data for each cluster, repeatedly fits the finite-size form, and estimates extrapolated values and errors.The procedure uses weighted least-squares fits and Mrs = 200 repetitions.
Methods: Finite-size scaling of spin gaps and charge correlations distinguishes the semimetallic 1/L behavior from antiferromagnetic 1/L2 behavior and charge-gap signatures.

Supplementary information

The supplementary information defines the half-filled Hubbard Hamiltonian and identifies hopping and on-site Coulomb interaction terms.

Supplementary information: The Hubbard model uses electron creation operators c†R,σ and densities nR,σ = c†R,σcR,σ for spin σ = ↑, ↓.
Supplementary information: The Hamiltonian sums over nearest-neighbor sites, with electron hopping in H0 and one-site Coulomb interaction in HI.

I. HONEYCOMB LATTICE

The supplementary construction uses a periodic honeycomb lattice with two sublattices per unit cell, symmetry-preserving clusters, and sizes chosen to include the Dirac-point momenta.

I. HONEYCOMB LATTICE: Each honeycomb unit cell contains two sites belonging to A and B sublattices, positioned using primitive vectors τ1 and τ2.
I. HONEYCOMB LATTICE: Periodic boundary conditions identify lattice points differing by the lattice vectors T1 and T2, producing N = 2L2 sites.The chosen finite lattices preserve the point-group symmetries of the infinite lattice.
I. HONEYCOMB LATTICE: The maximum-distance sites lie on the same sublattice at Lmax = L, with the distance defined using 120-degree rotation symmetry.
I. HONEYCOMB LATTICE: Only clusters with L a multiple of 3 contain the high-symmetry momenta K and K′ where the Dirac points occur.The study therefore chooses L values that are multiples of 3 for systematic finite-size scaling.

II. TRIAL WAVE FUNCTIONS

The simulations use auxiliary-field quantum Monte Carlo with symmetry-breaking trial wave functions and projection to the ground state. Separate singlet and triplet calculations avoid the sign problem and enable direct spin-gap estimation.

II. TRIAL WAVE FUNCTIONS: Auxiliary-field quantum Monte Carlo statistically evaluates projected expectation values using left and right trial wave functions.The exact ground-state expectation value is obtained in the limits τ → ∞ and ∆τ → 0.
II. TRIAL WAVE FUNCTIONS: The left trial function includes Gutzwiller projection and an antiferromagnetic mean-field order parameter while conserving spatial rotational symmetry.The parameters g and ∆ are optimized to accelerate convergence with projection time.
II. TRIAL WAVE FUNCTIONS: The right and left trial wave functions break different symmetries, improving convergence for symmetric operators through the larger singlet-sector gap.The left function breaks spin-rotation symmetry, while the right function breaks spatial symmetry.
II. TRIAL WAVE FUNCTIONS: Independent singlet and triplet simulations have no negative sign problem, allowing the spin gap to be estimated directly and accurately.This avoids relying on the asymptotic extraction of imaginary-time spin correlations.
II. TRIAL WAVE FUNCTIONS: The mixed-average estimator reduces energy statistical errors through the zero-variance principle when the left trial function is exact.Its statistical error becomes zero for an exact |ψL⟩.

III. CORRELATION FUNCTIONS AND SPIN GAP

The study analyzes long-distance spin correlations, magnetic structure factors, and the spin gap using finite-size extrapolations. Projection-time convergence and fit stability support the reported correlation estimates.

III. CORRELATION FUNCTIONS AND SPIN GAP: The maximum-distance spin correlation Cs(Lmax) averages spin products over symmetry-equivalent vectors connecting sites on the same sublattice.The corresponding distance is Lmax = |τmax|.
III. CORRELATION FUNCTIONS AND SPIN GAP: Both Cs(Lmax) and SAF converge by τt = L + 4 for L = 6, 18, and 36, including L = 36 near Uc/t ∼ 3.87.Their monotonic increase means finite projection time can at most underestimate the magnetic order parameter.
III. CORRELATION FUNCTIONS AND SPIN GAP: Spin structure factors and maximum-distance correlations are evaluated across system sizes and extrapolated to L → ∞ using cubic fits in 1/L.Three fitting choices test sensitivity to including the smallest and largest clusters.
III. CORRELATION FUNCTIONS AND SPIN GAP: The three cubic extrapolation fits give quantitatively consistent results within statistical errors, indicating stable thermodynamic estimates.Removing either the L = 36 or L = 6 data produces no material change within uncertainty.
III. CORRELATION FUNCTIONS AND SPIN GAP: The thermodynamic estimates from Cs(Lmax) and SAF are statistically consistent, although removing the smallest sizes slightly increases the extrapolated order parameter.The same trend is also observed with quadratic fits.
III. CORRELATION FUNCTIONS AND SPIN GAP: The spin gap is obtained from separate singlet and triplet ground-state energies, ∆s = E(S = 1) − E(S = 0).This definition uses the independently calculated energies of the two spin sectors.