The Hubbard model: A computational perspective

1. Introduction

• The Hubbard model combines lattice hopping with local interactions and supports diverse phases, making it a useful testbed as numerical methods increasingly yield consistent results.Its single-orbital form emphasizes local electron-correlation physics, while omitting non-local interactions, band structure, and multiple orbitals.Source: §1 · pp. 2–3

1.1. Purpose and structure of this article

• The review focuses on the single-band Hubbard model on the two-dimensional square lattice across weak, intermediate, and strong coupling, with connections to the t-J model and three dimensions.It defines these regimes as U/t ≲4, 2 ≲ U/t ≲6, and U/t ≳6, respectively, and reviews half-filling, doped systems, experimental probes, extensions, and open questions.Source: §1.1 · p. 4

1.2. ‘Solving’ the Hubbard model

• At U/t ∼8 and 1/8 doping, recent agreement among several methods finds stripe order but no dx2−y2 superconductivity, although the stripe-order magnitude still needs greater precision.The article defines ‘solving’ as obtaining a physically interesting observable with accuracy comparable to or better than comparison calculations or experiment; what that requires depends on the observable and regime.Source: §1.2 · p. 4
• Because every approach relies on approximate solutions and extrapolation, multi-method consensus—agreement within uncertainties across methods with different approximations—helps establish reliable results.Controlled methods can be difficult to extrapolate or used outside their small-parameter regimes, while uncontrolled methods may still be useful when residual errors are modest (152, 224).Source: §1.2 · pp. 4–5

1.3. Connection of the model to experiments

• Two-dimensional square-lattice simulations show cuprate-like phases and trends with doping and temperature, but precise correspondence between cuprates and the single-orbital model is difficult.Simulated phases include dx2−y2 superconductivity, antiferromagnetic correlations, stripes, pseudogaps, Fermi-liquid behavior, and bad metallic behavior.Source: §1.3 · p. 5

1.3.1. Cuprate physics.

• Fits associate cuprates with U/t ∼8 and t′/t from −0.1 to −0.4 for most hole-doped materials, but the fits are neither systematic nor controlled.A more accurate description may require multiple orbitals, non-local interactions, and additional non-local hopping terms.Source: §1.3.1 · p. 6

1.3.2. Quantum emulators and quantum computers.

• Ultracold-atom and quantum-dot emulators have yielded interesting time-dependent results but not equilibrium predictive power matching numerical simulations; current quantum-computer calculations lag classical capabilities.The authors expect Hubbard-model results in interesting regimes from quantum computers before quantum-computing solutions to materials or chemistry problems (272).Source: §1.3.2 · p. 6

Phase diagram and phase boundaries

• The precise intermediate-to-strong-coupling phase boundaries remain debated, although there is broad consensus on phases and their approximate locations.The weak-coupling ground-state phase diagram is fairly well understood (49, 266).Source: §2 · pp. 6–7
• Diverging correlation lengths, large distinct unit cells, tiny energy differences, and broad crossovers hinder precise boundaries; approximate methods can also overemphasize ordered phases.Finite-size methods may predict spurious nonzero-temperature magnetic phases, forbidden by the Mermin-Wagner theorem (176), which disappear only slowly at very large sizes (170).Source: §2 · pp. 6–7
• Phase diagrams such as Fig. 1 represent results within a particular approximation or are inspired by simple theories or experimental target systems.The precise shape of transition lines is often unknown.Source: §2 · p. 7
• Many orders and phenomena coexist or lie close in energy or parameter space, making phase locations and properties difficult to pin down; additional proposed phases include nematic phases, reviewed in Ref. (7).This situation is described as “intertwined orders” (66).Source: §2 · p. 7

3.1. The two-dimensional model at half filling

• At half filling on bipartite lattices, the fermion sign problem is absent, enabling efficient large-lattice QMC or quantum-embedded cluster simulations whose results can be extrapolated to the thermodynamic limit (99, 165, 152, 224).Other numerical methods also provide precise results across a wide range of physical quantities (152, 224).Source: §3.1 · p. 7

3.1.1. Static quantities: energies, entropies, and spatial correlation functions.

• The weak-coupling magnetic correlation length grows exponentially on cooling, contributing to the destruction of finite-temperature magnetic order and the development of a pseudogap near the AF ground state (27, 284, 181, 221, 211, 224; 176).Thermodynamic-limit ground-state energies, double occupancies, and AF order parameters were obtained through finite-size scaling with twist-averaged boundary conditions (205); short-range correlations and finite-temperature quantities were also studied (206; 259; 137; 154).Source: §3.1.1 · pp. 7–8

