Hubbard–Holstein Model Overview

• Hubbard–Holstein Model is a minimal lattice model coupling electron correlations with local phonon interactions to simulate Mott, Peierls, polaronic, and charge-order phenomena.
• It employs both quantum-phonon and static adiabatic limits, allowing detailed analysis through methods such as QMC, DMRG, and DMFT.
• Key insights include the interplay between electron-electron repulsion and electron-phonon coupling, leading to rich phase diagrams and emergent metallic behavior.

The Hubbard–Holstein model is a minimal lattice model for correlated electrons coupled to local lattice degrees of freedom. In its standard single-band form it combines nearest- or longer-range hopping, an on-site Hubbard repulsion $U$, and a local Holstein coupling between the electronic density and Einstein phonons or static lattice displacements, thereby unifying Mott, Peierls, polaronic, and charge-ordering tendencies within one framework (Karakuzu et al., 2022). Across the literature, it appears both in quantum-phonon form and in static adiabatic limits, and it has been used to study half-filled antiferromagnetism and charge order, doped metal–insulator transitions, stripe correlations, multiorbital orbital-selective regimes, and nonequilibrium thermalization (Kurdestany et al., 2017).

1. Canonical Hamiltonians and effective couplings

A standard quantum-phonon parameterization is

$$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$

with $c_{i\sigma}^\dagger$ and $c_{i\sigma}$ fermionic creation and annihilation operators, $n_i=n_{i\uparrow}+n_{i\downarrow}$, $a_i^\dagger$ and $a_i$ local phonon operators, and $\Omega$ the bare phonon frequency (Karakuzu et al., 2022). In cuprate-motivated two-dimensional studies this Hamiltonian is often implemented on a square lattice with $t$, $t'$, and sometimes $$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$0, while in other contexts it is specialized to nearest-neighbor hopping, single-band square or cubic lattices, or multiorbital chains (Mendl et al., 2017).

A common static adiabatic form replaces the phonon operators by classical distortions $$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$1,

$$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$2

which is the form used in the three-dimensional doped study of the Mott transition (Kurdestany et al., 2017). In this limit minimization with respect to $$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$3 yields $$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$4, so the lattice contributes an energy gain $$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$5 for homogeneous density.

Dimensionless electron–phonon couplings are convention-dependent. For quantum Einstein phonons a standard choice is

$$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$6

with $$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$7 the bare bandwidth; for a nearest-neighbor square lattice, $$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$8 (Karakuzu et al., 2022). In the static Holstein limit a common alternative is

