Hubbard-Extended Holstein Model Overview

• The Hubbard-Extended Holstein Model is a class of lattice Hamiltonians that combine local Coulomb repulsion with Holstein-type electron-phonon coupling to produce diverse magnetic, charge, and superconducting phases.
• Model variants extend interactions to include nonlocal density-density and phonon couplings, enabling studies of bipolaron formation, metallic regimes, and time-dependent nonequilibrium dynamics.
• Advanced computational methods such as quantum Monte Carlo, DMFT, and Lang-Firsov transformations are used to explore phase competition, thermalization, and lattice dynamics in both equilibrium and driven settings.

The Hubbard-Extended Holstein Model denotes a family of lattice Hamiltonians in which Holstein-type electron-phonon coupling is combined with Hubbard electron-electron repulsion. In the narrowest usage this is the Hubbard-Holstein model, while in broader usage it also includes nearest-neighbor density interactions, nonlocal electron-phonon form factors, multiorbital generalizations, doping away from half filling, and explicitly time-dependent nonequilibrium variants. Across these formulations, the defining problem is the competition between local Coulomb repulsion, phonon-mediated attraction, retardation, and lattice geometry, which yields antiferromagnetic, charge-ordered, bond-ordered, metallic, superconducting, bipolaronic, and nonequilibrium thermalizing regimes (Picano et al., 12 Apr 2026, Miyao, 2016, Weber, 2024).

1. Canonical Hamiltonians and model variants

In its standard static form, the model combines nearest-neighbor hopping, an onsite Hubbard repulsion, and a local Holstein coupling between the electron density and a local lattice displacement. A representative Hamiltonian is

$${\cal H} = \sum_{\langle ij \rangle \sigma} t_{ij}\left(c_{i\sigma}^{\dagger} c_{j\sigma} + H.c.\right) + U \sum_i n_{i\uparrow} n_{i\downarrow} + \sum_i \left( \frac{1}{2} K Q_i^2 - g Q_i n_i \right) - \mu \sum_{i\sigma} n_{i\sigma},$$

with $t_{ij}=-t$, local distortion $Q_i$, lattice stiffness $K$, and chemical potential $\mu$ (Kurdestany et al., 2017). In the square-lattice formulation used for variational Monte Carlo studies, the phonon coordinate is equivalently written as

$x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$

so that the Hamiltonian contains a harmonic phonon sector with mass $M$ and frequency $\Omega$ together with the local density-displacement term $g\sum_i x_i n_i$ (Ohgoe et al., 2017).

The “extended” qualifier is used in two distinct senses in the literature. One extends the electron-electron sector by adding nonlocal density-density interactions. A rigorous half-filled formulation on a finite cubic lattice is

$$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$

where $t_{ij}=-t$0 is a nearest-neighbor repulsion and $t_{ij}=-t$1 the Einstein-phonon frequency (Miyao, 2016). The other extends the electron-phonon sector itself, as in the one-dimensional extended Holstein-extended Hubbard model,

$t_{ij}=-t$2

where the electron on site $t_{ij}=-t$3 couples to phonons on two sites, $t_{ij}=-t$4 and $t_{ij}=-t$5 (Tam et al., 2013).

Nonequilibrium work commonly adopts a time-dependent Hubbard-Holstein Hamiltonian with local Einstein phonons and explicit quenches. In the weak-coupling nonequilibrium DMFT study, the Hamiltonian contains a time-dependent Holstein coupling $t_{ij}=-t$6, and at half filling the choice $t_{ij}=-t$7, $t_{ij}=-t$8, and $t_{ij}=-t$9 removes the Hartree term (Picano et al., 12 Apr 2026). This half-filled, local, dispersionless-phonon setting is often used as the minimal “Hubbard-extended Holstein” reference point.

2. Effective interactions, retardation, and transformed descriptions

A central organizing quantity is the phonon-renormalized onsite interaction. After a Lang-Firsov transformation in the extended Holstein-Hubbard model, the effective interaction becomes

$Q_i$0

so the phonons contribute an attractive term $Q_i$1 that can reduce or overcome the bare onsite repulsion (Miyao, 2016). In more general bipartite-lattice formulations,

$Q_i$2

which makes explicit that the phonon-mediated attraction is a matrix correction to the Coulomb interaction (Miyao, 2016).

In path-integral treatments, integrating out the phonons yields a retarded density-density interaction rather than a purely static one. For the one-dimensional half-filled Hubbard-Holstein model,

$Q_i$3

and only in the antiadiabatic limit does one recover the instantaneous interaction

$Q_i$4

This is the basis for the statement that the model is “reminiscent of the extended Hubbard model,” while the mechanism is distinct: the competition is not between onsite and nearest-neighbor instantaneous repulsion, but between a static repulsive part and a retarded attractive screening in the same local channel (Weber, 2024).

