Two-Dimensional Holstein Model:
Critical Temperature, Ising Universality, and Bipolaron Liquid
Fundamental open questions about the two-dimensional Holstein model of electrons coupled to quantum phonons are addressed by continuous-time quantum Monte Carlo simulations. The critical temperature of the charge-density-wave transition is determined by a finite-size scaling of a renormalization-group-invariant correlation ratio. is finite for any nonzero coupling for classical phonons, and suppressed by quantum lattice fluctuations. The phase transition—also detectable via the fidelity susceptibility and machine learning—is demonstrated to be in the universality class of the two-dimensional Ising model. We discuss the possibility of at weak coupling and present evidence for a spin-gapped, bipolaronic metal above .
Introduction.—Phase transitions in two-dimensional (2D) fermionic systems are a central topic of theoretical and experimental condensed matter physics. Correlated quasi-2D materials with rich phase diagrams include high-temperature superconductors  and transition-metal dichalcogenides . Dirac fermions in two dimensions can be investigated in graphene . Strongly correlated 2D fermions exhibit exotic phases  and phase transitions , and support long-range order at . While magnetic order originates from short-range Coulomb repulsion, the mechanism behind the numerous charge-density-wave (CDW) phases found experimentally is electron-phonon coupling. In addition to polaron effects, the latter leads to a phonon-mediated, retarded electron-electron interaction and an intricate interplay of spin, charge, and lattice fluctuations.
Quantum Monte Carlo (QMC) simulations are a key tool to understand the physics of correlated 2D quantum systems. Although simulations are significantly harder for fermions than for spins or bosons, QMC methods have been very successfully applied to fermionic models. However, whereas the phase diagram and critical behavior of, e.g., the 2D honeycomb Hubbard model is known in detail [7; 8; 9; 10], the same is not true even for the simplest electron-phonon models. Simulations with phonons are often severely restricted by long autocorrelation times also away from critical points . These challenges are the main reason why reliable results for critical temperatures and a convincing analysis of the critical behavior of the most widely studied Holstein molecular-crystal model are still missing. In fact, even the simpler 1D case had until recently been discussed controversially .
In this Letter, we use large-scale continuous-time QMC simulations to study the CDW transition in the 2D Holstein model. In contrast to previous work, we use finite-size scaling to determine . We demonstrate that the transition can be detected by the fidelity susceptibility and deep neural networks and provide evidence for Ising critical behavior as well as a metallic bipolaron phase. Finally, we point out the possibility of minimal CDW order (i.e., ) resulting from quantum fluctuations.
Model.—The Holstein Hamiltonian  reads
The first two terms describe free electrons and free phonons, respectively. Here, creates an electron with spin at lattice site and electrons hop with amplitude between nearest-neighbor sites on a square lattice. The phonons are of the Einstein type with frequency ; their displacements couple to local fluctuations of the electron occupation . We simulated lattices with periodic boundary conditions at half-filling (, chemical potential ). A useful dimensionless coupling parameter is with the free bandwidth . We set , , and the lattice constant to one and use as the energy unit.
The relative simplicity of Eq. (1) has led to numerous QMC investigations of CDW formation and superconductivity [14; 15; 16; 17; 18; 19; 20; 21; 22]. The model (1) also corresponds to the limit of the Holstein-Hubbard Hamiltonian with competing Mott and CDW ground states [23; 24; 25; 26; 27]. For Eq. (1), mean-field theory (exact for and ) predicts a CDW ground state with a checkerboard pattern for the lattice displacements and the charge density [ordering vector , see inset of Fig. 1] at half-filling [14; 15; 17]. Here, we explore the impact of quantum and thermal fluctuations.
Methods.—The determinant QMC (DetQMC) method for 2D fermions  used in almost all previous works is limited by autocorrelations . Instead, we employed the exact continuous-time interaction expansion (CT-INT) method  to simulate the fermionic action obtained from Eq. (1) by integrating out the phonons in the coherent-state path integral . CT-INT simulations amount to a stochastic but exact summation of a weak-coupling Dyson expansion in the retarded interaction around ; for reviews see Refs. [31; 32]. The simulation time scales as , where is the average expansion order and . DetQMC formally has a better scaling, but CT-INT benefits from small expansion orders at weak coupling, noncritical autocorrelation times 111We used 5000 single-vertex updates and 8 Ising spin flips per sweep for all results shown., and the absence of discretization errors. It outperforms DetQMC for the parameters considered despite being limited for by a sign problem. Whereas even the noninteracting case is challenging for DetQMC, CT-INT trivially gives exact results for and can in principle simulate the entire range of phonon frequencies, including the experimentally important adiabatic regime . Finally, we also simulated the classical limit using the method of Ref.  in combination with parallel tempering.
(with ) from the charge structure factor
either at fixed or at fixed . Being a renormalization-group invariant, becomes independent of at the critical point for sufficiently large . This gives a crossing of curves for different independent of critical exponents, as illustrated in Fig. 2(a) for and . For smaller , the crossing points drift and critical values can be obtained by extrapolating to [Fig. 2(a), inset]. The finite-size corrections are well described by the Ising form [and similar for ] with . In Fig. 2(a), the extrapolation gives . CT-INT is particularly well suited to study weak-coupling problems. In particular, for can be tracked down to very small values because the larger is compensated by a smaller so that the expansion order remains manageable. Finally, similar to previous detQMC results [14; 15; 17], we only find short-range superconducting correlations that are further suppressed at .
Given the current interest in unbiased diagnostics for phase transitions, we illustrate in Figs. 2(b) and 2(c) that the CDW transition can be detected using the fidelity susceptibility  or deep learning . The estimator for  gives rather large errors at low temperatures, but in Fig. 2(b) shows a peak for consistent with from Fig. 1. A simple way to use deep learning is to train a neural network at and with labeled QMC data for and then estimate the probability 222We used Keras to implement a network similar to Ref. : input neurons followed by a convolutional layer with 8 features, receptive fields, and sigmoid activation. After max-pooling, a second (identical) convolutional layer with ReLU activation and again max-pooling. Next, a dense layer with 256 ReLU neurons and dropout regularization. Finally, a single sigmoid neuron to produce . labeled QMC samples of the density correlator measured on a single time slice were used for supervised learning (20 epochs, batch size 20) at and , respectively. was averaged over 2000 samples to obtain Fig. 2(c).. The critical point then manifests itself as the crossover from 0 to 1 of visible in Fig. 2(c), with at the estimated .
Figure 1 shows as a function of for different , covering in particular the entire adiabatic range . For classical phonons (), is zero for but nonzero for . For , simulations are more restricted regarding the accessible temperatures. Quantum lattice fluctuations suppress at a given . For , deviates from the classical value only at very weak coupling in the regime . For larger , quantum fluctuation effects are visible over the entire parameter range shown. The systematic suppression of with increasing is perfectly consistent with the fact that for the attractive Hubbard model , to which the Holstein model maps in the limit . This connection and a possible metallic phase at low temperatures as a result of quantum fluctuations will be discussed below. A metallic phase at is naturally expected in the 2D Holstein-Hubbard model because the antiferromagnetic Mott state is confined to .
Critical behavior.—In the thermodynamic limit, the long-range CDW order at spontaneously breaks the sublattice symmetry. The two possible CDW patterns (cf. Fig. 1) imply the same critical behavior as the 2D Ising model. Accordingly, the charge structure factor (3) should obey the finite-size scaling relation 
with and 333At weak coupling, scaling with the 2D quantum Ising exponents is observed on similar system sizes, see arXiv:1709.01096v1. However, classical 2D Ising criticality is expected for any at sufficiently large .. This is confirmed in Fig. 3 for and . The best scaling collapse  over in Fig. 3(b) gives .
Phase diagram.—While Fig. 1 gives accurate critical values for the CDW transition at , we are unable to directly simulate the ground state. To gain further insight into the effect of quantum fluctuations on the low-temperature physics and the phase diagram we exploit the exact mapping to the 2D attractive Hubbard model with in the antiadiabatic limit , corresponding to maximally strong quantum fluctuations. At , the latter has coexisting CDW and superconducting order for any but minimal CDW order in the sense that .
Figure 1 suggests that the classical Holstein model () has a CDW ground state for any . Exact numerical results show that the 1D Holstein model with quantum phonons remains metallic at for , despite a nesting-related divergence of the noninteracting charge susceptibility . Its limit, the 1D attractive Hubbard model, has a metallic but spin-gapped Luther-Emery liquid  ground state and no long-range order. For a half-filled square lattice and nearest-neighbor hopping, due to the combined effect of nesting and van Hove singularities [47; 15]. For the Hubbard model, this stronger divergence gives rise to a weak coupling instability and long-range order for any . For the Holstein model, it suggests a CDW ground state for any . Finally, an extended metallic phase is found within dynamical mean-field theory [48; 49; 50] where van Hove singularities are absent.
To explore the impact of quantum fluctuations for the 2D Holstein model, we calculated the susceptibility
for and phonon frequencies ranging from to . To avoid spurious finite-size effects, we set [51; 15]. The results in Fig. 4(a) reveal a dependence for the value close to the classical limit, but a significant suppression upon increasing fluctuations via larger values of . As expected, the results for the Holstein model approach those for the attractive Hubbard model, but for any remains larger than for . Figure 4(b) shows a finite-size scaling of the charge structure factor (3) at fixed and , corresponding to for all (cf. Fig. 1). At small distances, CDW correlations are stronger for the Holstein model than for the attractive Hubbard model. However, even the data for are essentially identical to those for at large and vanish faster than , consistent with .
Before making the connection with our results, we note that there are two distinct scenarios for the shape of the phase boundary defined by , as illustrated in Fig. 5. As a function of , the Holstein model interpolates between (classical case, CDW order for ) and (attractive Hubbard model, CDW order at only). In scenario (I), for any if , whereas for any if . In contrast, in scenario (II), for and for . Case (II) can further be divided into (IIa) where CDW order exists at for any , and (IIb) with a disordered phase at below . In scenario (I), the adiabatic (classical) fixed point determines the behavior for any finite . By contrast, in scenario (IIa), the physics is determined by the antiadiabatic fixed point for and by the adiabatic fixed point for . Finally, scenario (IIb) is realized in the 1D Holstein model (where for all and ) .
At a given temperature, the phase boundary in Fig. 1 undergoes an increasingly strong shift to larger with increasing . For , scenario (I) would hence imply an onset at weak coupling with a very small exponent. Support for scenario (II) also comes from the quasi-degeneracy of results for and in Fig. 4(b). The fact that the susceptibility and CDW correlations for are at least as strong as for (Fig. 4) favor scenario (iia) with CDW order at . On the other hand, recent approximate variational QMC studies report an extended metallic  or superconducting  phase at . The exact ground-state phase diagram hence remains an important open issue.
Bipolaron liquid.—A final interesting point is the nature of the metallic phase at . In the CDW phase, spin, charge, and hence also single-particle excitations are gapped. For 1D electron-phonon models, the spin gap persists in the metallic phase  and the CDW transition occurs at the two-particle level via the ordering of preformed pairs (singlet bipolarons) and the opening of a charge gap. The same is true for the attractive Hubbard model for which the spin gap can in fact be made arbitrarily large by increasing while keeping . Hence, the disordered phase at low but finite temperatures is not a Fermi liquid but a metal with a gap for single-particle and spin excitations [53; 54]. Such a phase is the 2D analog of a Luther-Emery liquid . Singlet bipolarons in principle also form for any in the 2D Holstein model, although their binding energy () can be small . Nevertheless, we expect a spin-gapped metallic phase to exist for suitable parameters at half-filling. A potential metallic phase  with a spin gap at would realize a Bose metal  in a fermionic system. Finally, at sufficiently high temperatures, bipolarons undergo thermal dissociation .
To detect signatures of a spin-gapped metal, we consider the static charge and spin susceptibilities
with , . Figure 6(a) shows results for and . Whereas diverges with decreasing temperature in a Fermi liquid (open symbols), it is strongly suppressed as by the spin gap. The charge susceptibility is also suppressed at very low , but is expected to remain finite at since the attractive Hubbard model has a vanishing charge gap. The distinct temperature scales reflected by the maxima of and reveal the spin-gapped metallic phase at . For the Holstein model, is cut off by the spin gap, whereas is cut off by the charge gap that opens at the CDW transition at . The distinct maxima visible even in the adiabatic regime [Figs. 6(b) and 6(c)] are consistent with a spin-gapped phase at . The extent of the latter appears to decrease with decreasing and the phase is expected to be absent in the classical or mean-field limit () where charge and spin gaps are identical. An immediate and important corollary of the existence of a spin-gapped metal of bipolarons above would be that, contrary to expectations in previous work [20; 21; 26], a single-particle and/or spin gap does not imply a CDW phase.
Conclusions.—Using exact QMC simulations, we determined critical temperatures and confirmed the Ising universality of the CDW transition of the 2D Holstein model. Our results suggest long-range CDW order at for any nonzero coupling, potentially with due to quantum fluctuations, and a spin-gapped metallic phase of bipolarons above .
We thank F. Assaad, P. Bröcker, and T. Lang for helpful discussions and acknowledge support by the DFG through SFB 1170 ToCoTronics and FOR 1807. We gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA  at the Jülich Supercomputing Centre.
- Lee et al.  P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
- Manzeli et al.  S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nat. Rev. Materials 2, 17033 (2017).
- Neto et al.  A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- Zhou et al.  Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
- Senthil et al.  T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Science 303, 1490 (2004).
- Mermin and Wagner  N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
- Sorella et al.  S. Sorella, Y. Otsuka, and S. Yunoki, Sci. Rep. 2, 992 (2012).
- Assaad and Herbut  F. F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013).
- Parisen Toldin et al.  F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Phys. Rev. B 91, 165108 (2015).
- Otsuka et al.  Y. Otsuka, S. Yunoki, and S. Sorella, Phys. Rev. X 6, 011029 (2016).
- Hohenadler and Lang  M. Hohenadler and T. C. Lang, in Computational Many-Particle Physics (Springer Berlin Heidelberg, 2008).
- Hohenadler and Fehske  M. Hohenadler and H. Fehske, arXiv:1706.00470 (2017).
- Holstein  T. Holstein, Annal. Phys. 8, 325 (1959).
- Scalettar et al.  R. T. Scalettar, N. E. Bickers, and D. J. Scalapino, Phys. Rev. B 40, 197 (1989).
- Marsiglio  F. Marsiglio, Phys. Rev. B 42, 2416 (1990).
- Levine and Su  G. Levine and W. P. Su, Phys. Rev. B 42, 4143 (1990).
- Noack et al.  R. M. Noack, D. J. Scalapino, and R. T. Scalettar, Phys. Rev. Lett. 66, 778 (1991).
- Levine and Su  G. Levine and W. P. Su, Phys. Rev. B 43, 10413 (1991).
- Vekić et al.  M. Vekić, R. M. Noack, and S. R. White, Phys. Rev. B 46, 271 (1992).
- Vekić and White  M. Vekić and S. R. White, Phys. Rev. B 48, 7643 (1993).
- Niyaz et al.  P. Niyaz, J. E. Gubernatis, R. T. Scalettar, and C. Y. Fong, Phys. Rev. B 48, 16011 (1993).
- Zheng and Zhu  H. Zheng and S. Y. Zhu, Phys. Rev. B 55, 3803 (1997).
- Berger et al.  E. Berger, P. Valášek, and W. von der Linden, Phys. Rev. B 52, 4806 (1995).
- Nowadnick et al.  E. A. Nowadnick, S. Johnston, B. Moritz, R. T. Scalettar, and T. P. Devereaux, Phys. Rev. Lett. 109, 246404 (2012).
- Johnston et al.  S. Johnston, E. A. Nowadnick, Y. F. Kung, B. Moritz, R. T. Scalettar, and T. P. Devereaux, Phys. Rev. B 87, 235133 (2013).
- Nowadnick et al.  E. A. Nowadnick, S. Johnston, B. Moritz, and T. P. Devereaux, Phys. Rev. B 91, 165127 (2015).
- Ohgoe and Imada  T. Ohgoe and M. Imada, arXiv:1703.08899 (2017).
- Blankenbecler et al.  R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
- Rubtsov et al.  A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005).
- Assaad and Lang  F. F. Assaad and T. C. Lang, Phys. Rev. B 76, 035116 (2007).
- Gull et al.  E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
- Assaad  F. F. Assaad, “DMFT at 25: Infinite Dimensions: Lect. Notes of the Autumn School on Correlated Electrons,” (Verlag des Forschungszentrum Jülich, Jülich, 2014).
-  We used 5000 single-vertex updates and 8 Ising spin flips per sweep for all results shown.
- Michielsen and de Raedt  K. Michielsen and H. de Raedt, Mod. Phys. Lett. B 10, 467 (1996).
- Binder  K. Binder, Z. Phys. B Con. Mat. 43, 119 (1981).
- Hesselmann and Wessel  S. Hesselmann and S. Wessel, Phys. Rev. B 93, 155157 (2016).
- Wang et al.  L. Wang, Y.-H. Liu, J. Imriška, P. N. Ma, and M. Troyer, Phys. Rev. X 5, 031007 (2015).
- LeCun et al.  Y. LeCun, Y. Bengio, and G. Hinton, Nature 521, 436 (2015).
- Weber et al.  M. Weber, F. F. Assaad, and M. Hohenadler, Phys. Rev. B 94, 245138 (2016).
-  We used Keras to implement a network similar to Ref. : input neurons followed by a convolutional layer with 8 features, receptive fields, and sigmoid activation. After max-pooling, a second (identical) convolutional layer with ReLU activation and again max-pooling. Next, a dense layer with 256 ReLU neurons and dropout regularization. Finally, a single sigmoid neuron to produce . labeled QMC samples of the density correlator measured on a single time slice were used for supervised learning (20 epochs, batch size 20) at and , respectively. was averaged over 2000 samples to obtain Fig. 2(c).
- Carrasquilla and Melko  J. Carrasquilla and R. G. Melko, Nat. Phys. 13, 431 (2017).
- Hirsch  J. E. Hirsch, Phys. Rev. B 31, 4403 (1985).
- Hirsch and Fradkin  J. E. Hirsch and E. Fradkin, Phys. Rev. B 27, 4302 (1983).
-  At weak coupling, scaling with the 2D quantum Ising exponents is observed on similar system sizes, see arXiv:1709.01096v1. However, classical 2D Ising criticality is expected for any at sufficiently large .
- Melchert  O. Melchert, arXiv:0910.5403 (2009).
- Luther and Emery  A. Luther and V. J. Emery, Phys. Rev. Lett. 33, 589 (1974).
- Hirsch and Scalapino  J. E. Hirsch and D. J. Scalapino, Phys. Rev. Lett. 56, 2732 (1986).
- Koller et al.  W. Koller, D. Meyer, Y. Ono, and A. Hewson, Euro. Phys. Lett. 66, 559 (2004).
- Jeon et al.  G. S. Jeon, T.-H. Park, J. H. Han, H. C. Lee, and H.-Y. Choi, Phys. Rev. B 70, 125114 (2004).
- Murakami et al.  Y. Murakami, P. Werner, N. Tsuji, and H. Aoki, Phys. Rev. B 88, 125126 (2013).
- Hirsch and Scalapino  J. E. Hirsch and D. J. Scalapino, Phys. Rev. B 32, 117 (1985).
- Karakuzu et al.  S. Karakuzu, L. F. Tocchio, S. Sorella, and F. Becca, arXiv:1709.00278 (2017).
- dos Santos  R. R. dos Santos, Phys. Rev. B 50, 635 (1994).
- Randeria et al.  M. Randeria, N. Trivedi, A. Moreo, and R. T. Scalettar, Phys. Rev. Lett. 69, 2001 (1992).
- Macridin et al.  A. Macridin, G. A. Sawatzky, and M. Jarrell, Phys. Rev. B 69, 245111 (2004).
- Das and Doniach  D. Das and S. Doniach, Phys. Rev. B 60, 1261 (1999).
- Hohenadler and von der Linden  M. Hohenadler and W. von der Linden, Phys. Rev. B 71, 184309 (2005).
- Jülich Supercomputing Centre  Jülich Supercomputing Centre, J. Large-Scale Res. Facilities 2, A62 (2016).
- Broecker et al.  P. Broecker, J. Carrasquilla, R. G. Melko, and S. Trebst, arXiv:1608.07848 (2016).