Introduction

Metal-insulator transitions (MITs) driven by the electron-phonon (EP) coupling have been the focus of solid-state physics studies for decades. The Peierls instability1,2 is perhaps the most prominent and fundamental example. Acting in the static (frozen-phonon) limit, it establishes an insulating—charge-density-wave (CDW)—broken-symmetry state, related to a structural distortion, as observed in most one-dimensional (1D) inorganic and organic conductors3,4. Quantum phonon fluctuations, on the other hand, become increasingly important in low dimensions, and counteract any development of long-range order5,6,7.

While the basic mechanisms promoting or hindering a metal-insulator quantum-phase transition are well known8, their detailed understanding within the framework of minimal theoretical models is challenging. Actually there are only a very limited number of microscopic model Hamiltonians for which such an MIT could be rigorously proven. Holstein-like models9 provide a paradigm in this respect, both in the spinless and the spinful case10,11,12 By using involved numerical techniques, such as exact diagonalization and kernel polynomial13,14,15,16, diagrammatic and quantum Monte Carlo17,18,19, or density-matrix renormalization group (DMRG)5,20 methods, this type of model could be solved numerically exact in the last years. This concerns the ground-state, spectral, transport and thermodynamic properties in 1D, both for a few particles and for the case of the half-filled band. As a result, in the latter case, the ground-state phase diagram of, e.g., the spinless fermion Holstein model was determined, and it has been shown that an MIT from a (repulsive) Tomonaga-Luttinger-liquid (TLL)21,22 to a CDW occurs when the EP coupling is increased at finite phonon frequency23,24,25,26,27,28.

The original Holstein model considers a spatially localized EP coupling to a dispersionless phonon mode9. While extensions regarding the range and kind of the EP coupling and the electron hopping have been discussed for some time29,30,31, the influence of the missing phonon dispersion was recently questioned, but largely only for the few particle (polaron and bipolaron) problem32,33,34,35,36,37. Interestingly, it turned out that the phonon dispersion had a profound effect on the transport properties of Holstein and Edwards polarons in 1D. In two dimensions, it has been demonstrated for the spinful Holstein model that the competition between pairing and charge order can be tuned by even a weak bare phonon dispersion7. Against this background, it seems necessary to re-examine the influence of phonon dispersion on the TLL-CDW MIT of the spinless fermion Holstein model as well. That is the main purpose of this work. To this end, we extend the model accordingly and determine the ground state and transport properties using unbiased numerical techniques for the 1D infinite half-filled Holstein system. The MIT phase boundary is compared with that of an effective electronic Hamiltonian, derived for weak phonon dispersion in Appendix A. Results for the related Edwards fermion-boson transport model 38,39,40,41,42 are presented and discussed in Appendix B.

Model and methods

The modified Holstein model under consideration is

$$\begin{aligned} \hat{H}&= -t_f \sum _j ( \hat{f}_j^{\dagger } \hat{f}_{j+1} + \text {H.c.}) - g \omega _0 \sum _j \hat{n}_j (\hat{b}_j^{\dagger } + \hat{b}_j ) + \omega _0 \sum _j \hat{b}_j^{\dagger } \hat{b}_j + t_\omega \sum _j ( \hat{b}_j^{\dagger } \hat{b}_{j+1} + \text {H.c.}) , \end{aligned}$$
(1)

where \(\hat{f}_j\) are fermion annihilation operators describing the conduction electrons, \(\hat{b}_j\) are boson annihilation operators for the phonons, and \(\hat{n}_j = \hat{f}_j^{\dagger } \hat{f}_j\). Compared with the regular Holstein model, there is an additional nearest-nearest neighbor hopping of the phonons that changes their dispersion to \(\omega (k) = \omega _0 + 2 t_\omega \cos (k)\) (note the different sign convention compared with the electron hopping). To simplify the notation, we will use \(t_f = 1\) as the unit of energy throughout this work. We will moreover restrict ourselves to the half-filled band case with the number of electrons equal to half the number of sites.