The same structure appears in the square-lattice Holstein-Hubbard model treated by variational Monte Carlo, where integrating out phonons gives

$Q_i$5

and in the antiadiabatic limit

$Q_i$6

This language is useful for locating the regime $Q_i$7, but the finite-frequency problem remains genuinely dynamical (Ohgoe et al., 2017). Nonequilibrium studies make this limitation sharper: after an electron-phonon coupling quench, the generalized Lang-Firsov construction yields a time-dependent effective interaction $Q_i$8, so a static screened-$Q_i$9 picture is insufficient once coherent phonon dynamics are excited (Werner et al., 2013).

Strong-coupling expansions recast the problem into effective low-energy Hamiltonians whose form depends on the sign of

$K$0

For $K$1, the low-energy sector contains no doubly occupied sites and maps to a generalized $K$2-$K$3-$K$4 model; for $K$5, the low-energy sector is one of onsite bipolarons and maps to a hard-core boson model (Han et al., 2020). Retardation, encoded by the phonon frequency relative to interaction scales, controls how strongly hopping is suppressed by Franck-Condon factors and therefore how readily the system enters cluster, bipolaronic, or charge-ordered regimes.

3. Half-filled phase competition

At half filling, the model is a canonical setting for competition among spin-density-wave, charge-density-wave, bond-ordered, metallic, and pairing tendencies. In one dimension, exact directed-loop quantum Monte Carlo for retarded interactions finds four regimes at $K$6: SDW, Luther-Emery liquid, BOW, and CDW. For $K$7, the SDW-to-BOW boundary is consistent with a Berezinskii-Kosterlitz-Thouless transition at $K$8, while the BOW-to-CDW transition is continuous and second-order at $K$9, with finite-size extrapolation giving $\mu$0. At stronger coupling, the BOW phase shrinks; at $\mu$1, the intermediate BOW phase is gone and the SDW-CDW transition is strongly first-order (Weber, 2024).

Two-dimensional studies reveal a different ordering balance. Variational Monte Carlo on the square-lattice Holstein-Hubbard model finds an AF insulating phase at large $\mu$2 and small $\mu$3, a CO insulating phase at large $\mu$4 and small $\mu$5, and an “extended intermediate metallic or weakly superconducting phase” near $\mu$6. For the example $\mu$7, $\mu$8, the AF-to-intermediate and intermediate-to-CO boundaries are marked by energy crossings near $\mu$9 and $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$0, respectively, consistent with first-order-like transitions (Ohgoe et al., 2017). In a complementary finite-temperature CT-INT study of the half-filled two-dimensional model, the CDW transition is in the universality class of the two-dimensional Ising model. At $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$1, $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$2, and $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$3, the critical temperature is $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$4, and increasing $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$5 to $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$6 reduces $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$7 by about $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$8. Above $x_i = \sqrt{\frac{1}{2M\Omega}}(b_i+b_i^\dagger),$9, the data are consistent with a spin-gapped bipolaronic metal (Weber et al., 2017).

Phonon dynamics sharpen the distinction between Peierls, metallic, and Mott regimes. In the one-dimensional half-filled Holstein-Hubbard chain, CT-INT calculations of the phonon spectral function find behavior consistent with a soft-mode Peierls transition in the adiabatic regime, a central peak related to long-range order in the Peierls phase, substantial phonon renormalization in the metallic regime, and only weakly modified dispersion in the Mott phase (Weber et al., 2015). The result ties the lattice sector directly to charge fluctuations through exact relations between the phonon propagator and the dynamic charge structure factor.

In three and higher dimensions, rigorous results confirm that the half-filled extended Holstein-Hubbard model supports genuine long-range charge order when the electron-phonon side is sufficiently strong. Under the condition $M$0 and $M$1, there is staggered long-range charge order at sufficiently low temperature; in the opposite regime $M$2, earlier results imply no long-range charge order and an antiferromagnetic ground state (Miyao, 2016). This provides a mathematically controlled counterpart to the phase-competition picture found numerically in lower-dimensional formulations.

4. Doping, nonlocal couplings, multiorbital structure, and bipolaron physics

Away from half filling, the balance between Mottness and lattice-driven localization changes qualitatively. In a Hartree-Fock treatment of the doped Hubbard-Holstein model on the cubic lattice, the ground-state energy $M$3 is not everywhere convex, so the uniform state is unstable to coexistence between an insulating AF phase at $M$4 and a carrier-rich metallic phase at $M$5. The metallic volume fraction is $M$6, and conduction begins when $M$7 exceeds the percolation threshold $M$8, giving $M$9. Finite electron-lattice coupling $\Omega$0 does not change the ordering of homogeneous AF, F, P, or spiral phases at fixed density, but it increases $\Omega$1 and hence $\Omega$2, stabilizing the insulating side (Kurdestany et al., 2017).