3.1.2. Dynamical quantities:

• The single-particle gap persists at all interactions, decaying exponentially at weak coupling and scaling linearly at strong coupling; tensor-network spectra agree with earlier QMC results at higher resolution (263; 21, 202, 201, 76).Imaginary-time or Matsubara methods require ill-conditioned analytic continuation, while tensor-network momentum-resolved spectra remain difficult to obtain across the full Brillouin zone (122; 99; 161; 141).Source: §3.1.2 · p. 8

3.1.3. Nature of the gap.

• The square-lattice ground state is insulating with AF Néel order at all U; its gap changes from Slater-like at weak coupling to Mott-like at strong coupling, crossing over near U/t ≈4 without a phase transition (99, 284, 221, 263, 265; 101, 259; 87, 211, 137; 263, 137).On cooling, the system crosses from a high-temperature incoherent regime through a coherent metallic regime to a low-temperature insulating AF pseudogap, while the inverse quasiparticle lifetime has a temperature minimum (221, 265, 137; 211, 223).Source: §3.1.2 · p. 8, §3.1.3 · pp. 8–9

3.2. The half-filled model in 3D

• In 3D, the half-filled model transitions from a high-temperature paramagnet to a low-temperature AF insulator, with maximum TN/t∼0.35 near U/t∼8 and a high-temperature metal–insulator crossover near U/t ∼6 −8 (220, 242, 136, 212, 210; 242, 69).The highest critical entropy is S ∼0.65 near U ∼8t; experiments reached entropy 0.77 and observed AF correlations (68; 129; 75, 91).Source: §3.2 · p. 9

4.1. Magnetic and superconducting properties

• Away from half-filling, studies find incommensurate magnetic order, and fRG indicates it can persist to fairly high doping when superconductivity is suppressed.Near the van Hove singularity, spin-density-wave, pairing, and ferromagnetic instabilities compete; at finite temperature, a superconducting dome centered around optimal doping has been found.Source: §4.1 · p. 10
• D-wave superconductivity is supported by perturbative, DiagMC, and fRG studies; at finite next-nearest-neighbor hopping it can coexist with antiferromagnetism across a broad doping range.fRG finds a sizable gap, while DMF2RG finds robust d-wave pairing at strong coupling, driven by magnetic correlations at the edge of the antiferromagnetic regime.Source: §4.1 · p. 10

4.2. Weak-coupling pseudogap

• At weak to intermediate coupling, long-range antiferromagnetic correlations induce a pseudogap when the correlation length exceeds the thermal de Broglie wavelength, ξ ≫ vF/(πT).The gap is predicted near hot spots, opening first at the antinode as temperature falls and spreading toward the node.Source: §4.2 · p. 10
• The pseudogap occurs only for a hole-like Fermi surface, where low-energy quasiparticle coherence is suppressed at hot spots; electron-like topology restores coherence along the Fermi surface.The Lifshitz transition from hole-like to electron-like topology controls both the location of self-energy singularities and the Fermi-surface topological transition.Source: §4.2 · p. 11

5.1. Competition of low-energy ground states: uniform vs stripe states

• Uniform d-wave states and charge- or spin-modulated stripes are near-degenerate competitors; a PDW has a π phase shift between stripes and was proposed to explain suppressed superconductivity in LBCO near δ = 1/8.Stripes can coexist with d-wave superconductivity, while early Hartree-Fock studies found insulating stripes with one hole per unit length.Source: §5.1 · p. 11

5. The doped 2D Hubbard model at intermediate-to-strong coupling

• At U/t = 8 and δ = 1/8, the period-8 stripe without d-wave superconductivity is the ground state; the uniform d-wave state is ≈0.01t higher in energy.Four methods estimate the period-8 stripe energy at −0.767±0.004t, and subsequent studies confirmed the stripe result; the earlier VMC discrepancy arose from testing a higher-energy period-4 stripe.Source: §5.1 · pp. 11–12, §5.1.1 · p. 13

5.1. Competition of low-energy ground states: uniform vs stripe states