All our numerical simulations are based on the matrix-product-state (MPS) formalism43. For ground-state calculations with finite system sizes, we employ the regular DMRG algorithm44. Is is often more efficient, however, to work directly in the thermodynamic limit by using infinite matrix-product states (iMPS). For those simulations we apply the variational uniform MPS algorithm (VUMPS)45 to obtain a ground-state approximation, and the time-dependent variational principle (TDVP)46 with bond expansion47,48 to carry out time evolutions. Since we consider perturbations that only affect the state in a finite region, the latter amounts to the TDVP method for finite systems with infinite boundary conditions49,50,51. The maximum bond dimension in our simulations was 400, except for systems with periodic boundary conditions, in which case it was 1600.

The infinite-dimensional local Hilbert spaces of the phonons must be truncated to a finite dimension \(D_{b}\) for the numerical simulations. Depending on the model parameters, one may need a relatively large \(D_b\) to accurately represent the low-energy states, particularly in the CDW phase. Several methods have been developed to make MPS simulations more efficient for systems with large local Hilbert spaces52. Here, we use the pseudo-site approach for the finite-system DMRG calculations53. For the VUMPS and TDVP simulations in the TLL phase, we did not utilize any such techniques because we found a relatively small boson dimension \(D_b=8\) to be sufficient.

Results

Phase diagram

To gain insight into the general ground-state phase diagram of the model (1), it is helpful to apply a Lang-Firsov transformation that changes the fermion and boson operators while leaving the fermion density operators \(\hat{n}_j\) invariant54. As shown in Appendix A, this transformation reveals an effective electron-electron interaction that is not present in the regular Holstein model with dispersionless phonons. Since the interaction is exponentially decaying with an exponent \(\text {arcosh}(\omega _0 / |2 t_\omega |)\), it suffices to consider only the nearest-neighbor part if the dispersion is small. Applying standard perturbation theory10,55 then yields the following purely electronic Hamiltonian that is valid for \(t_f/g^2, |t_\omega | \ll \omega _0\):

$$\begin{aligned} \hat{H}_{\text {eff}}&= -\tilde{t}_f \sum _j \left( \hat{f}_j^{\dagger } \hat{f}_{j+1}+ \text {H.c.} \right) + \Big [\frac{\tilde{t}_f^2}{\omega _0} \big (4 g_1 + 2 g_2 \big ) + 2 t_\omega g^2 \Big ] \sum _j \hat{n}_j \hat{n}_{j+1} + \frac{\tilde{t}_f^2 }{\omega _0} g_1 \sum _j \big [ \hat{f}_{j-1}^{\dagger } (2 \hat{n}_j -1) \hat{f}_{j+1} + \text {H.c.} \big ] \,, \end{aligned}$$
(2)

where \(\tilde{t}_f=t_f e^{-a^2}\), \(g_1 = \sum _{n=1}^\infty \frac{a^{2n}}{n!n}\), \(g_2 = \sum _{n=1}^\infty \sum _{m=1}^\infty \frac{a^{2(n+m)}}{n!m!(n+m)}\), and \(a = g(1 - t_\omega / \omega _0)\). The main difference compared with the effective model for the regular Holstein model is the additional contribution \(\propto t_\omega\) to the interaction term. From Eq. (2) one can draw some qualitative conclusions about the phase diagram. For \(t_\omega = 0\), the interaction is always repulsive, and since it decreases slower with g than the effective hopping strength, it causes a transition from a TLL to a CDW as g is increased10. When \(t_\omega\) is finite, however, there is an additional term that grows quadratically with g and therefore becomes decisive at strong EP coupling. In particular, it leads to an attractive interaction if \(t_\omega < 0\), which indicates the possibility of an attractive TLL, and even phase separation, both of which do not occur in the regular Holstein model. For \(t_\omega > 0\), on the other hand, the phonon dispersion should simply move the CDW-TLL transition to lower EP couplings. Since it is caused by a strong nearest-neighbor attraction, the phase-separated state considered here consists of two subsystems that are empty and occupied, respectively, with small particle-number fluctations at the interface. For open boundary conditions, the electrons will accumulate on either the left or the right side of the system, so that the number of bonds where both sites are empty or both sites are occupied is maximized.

