Mean-field analysis of a Hubbard interaction on Bernal Bilayer Graphene

We perform unrestricted Hartree-Fock calculations on the 2D Hubbard model on a honeycomb and bilayer honeycomb lattice at both zero and finite temperatures. Finite size real space calculations are supplemented with RPA calculations in the thermodynamic limit. Our motivation comes from high doping levels achieved in graphene and Bernal bilayer graphene by interacalation. We present phase diagrams in doping and temperature for a moderate Hubbard interaction. The magnetic states we find are classified systematically based on the dominant Fourier components of their spin patterns, their average magnetization and spin incommensurabilities. The dominant spin patterns are Néel order and various types of stripes. Around Van Hove filling, we resolve the competition between stripe and chiral spin density waves in the symmetry-broken regime. We also investigate the effect of an applied external displacement field on the spin patterns of BBG.

💡 Research Summary

In this work the authors investigate the magnetic and charge‑ordered phases of doped single‑layer graphene and Bernal‑stacked bilayer graphene (BBG) within the two‑dimensional Hubbard model. The study combines unrestricted Hartree‑Fock (HF) calculations performed in real space on finite lattices with a Random Phase Approximation (RPA) analysis carried out directly in the thermodynamic limit. The motivation stems from recent experimental advances that enable very high carrier densities in graphene and BBG through intercalation, as well as from the possibility to apply a perpendicular electric displacement field that creates an asymmetry between the two layers of BBG.

Model and computational approach
The Hamiltonian consists of a kinetic part H₀ that includes nearest‑neighbour hopping t for the monolayer and, for BBG, additional inter‑layer hopping t⊥ together with skew hoppings t₃ and t₄. An on‑site repulsive interaction U>0 is added, and a displacement field D is modeled by a term (D/2)∑_j ε_j n_j with ε_j = +1 on the upper layer and –1 on the lower layer. The authors explore a broad range of electron fillings n (from empty to fully filled bands) and temperatures T, focusing on a moderate interaction strength (U comparable to the bandwidth).

In the HF part, the interaction is decoupled in both Hartree and Fock channels, yielding site‑dependent order parameters Δ_{jσ}=U⟨n_{jσ}⟩ and Δ_{j+}=−U⟨c†{j↑}c{j↓}⟩. Self‑consistency is achieved on lattices of 18×18 unit cells for the monolayer (648 sites) and 12×12 unit cells for BBG (576 sites) with periodic boundary conditions. The chemical potential μ is adjusted each iteration to keep the average density fixed. Multiple random initial configurations and continuation from neighboring points in parameter space are used to avoid trapping in metastable minima. After convergence, the free energy per site is evaluated and the lowest‑energy solution is identified.

The RPA analysis computes the non‑interacting Green’s function G₀(k) and the static susceptibility χ₀(q)=∑_k G₀(k+q)G₀(k). The interacting susceptibility χ(q)=χ₀(q)