Nonlocal electron-phonon coupling can qualitatively alter the half-filled low-energy instability. In the one-dimensional extended Holstein-extended Hubbard model, the momentum dependence of the phonon-mediated interaction contains the factor $\Omega$3. At half filling, $\Omega$4 implies $\Omega$5, so the nonlinear $\Omega$6 electron-phonon term vanishes, spin-charge coupling is weakened, Umklapp is suppressed, and dominant singlet superconducting fluctuations appear in a substantial weak-repulsion region of the $\Omega$7-$\Omega$8 plane (Tam et al., 2013). This is a sharp contrast with the purely local Holstein limit, where half-filled density-wave physics is typically stronger.

Extended electron-phonon range also changes bipolaron binding. In one dimension, the local Holstein-Hubbard model shows no bipolaron-bipolaron attraction, whereas the extended-Holstein-Hubbard $\Omega$9 model exhibits clear bipolaron-bipolaron attraction above $g\sum_i x_i n_i$0 at $g\sum_i x_i n_i$1. The threshold for a single spin-parallel bipolaron is $g\sum_i x_i n_i$2, and above this value the resulting composite survives even as $g\sum_i x_i n_i$3 (Chakraborty et al., 2013). In three-dimensional BCC and FCC lattices, continuous-time path-integral QMC finds that purely local Holstein coupling favors only onsite $g\sum_i x_i n_i$4 bipolarons, while extended coupling stabilizes intersite $g\sum_i x_i n_i$5 bipolarons. In the extended model, there is a region of light pairing in both lattices, and on the FCC lattice at large phonon frequency and large $g\sum_i x_i n_i$6, intersite bipolarons become superlight because they move by first-order hopping (Adebanjo et al., 23 Jul 2025).

Multiorbital generalizations add another layer of competition. Determinant quantum Monte Carlo for a one-dimensional three-orbital Hubbard-Holstein model at filling $g\sum_i x_i n_i$7 finds a metallic phase at weak couplings, an orbital-selective Mott phase at larger $g\sum_i x_i n_i$8 and small $g\sum_i x_i n_i$9, a multicomponent CDW insulator at larger $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$0 and small $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$1, and, when both couplings are large but comparable, an orbitally correlated insulating phase with strong short-range orbital correlations but without a large charge susceptibility. In the intermediate-coupling regime, the metallic region sits roughly where $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$2 (Li et al., 2018). This indicates that the Hubbard-Extended Holstein framework remains nontrivial even when the local charge channel is embedded in multiorbital Hund-coupled physics.

5. Nonequilibrium dynamics and dynamical diagnostics

Nonequilibrium formulations treat the model as a laboratory for thermalization, screening dynamics, and coupled electron-phonon memory effects. In the weak-coupling half-filled Hubbard-Holstein model, a sudden quench $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$3 from an initial equilibrium state at $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$4 and $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$5 produces a crossover between electron-dominated and phonon-dominated relaxation. For the parameters studied, the crossover occurs around $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$6. Step-by-Step DMFT reveals a sharp thermalization front in the plane of real time and DMFT iteration number; electronic observables show a clear front already for weak quenches, while the local dispersionless phonons exhibit a delayed front that becomes visible only near and beyond the crossover. Whenever both fronts are resolved, they propagate with the same velocity, indicating coherent thermalization of the coupled electron-phonon system. Moderate Hubbard repulsion $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$7 changes rates quantitatively but not the overall mechanism (Picano et al., 12 Apr 2026).

In the strong-coupling nonequilibrium regime, interaction quenches and electron-phonon coupling quenches produce different but related effects. Nonequilibrium DMFT with a generalized Lang-Firsov transformation and NCA/OCA impurity solvers shows that a rapid $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$8-pulse creates doublon-hole excitations whose decay is phonon-assisted and strongly enhanced when the Mott gap is commensurate with an integer multiple of $$H_A=\sum_{(x;y)\,\circ=1} -t\,(c_x^*c_y+c_y^*c_x) +U\sum_{x\in A}(n_x-1)^2 +V\sum_{(x;y)}(n_x-1)(n_y-1) +g\sum_{x\in A}(n_x-1)(b_x+b_x^*) +\omega\sum_{x\in A} b_x^*b_x,$$9. A quench in $t_{ij}=-t$00 generates a time-dependent effective interaction

$t_{ij}=-t$01

persistent oscillations of the phonons, and phonon-enhanced doublon production (Werner et al., 2013). This establishes that nonequilibrium phonons act not only as a dissipative channel but also as a coherent drive.

