Langevin dynamics simulations of the two-dimensional Su-Schrieffer-Heeger model
Anika Goetz, Stefan Beyl, Martin Hohenadler, Fakher F. Assaad
LLangevin dynamics simulations of the two-dimensional Su-Schrieffer-Heeger model
A. G¨otz, S. Beyl, M. Hohenadler, and F. F. Assaad
1, 2 Institut f¨ur Theoretische Physik und Astrophysik, Universit¨at W¨urzburg, 97074 W¨urzburg, Germany W¨urzburg-Dresden Cluster of Excellence ct.qmat, Am Hubland, 97074 W¨urzburg, Germany
We use Langevin sampling methods within the auxiliary-field quantum Monte Carlo algorithm toinvestigate the phases of the Su-Schrieffer-Heeger model on the square lattice at the O(4) symmetricpoint. Based on an explicit determination of the density of zeros of the fermion determinant, weargue that this method is efficient in the adiabatic limit. By analyzing dynamical and static quanti-ties of the model, we demonstrate that a ( π, π ) valence bond solid gives way to an antiferromagneticphase with increasing phonon frequency.
I. INTRODUCTION
Phonons in systems of electrons or spins can generatemany fascinating states of matter. Apart from supercon-ductivity [1], this also includes phases that break latticesymmetries such as charge-density wave states (CDW)and various flavors of valence-bond solid (VBS) states[2]. In spin systems, the coupling to phonons can gen-erate quantum phase transitions between antiferromag-netic (AFM) and VBS states [3]. Such interactions henceoffer the possibility of investigating quantum phase tran-sitions between various states of matter all characterizedby a local order parameter. A particularly intriguingaspect is the possibility of realizing quantum phase tran-sitions beyond the Landau-Ginzburg-Wilson paradigm inmodels relevant for materials.The Debye frequency ω D is typically much smaller thanthe Fermi energy (cid:15) F . This separation of energy scalesunderlies Migdal’s theorem [1, 4], which provides a smallparameter, (cid:126) ω D /(cid:15) F , to justify perturbative approachesto the electron-phonon problem. Fermion Monte Carlosimulations offer the possibility to take a step beyondperturbative approaches and thereby investigate compet-ing instabilities [5]. In fact, it is known that the genericelectron-phonon problem does not suffer from a negativesign problem, irrespective of the lattice geometry andband filling. In particular, for each space-time configu-ration of phonon fields, time-reversal symmetry ensuresthat the eigenvalues of the fermion determinant comein complex conjugate pairs [6]. This argument for theabsence of a negative sign problem is valid only if thephonons are not integrated out, as done in Refs. [7–9]. Assuch, it should be a technically simple problem. However,this is not the case, since the key challenge is to find anadequate sampling scheme that deals with the separationof energy scales. Adopting a local updating scheme—ascommonly used in Monte Carlo simulations of fermions—in which the phonon field is updated on a time scale setby the electron motion, leads to prohibitively long au-tocorrelation times [10]. Global updates of the phononfields on the imaginary time scale of the inverse Debyetemperature are highly desirable and have been achievedusing, for example, self-learning methods [11].In the one-dimensional case, there has been a numberof simulations of the Su-Schrieffer-Heeger (SSH) model FIG. 1. Schematic phase diagram of the O(4) symmetric SSHmodel on the square lattice as a function of phonon frequency. [2, 12]. The model exhibits many fascinating quantumphase transitions, including instances of deconfined quan-tum criticality in one dimension [13] between VBS andCDW states [14]. Simulations in two dimensions are stillin their infancy. In Ref. [15], the authors show that themodel sustains a VBS state. A parallel work to the re-sults presented here shows that the model also supportsan AFM phase [16].The present work uses Langevin updates, as success-fully applied before to the Holstein model [17]. Here,we will use an identical algorithm for the SSH modelwith Einstein phonons [18] on a square lattice at theO(4) symmetric point. In the very same way as the hy-brid Monte-Carlo (HMC) approach discussed in Ref. [19],Langevin and HMC updates are global moves for whichthe acceptance rate is unity or can be made arbitrarilylarge. The key point for the success of Langevin andHMC approaches is the absence of singularities in theaction. Since the action contains the logarithm of thefermion determinant, this requires that the latter doesnot vanish. In special cases where the action has nosingularities, this class of updating schemes works verywell. For a specific example, we refer to previous work onthe one-dimensional Hubbard model with open boundaryconditions [20].For the electron-phonon problem in the adiabatic limit,the phonon fields are frozen in imaginary time. In thislimiting case, the fermion determinant is strictly positive.A small phonon frequencies, we will show that the den-sity of zeros of the fermion determinant is small so that a r X i v : . [ c ond - m a t . s t r- e l ] F e b Langevin dynamics performs well. As the frequency ofthe phonons is increased, the density of zeros grows andthe Langevin approach becomes unstable.The main results and organization of the article are asfollows. In Sec. II, we define the SSH model considered.We will concentrate on a high-symmetry point with par-tial particle-hole symmetry that exhibits an O(2 N ) sym-metry for the general case of fermions with N flavors. InSec. III, we discuss the numerical implementation, whichmakes use of a standard auxiliary-field quantum MonteCarlo (AFQMC) [21] formulation that allows us to usethe ALF-2.0 package [20]. For the updating scheme, wehave followed the work of Ref. [17] and implemented theLangevin updating scheme with Fourier acceleration. ForO(2 N ) symmetric problems, the fermion determinant fora single fermion flavor can be written in terms of thesquare of a Pfaffian [22]. The sign of the Pfaffian is notpinned to unity. Configurations at which the sign changescorrespond to zeros of the fermion determinant. As such,the average sign of the Pfaffian provides a measure for thedensity of zeros. We will show that the average sign ofthe Pfaffian decreases with increasing phonon frequency,thus rendering Langevin simulations prohibitively hardand unstable in the high frequency limit. In Sec. IV, wepresent our results for the adiabatic limit and nonzerophonon frequencies. Within our Langevin simulations,we observe two phases. A ( π, π ) VBS phase that sponta-neously breaks the C symmetry of the lattice is observedat small phonon frequencies. At higher frequencies, theVBS state gives way to an AFM phase (see Fig 1). Thisphase spontaneously breaks the O(4) symmetry down toSU(2). We conclude in Sec. V.Aspects of this work, and in particular the transitionfrom VBS to AFM, were already published in Ref. [23]. II. MODEL AND SYMMETRIESA. Hamiltonian
We consider an SSH model with optical phonons, de-fined by the Hamiltonianˆ H el = − t (cid:88) (cid:104) i , j (cid:105) ˆ K b + g (cid:88) (cid:104) i , j (cid:105) ˆ Q b ˆ K b (1)+ (cid:88) b (cid:34) ˆ P b m + k Q b (cid:35) . The operator ˆ K b = (cid:80) Nσ =1 (cid:0) ˆ c † i ,σ ˆ c j ,σ + h.c. (cid:1) describesfermion hopping on a bond b = (cid:104) i , j (cid:105) connecting twonearest-neighbor sites i , j with hopping amplitude t . Theoperator ˆ c † i ,σ creates an electron centered at site i andwith z -component of spin σ that runs over N flavors. Weuse anti-periodic boundary conditions ˆ c † i + L a ,σ = − ˆ c † i ,σ in the direction of the primitive vector a of the latticeand periodic boundary conditions ˆ c † i + L a ,σ = ˆ c † i ,σ in the direction of a . The phonons are represented by har-monic oscillators that reside on the bonds. They aredescribed by the momentum and position operators ˆ P b and ˆ Q b as well as the frequency ω = km . The hoppingof the electrons is modulated by the phonon coordinateˆ Q b on the respective bond b with coupling strength g . B. Symmetries
The SSH model at half-filling and on a bipartite lattice,such as the square lattice considered here, is invariantunder the partial particle-hole transformationˆ P − σ ˆ c † i ,σ (cid:48) ˆ P σ = δ σ,σ (cid:48) e i Q · i ˆ c i ,σ (cid:48) + (1 − δ σ,σ (cid:48) ) ˆ c † i ,σ (2)where Q = ( π, π ) for the square lattice. We can definea corresponding Z order parameter, the fermion parityon site i [24] ˆ p i = N σ (cid:89) σ =1 (1 − n i ,σ ) . (3)This Ising-like order parameter supports order at finitetemperature. Since it changes sign under the partialparticle-hole transformation, it indicates a spontaneousbreaking of this symmetry.In addition to the apparent global SU( N ) spin rota-tional symmetry, the model possesses an enlarged O(2 N )symmetry on a bipartite lattice. To prove this, we refor-mulate the Hamiltonian using Majorana fermions [19, 24]ˆ c † i ,σ = 12 (ˆ γ i ,σ, − iˆ γ i ,σ, ) . (4)After a canonical transformation ˆ c † i → iˆ c † i on one sublat-tice, the hopping operator can be written asˆ K b = i2 (cid:88) σ (cid:88) α =1 ˆ γ i ,σ,α ˆ γ j ,σ,α (5)and the O(2 N ) symmetry becomes obvious. Because ofthis symmetry, the model is free of a sign problem for oddvalues of N [25]. For even N , time-reversal symmetry issufficient to show the absence of a sign problem [6]. In thecase of N = 2 considered here, the spin operators and theAnderson pseudospin operators [26] are the infinitesimalgenerators of the SO(4) symmetry. They are defined byˆ S i = (cid:80) σ,σ (cid:48) ˆ c † i ,σ σ σ,σ (cid:48) ˆ c i ,σ (cid:48) , ˆ η i = ˆ P − ↑ ˆ S i ˆ P ↑ , (6)where the vector σ contains the three Pauli matrices.Both, the spin and the pseudospin components ( l , m , n )fulfill the Lie algebra of the SU(2) group (cid:2) ˆ S i ,l , ˆ S j ,m (cid:3) =i δ i , j (cid:80) n ε lmn ˆ S i ,n and commute among each other. Here, ε lmn is the Levi-Civita symbol. The Lie algebra of theglobal O(4) symmetry can be interpreted as O(4) =SU(2) × SU(2) × Z , where the additional Z symmetrycorresponds to the partial particle-hole symmetry [24].Hence, an AFM phase is degenerate with a charge-densitywave (CDW) and an s-wave superconductor (SC). If theparity ˆ p i orders and the particle-hole symmetry is spon-taneously broken, either the spin or charge sector is ex-plicitly chosen. C. Limiting cases
In the adiabatic limit ω →
0, imaginary-time fluctu-ations of the phonon fields are exponentially suppressed.The phonon displacements can be treated classically andthe Hamiltonian can be written asˆ H = (cid:88) b ( − t + gq b ) ˆ K b + (cid:88) b q b . (7)Here, ˆ Q b | q (cid:105) = q b | q (cid:105) . The Hamiltonian consists solely ofa modulated hopping of the electrons and the potentialenergy of the phonon fields.For ω > S eff = − g k (cid:82) β (cid:82) β d τ d τ (cid:48) (cid:80) b ˆ K b ( τ ) D ( τ − τ (cid:48) ) ˆ K b ( τ (cid:48) ) (8)with the local and retarded interaction D ( τ ) = ω − ω | τ | + e − ω ( β −| τ | ) − e − ω β . (9)By taking the antiadiabatic limit ω → ∞ the interactionbecomes instantaneous [24],lim ω →∞ D ( τ ) = δ ( τ ) . (10)The effective Hamiltonian of the SSH model in the an-tiadiabatic limit is given byˆ H eff = − t (cid:88) (cid:104) i , j (cid:105) ˆ K b − g k (cid:88) (cid:104) i , j (cid:105) ˆ K b . (11)For N = 1, this expression is equivalent to the Hamil-tonian of the t - V model if we identify the interactionstrength with g = √ kV . For two fermion flavors, N = 2,we can rewrite the interaction term as −
14 ˆ K b = ˆ S i ˆ S j + ˆ η i ˆ η j . (12)In this form, the O(4) symmetry of the system is obvi-ous. In the high-frequency limit ω → ∞ , we expect aAFM/CDW/SC ground state. III. METHODSA. Langevin dynamics
Using the real-space formulation of the path integralfor the phonon degrees of freedom with the eigenstates of the position operator ˆ Q b | q (cid:105) = q b | q (cid:105) , the partition func-tion of the model can be written as Z = (cid:90) (cid:89) b,τ d q b,τ e − S , (13) S = S + S F = S − N ln det [ + B ( β, B ( τ , τ ) = τ (cid:89) τ = τ +∆ τ (cid:89) b (cid:0) e − ∆ τgq b,τ K b (cid:1) e ∆ τ (cid:80) b t b K b (14)and ˆ K b = (cid:80) i , j ,σ ˆ c † i ,σ ( K b ) i , j ˆ c † j ,σ . We have used theformula by Blankenbecler, Scalapino, and Sugar [21]to rewrite the trace over the fermionic degrees of free-dom as a determinant. Furthermore, τ runs from1 · · · L Trot where the inverse temperature β = L Trot ∆ τ .The phonon degrees of freedom have their own intrinsicimaginary-time dynamics governed by the action S = ∆ τ (cid:88) b,τ (cid:32) ω (cid:20) q b,τ +1 − q b,τ ∆ τ (cid:21) + q b,τ (cid:33) . (15)Since a harmonic oscillator has only one free parameter,namely the frequency ω , we set the force constant to k = 2.To update the bosonic fields q = { q b,τ } we useLangevin dynamics. The corresponding Langevin equa-tion is a stochastic differential equation for the fields[17, 20], d q ( t l )d t l = − Q ∂S ( q ( t l )) ∂ q ( t l ) + (cid:112) Q η ( t l ) , (16)with an additional Langevin time t l . The independentGaussian random variables η satisfy (cid:104) η b,τ ( t l ) (cid:105) = 0 , (cid:104) η b,τ ( t l ) η b (cid:48) ,τ (cid:48) ( t (cid:48) l ) (cid:105) = δ b,b (cid:48) δ τ,τ (cid:48) δ ( t l − t (cid:48) l )(17)with δ being the Kronecker delta for the discrete indicesand the Dirac delta function for the continuous t l . Thematrix Q is an arbitrary positive-definite matrix. In or-der to use the Langevin equation in our AFQMC codewe discretize the Langevin time t l with a finite time step δt l . Using the Euler method, the discretized equation isgiven by [28] q ( t l + δt l ) = q ( t l ) − Q ∂S ( q ( t l )) ∂ q ( t l ) δt l + (cid:112) δt l Q η ( t l ) . (18)For the random variables η we have to make the replace-ment δ ( t l − t (cid:48) l ) → δ t l ,t (cid:48) l . The systematic error introducedby discretizing the Langevin time is of linear order in δt l [17, 20]. By transforming the Langevin equation into theform of a Fokker-Planck equation one can show that thestationary probability distribution of finding the systemin state q is given by [29] P ( q ) = e − S ( q ) (cid:82) D q e − S ( q ) . (19)A major aspect of Langevin dynamics are the forces,their computation, and characteristics. For the SSHmodel, using Eq. (13), the forces read ∂S∂q b,τ = ∆ τ kq b,τ + m ∆ τ (2 q b,τ − q b,τ +1 − q b,τ − ) (20)+ N g ∆ τ Tr { K b (1 − G ( b, τ )) } with the Green function G i , j ( b, τ ) = Tr (cid:104) ˆ U < ( b, τ )ˆ c i ˆ c † j ˆ U > ( b, τ ) (cid:105) Tr (cid:104) ˆ U ( β, (cid:105) (21)and the propagatorsˆ U < ( b (cid:48) , τ ) = ˆ U ( β, τ ) N b (cid:89) b = b (cid:48) e − ∆ τgq b,τ ˆ c † K b ˆ c , (22)ˆ U > ( b (cid:48) , τ ) = b (cid:48) − (cid:89) b =1 e − ∆ τgq b,τ ˆ c † K b ˆ c e ∆ τ (cid:80) b t b ˆ c † K b ˆ c ˆ U ( τ − ∆ τ, . From Eq. (13) we see that the action has logarithmic di-vergences if the determinant vanishes. The O(2 N ) sym-metry of the model only guarantees that the determinantis non-negative. B. Calculation of the Pfaffian
To consider possible divergences of the action we an-alyze the determinant that appears in the action (13).First, we use the O(2 N ) symmetry of the model to ex-press the determinant as a square of a trace over one ofthe two Majorana fermions: Z γ = det [ + B ( β, Z γ = Tr (cid:34)(cid:89) τ (cid:89) b (cid:16) e − i2 ∆ τgq b,τ ˆ γ i ˆ γ j (cid:17) e i2 ∆ τ (cid:80) b t b ˆ γ i ˆ γ j (cid:35) . Here and in the rest of this section, we drop the spinand Majorana kind indices since none of the quantitiesconsidered explicitly depend on them. One can show witha canonical transformation of the Majorana fermions ononly one sublattice, ˆ γ i → − ˆ γ i , that Z γ is real [22] andits square is non-negative.The quantity Z γ can either be positive or negative indistinct regions of the configuration space. Since Z γ isan entire function, it necessarily has to vanish betweenthese regions. Hence, the average sign of Z γ serves as anestimate of the number of zeros of the determinant. Ifthe average sign is close to plus or minus unity, it is moreimprobable to cross a boundary between two regions inwhich Z γ has different signs. In contrast, for a vanishingaverage sign, there are more zeros.To measure the sign of Z γ we reformulate it as a Pfaf-fian. First, we use an alternative Trotter decomposition and rewrite the exponentials as hyperbolic functions byusing (ˆ γ i ˆ γ j ) = − Z γ = Tr (cid:34)(cid:89) t e i y t ˆ γ i ˆ γ j (cid:35) (24)= (cid:89) t (cosh y t ) Tr (cid:34)(cid:89) t (1 + iˆ γ i ˆ γ j tanh y t ) (cid:35) . (25)The tuple t = ( b, τ ) combines the bond index and theimaginary time slice into a new index that is ordered ac-cording to its position in the product (cid:81) τ (cid:81) b . To lightenthe notation, we used y b,τ = ∆ τ ( t b − gq b,τ ). Next, weintroduce Grassmann variables ξ i / j ,τ on every site andimaginary time slice, where i is on sublattice A and j on sublattice B, and use the following formula for even n [30]: C ± n (cid:89) t =1 (cid:112) a ( t ) = (cid:90) [d ξ ] e ± (cid:80) t 81 0 1 2 3 4 (a) s g n ( P f A ) VL = 4 , βt = 4 , δt l = 0 . . . . . 981 0 1 2 3 4 5 6 (b) s g n ( P f A ) ω a-p, L = 4p-p, L = 4a-p, L = 6p-p, L = 6 FIG. 2. (a) Average sign of the Pfaffian for the t - V modelas a function of V /t . Here, we consider a square lattice witha π -flux threading each plaquette. This model supports aGross-Neveu transition at V c = 1 . t [32]. We observea strong decrease of the average sign of the Pfaffian alreadyat high temperatures for increasing coupling strength. (b)Average sign of the Pfaffian for the SSH model with differentboundary conditions at βt = 5 . 0, ∆ τ t = 0 . 1, and δt l = 0 . time step, the fermionic forces for every b and τ are com-pared to a preset maximal force F max . If the maximalcomputed force max (cid:0) ∂S F ∂q b,τ (cid:1) exceeds F max the Langevintime step is decreased by the ratio of the two forces,¯ δt l = F max max (cid:16) ∂S F ∂q b,τ (cid:17) δt l . (28)The variations of the Langevin time step have to be takeninto account when measuring observables, (cid:68) ˆ O (cid:69) = 1 N m (cid:80) N m α =1 (cid:0) ¯ δt l (cid:1) α (cid:68)(cid:68) ˆ O (cid:69)(cid:69) α (cid:80) N m α =1 (cid:0) ¯ δt l (cid:1) α . (29)Here, N m is the total number of measurements and (cid:104)(cid:104) ˆ O (cid:105)(cid:105) α denotes the value of the observable ˆ O for con-figuration C α of the phonon fields. C. Fourier acceleration Following Refs. [17, 29] we used Fourier acceleration(FA) to decrease autocorrelations. The main idea of FAis to adapt the Langevin time step δt l according to themomentum of the phononic modes. Slow modes are mul-tiplied with a larger step size than fast modes by using anadequate choice of the matrix Q in the Langevin equation[33].As a foundation for the choice of Q we consider thenon-interacting case ( g = 0). We carry out a Fouriertransformation of the force in imaginary time,ˆ F (cid:20) d S d q b,τ (cid:21) = (cid:104) ∆ τ k + m ∆ τ (cid:16) − cos (cid:16) πk τ L τ (cid:17)(cid:17)(cid:105) ˜ q b,k τ , (30)with ˆ F [ f ( τ )] = ˜ f ( k τ ) = L τ (cid:80) L τ τ =1 e i πLτ k τ τ f ( τ ) . (31)The momentum takes the values k τ = − L τ + 1 , − L τ +2 , ..., L τ . As noted in Ref. [17], the ratio of the slowest . . . . . . 001 0 . 01 0 . (a) S S ( q = ( π , π ) , ) /t l no FAFA − . . . . . 81 0 50 100 150 200 (b) C S ( t l ) t l no FAFA FIG. 3. (a) Equal-time spin correlation function S S ( q , C S ( t l ) at wavevector q = ( π, π ) as a function of Langevin time t l with andwithout the use of Fourier acceleration (FA). Here, ω = 0 . L = 6, β = 7 . 0, and δt l = 0 . to fastest mode is (∆ τ ) k (∆ τ ) k + 4 m (cid:28) . (32)Especially in the limit of small phonon frequencies ω theratio is close to zero. We choose the factor Q in Fourierspace such that the prefactor of ˜ q b,k τ in Eq. (30) becomesindependent of the momentum k τ ,˜ Q ( k τ ) = ∆ τ k + m ∆ τ ∆ τ k + m ∆ τ (cid:16) − (cid:16) πk τ L τ (cid:17)(cid:17) . (33)Although this choice is guided by the non-interactingcase g = 0, we also use it for g > q ( t l + δt l ) = q ( t l ) (34) − ˆ F − (cid:20) δt l ˜ Q ( k τ ) ˆ F (cid:20) d S d q ( t l ) (cid:21) − (cid:112) δt l (cid:113) ˜ Q ( k τ ) ˆ F [ η ( t l )] (cid:21) . To see the effect of FA on autocorrelation times, wemeasured the equal-time spin correlation function S S ( q , τ = 0) = 1 L (cid:88) i , j e − i q ( i − j ) (35) × (cid:16)(cid:68) ˆ S i ,z ˆ S j ,z (cid:69) − (cid:68) ˆ S i ,z (cid:69) (cid:68) ˆ S j ,z (cid:69)(cid:17) with FA [using Eq. (34) to update the fields] and withoutFA (by setting Q = 1). The spin correlation function atwave vector q = ( π, π ) is shown in Fig. 3(a) as a functionof the inverse Langevin time. The equilibration time isobviously reduced by FA and, as expected, the results ofboth methods agree at sufficiently long times.We also consider the autocorrelation function [20] C ˆ O ( t l ) = T l − t l (cid:88) t (cid:48) l =0 (cid:16) O ( t (cid:48) l ) − (cid:68) ˆ O (cid:69)(cid:17) (cid:16) O ( t (cid:48) l + t l ) − (cid:68) ˆ O (cid:69)(cid:17)(cid:16) O ( t (cid:48) l ) − (cid:68) ˆ O (cid:69)(cid:17) . The observable evaluated at time t l is given by O ( t l )and T l is the maximal time at which measurements aretaken. For uncorrelated measurements, the autocorrela-tion function drops to zero. The shorter the autocorrela-tion time, the faster the autocorrelation function shoulddecay to zero. A decrease of autocorrelations by FA isclearly visible from Fig. 3(b). IV. RESULTS All simulations were carried out using the ALF package[20]. In our simulations, we set k = 2, t = 1, g = 1 . τ = 0 . 1. We used a Langevin time step of δt l = 0 . 01 inthe adiabatic limit and δt l = 0 . A. Adiabatic limit With the aid of a mean-field approach we study thepattern of the fields q b in the ground state. The nestingof the Fermi surface of the non-interacting model givesrise to a log divergence in the Q = ( π, π ) susceptibilityat low temperatures. Hence, we expect an instability toa ( π, π ) VBS. In contrast to one dimension, where theordering pattern is unique, in two dimensions there is amultitude of generalizations including staircase, colum-nar, staggered, and plaquette arrangements [34, 35]. Inour calculation, we use an enlarged 2 × π, π ) and (0 , π ) VBSorders. After minimization of the energy, the pattern weobtain is the ( π, π ) staggered VBS state illustrated inFig. 4. q q q q q q q q ……… … (a) q q q q q q q q ……… … (a) · · · · · · ...... · · · · · · ...... · · · · · · ...... · · · · · · ...... Figure 6.6: ( π, π ) VBS pattern according to our mean field approach. Strong bonds are coloredwhile weak bonds are represented by solid lines.In order to check our implementation we compared the energy of the free system and predefinedpatterns with results from Ref. [55]. Here we took into account that the authors also neglectedthe φ b terms in this specific publication and reproduced the benchmark data. Contrary to theexpectations in Refs. [66, 67], we achieve the lowest energy according to the ( π, π ) pattern ofFig. 6.4b for including the φ b terms. For more detailed discussions about the patterns of Fig. 6.5we refer to the literature, see e.g. Refs. [57, 58].From now on we do not restrict ourselves to predefined patterns. Instead we are using all degreesof freedom to find the ground state. Due to the C lattice symmetry, the global minimum ofthe energy function in Eq. (6.9) is four-fold degenerate. Minimizing the mean field energy for agiven set of parameters reveals a φ field configuration that can be interpreted as a VBS pattern.Table 6.1 shows a sample ground state configuration of φ . φ φ φ φ φ φ φ φ -1.26 -0.32 -0.32 -1.26 -0.32 -0.32 -0.32 -0.32 Table 6.1: A sample mean field solution for t = 1 . , ˜ g = 1 . and N x = N y = 50 .In general, the phonons strengthen the electronic hopping on all bonds. However, two bondsdominate. The values on the first and fourth bond are much larger than the others. Every latticesite has only one strong bond and three relatively weak bonds, see Fig. 6.6. The pattern ofstrong bonds corresponds to a ( π, π ) ordering as expected from the shape of the Fermi surface inFig. 6.3.Proceeding from this mean field observation we will show how the ground state evolves at finitephonon frequencies ω using the HQMC method.60 (b) FIG. 4. (a) Unit cell used in the mean-field ansatz and(b) ( π, π ) VBS pattern suggested by our mean-field analysis.Strong bonds are colored while weak bonds are representedby black lines. The phonons enhance the hopping amplitude on allbonds and effectively renormalize the bandwidth of theelectronic band structure. Furthermore, they modulatethe hopping in a ( π, π ) pattern (see Fig. 4(b)) that leadsto a finite gap at the Fermi surface. Since the VBS order-ing breaks the discrete C symmetry of the lattice, theordering can survive at finite temperatures.In addition to the mean-field analysis, we used a clas-sical Monte Carlo approach to study the adiabatic limit.In this way, we can observe the melting of VBS order asa function of the temperature T = β − . The partition function simplifies to Z = (cid:90) (cid:89) b d q b e − S , (36) S = β (cid:88) b k q b − N σ ln det (cid:104) + e − β (cid:80) b ( − t + gq b ) K b (cid:105) . The determinant is strictly positive since its argumentis a symmetric matrix. Therefore, we can use Langevindynamics without divergences in the forces to update thephonon fields via ∂S∂q b = βkq b + βgN Tr { K b (1 − G ) } , (37) G i , j = Tr (cid:104) e − β (cid:80) b ( − t b + gq b ) ˆ K b ˆ c i ˆ c † j (cid:105) Tr (cid:104) e − β (cid:80) b ( − t b + gq b ) ˆ K b (cid:105) . To capture VBS order we measured the bond kinetic sus-ceptibility χ δ,δ (cid:48) K ( q ) = (cid:90) β d τ S δ,δ (cid:48) K ( q , τ ) (38)with the imaginary-time-displaced correlation function S δ,δ (cid:48) K ( q , τ ) = (cid:68) ˆ K δ ( q , τ ) ˆ K δ (cid:48) ( − q ) (cid:69) − (cid:68) ˆ K δ ( q ) (cid:69) (cid:68) ˆ K δ (cid:48) ( − q ) (cid:69) (39)and ˆ K δ ( q ) = 1 √ N (cid:88) i ,σ e i q · i (cid:16) ˆ c † i ,σ ˆ c i + a δ ,σ + H.c. (cid:17) . (40)Figure 5 shows the bond susceptibility as a function oftemperature and system size at the ordering wave vector Q = ( π, π ). . 02 0 . 06 0 . . 14 0 . (a) T r χ K ( q = ( π , π )) T L = 4 L = 6 L = 8 L = 10 L = 12 − − − − . 02 0 . 06 0 . . 14 0 . (b) h ˆ H i T FIG. 5. (a) Bond susceptibility as a function of temperatureand system size. (b) Energy as a function of temperature for L = 12. Simulations were done at δt l = 0 . 01 and startingfrom the mean-field configuration. The simulations were started in the mean-field config-uration to reduce warm-up times. At low temperatures,the susceptibility grows in accordance with a long-ranged( π, π ) VBS state. On our largest lattice size ( L = 12),we observe a sudden drop of the signal at T ≈ . 06 (seeFig. 5(a)). The energy (cid:104) ˆ H (cid:105) shows a kink at the sametemperature [Fig. 5(b)]. Above this temperature, ther-mal fluctuations destroy the long-range order. B. Finite phonon frequencies 1. Equal time and static quantities To map out the phases as a function of phonon fre-quency we computed spin-spin correlations, S S ( q , τ ) = (cid:68) ˆ S z ( q , τ ) ˆ S z ( − q ) (cid:69) − (cid:68) ˆ S z ( q , τ ) (cid:69) (cid:68) ˆ S z ( − q ) (cid:69) , (41)as well as the imaginary-time-displaced correlations ofthe bond kinetic energy defined in Eq. (39). Here,ˆ S z ( q ) = 1 √ N (cid:88) i e i q · i (ˆ n i , ↑ − ˆ n i , ↓ ) . (42)We also consider the bond-kinetic susceptibility ofEq. (38) and the equivalent form of the spin suscepti-bility, χ S ( q ). As discussed in Sec. II, our Hamiltonianhas an O(4) symmetry. Hence, the three components ofthe spin-spin correlations are degenerate with the charge-density wave (CDW) and s-wave superconducting (SC)correlations. As argued in Ref. [36], introducing an in-finitesimally small positive (negative) Hubbard interac-tion will select the spin (CDW/SC) correlations. Here,we will discuss the results from the point of view of spin-spin correlations. . 08 0 . 12 0 . 16 0 . . (a) T r χ K ( q = ( π , π )) /L ω = 0 . ω = 0 . ω = 0 . ω = 0 . ω = 1 . ω = 2 . . 02 0 . 06 0 . . 14 0 . (b) T r S K ( q = ( π , π ) , ) T . 02 0 . 06 0 . . 14 0 . (c) T r χ K ( q = ( π , π )) T FIG. 6. (a) Size scans of the bond-kinetic susceptibility at βt = 40. (b) Temperature scans at L = 12 for the bond-kinetic structure factor Tr S K ( q = ( π, π ) , τ = 0). (c) Temper-ature scans of the bond-kinetic susceptibility for L = 12. In Fig. 6 we present the dependence of the bond-kineticsusceptibility and structure factor on lattice size, temper-ature, and phonon frequency. As apparent, at the low-est frequency considered ( ω = 0 . χ K ( q = ( π, π ))grows as a function of size and inverse temperature, sug-gesting long-range ( π, π ) VBS order as in the mean-field approximation. Note that cooling down our sys-tem to βt = 40 was not sufficient to achieve conver-gence of Tr χ K ( q = ( π, π )) on an L = 12 lattice. Asthe phonon frequency grows, we observe a rapid drop . 08 0 . 12 0 . 16 0 . . (a) χ S ( q = ( π , π )) /L ω = 0 . ω = 0 . ω = 1 . ω = 1 . ω = 1 . ω = 2 . . . . . . . . 02 0 . 06 0 . . 14 0 . (b) S S ( q = ( π , π ) , ) T . . . . . 02 0 . 06 0 . . 14 0 . (c) χ S ( q = ( π , π )) T FIG. 7. (a) Size scans of the spin susceptibility at βt = 40.(b) Temperature scans of the spin structure factor S S ( q =( π, π ) , τ = 0) for L = 10. (c) Temperature scans of the spinsusceptibility for L = 10. in Tr χ K ( q = ( π, π )) that indicates that the ( π, π ) VBSstate gives way to another phase. . . . . . . . . . . . . . . . R χ , S ω L = 4 L = 6 L = 8 L = 10 L = 12 FIG. 8. Correlation ratio based on the spin susceptibility asdefined in Eq. (43). Here, βt = 40, which is representative ofthe ground state for this quantity. In Fig. 7(a), we show results for the spin degrees of free-dom. At low temperatures, the size-dependence of theantiferromagnetic spin susceptibility shows a marked in-crease at large phonon frequencies. In Figs. 7(b) and (c),the temperature dependence at fixed lattice size showsthat we are able to achieve convergence with respect totemperature. This allows us to compute the correlationratio R χ,S = 1 − χ S ( q = ( π, π ) + ∆ q ) χ S ( q = ( π, π )) (43)where | ∆ q | = 2 π/L . This renormalization group invari-ant quantity takes the value of unity (zero) in the ordered(disordered) phase. At T = 0 and for a continuous tran-sition, it scales as R χ,S = f (cid:16) [ ω − ω c ] L /ν (cid:17) . (44)In Fig. 8 we plot this quantity as a function of systemsize for the lowest temperature available (representativeof the ground state). Although corrections to scaling,not included in Eq. (44), lead to a meandering of thecrossing point, the results suggest a critical phonon fre-quency ω c (cid:39) . . . . . . . . . . . . . . . χ p ( q = ( , )) ω L = 4 L = 6 L = 8 L = 10 L = 12 FIG. 9. Parity susceptibility as defined in Eq. 45 at βt = 10as a function of phonon frequency ω . Being a modulation of the bond kinetic energy, the( π, π ) VBS state does not break the underlying O(4) sym-metry of the lattice. However, it does break translationand rotational symmetries. On the other hand, the AFMphase breaks the O(4) symmetry. To document this sym-metry breaking, we have computed the susceptibility ofthe parity operator, χ p ( q ) = (cid:90) β d τ (cid:88) r e i q · r (cid:104) ˆ p r ( τ )ˆ p (cid:105) . (45)Since ˆ p i is an Ising variable that changes sign under an O (4) transformation M with det M = − 1, we expect χ p ( q = ) to diverge at a critical temperature corre-sponding to a 2D Ising transition. Being an 8-point cor-relation function, χ p ( q ) becomes very noisy at low tem-peratures and we are restricted to βt = 10. Our data as afunction of system size and phonon frequency are plottedin Fig. 9. For ω = 2, χ p ( q = ) grows as a function ofsystem size, suggesting that for this frequency the Isingtemperature is below T = 0 . t . On the other hand, for ω = 1 (still in the AFM phase) our temperature is cer-tainly too high to capture the Ising transition. Based onthese data, we conclude that the AFM phase breaks theO(4) symmetry down to SO(4) at a finite-temperatureIsing transition occurring at T Ic . A natural conjecture isthat T Ic vanishes at ω c . 2. Dynamical quantities We have used the ALF implementation [20] of thestochastic maximum entropy algorithm [37, 38] to carryout the Wick rotation from imaginary to real time. Thesingle-particle spectral function, A ( k, ω ) is related to theimaginary-time Green function via (cid:104) ˆ c k .σ ( τ )ˆ c † k ,σ (0) (cid:105) = 1 π (cid:90) d ω e − τω e − βω A ( k , ω ) . (46) -10-5 0 5 10(0,0) (0, π ) ( π , π ) (0,0) -10-5 0 5 10(0,0) (0, π ) ( π , π ) (0,0) -10-5 0 5 10(0,0) (0, π ) ( π , π ) (0,0) -10-5 0 5 10(0,0) (0, π ) ( π , π ) (0,0) π ) ( π , π ) (0,0) ω = 0.4 0.01 0.1 1 10 FIG. 10. Single particle spectral function at (a) ω = 0 . ω = 1, (c) ω = 1 . 5, (d) ω = 2 . 0. Here, L = 12, βt = 40. Figure 10(a) shows the single-particle spectral func-tion at βt = 40 for L = 12. The coupling of the Einsteinphonon mode to the electrons breaks the ˆ Q b → − ˆ Q b symmetry. Consequently, N (cid:80) b (cid:10) ˆ Q b (cid:11) acquires a non-zero expectation value that renormalizes the bandwidth.At ω = 0 . 4, we measure N (cid:80) b (cid:10) ˆ Q b (cid:11) = − . g = 1 . t eff = 1 . 85. This explains the observed range ofthe band from − t eff at k = 0 to 4 t eff at k = ( π, π ). Atthis phonon frequency, the ( π, π ) modulation of the hop-ping is at the origin of the gap that opens along the non-interacting Fermi surface ( k = (0 , π ) and k = ( π/ , π/ t eff are features that can be accounted for at themean-field level. However, the spectral function exhibitslow-lying spectral weight that extends over the consid-ered path in the Brillouin zone. In analogy with theone-dimensional case [39], we attribute this low-energyweight to polaron formation.As the phonon frequency is enhanced to ω = 2 (seeFig. 10(b)-(d)), N (cid:80) b (cid:10) ˆ Q b (cid:11) remains almost constantand, consequently, the width of the cosine band doesnot change substantially. Upon inspection, we see thatweight is transferred to the polaron. The origin of the gapat ω > ω c is to be found in antiferromagnetic ordering.On a general basis, the O (4) symmetry of the modelleads to the relation A ( k , ω ) = A ( k + Q , − ω ) with Q =( π, π ). Hence, in the absence of symmetry breaking, thepolaron band is nested and should show instabilities toantiferromagnetic order or to a ( π, π ) VBS order.Figure 11 shows the VBS dynamical structure factor S K ( q , ω ) at four different phonon frequencies. Using π ) ( π , π ) (0,0) ω = 0.4 0.01 0.1 1 10 π ) ( π , π ) (0,0) π ) ( π , π ) (0,0) π ) ( π , π ) (0,0) π ) ( π , π ) (0,0) FIG. 11. VBS dynamical structure factor at (a) ω = 0 . ω = 1, (c) ω = 1 . 5, (d) ω = 2 . 0. Here, L = 12, βt = 40. the maximum entropy method, we compute the imag-inary part of the VBS dynamical spin susceptibility,Tr χ (cid:48)(cid:48) ( q, ω ), usingTr S K ( q, τ ) = 1 π (cid:90) d ω e − τω − e − βω Tr χ (cid:48)(cid:48) ( q, ω ) . (47)The dynamical VBS structure factor then reads S K ( q , ω ) = Tr χ (cid:48)(cid:48) ( q, ω )(1 − e − βω ) . (48)Since the phonons couple to the bond kinetic energy,we expect S K ( q , ω ) to reveal both the phonon dynam-ics and the particle-hole continuum. At ω = 0 . Q = ( π, π ) in accordancewith our expectations. However, we cannot resolve thedispersion relation. In comparison to the bandwidth, therenormalized phonon modes are very slow and are at theorigin of long autocorrelation times. As the frequencygrows and in the antiferromagnetic phase (Fig. 11(b)-(c)), the phonon mode acquires a gap. π ) ( π , π ) (0,0) ω = 0.4 0.01 0.1 1 10 π ) ( π , π ) (0,0) π ) ( π , π ) (0,0) π ) ( π , π ) (0,0) π ) ( π , π ) (0,0) FIG. 12. Spin dynamical structure factor at (a) ω = 0 . ω = 1, (c) ω = 1 . 5, (d) ω = 2 . 0. Here, L = 12, βt = 40. Finally, in Fig. 12, we show the dynamical spin struc-ture factor, S S ( q , ω ). Because phonons do not carry a spin, they are not visible in spin-flip scattering processes.The data in Fig. 12 (a) show that S S ( q , ω ) is dominatedby the particle-hole continuum. Importantly, at low en-ergies, and as we enter the antiferromagnetic phase (seeFig. 12(b)-(d)), spectral weight develops at Q = ( π, π )and reflects the onset of antiferromagnetic ordering. V. DISCUSSION AND CONCLUSIONS We have used a Langevin dynamics updating schemewith Fourier acceleration [17] in the framework of theauxiliary-field quantum Monte Carlo method to studythe physics of the two-dimensional SSH model as a func-tion of phonon frequency.In contrast to Ref. [19], we computed forces exactlyfor a given field configuration. A comparison betweenstochastic and deterministic calculations of forces can befound in Ref. [40]. Although the CPU time per sweepis longer and scales as L β , fluctuations, especially fortime-displaced correlation functions, are smaller. Oneof the key difficulties encountered in Langevin dynamicsare zeros of the determinant that lead to logarithmic sin-gularities of the action. Using the O(4) symmetry of theSSH model, the determinant can be written as the squareof a Pfaffian. The average sign of the Pfaffian provides ameasure for the density of zeros. We demonstrated thatfor small phonon frequencies the density of zeros is smallso that Langevin simulations can be stabilized using anadaptive time step scheme. Nevertheless, simulationsoccasionally suffer from spikes in observables when thestochastic walk approaches a zero. Obviously, such con-figurations have very small weight and a hybrid molecu-lar dynamics update may be more efficient. However, asshown in Ref. [19], this can lead to ergodicity issues thatcan be overcome by a complexification of the fields. Forthe Hubbard model, this is possible since the decouplingof the interaction can be done in various channels. In thecase of phonons, we do not have such liberty. We believethat global Langevin updates are a good choice in theadiabatic limit where local moves fail. As the phonon fre-quency grows, Langevin dynamics becomes increasinglychallenging but local updates, as discussed for examplein Ref. [41], become increasingly attractive.We have computed static and dynamical quantities toelucidate the physics of the SSH model at a fixed electron-phonon coupling as a function of the phonon frequency.The phase diagram is tied to the O(4) symmetry of themodel. In particular, the single-particle spectral functionsatisfies the condition A ( k , ω ) = A ( k + Q , − ω ) with Q =( π, π ). Hence, any Fermi liquid state that does not breakthis symmetry will ultimately be unstable to orders thatcan open up a gap. This includes the ( π, π )-VBS phaseas well as antiferromagnetism. In the adiabatic limit,the problem simplifies since the phonons become classi-cal. The ground state is a ( π, π ) VBS phase. At nonzerophonon frequencies, we can integrate out the phonons infavor of a retarded interaction. At equal times, the latter0reduces to g k (cid:80) (cid:104) i , j (cid:105) ˆ S i ˆ S j +ˆ η z i ˆ η z j − ˆ η y i ˆ η y j − ˆ η x i ˆ η x j . Note thatfor non-frustrating interactions, the sign of the transverseAnderson ˆ η -operators can be reversed via a canonicaltransformation, so as to yield Eq. (12). This interac-tion triggers AFM order that is degenerate with CDWand SC states. The ( π, π ) VBS breaks lattice symmetriesbut not the O(4) symmetry. On the other hand, AFMor SC/CDW phases break the O(4) symmetry down toSU(2). This symmetry reduction occurs in two steps. Atfinite temperature, we expect an Ising transition corre-sponding to the spontaneous ordering of the parity opera-tor. Even parity corresponds to the SC/CDW phase andodd parity to the AFM phase. At T = 0, the SU(2) spin(pseudospin) symmetry is spontaneously broken, leavingthe SU(2) pseudospin (spin) symmetry unbroken. Thesingle-particle spectral function supports the picture ofa narrow polaronic band undergoing a transition from a( π, π ) VBS to AFM/CDW/SC. In the particle-hole chan-nel, the dynamical VBS correlation function reveals thephonon dynamics as a function of decreasing phonon fre-quency, including a softening at ( π, π ). On the otherhand, as ω grows, we observe enhanced spectral weightat low energies and at Q = ( π, π ) in the spin channel. (a) ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... (b) ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... (c) ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... ··· ··· ...... Figure 6.18: Defects of the ( π, ordered cVBS states (a) and the ( π, π ) ordered sVBS patterns,(b) and (c).of Ref. [77] investigate the low-energy field theory of sVBS phase and see large differences to thecase of cVBS.Unfortunately, at the end of this chapter we remain without a clear interpretation of the frequency-dependent phase transition in the two-dimensional half filled SSH model. After all, a weak firstorder transition can not be excluded for sure. A standard DQCP mechanism does not apply.Furthermore, a fine tuned Ginzburg-Landau theory for ˜ g = 1 . seems to be unlikely too.At the end several scenarios seem to be conceivable to explain our observations. A metallic pointexactly at the critical point separating the AFM and VBS phase could be one interpretationconsistent with the nesting of the Fermi surface and explain the suppression/closing of the singleparticle gap. This could be some kind of special fine tuning that occurs for every strength of thecoupling. Considering the nesting of the Fermi surface we do also not expect a disordered phasebetween the VBS and AFM state. However, we can not finally exclude a small regime of coexistingphases. A scenario where the AFM starts to form while the VBS is weakening for increasing thephononic frequency could be an option that fits in the Ginzburg-Landau theory. Finally, alsoa different mechanism we do not know yet could be responsible for the phase transition. Asconcluded in Ref. [77], some field theoretical aspects of the staggered VBS phase are still underinvestigation and may hold some surprises.In this chapter we provided some first insights into the physics of the half filled two-dimensionalSSH model using a Monte Carlo method. We discussed the VBS pattern in a mean field approachas well as for the results of the HQMC method and achieved good agreement. We have shownthe O(4) symmetric AFM state in the simulation results of our Monte Carlo method and we75 FIG. 13. A vortex of the C symmetry breaking VBS statewith a featureless core. Generically the O(4) symmetry will be broken downto SU(2) by, for example, adding a next-nearest-neighbor hopping. In this case, we expect that the phase diagramwill be dominated by an SC phase. This conjecture isbased on the fact that the Cooper instability is insen-sitive to the shape of the Fermi surface. The stabilityof the VBS phase as a function of an O(4) symmetrybreaking interaction such as a chemical potential or anext-nearest-neighbor hopping is certainly an interestingtopic for future investigations.The nature of the transition remains elusive. We can-not understand it in terms of a deconfined quantum crit-ical point (DQCP) [42, 43]. To see this, we can adoptthe Dirac fermion understanding of DQCP put forwardin Refs. [44–46] in terms of five anti-commuting AFMand VBS mass terms of an 8-component Dirac metal.The algebra of the mass terms guarantees that the vor-tex core of the VBS order carries a spin-1/2 excitation.In this framework, the VBS order has momentum (0 , π )and ( π, π, π ) VBS observed heredoes not correspond to a Dirac mass term. This point ofview is substantiated by noticing that a C vortex of the( π, π ) VBS can be trivial, as shown in Fig. 13.During the preparation of this manuscript, we becameaware of the work of Ref. [16]. Our results are in agree-ment with their findings. ACKNOWLEDGMENTS [1] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory ofsuperconductivity, Phys. Rev. , 1175 (1957).[2] M. Hohenadler and H. Fehske, Density waves in stronglycorrelated quantum chains, The European Physical Jour-nal B , 204 (2018).[3] M. Weber, Valence bond order in a honeycomb antiferro-magnet coupled to quantum phonons, Phys. Rev. B ,L041105 (2021).[4] A. Migdal, Interaction between electrons and lattice vi-brations in a normal metal, JETP , 996 (1958).[5] I. Esterlis, B. Nosarzewski, E. W. Huang, B. Moritz, T. P.Devereaux, D. J. Scalapino, and S. A. Kivelson, Break-down of the migdal-eliashberg theory: A determinantquantum monte carlo study, Phys. Rev. B , 140501 (2018).[6] C. Wu and S.-C. Zhang, Sufficient condition for absenceof the sign problem in the fermionic quantum monte carloalgorithm, Phys. Rev. B , 155115 (2005).[7] F. F. Assaad and T. C. Lang, Diagrammatic determinan-tal quantum monte carlo methods: Projective schemesand applications to the hubbard-holstein model, Phys.Rev. B , 035116 (2007).[8] C. Chen, X. Y. Xu, Z. Y. Meng, and M. Hohenadler,Charge-density-wave transitions of dirac fermions cou-pled to phonons, Phys. Rev. Lett. , 077601 (2019).[9] S. Karakuzu, K. Seki, and S. Sorella, Solution of the signproblem for the half-filled hubbard-holstein model, Phys.Rev. B , 201108 (2018). [10] M. Hohenadler and T. C. Lang, Autocorrelations inquantum monte carlo simulations of electron-phononmodels, in Computational Many-Particle Physics , editedby H. Fehske, R. Schneider, and A. Weisse (Springer-Verlag Berlin Heidelberg, 2008) pp. 357–366.[11] C. Chen, X. Y. Xu, J. Liu, G. Batrouni, R. Scalettar,and Z. Y. Meng, Symmetry-enforced self-learning montecarlo method applied to the holstein model, Phys. Rev.B , 041102 (2018).[12] E. Fradkin and J. E. Hirsch, Phase diagram of one-dimensional electron-phonon systems. i. the su-schrieffer-heeger model, Phys. Rev. B , 1680 (1983).[13] C. Mudry, A. Furusaki, T. Morimoto, and T. Hikihara,Quantum phase transitions beyond landau-ginzburg the-ory in one-dimensional space revisited, Phys. Rev. B ,205153 (2019).[14] M. Weber, F. Parisen Toldin, and M. Hohenadler, Com-peting orders and unconventional criticality in the su-schrieffer-heeger model, Phys. Rev. Research , 023013(2020).[15] B. Xing, W.-T. Chiu, D. Poletti, R. T. Scalettar, andG. Batrouni, Quantum monte carlo simulations of the 2dsu-schrieffer-heeger model, Phys. Rev. Lett. , 017601(2021).[16] X. Cai, Z.-X. Li, and H. Yao, Antiferromagnetisminduced by electron-phonon-coupling, arXiv:2102.05060(2021), arXiv:https://arxiv.org/abs/2102.05060 [cond-mat.str-el].[17] G. G. Batrouni and R. T. Scalettar, Langevin simulationsof a long-range electron-phonon model, Phys. Rev. B ,035114 (2019).[18] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Soliton ex-citations in polyacetylene, Phys. Rev. B , 2099 (1980).[19] S. Beyl, F. Goth, and F. F. Assaad, Revisiting the hybridquantum monte carlo method for hubbard and electron-phonon models, Phys. Rev. B , 085144 (2018).[20] ALF Collaboration, F. F. Assaad, M. Bercx, F. Goth,A. G¨otz, J. S. Hofmann, E. Huffman, Z. Liu,F. Parisen Toldin, J. S. E. Portela, and J. Schwab,The ALF (Algorithms for Lattice Fermions) project re-lease 2.0. Documentation for the auxiliary-field quan-tum Monte Carlo code, arXiv:2012.11914 (2021),https://arxiv.org/abs/2012.11914 [cond-mat.str-el].[21] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar,Monte carlo calculations of coupled boson-fermion sys-tems., Phys. Rev. D , 2278 (1981).[22] Z.-X. Li, Y.-F. Jiang, and H. Yao, Solving the fermionsign problem in quantum monte carlo simulations by ma-jorana representation, Phys. Rev. B , 241117 (2015).[23] S. Beyl, Hybrid Quantum Monte Carlo for CondensedMatter Models , Doctoral thesis, Universit¨at W¨urzburg(2020).[24] F. F. Assaad and T. Grover, Simple fermionic model ofdeconfined phases and phase transitions, Phys. Rev. X ,041049 (2016).[25] Z.-X. Li, Y.-F. Jiang, and H. Yao, Majorana-time-reversal symmetries: A fundamental principle for sign-problem-free quantum monte carlo simulations, Phys.Rev. Lett. , 267002 (2016).[26] P. W. Anderson, Random-phase approximation in thetheory of superconductivity, Phys. Rev. , 1900 (1958).[27] J. W. Negele and H. Orland, Quantum Many-particle Systems (Westview Press, 1998).[28] C. Gardiner, Stochastic Methods (Springer, Berlin, Hei-delberg, 2009).[29] G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage,B. Svetitsky, and K. G. Wilson, Langevin simulations oflattice field theories, Phys. Rev. D , 2736 (1985).[30] E. Huffman, Fermion Bag Approach for Hamiltonian Lat-tice Field Theories , Ph.D. thesis, Duke University (2018).[31] M. Wimmer, Algorithm 923, ACM Transactions onMathematical Software , 30 (2012).[32] E. Huffman and S. Chandrasekharan, Fermion-bag in-spired hamiltonian lattice field theory for fermionic quan-tum criticality, Phys. Rev. D , 074501 (2020).[33] C. Davies, G. Batrouni, G. Katz, A. Kronfeld, P. Lepage,P. Rossi, B. Svetitsky, and K. Wilson, Langevin simula-tions of lattice field theories using fourier acceleration,Journal of Statistical Physics , 1073 (1986).[34] S. Tang and J. E. Hirsch, Peierls instability in the two-dimensional half-filled hubbard model, Phys. Rev. B ,9546 (1988).[35] S. Mazumdar, Valence-bond approach to two-dimensional broken symmetries: Application to la cuo ,Phys. Rev. B , 7190 (1987).[36] M. Feldbacher, F. F. Assaad, F. H´ebert, and G. G. Ba-trouni, Coexistence of s -wave superconductivity and an-tiferromagnetism, Phys. Rev. Lett. , 056401 (2003).[37] A. Sandvik, Stochastic method for analytic continuationof quantum monte carlo data, Phys. Rev. B , 10287(1998).[38] K. S. D. Beach, Identifying the maximum entropymethod as a special limit of stochastic analytic contin-uation, eprint arXiv:cond-mat/0403055 (2004), cond-mat/0403055.[39] F. F. Assaad, Spin, charge, and single-particle spectralfunctions of the one-dimensional quarter filled holsteinmodel, Phys. Rev. B , 155124 (2008).[40] M. Ulybyshev, N. Kintscher, K. Kahl, and P. Buivi-dovich, Schur complement solver for quantum monte-carlo simulations of strongly interacting fermions, Com-puter Physics Communications , 118 (2019).[41] M. Hohenadler, F. Parisen Toldin, I. F. Herbut, and F. F.Assaad, Phase diagram of the kane-mele-coulomb model,Phys. Rev. B , 085146 (2014).[42] T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, andM. P. A. Fisher, Quantum criticality beyond the landau-ginzburg-wilson paradigm, Phys. Rev. B , 144407(2004).[43] T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, andM. P. A. Fisher, Deconfined quantum critical points, Sci-ence , 1490 (2004).[44] A. Tanaka and X. Hu, Many-body spin berry phasesemerging from the π -flux state: Competition between an-tiferromagnetism and the valence-bond-solid state, Phys.Rev. Lett. , 036402 (2005).[45] T. Senthil and M. P. A. Fisher, Competing orders, non-linear sigma models, and topological terms in quantummagnets, Phys. Rev. B , 064405 (2006).[46] Y. Liu, Z. Wang, T. Sato, M. Hohenadler, C. Wang,W. Guo, and F. F. Assaad, Superconductivity from thecondensation of topological defects in a quantum spin-hall insulator, Nature Communications10