• A period-8 stripe remains stable for U/t ∼ 6–12, with four methods estimating its energy in the range −0.767±0.004t; prior VMC’s period-4 stripe was substantially higher than period-5–8 stripes.Source: §5.1 · p. 12
• Away from δ = 1/8, stripe periods tend to decrease as doping rises, while the presence of superconductivity and low-doping phase separation remain method- and regime-dependent.For U/t = 8, t′ = 0, superconductivity in stripe ground states was reported absent over δ ∼0.1–0.2 and U/t ∼6–8; phase separation was found near δ ≈0.08, but still awaits confirmation with tensor-network methods or CP-AFQMC.Source: §5.1 · p. 12
• Uniform d-wave stability is predicted at higher doping in the simplest Hubbard model, while long-range Coulomb interactions in an Hg-based effective model extend it to δ ≳0.1.Predicted uniform-state regions include 0.20 ≲δ ≲0.27 for VMC and 0.19 ≲δ ≲0.22 for variational AFQMC at U/t = 8; the differing phase boundaries motivate more accurate models and methods.Source: §5.1 · p. 12

5.1.1. Competition at intermediate U/t.

• At intermediate coupling U/t = 4, DMRG and CP-AFQMC find weaker stripe tendencies, while DiagMC predicts uniform d-wave order only for δ ≥0.3.Smaller dopings were inaccessible to DiagMC because of convergence problems, leaving the crossover between uniform and stripe states unresolved.Source: §5.1.1 · pp. 12–13

5.1.2. Stripes at finite t′.

• With negative t′, stripe periods shorten as |t′|/t increases; at δ = 1/8, period-4 stripes occur for 0.16(4) < |t′|/t < 0.423(10), whereas sufficiently positive t′ favors uniform d-wave order.iPEPS and VMC find period-4 stripes without coexisting d-wave order at 1/8 doping, but superconducting coexistence predictions differ for other stripe periods and dopings; width-4 cylinder pairing may be plaquette rather than ordinary d-wave.Source: §5.1.2 · pp. 13–14
• Three-band cuprate models also show competition between uniform d-wave states and stripes, with stripe order depending on charge-transfer energy and superconducting coexistence still unconfirmed by newer methods.CP-AFQMC found charge and spin order at Δ = 4.4 eV, but weaker charge order and subtle competition between spin orders at Δ = 2.5 eV.Source: §5.1.2 · p. 14

5.2. Pseudogap

• In cuprates, below T*∼300K and below optimal doping, the pseudogap suppresses antinodal but not nodal density of states; experimental correlation lengths are a few lattice spacings.It was originally identified through NMR Knight-shift suppression (2); signatures occur in varied probes (111), particularly angle-resolved photoemission spectroscopy (131, 229).Source: §5.2 · p. 14

5.1.3. Competition in 3-band models of the cuprates.

• DCA and CDMFT capture short-range correlations and spectral functions; calculations show these correlations suffice for temperature-driven antinodal suppression, corroborated by 2×2 CDMFT across solvers and parameters.For t′/t < 0, Fermi-surface curvature becomes more hole-like toward the insulating state, whereas t′/t > 0 drives the system toward nesting (37).Source: §5.1.3 · pp. 14–15

5.2.2. Organizing principles: momentum-sector-selective metal-insulator transition, Widom line and Fermi-surface topology.

• Doping toward half-filling produces a sector-selective transition: an antinodal gap opens while the node remains metallic, followed by a larger, isotropic Mott gap (86).At higher doping, the system is nearly an isotropic Fermi liquid; approaching half-filling first differentiates nodal and antinodal quasiparticle lifetimes. On the electron-doped side, the sector-selective regime is smaller or absent for sufficiently large t′, and momentum differentiation is weaker.Source: §5.2.2 · p. 15
• In 2×2 CDMFT, multiple thermodynamic and dynamic criteria for T* collapse onto a Widom line of anomalies emanating from the finite-doping Mott critical endpoint (236).Source: §5.2.2 · p. 15
• At strong coupling, van Hove effects on the pseudogap remain unresolved: 8-site DCA finds no substantial effect (86), whereas DCA and DQMC report substantial effects (288, 289).Source: §5.2.2 · p. 15

5.2.3. Insights into the origin of the pseudogap via fluctuation diagnostic approaches.

• DCA fluctuation diagnostics finds a peaked Q=(π, π), zero-frequency contribution only in the spin channel, implicating long-lived short-range magnetic fluctuations in pseudogap opening.Charge and particle-particle contributions are distributed uniformly over the momenta considered, unlike the spin contribution.Source: §5.2.3 · p. 16

5.3. Superconductivity