Several diagnostics expose the local lattice sector more directly than conventional spectra. NRG combined with DMFT yields the full probability distribution $t_{ij}=-t$02 of the local phonon displacement and an effective local potential $t_{ij}=-t$03. In the impurity problem and the infinite-dimensional lattice model, increasing electron-phonon coupling broadens $t_{ij}=-t$04 and can produce a double-peak structure or a double-well effective potential, while finite $t_{ij}=-t$05 delays these features; in the lattice problem the normal, antiferromagnetic, and charge-ordered phases display distinct displacement distributions (Hewson et al., 2010). Exact nonequilibrium spectral moment sum rules for the retarded Green’s function and self-energy further constrain any time-dependent solution. For the fully time-dependent Holstein-Hubbard Hamiltonian, zeroth through third moments were derived in real and momentum space, providing exact benchmarks for time-resolved many-body calculations and for the interpretation of ultrafast spectra (Najafi et al., 2021).

The phonon spectral function provides an additional dynamical window. In the one-dimensional half-filled model, the exact relation

$t_{ij}=-t$06

shows that the phonon propagator is directly tied to charge fluctuations. This explains the small-$t_{ij}=-t$07 anomaly as a hybridization of charge and phonon excitations and clarifies why the central peak in the Peierls phase is suppressed when Hubbard repulsion drives the system toward a metallic or Mott regime (Weber et al., 2015).

6. Rigorous results, methods, and realizations

The Hubbard-Extended Holstein family has an unusually strong rigorous foundation. For the extended Holstein-Hubbard model on a connected bipartite lattice at half filling, if $t_{ij}=-t$08 is positive definite, then for each spin sector $t_{ij}=-t$09 the ground state of $t_{ij}=-t$10 is unique. The same work proves an upper bound on the charge susceptibility,

$t_{ij}=-t$11

for momenta with $t_{ij}=-t$12, implying absence of charge long-range order when $t_{ij}=-t$13 is uniformly positive (Miyao, 2014). A related theorem for the Holstein-Hubbard model on connected bipartite lattices states that if $t_{ij}=-t$14 is positive definite, the ground state has Lieb-type ferrimagnetic total spin $t_{ij}=-t$15, is unique up to $t_{ij}=-t$16-fold spin degeneracy, exhibits antiferromagnetic long-range order when the sublattice imbalance is extensive, and shows no long-range charge order when $t_{ij}=-t$17 (Miyao, 2016). In the single-hole, $t_{ij}=-t$18 limit, Nagaoka ferromagnetism remains stable after adding Holstein coupling, so the projected ground state still has maximal total spin $t_{ij}=-t$19 under the connectivity condition (Miyao, 2016).

The model has therefore been attacked with a broad methodological arsenal rather than a single dominant technique. Weak-coupling nonequilibrium studies use conserving nonequilibrium DMFT with self-consistent Migdal and self-consistent second-order perturbation theory (Picano et al., 12 Apr 2026). Strong-coupling nonequilibrium work uses generalized Lang-Firsov transformations together with NCA and OCA impurity solvers (Werner et al., 2013). Equilibrium lattice problems have been studied with CT-INT QMC (Weber et al., 2015, Weber et al., 2017), exact directed-loop QMC for retarded interactions (Weber, 2024), many-variable variational Monte Carlo (Ohgoe et al., 2017), determinant QMC (Li et al., 2018), NRG+DMFT (Hewson et al., 2010), and continuous-time path-integral QMC for bipolaron formation (Adebanjo et al., 23 Jul 2025). This methodological diversity reflects the fact that the model interpolates between weak-coupling Fermi-surface physics, intermediate-coupling competing orders, and strong-coupling polaronic or bipolaronic sectors.

The same breadth appears in proposed realizations. A cold-atom quantum simulator based on two dressed Rydberg species in a monolayer with species-selective painted optical potentials was proposed as an analogue of the extended Hubbard-Holstein problem. The construction is designed to realize hopping, onsite Hubbard interaction, phonon modes, and tunable phonon-mediated onsite and offsite interactions. It was argued that both boson-mediated preformed pairing and BKT transition temperatures are experimentally accessible, with $t_{ij}=-t$20 and $t_{ij}=-t$21 in representative regimes (Hague et al., 2020).

Taken together, these results define the Hubbard-Extended Holstein Model not as a single Hamiltonian but as a controlled hierarchy of local and extended electron-phonon Hubbard systems. Its unifying theme is that a local repulsion $t_{ij}=-t$22, a Holstein coupling $t_{ij}=-t$23, and, when present, nonlocal interactions or nonequilibrium driving do not merely compete quantitatively; they reorganize the accessible low-energy manifold, the nature of order parameters, and the route to thermalization.