Figure 1
figure 1

Ground-state phase diagram for different values of the phonon hopping \(t_\omega\). Circles mark the numerically calculated phase boundaries, with the solid lines as a guide to the eye. For comparison, the dashed line in (a) shows the TLL-CDW boundary in the dispersionless Holstein model28. The dash-dotted lines in (b) indicate the transitions according to the effective model (2) when neglecting the correlated next-nearest-neighbor hopping.

Full size image

Since the effective Hamiltonian (2) has a limited region of validity, it is important to check the above predictions with numerical calculations. Let us first discuss the resulting phase diagram shown in Fig. 1. For \(t_\omega = 0.1\), there is indeed a significant shift of the TLL-CDW transition to lower values of g for all \(\omega _0\). The stronger effect in the adiabatic regime is likely related to the fact that we consider a fixed \(t_\omega\), so that the effective interaction becomes longer ranged at small \(\omega _0\). Going to a small negative \(t_\omega = -0.0025\), a region with phase separation appears at large g. It is mostly located above the CDW phase, an exception being the anti-adiabatic regime, where the CDW phase disappears and there is a direct transition from a TLL to phase separation. The phase-transition points at strong-coupling can also be estimated by dropping the last term in Eq. (2), which leaves only a t-V Hamiltonian whose phase diagram is known exactly56. As demonstrated in Fig. 1b, the phase-separation boundary obtained this way agrees quite well with that of the full model. Interestingly, the competition between the effective repulsive interaction from perturbation theory and the attractive interaction due to the phonon dispersion causes a reentrant CDW-TLL transition that is visible at \(\omega _0 = 5\). When \(\omega _0\) is lowered, this second TLL region becomes smaller and effectively disappears. Already at moderate negative hopping \(t_\omega = -0.1\), we no longer observe a CDW phase, i.e., the system appears to always be in an attractive TLL or phase-separated state. We would like to point out that, in contrast, the half-filled Edwards model with dispersive bosons neither forms an attractive TLL nor a phase-separated state for negative values of \(t_\omega\), see Appendix B.

Within the TLL phase, the system is characterized by the TLL parameter K, which can be used to distinguish between effective repulsive \((K < 1)\) and attractive interactions \((K > 1)\). Furthermore, \(K=1/2\) signals a transition to the CDW phase, while \(K \rightarrow \infty\) indicates the onset of phase separation. There are several ways to numerically determine the TLL parameter K. Here, we use its relation to the linear conductance57,58 and the charge structure factor59,60,61. Figure 2 displays the results for an intermediate phonon energy \(\omega _0=1\). At weak EP coupling, the TLL parameter is close to the value \(K=1\) of a non-interacting chain of fermions regardless of the phonon parameters. As the EP coupling g increases, however, the effect of the phonon dispersion on K becomes significant. For \(t_\omega =0.1\) and \(t_\omega =-0.0025\), K decreases with g up to the TLL-CDW transition where \(K=1/2\), indicating that the TLL for these parameters is repulsive as in the regular Holstein model. In contrast, setting \(t_\omega = -0.1\) clearly leads to an attractive TLL with \(K > 1\), similar to a sufficiently strong longer-ranged EP coupling62. The rapid increase of K near \(g \approx 1.5\) also suggests that a phase-separation transition is approached. However, while calculating K allows to accurately locate the TLL-CDW transition, it is difficult to prove the occurrence of phase separation in this manner. We therefore also compute the inverse compressibility,