• Weak coupling supports superconductivity, but near U/t∼8 and moderate doping stripes preempt it; superconductivity remains close in energy and modest model changes may favor it.Its relationship to pseudogap, charge-ordered, and metallic states remains an active research topic.Source: §5.3 · p. 16
• The superconducting phase boundary can be located by tracking susceptibility divergence while cooling from the disordered phase or order-parameter disappearance while heating from low temperature.Source: §5.3 · p. 16
• Finite-cluster QMC favored dx2−y2 pairing at moderate doping (283), while DCA enabled extrapolated large-cluster critical-temperature determination (170) and exploration of favorable regimes (30).Source: §5.3 · p. 16
• Four-site (158, 167, 93, 36, 236) and eight-site (85, 80) cluster DMFT studies ordered superconductivity, its energetics, spectra, response functions, and phase boundary.Finite-size effects remain, but generic energetic and response-function trends coincide with other methods and cuprate observations.Source: §5.3 · pp. 16–17
• AF fluctuations act as pairing glue (180, 216), including in DCA studies (171, 168, 169, 172, 82), while particle-particle scattering is a main Tc suppressor in DΓA analysis (139).The DΓA study analyzed how dynamical vertex structure influences the superconducting transition temperature.Source: §5.3 · p. 17

5.4. Bad metal

• At high temperature and intermediate doping, the model is metallic, but its transport properties are inconsistent with Fermi-liquid behavior.Resistivity is believed to be linear in temperature over a wide range, and transport in this regime has been studied using linked cluster expansion (195), dynamical mean field theory (204, 48, 195, 26), and lattice Monte Carlo (109).Source: §5.4 · p. 17
• Matsubara-axis resistivity results are uncertain because analytical continuation is unreliable at high temperature and Matsubara-point spacing grows linearly with temperature.Real-frequency DMFT calculations using NRG (48) or ED (26) impurity solvers neglect vertex corrections (25).Source: §5.4 · p. 17
• The very-high-temperature regime is well understood, but no reliable calculations currently exist for intermediate temperatures; real-frequency DiagMC methods (245) promise a resolution.Source: §5.4 · p. 17

6.1. Angle-resolved photoemission spectroscopy

• Single- and two-particle excitations remain less understood at intermediate-to-strong coupling; analytic continuation makes spectral-function values uncertain, while imaginary-time and Matsubara methods provide alternatives (122, 61; 255; 79, 71).Continuation noise is amplified, especially at high temperature, away from zero frequency, and for bosonic quantities; gap existence and quasiparticle peaks are typically more robust than their sizes or heights.Source: §6.1 · pp. 17–18

6.2. Optical conductivities

• Cluster methods find a large undoped insulating gap, doping-induced Drude and mid-IR features, and a metallic Drude peak (161); c-axis response lacks vertex corrections, unlike in-plane response (65, 161; 160, 161).Vertex corrections also disappear in single-site (121) and four-site CDMFT (92) by symmetry; conductivity results span cluster DMFT at low temperature and several methods at high temperature.Source: §6.2 · p. 18

6.3. Raman spectroscopy

• Simulated 2D-model Raman results from cluster DMFT, with (162) and without (161, 214) vertex corrections and in the superconducting state (81), show an insulating two-magnon peak, pseudogap signatures upon doping, and temperature- and doping-evolution generally consistent with experiment.Source: §6.3 · p. 18

6.4. Magnetic response

• Quantum cluster methods reproduce pseudogap suppression of the NMR Knight shift and distinguish oxygen and copper relaxation rates (31).Full momentum- and energy-dependent magnetic susceptibility remains an open simulation problem because lattice and quantum cluster methods lack sufficient momentum resolution; approximate dual-fermion results can be compared with neutron spectroscopy (153, 156; 38).Source: §6.4 · pp. 18–19

7.1. The 2D Hubbard model on other lattices

• At half filling, the honeycomb model has a direct Dirac-semimetal-to-AF-Mott-insulator transition at Uc/t ≈3.8, with no intermediate spin liquid; at quarter doping, chiral d+id superconductivity was found.The transition is in the Gross-Neveu-Yukawa universality class (9, 189); quarter-doped chiral d+id order was found by several groups (184, 269, 78).Source: §7.1 · p. 19

Triangular lattice.