$$H = -\sum_{ij,\sigma} t_{ij}\left(c^{\dagger}_{i\sigma} c_{j\sigma} + \text{H.c.}\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + g \sum_i (n_i-\bar n)(a_i+a_i^\dagger) + \Omega \sum_i \left(a_i^\dagger a_i + \tfrac12\right) - \mu \sum_i n_i,$$9

again expressed relative to the electronic bandwidth (Kurdestany et al., 2017). In the antiadiabatic limit, integrating out phonons gives a local instantaneous attraction and hence

$c_{i\sigma}^\dagger$0

or, in some bandwidth-based conventions, $c_{i\sigma}^\dagger$1 (Karakuzu et al., 2022). The extended Holstein–Hubbard model supplements these local terms by a nearest-neighbor density interaction $c_{i\sigma}^\dagger$2, which is central to rigorous charge-order results (Miyao, 2016).

2. Half-filling: antiferromagnetism, charge order, and metallic windows

At half filling, the Hubbard–Holstein model is organized by the competition between a repulsion-driven antiferromagnetic or Mott tendency and a phonon-driven charge-density-wave or Peierls tendency. Determinant QMC on the two-dimensional half-filled model finds a strong competition between AFM and CDW order and shows that a static effective-$c_{i\sigma}^\dagger$3 picture requires significant refinement: retardation effects slow the onset of charge order, so CDW order remains absent even when the effective $c_{i\sigma}^\dagger$4 is negative, opening a window where neither AFM nor CDW order is well established and where there are signatures of a possible metallic phase (Nowadnick et al., 2012). A static auxiliary-field study at half filling on a two-dimensional square lattice likewise finds AFI and COI/BPI regimes, and reports that when both interactions are weak, mutual competition leads to a metallic phase in an otherwise insulator-dominated phase diagram, while spatial correlations induced by thermal fluctuations lead to pseudogap features at intermediate coupling (Pradhan et al., 2015).

The strong-coupling and low-dimensional limits sharpen this competition in different ways. In the one-dimensional half-filled model, the phonon spectral function is consistent with a soft-mode Peierls transition in the adiabatic regime, together with a central peak related to long-range order in the Peierls phase; tuning the system from a Peierls to a metallic phase with nonzero $c_{i\sigma}^\dagger$5 suppresses the central peak, whereas the dispersion is only weakly modified in the Mott phase (Weber et al., 2015). In DMET studies of the one-dimensional model, the adiabatic regime displays a direct Peierls-insulator to Mott-insulator transition, while the antiadiabatic regime exhibits a rather large intervening metallic phase; the Born–Oppenheimer approximation is qualitatively similar only in the adiabatic regime and fails entirely in the antiadiabatic regime, where it predicts a sharp direct transition that the fully quantum treatment replaces by a broad metallic region (Reinhard et al., 2018).

Rigorous results delimit part of this landscape. For the half-filled Holstein–Hubbard model on bipartite lattices, positive-definite effective interactions $c_{i\sigma}^\dagger$6 imply ferrimagnetism whenever the electron–phonon interaction is not so strong; antiferromagnetic long-range order is then present and charge long-range order is absent (Miyao, 2016). For the extended Holstein–Hubbard model in $c_{i\sigma}^\dagger$7, the opposite inequality $c_{i\sigma}^\dagger$8, with $c_{i\sigma}^\dagger$9, yields rigorous staggered long-range charge order at sufficiently low temperature, thereby justifying the phase competition between antiferromagnetic and charge orders (Miyao, 2016).

3. Doping, phase separation, and density-driven metal–insulator transitions

Away from half filling, the Hubbard–Holstein model supports inhomogeneous and percolative metal–insulator scenarios that do not arise naturally in the simplest homogeneous Hubbard description. In a three-dimensional single-band model on a simple cubic lattice with static Holstein phonons, Hartree–Fock analysis finds that $c_{i\sigma}$0 is not convex near half filling, so the homogeneous state is unstable to phase separation into an undoped AF insulating region and a hole-rich metallic region with concentration $c_{i\sigma}$1 (Kurdestany et al., 2017). Standard Maxwell construction then fixes $c_{i\sigma}$2, and if the metallic and insulating volume fractions are $c_{i\sigma}$3 and $c_{i\sigma}$4, one obtains

$c_{i\sigma}$5

With a representative percolation threshold $c_{i\sigma}$6, global conduction emerges only when $c_{i\sigma}$7, giving a critical concentration

$c_{i\sigma}$8

In that formulation, the Mott transition is density-driven and percolative rather than a uniform collapse of the gap, and increasing Holstein coupling increases $c_{i\sigma}$9 and therefore raises $n_i=n_{i\uparrow}+n_{i\downarrow}$0 while leaving the magnetic phase boundaries of homogeneous phases unchanged (Kurdestany et al., 2017).

Determinant QMC on the hole-doped two-dimensional model yields a complementary asymmetry between doped Mott and doped Peierls regimes. In the Mott-insulating regime at large $n_i=n_{i\uparrow}+n_{i\downarrow}$1 and small $n_i=n_{i\uparrow}+n_{i\downarrow}$2, doping rapidly suppresses the insulator by transferring spectral weight from the upper Hubbard band, shifting the lower Hubbard band toward the Fermi level, and generating emergent quasiparticles at the Fermi level (Mendl et al., 2017). By contrast, in the Peierls-insulating regime at large $n_i=n_{i\uparrow}+n_{i\downarrow}$3 and small $n_i=n_{i\uparrow}+n_{i\downarrow}$4, the CDW gap remains robust out to relatively high doping. The same study reports that the $n_i=n_{i\uparrow}+n_{i\downarrow}$5-wave superconducting susceptibility increases with lowering temperature in a regime of intermediate $n_i=n_{i\uparrow}+n_{i\downarrow}$6 and $n_i=n_{i\uparrow}+n_{i\downarrow}$7, which suggests that the doped model can support a metallic state in which strong e–e and e–ph interactions cooperate rather than simply cancel (Mendl et al., 2017).

4. One-dimensional strong-coupling structures and quantum-phonon effects

The one-dimensional model is a particularly sharp laboratory for the role of phonon quantum fluctuations. In the phonon spectral function, adiabatic Peierls order appears through the softening of the phonon mode at $n_i=n_{i\uparrow}+n_{i\downarrow}$8, the emergence of a central peak in the ordered phase, and hybridization between charge and phonon excitations at small momenta; a nonzero Hubbard interaction can suppress the central peak while leaving a sizable renormalization of the phonon dispersion, whereas the Mott phase leaves the phonon dispersion only weakly modified (Weber et al., 2015). This identifies the phonon sector itself as a diagnostic of whether the system is in a Peierls, metallic, or Mott regime.

A different strong-coupling limit, relevant for $n_i=n_{i\uparrow}+n_{i\downarrow}$9 and strong e–e and e–ph couplings, produces a correlated nearest-neighbor singlet phase. Starting from an effective Hamiltonian for the one-dimensional Hubbard–Holstein model, the singlet sector maps onto a hard-core-boson $a_i^\dagger$0–$a_i^\dagger$1 model at a transformed filling, and in that description superfluidity and CDW occur mutually exclusively, with diagonal long-range order appearing only at one-third filling (Reja et al., 2011). The same analysis shows that the BEC occupation number $a_i^\dagger$2 scales as $a_i^\dagger$3, as in the HCB tight-binding problem, but with a smaller prefactor in the interacting singlet phase (Reja et al., 2011).

DMET and DMRG refine this picture by showing that quantum phonons can qualitatively change the phase topology. In the adiabatic regime, the Born–Oppenheimer approximation gives qualitatively similar results to the fully quantum treatment, but in the antiadiabatic regime it fails by predicting a sharp direct Mott–Peierls transition rather than the large intervening metallic phase found by DMET and benchmarked against DMRG (Reinhard et al., 2018). A plausible implication is that metallicity in one dimension is stabilized not merely by the cancellation of $a_i^\dagger$4 against a phonon-mediated attraction, but by phonon quantum fluctuations themselves.

5. Stripe order, multiorbital variants, and orbital selectivity

In the doped two-dimensional single-band model, local electron–phonon coupling strongly reshapes inhomogeneous ordering tendencies. Zero-temperature VMC together with finite-temperature DQMC shows that the lattice couples more strongly to the charge component of stripes than to the spin component, so the e–ph interaction can either enhance or suppress stripe correlations depending on parameters such as $a_i^\dagger$5 and $a_i^\dagger$6 (Karakuzu et al., 2022). In that setting, the Hubbard–Holstein model is used as a minimal description of a CuO$a_i^\dagger$7-like plane in which stripe order and $a_i^\dagger$8-wave superconductivity compete or coexist. The central result is parameter dependent: intermediate e–ph coupling can amplify charge modulations associated with stripes, while stronger coupling can instead favor more uniform or differently ordered states and thereby alter the balance between stripe and superconducting correlations (Karakuzu et al., 2022).

Multiorbital generalizations introduce an additional hierarchy of bandwidths, Hund couplings, and orbital fillings. In a one-dimensional three-orbital Hubbard–Holstein model with average filling $a_i^\dagger$9, determinant QMC finds competition between an orbital-selective Mott phase and a multicomponent CDW insulating phase, with an intermediate metallic phase at weak e–e and e–ph couplings (Li et al., 2018). For larger and comparable couplings, that metallic phase develops short-range orbital correlations and becomes insulating, producing an orbital-correlated phase in addition to the metallic, CDW, and OSMP regimes. The OSMP is identified by the pinning of one orbital near integer occupancy and the suppression of its low-energy spectral weight, while the CDW phase shows uniform orbital occupations near $a_i$0, enhanced double occupancies, and charge susceptibility peaks at $a_i$1 together with secondary features near $a_i$2 at lower temperature (Li et al., 2018). This suggests that in multiorbital settings the Holstein channel can compete not only with Mott localization but also with orbital selectivity itself.

6. Computational frameworks, density-functional formulations, and nonequilibrium dynamics

Much of the modern understanding of the Hubbard–Holstein model is tied to nonperturbative numerical methods and to reformulations that treat electronic and bosonic sectors on equal footing. Determinant QMC has been developed in detail for the two-dimensional single-band model, with explicit treatment of the lattice degrees of freedom. A notable methodological result is that, although Holstein coupling introduces a sign problem even at half filling, there exists a region of large Holstein coupling where the fermion sign recovers despite large values of the Hubbard interaction, indicating that studies of correlated polarons at finite carrier concentrations are likely accessible to DQMC simulations (Johnston et al., 2013).

Density-functional reformulations recast the model in terms of site densities and phonon coordinates. For the inhomogeneous Hubbard–Holstein model with local e–e and e–ph interactions, a lattice DFT has been constructed in which exchange-correlation potentials are derived analytically for an isolated site and via DMFT for the infinite-dimensional Bethe lattice (Boström et al., 2019). In that formulation the electronic XC potential displays discontinuities as a function of density, with their location and character controlled by the screened interaction $a_i$3, while the phonon XC potential vanishes for the strictly local harmonic coupling considered there. Benchmarks against exact calculations show that this DFT gives a good description of linear conductance and real-time dynamics (Boström et al., 2019).

Nonequilibrium DMFT extends the same model family into real time. After a sudden switch-on of the electron–phonon interaction in the weak-coupling Hubbard–Holstein model, self-consistent Migdal plus second-order Hubbard perturbation theory yields a crossover between electron-dominated and phonon-dominated relaxation regimes (Picano et al., 12 Apr 2026). In the Step-by-Step DMFT construction, thermalization appears in the plane of real time and DMFT iteration number as a sharp propagating front. Electronic observables show such a front already for weak quenches within the simulated time window, whereas the phononic sector exhibits a visible front only at sufficiently strong coupling; whenever both fronts are resolved, they propagate with the same velocity, indicating coherent thermalization of the coupled electron–phonon system (Picano et al., 12 Apr 2026).

Across these formulations, several open directions recur. Numerical studies emphasize momentum-dependent e–ph couplings, long-range Coulomb interactions, and more realistic phonon spectra as natural extensions beyond local Einstein phonons, while the rigorous literature identifies away-from-half-filling behavior, precise phase boundaries, and superconductivity as outstanding problems (Karakuzu et al., 2022). The Hubbard–Holstein model remains minimal in the technical sense used throughout this literature: simple enough to expose the interplay of local repulsion and lattice coupling, but rich enough to support Mott, Peierls, stripe, percolative, orbital-selective, and nonequilibrium phenomena within a unified framework (Miyao, 2016).