$$\begin{aligned} \kappa ^{-1}(L)&= \frac{L}{2}\left[ E_0(L/2+1) + E_0(L/2-1) - 2 E_0(L/2)\right] , \end{aligned}$$
(3)

where \(E_0(N)\) is the ground-state energy for a system with N electrons, and L is the number of sites63,64. For \(\omega _0 = 1\) and \(t_\omega = -0.1\), \(\kappa ^{-1} = \lim _{L \rightarrow \infty }\kappa ^{-1}(L)\) vanishes around \(g \approx 1.74\), which indicates that the system indeed becomes unstable towards phase separation at this point. Furthermore, we found no evidence for electron pairing, i.e., defining \(\kappa ^{-1}(L)\) in terms of the two-particle gap instead of the single-particle gap as in Eq. (3) does not lead to a smaller \(\kappa ^{-1}\) for the parameters considered here. As an alternative to the inverse compressibility, one can also look at the ground-state energy per site \(\epsilon _0\) from an iMPS simulation and extrapolate where it becomes larger than that of a phase-separated state with an empty and a fully occupied subsystem, i.e., where \(\Delta \epsilon = - \frac{1}{2} g^2 \omega _0^2 / (\omega _0 + 2 t_\omega ) - \epsilon _{0} = 0\). The transition point obtained this way agrees reasonably well with that from the inverse compressibility, which shows that phase separation indeed occurs between an empty and an occupied phase. We therefore use this simpler method to determine the remaining phase-separation boundaries.

Figure 2
figure 2

(a) TLL parameter K as a function of the EP coupling g for \(\omega _0=1\) and different phonon dispersions. Dotted lines mark the locations of phase transitions. (b) Inverse compressibility \(\kappa ^{-1}\) and ground-state energy per site \(\Delta \varepsilon\) relative to a phase-separated state for \(\omega _0=1\) and \(t_\omega =-0.1\). The inset shows the extrapolation from finite-system DMRG calculations of \(\kappa ^{-1}(L)\) with open boundary conditions. (c) Charge gap in the CDW phase for \(\omega _0=1\) and \(t_\omega = -0.1\). The dashed line shows the charge gap when only the nearest-neighbor-repulsion term is considered in Eq. (2). The calculations were done for finite systems with open boundary conditions and then extrapolated to the thermodynamic limit using a second-order polynomial. For parameters near the TLL-CDW transition the finite-size extrapolation is presented in the inset.

Full size image

The TLL-CDW transition can also be determined by locating the breakdown of the CDW state, which is signaled by a closing of the charge gap \(\Delta _c = E_0(L/2+1) + E_0(L/2-1) - 2 E_0(L/2)\). Figure 2c shows \(\Delta _c\) in the CDW phase near the phase transitions for model parameters \(\omega _0 = 1\) and \(t_\omega =-0.0025\). For EP couplings g close to the lower phase boundary, the results are clearly consistent with those for the TLL parameter K, although accurately determining the transition point this way is difficult because of the exponential closing of the gap at the Kosterlitz-Thouless transition. Near the second transition, on the other hand, the electron motion is almost completely frozen because of the exponential decrease of the effective hopping amplitude \(\tilde{t}_f\) with the EP coupling g, so that the charge gap in this region is approximately proportional to the coefficient of the nearest-neighbor interaction in Eq. (2). The closing of the gap coincides with \(\Delta \varepsilon\) becoming zero, which suggests that only a vanishingly small intermediate TLL phase exists.

Figure 3
figure 3

(a,b) Examples for the calculation of K by means of the charge conductance for \(\omega _0=1\). The inset of shows an alternative method based on the charge structure factor in finite systems with periodic boundary conditions. In the main panels, the extrapolated values of K are displayed as dashed lines.

Full size image

Calculation of the Tomonaga-Luttinger parameter K