• On the triangular lattice, a chiral spin liquid between metallic and 120-degree ordered phases was found by DMRG (244), but VMC (252) found none except with sufficiently large t′.The phase diagram remains under debate; at intermediate coupling, increased chiral correlations coexist with stripy antiferromagnetic correlations (286), and ladder dual fermion work advanced frequency-dependent susceptibility calculations (156).Source: Triangular lattice · p. 19

Kagome lattice.

• At half filling, the Kagome model has no magnetic order in the strongly interacting limit, with growing consensus for a quantum spin liquid; cluster DMFT found a Mott-transition critical coupling Uc/t ∼8.22 (186).DMRG found a Wigner crystal in the lightly doped t-J model (123); larger doping remains under study, and DQMC, DMFT, and DΓA were used to investigate magnetic correlations across the Mott transition (135).Source: Kagome lattice · pp. 19–20

7.2. The attractive model

• Away from half filling, the attractive Hubbard model retains onsite pairing and has a finite-temperature KT transition whose temperature vanishes on approaching half filling (219, 181).At half filling, onsite s-wave pairing and (π, π) charge-density-wave order coexist (99); superconductivity can be enhanced near a logarithmic Van Hove singularity (100), and Rashba spin-orbit coupling does not introduce a sign problem (230).Source: §7.2 · p. 20

7.3. Multi-orbital models

• Multi-orbital Hubbard generalizations model material-specific systems, including cuprate-inspired Emery and bilayer models and oxide perovskites with Slater-Kanamori interactions U, U′, and J.Single-site multi-orbital DMFT has yielded Hunds’ metals (72), while combining band downfolding with DMFT (DFT+DMFT) is frequently used to describe realistic materials (144).Source: §7.3 · p. 20

7.4. Models with non-local interactions

• The extended Hubbard model adds repulsive nearest-neighbor density interaction V to local U, driving charge order and an interplay among charge order, antiferromagnetism, metallicity, and Mott insulating behavior.DCA treats these interactions explicitly (247, 248, 126, 192), whereas the dual boson approach treats them through a diagrammatic expansion around DMFT (210, 243).Source: §7.4 · pp. 20–21

7.5. Non-equilibrium and driven system

• Reliable, generic numerical tools for nonequilibrium Hubbard setups are mostly absent, despite their promise and the range of driven and transient questions they raise.These include quenches, external-field driving, voltage-induced currents, thermalization, and the destruction or reestablishment of order; single-site and small-cluster DMFT studies and cold-atom Floquet studies have been performed (6, 177, 215).Source: §7.5 · p. 21

Conclusions

• Stripe states are widely accepted as ground states across an extended doping, interaction, and t′ range, while strong short-range antiferromagnetic fluctuations cause the strong-coupling pseudogap.These are among the aspects of the model on which computational studies have established a consensus.Source: §8 · p. 21
• Accurately locating phase boundaries, especially d-wave superconductivity, and resolving phase competition remain major open challenges.These include studying low-doping phase separation with more methods; the interplay of superconductivity and stripe order in partially filled stripes; and finite-temperature competition between stripe and uniform phases. More accurate beyond-single-site DMFT transport calculations, controlled access to real-time evolution and real-frequency observables, and reliable methods for multiorbital or longer-range-interaction models are also needed to connect single-band calculations to multiband Hubbard models and electronic-structure Hamiltonians without crude downfolding.Source: §8 · pp. 21–22

Perspectives

• Tensor-network, dynamical variational Monte Carlo, and two-particle-vertex diagrammatic approaches are advancing finite-temperature, dynamical, and real-frequency studies.Two-dimensional tensor networks have been extended to finite temperature and used to study stripe correlations at δ = 1/16 (285); excitation spectra and results in other parameter regimes are anticipated. Dynamical variational Monte Carlo has been applied to the 2D Hubbard model (62, 28), and two-particle-vertex diagrammatic approaches are being extended to real frequencies (147).Source: §8 · p. 22
• Finite-temperature methods have advanced through large-cluster dynamical mean field studies extrapolatable to the thermodynamic limit, more accurate DMFT extensions, and more powerful continuous-time and DiagMC methods.The passage characterizes the latter methods as orders of magnitude more powerful than previous techniques.Source: §8 · p. 22
• Developing generally applicable numerical techniques that overcome existing methods’ shortcomings remains the most urgent need for addressing open issues.Most advances over the last 15 years can be traced directly to algorithmic and numerical progress.Source: §8 · p. 22

• The authors are unaware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the review’s objectivity.Source: DISCLOSURE STATEMENT · p. 22