In the TLL phase the linear conductance at zero temperature is given by \(K / (2 \pi )\)57. Accordingly, one can determine the TLL parameter from the charge-current response to a small voltage V as

$$\begin{aligned} K&= \lim _{t \rightarrow \infty } \lim _{V \rightarrow 0} \frac{2 \pi J(t)}{V} . \end{aligned}$$
(4)

The charge current at time t for an arbitrary bond j is \(J_j(t) = -it_f \langle \hat{f}_j^\dagger \hat{f}_{j+1} - \hat{f}_{j+1}^\dagger \hat{f}_j \rangle _t\). Specifically, we apply a constant potential gradient over \(L_V = 10\) sites and calculate the average current J(t) in that region. The corresponding perturbation to the Hamiltonian for times \(t>0\) is \(\sum _j V_j \hat{n}_j\), with \(V_j = \min (\max (-\tfrac{V}{2} + \tfrac{jV}{L_V+1},-\tfrac{V}{2}),\tfrac{V}{2})\). Figure 3 shows the obtained current J(t) for parameters near the phase transitions for \(\omega _0 = 1\) and \(t_\omega = \pm 0.1\). Except for small fluctuations, J(t) saturates with time, which allows to estimate the conductance and thereby K. Only for \(t_\omega = 0.1\) and \(g = 1.2\), the current J(t) appears to decrease at long times, likely because the system is already slightly in the insulating CDW regime.

Another way to determine the TLL parameter K is the relation60

$$\begin{aligned} K&= \lim _{q \rightarrow 0 } \frac{2 \pi S(q)}{q} \,, \end{aligned}$$
(5)

which connects K to the static structure factor \(S(q)= \frac{1}{L}\sum _{j \ell } e^{i q(j - \ell )} \langle \hat{n}_j \hat{n}_\ell \rangle\). Although S(q) can be calculated from the iMPS approximation of the ground state, we found that using finite systems with periodic boundary conditions and \(K = \lim _{L \rightarrow \infty } L S(2 \pi /L)\) leads to a weaker momentum dependence and thus a clearer extrapolation. As demonstrated in Fig. 3, the TLL parameters obtained this way are consistent with the charge current at long times. We also used the approach based on the structure factor for the 3 points in Fig. 1b around \(\omega _0=5\), where we encountered convergence issues in the iMPS method.

Figure 4
figure 4

(a) TLL parameter in the Edwards model with dispersive bosons and (b) the corresponding ground-state phase diagram.

Full size image

Conclusions

We used MPS techniques to numerically investigate the effect of a finite phonon dispersion on the phase diagram of the one-dimensional Holstein model, focusing in particular on the TLL parameter K, which we have extracted from both static and dynamic quantities. In agreement with the derived effective strong-coupling Hamiltonian, our results demonstrate that a downward dispersion increases the tendency to CDW order, while an upward dispersion leads to an effective electron-electron attraction. As the EP coupling is increased, the attractive interaction manifests itself first as a TLL with \(K > 1\), and then as phase separation between an empty and a full subsystem. This is markedly different from the phase diagram of the Holstein model without phonon dispersion that consists only of a repulsive TLL and a CDW, but resembles the situation in models with longer-ranged electron-phonon coupling62. As shown in Appendix B, the corresponding Edwards fermion-boson model also shows no attractive TLL or phase separation.

While we studied the half-filled case in this work, it would also be interesting to investigate the phase diagram at other densities. For example, a downward phonon dispersion and strong EP coupling should lead to phase separation between a CDW and an empty region, since the effective interaction in that case favours a CDW with wave number \(\pi\). Lastly, it might be worthwhile to consider extensions of the model with other dispersion or more general types of EP coupling. A natural question in this regard is, whether the different effective electron-electron interactions can lead to more exotic phases, such as TLLs with electron pairing, which were observed in purely fermionic chains with specific longer-ranged interactions or pair-hopping terms65.