Charge-Density Wave Order on a π -flux Square Lattice
CCharge-Density Wave Order on a π -flux Square Lattice Y.-X. Zhang, H.-M. Guo, and R. T. Scalettar Department of Physics, University of California, Davis, CA 95616,USA Department of Physics, Key Laboratory of Micro-nano Measurement-Manipulationand Physics (Ministry of Education), Beihang University, Beijing, 100191, China (Dated: May 28, 2020)The effect of electron-phonon coupling (EPC) on Dirac fermions has recently been explorednumerically on a honeycomb lattice, leading to precise quantitative values for the finite temperatureand quantum critical points. In this paper, we use the unbiased determinant Quantum Monte Carlo(DQMC) method to study the Holstein model on a half-filled staggered-flux square lattice, andcompare with the honeycomb lattice geometry, presenting results for a range of phonon frequencies0 . (cid:54) ω (cid:54) .
0. We find that the interactions give rise to charge-density wave (CDW) order, butonly above a finite coupling strength λ crit . The transition temperature is evaluated and presentedin a T c - λ phase diagram. An accompanying mean-field theory (MFT) calculation also predicts theexistence of quantum phase transition (QPT), but at a substantially smaller coupling strength. I. INTRODUCTION
The physics of massless Dirac points, as exhibited inthe band structure of the honeycomb lattice of graphene,has driven intense study[1–4]. The square latticewith π -flux per plaquette is an alternate tight-bindingHamiltonian which also contains Dirac points in its bandstructure. Initial investigations of the π -flux modelfocused on the non-interacting limit[5], but, as with thehoneycomb lattice, considerable subsequent effort hasgone into extending this understanding to incorporatethe effect of electron-electron interactions. Numericalsimulations of the Hubbard Hamiltonian with an on-siterepulsion U between spin up and spin down fermions,including Exact Diagonalization [6] and Quantum MonteCarlo (QMC)[7–14] revealed a quantum phase transitionat U c ∼ . t into a Mott antiferromagnetic (AF)phase in the chiral Heisenberg Gross-Neveu universalityclass. For a spinless fermion system with near-neighborinteractions a chiral Ising Gross-Neveu universality classis suggested[15]. These results have been contrasted withthose on a honeycomb lattice, which has a similar Diracpoint structure, though at a smaller critical interaction U c ∼ . t [11].In the case of the repulsive Hubbard Hamiltonian,there were two motivations for studying both thehoneycomb and the π -flux geometries. The first wasto verify that the quantum critical transitions to AForder as the on-site repulsion U increases share the sameuniversality class, that of the Gross-Neveu model. Thesecond was to confirm that an intermediate spin-liquid(SL) phase between the semi-metal and AF phases[16],which had been shown not to be present on a honeycomblattice[17], was also absent on the π -flux geometry.Studies of the SU(2) π -flux Hubbard model have alsobeen extended to SU(4), using projector QMC[18], andto staggered flux where ± π hopping phases alternateon the lattice[19]. In the former case, the semi-metalto AF order transition was shown to be replaced by a semi-metal to valence bond solid transition characterizedby breaking of a Z symmetry. In the latter work, anintermediate phase with power-law decaying spin-spincorrelations was suggested to exist between the semi-metal and AF.A largely open question is how this physics is affectedin the presence of electron-phonon rather than electron-electron interactions. A fundamental Hamiltonian,proposed by Holstein[20], includes an on-site couplingof electron density to the linear displacement of thephonon field. In the low density limit, extensivenumerical work has quantified polaron and bipolaronformation, in which electrons are “dressed’ by anaccompanying lattice distortion [21–28]. At sufficientlylarge coupling, electrons or pairs of electrons canbecome ‘self-trapped’ (localized). One of the mostessential features of the Holstein model is that thelattice distortion of one electron creates an energeticallyfavorable landscape for other electrons, so that there isan effective attraction mediated by the phonons. Athigher densities, collective phenomena such as Charge-Density Wave (CDW) phases, and superconductivity(SC) have been widely studied [24, 29–36]. CDW isespecially favored on bipartite lattices and at fillingswhich correspond to double occupation of one of the twosublattices. SC tends to occur when one dopes away fromthese commensurate fillings.Recent work on the Holstein model on the honeycomblattice suggested a quantum phase transition from semi-metal to gapped CDW order [37, 38] similar to theresults for the Hubbard Hamiltonian. However, a keydifference between the Hubbard and Holstein modelsis the absence of the SU(2) symmetry of the orderparameter in the latter case. Thus, while long-range AForder arising from electron-electron interaction occursonly at zero temperature in 2D, the CDW phasetransition induced by electron-phonon coupling can occurat finite temperature- the symmetry being broken is thatassociated with two discrete sub-lattices. For classical a r X i v : . [ c ond - m a t . s t r- e l ] M a y phonons ( ω = 0), the electron-phonon coupling becomesan on-site energy in the mean-field approximation. Inthe anti-adiabatic limit where phonon frequencies are setto infinity, the Holstein model maps onto the attractiveHubbard model.Here we extend the existing work on the effect ofEPC on Dirac fermions from the honeycomb geometryto the π -flux lattice. The π -flux state is realized bythreading half of a magnetic flux quantum througheach plaquette of a square lattice[39]. Recently it hasbeen experimentally realized in optical lattices usingRaman assisted hopping[40]. There are also theoreticalsuggestions that the π -flux lattice might be engineeredby the proximity of an Abrikosov lattice of vortices of atype-II superconductor, or via spontaneously generatinga π -flux by coupling fermions to a Z gauge theory in(2+1) dimensions[41]. The π -flux hopping configurationhas an additional interesting feature motivating ourcurrent work: it is the unique magnetic field value whichminimizes the ground state energy for non-interactingfermions at half-filled on a bipartite lattice. Indeed,Lieb has shown that this theorem is also true at finitetemperature, and furthermore holds in the presenceof Hubbard inteactions[42]. Here we consider thethermodynamics of the π -flux lattice with EPC.This paper is organized as follows: in the next section,we describe the Holstein model and the π -flux squarelattice. Section III presents, briefly, a mean-field theory(MFT) for the model. Section IV reviews our primarymethod, DQMC. Section V contains results from theDQMC simulations, detailing the nature of the CDWphase transition, both the finite temperature transitionat fixed EPC, and the QPT which occurs at T = 0 withvarying EPC. Section VI contains our conclusions. II. MODEL
The Holstein model [20] describes conduction electronslocally coupled to phonon degrees of freedom,ˆ H = − (cid:88) (cid:104) i , j (cid:105) ,σ (cid:0) t i , j ˆ d † i σ ˆ d j σ + h . c . (cid:1) − µ (cid:88) i ,σ ˆ n i ,σ + 12 M (cid:88) i ˆ P i + ω (cid:88) i ˆ X i + λ (cid:88) i ,σ ˆ n i ,σ ˆ X i . (1)The sums on i and σ run over all lattice sites and spins σ = ↑ , ↓ . (cid:104) i , j (cid:105) denotes nearest neighbors. ˆ d † i σ and ˆ d i σ are creation and annihilation operators of electrons withspin σ on a given site i ; ˆ n i ,σ = ˆ d † i σ ˆ d i σ is the numberoperator. The first term of Eq. (1) corresponds to thehopping of electrons K el , with chemical potential µ . Thenext line of the Hamiltonian describes optical phonons,local quantum harmonic oscillators of frequency ω andphonon position and momentum operators, ˆ X i and ˆ P i respectively. The phonons are dispersionless since there A B FIG. 1. π -flux phase on a 6 × t (cid:48) = − t , as opposite to black lines withhopping t . Arrows represent the basis vectors. are no terms connecting ˆ X i on different sites of thelattice. The phonon mass M is set to unity. The electron-phonon coupling is included in the last term. We sethopping | t i , j | = t = 1 as the energy scale and focus onhalf-filling, ( (cid:104) ˆ n (cid:105) = 1), which can be achieved by setting µ = − λ /ω . It is useful to present results in terms of thedimensionless coupling λ D = λ / ( ω W ) which representsthe ratio of the effective electron-electron interactionobtained after integrating out the phonon degrees offreedom, and W is the kinetic energy bandwidth.The two dimensional π -flux phase on a square latticeis schematically shown in Fig. 1. All hopping in the x direction are t , while half of the hoppings along the y -direction are set to t (cid:48) = t e iπ = − t , where thephase π in the hopping amplitude arises from the Peierlsprescription for the vector potential of the magneticfield. As a consequence, an electron hopping on acontour around each plaquette picks up a total phase π , corresponding to one half of a magnetic flux quantumΦ = hc/e per plaquette. The lattice is bipartite, withtwo sublattices A and B . Each unit cell consists of twosites. In reciprocal space, with the reduced Brillouinzone ( | k x | ≤ π, | k y | ≤ | k x | ), the non-interacting part ofHamiltonian Eq.(1) can be written as,ˆ H = (cid:88) k σ ˆ ψ † k σ H ( k ) ˆ ψ k σ , (2)where ˆ ψ k σ = (cid:0) ˆ d Aσ ˆ d Bσ (cid:1) T , (3)and the noninteracting Hamiltonian matrix H ( k ) = (cid:18) t cos k x + 2 i t sin k y t cos k x − i t sin k y (cid:19) . (4) k x - /2 0 /2 k y - /2 0 /2 E k / t -3-2-10123 | E k |/ t FIG. 2. The dispersion relation E k for π -flux phase on asquare lattice. There are two Dirac points at ( k x , k y ) =( ± π/ , π -flux model is W = 4 √ t . -3 -2 -1 0 1 2 3 E k / t D O S -fluxhoneycomb FIG. 3. The density of states for the π -flux phase squarelattice and the honeycomb lattice. The bandwidths are nearlyidentical, but the honeycomb lattice has a substantially largerslope of the linear increase of the DOS. The energy spectrum E k = ± t (cid:113) cos k x + sin k y describes a semi-metal with two inequivalent Dirac pointsat K ± = ( ± π/ , E k = 0,as shown in Fig. 3. The bandwidth of the π -flux phase is W = 4 √ t . In Fig. 3 the DOS of the honeycomb latticeis shown for comparison. The Dirac Fermi velocity is v F = 2 t (1 . t ) for the π -flux (honeycomb) lattice. Nearthe Dirac point, the DOS ρ ( ω ) ∼ | ω | /v F , and the π -fluxmodel has a smaller slope. III. MEAN-FIELD THEORY
In this section, we present a mean-field theoryapproach to solve the Holstein model. Semi-metal tosuperfluid transitions have previously been investigatedwith MFT in 2D and 3D [43, 44]. Here we focus onthe semi-metal to CDW transition. In the mean-fieldapproximation, the phonon displacement at site i isreplaced by its average value, modulated by a term whichhas opposite sign on the two sublattices, (cid:104) X i (cid:105) = X ± X mf ( − i . (5)Here X = − λ/ω is the “equilibrium position” at half-filling and X mf is the mean-field order parameter. When X mf = 0, phonons on all sites have the same averagedisplacement, indicating the system remains in semi-metal phase, whereas when X mf (cid:54) = 0, the last termin the Hamiltonian Eq. (1), i.e., λ (cid:80) i ,σ ˆ n i ,σ ˆ X i , becomesan on-site staggered potential, which corresponds to theCDW phase. The phonon kinetic energy term is zero asa result of the static field. The resulting static mean-field Hamiltonian is quadratic in the fermion operators.Diagonalizing gives energy eigenvalues (cid:15) n ( X mf ). The freeenergy F can then be directly obtained by, F ( β, X mf ) = − β (cid:88) n ln(1 + e − β(cid:15) n )+ N ω X + X ) , (6)Minimizing the free energy with respect to X mf (orequivalently, a self-consistent calculation) will determinethe order parameter. X mf is found to be zero at hightemperatures: the energy cost of the second term in Eq. 6exceeds the energy decrease in the first term associatedwith opening of a gap in the spectrum (cid:15) n . X mf becomesnonzero below a critical temperature T c . T c for the π -flux lattice is shown in Fig. 4, alongwith the result of analogous MFT calculations for thehoneycomb and (zero flux) square geometries. The latticesize L = 180 is chosen for all three models. This issufficiently large so that finite size effects are smaller thanthe statistical sampling error bars. At zero temperature,the CDW order exhibits a critical EPC for the π -flux andthe honeycomb lattices. This QCP arises from the Diracfermion dispersion, which has a vanishing DOS at theFermi energy. The honeycomb lattice QCP has a smallercritical value. However, when measured in units of theFermi velocity, the ratios λ D, crit /v F = 0 .
13 and 0 . π -flux geometriesrespectively. We will see this is also the case for theexact DQMC calculations. For the square lattice, on theother hand, the DOS has a Van-Hove singularity at theFermi energy, and the CDW develops at arbitrarily smallcoupling strength. D T C D W / t square-fluxhoneycomb FIG. 4. MFT T c for CDW phase transition as a functionof dimensionless coupling λ D for the square lattice withno magnetic flux, the π -flux phase square lattice, and thehoneycomb lattice. For the geometries with a Dirac spectrumMFT captures the existence for a QCP, a critical value of λ D below which there is no CDW order even at T = 0, and theabsence of a QCP for the conventional square lattice. Another feature of the MFT phase diagram is that,as the coupling increases, T c increases monotonically.This is in contrast to the exact DQMC results, where T c decreases at large coupling strengths (Fig. 13). Asimilar failure of MFT is well known for the HubbardHamiltonian where the formation of AF ordering isrelated to two factors: the local moment m z = ( n ↑ − n ↓ ) = 1 − (cid:104) n ↑ n ↓ (cid:105) and the exchange coupling J ∼ t /U . The double occupancy (cid:104) n ↑ n ↓ (cid:105) is suppressed bythe interaction, resulting in the growth of the localmoment. Thus upon cooling, the Hubbard model hastwo characteristic temperatures: the temperature of localmoment formation, which increases monotonically with U , and further the AF ordering scale, which falls as J . Since the interaction is simply decoupled locally andthe exchange coupling is not addressed, within MFT theformation of the local moments, and their ordering, occursimultaneously. MFT thus predicts a monotonicallyincreasing T c with U . IV. DQMC METHODOLOGY
We next describe the DQMC method[45, 46]. Inevaluating the partition function Z , the inversetemperature β is discretized as β = L τ ∆ τ , and completesets of phonon position eigenstates are introducedbetween each e − ∆ τ ˆ H . The phonon coordinates acquirean “imaginary time” index, converting the 2-dimensionalquantum system to a (2+1) dimensional classicalproblem. After tracing out the fermion degrees of D | K E e l |/ t =10/ t =8/ t =6/ t D D FIG. 5. Left: The magnitude of electron kinetic energy |K el | as a function of EPC strength λ D . Simulations are performedon a L = 10 lattice at inverse temperatures β = 6 /t, /t, /t and fixed ω = 1 . t . Right: Double occupancy D as a functionof EPC strength λ D . freedom, which appear only quadratically in the HolsteinHamiltonian, the partition function becomes Z = (cid:90) D x i ,l e −S ph [det M ( x i ,l )] , (7)where the “phonon action” is S ph = ∆ τ (cid:34) ω (cid:88) i x i ,l + 12 M (cid:88) i (cid:18) x i ,l +1 − x i ,l ∆ τ (cid:19) (cid:35) . (8)Because the spin up and spin down fermions havean identical coupling to the phonon field, the fermiondeterminants which result from the trace are the same,and the determinant is squared in Eq. 7. Thus thereis no fermion sign problem[47]. We use ∆ τ = 0 . /t ,small enough so that Trotter errors associated with thediscretization of β are of the same order of magnitudeas the statistical uncertainty from the Monte Carlosampling. V. DQMC RESULTSDouble occupancy and Kinetic Energy
We first show data for several local observables, theelectron kinetic energy |K el | = | (cid:80) (cid:104) i , j (cid:105) ,σ (cid:0) t i , j ˆ d † i σ ˆ d j σ +h . c . (cid:1) | and double occupancy D = (cid:104) n i ↑ n i ↓ (cid:105) . For a tight-binding model on a bipartite lattice at half-filling, Liebhas shown that the energy-minimizing magnetic flux is π per plaquette, both in for noninteracting fermions andin the presence of a Hubbard U [42]. Here we show |K el | for the Holstein model, a case not hitherto considered.Figure 5 shows |K el | (left panel) and D (right panel)as functions of the dimensionless EPC λ D for β = S C D W =0.1 t (a) D =0.177 D =0.398 D =0.707 D =1.105 D =1.591 =0.5 t (b) D =0.346 D =0.453 D =0.573 D =0.707 D =0.856 D =1.018 t S C D W =1.0 t (c) D =0.346 D =0.453 D =0.573 D =0.856 D =1.195 D =1.591 t =2.0 t (d) D =0.453 D =0.573 D =0.638 D =0.78 D =0.935 D =1.018 FIG. 6. The CDW structure factor of the π -flux phaseHolstein model as a function of inverse temperature β . Thephonon frequencies ω are (a), 0 . t ; (b), 0 . t ; (c), 1 . t ; (d),2 . t in the four panels. The lattice size L = 6. /t, /t, /t . There is little temperature dependencefor these local quantities. The magnitude of the kineticenergy |K el | decreases as λ D grows, reflecting the graduallocalization of the dressed electrons (“polarons”).At the same time, the double occupancy D evolvesfrom its noninteracting value D = (cid:104) n i ↑ n i ↓ (cid:105) = (cid:104) n i ↑ (cid:105) (cid:104) n i ↓ (cid:105) = 1 / D = 1 / λ D . In the strong coupling regime, we expect robust pairformation, so that half of the lattice sites will be emptyand half will be doubly occupied.The evolution of D and |K el | have largest slope at λ D ∼ .
42 which, as will be seen, coincides with the location ofthe QCP between the semi-metal and CDW phases.
Existence of Long-Range CDW Order
The structure factor S ( Q ) is the Fourier transform ofthe real-space spin-spin correlation function c ( r ), S ( Q ) = (cid:88) r e i Q · r c ( r ) ,c ( r ) = (cid:10) (cid:0) n i ↑ + n i ↓ (cid:1)(cid:0) n i + r ↑ + n i + r ↓ (cid:1) (cid:11) , (9)and characterizes the charge ordering. In a disorderedphase c ( r ) is short-ranged and S ( Q ) is independent oflattice size. In an ordered phase, c ( r ) remains largeout to long distances, and the structure factor will beproportional to the number of sites, at the appropriateordering wave vector Q . At half-filling S ( Q ) is largest at Q = ( π, π ). We define S cdw ≡ S ( π, π ). Figure 6 displays S cdw as a function of inverse temperature β at different S C D W = 4/ t = 5/ t = 6/ t = 8/ t = 10/ t / t FIG. 7. S cdw (a) as a function of λ at fixed ω = 1 . t ; and(b) as a function of ω at fixed λ =3.0, at different inversetemperatures β . Lattice size L = 6 is used in this figure. phonon frequencies ω and coupling strengths λ D . Thelinear lattice size L = 6. At fixed ω and strong coupling, S cdw grows as temperature is lowered, and saturates to S cdw ∼ N, indicating the development of long-rangeorder (LRO), i.e. the phase transition into CDW phase.Note that β = 10 /t is always in the plateau region,suggesting the correlation length has become larger thanthe lattice size, and the ground state has been reached. Inthe following, we use β = 10 /t to represent the propertiesat T → λ D is decreased sufficiently, S cdw eventually shows no signal of LRO even at large β ,providing an indication that there is a QCP, with CDWorder only occurring above a finite λ D value. Figure 6also suggests that the critical temperature T c is non-monotonic with increasing λ D . The values of β atwhich S cdw grows first shift downward, but then becomelarger again. This non-monotonicity agrees with previousstudies of Dirac fermions on the honeycomb lattice [37,38]. We can estimate the maximum T c to occur at λ D ≈ . , . , .
86 and 0 .
78 for ω = 0 . t, . t, . t, . t respectively. In the anti-adiabatic limit ω → ∞ , theHolstein model maps onto the attractive Hubbard model,and T c = 0 owing to the degeneracy of CDW andsuperconducting correlations[29]. (The order parameterhas a continuous symmetry.) A recent study[48] hasshown that ω (cid:38) t is required to achieve the − U Hubbard model limit, a surprisingly large value.Figure 7(a) shows S cdw as a function of λ at fixed ω =1 . t . At the highest temperature shown, β = 4 /t , S cdw reaches maximum at intermediate coupling λ ∼ .
0, thendecreases as λ gets larger. The region for which S cdw islarge is a measure of the range of λ for which the CDWordering temperature T c exceeds β − . As β increases,this range is enlarged. Figure 7(b) is an analogous plotof S cdw as a function of ω at fixed λ = 3 .
0. The twoplots appear as mirror images of each other since the D S C D W = 5/ t = 3.0 = 1.0 t D = 8/ t = 3.0 = 1.0 t FIG. 8. Comparison of the evolution of S cdw with couplingstrength by changing λ or changing ω . Data are taken fromFig. 7(a,b), for β = 5 /t (left) and β = 8 /t (right). Thedifference is negligible at λ D > . . < λ D < . D S C D W (a) = 6/ t = 8/ t = 10/ t D t = 16/ t = 20/ t FIG. 9. S cdw as a function of λ D for π -flux phase square lattice(left) and honeycomb model (right). The lattice size L = 6 isused for both geometries. λ D is varied by changing λ at fixed ω = 1 . t . S cdw does not change for the lowest temperatures,indicating that the ground state has been reached for thisfinite lattice size. dimensionless EPC λ D = λ / ( ω W ) increases with λ ,but decreases with ω .It is interesting to ascertain the extent to which thephysics of the Holstein Hamiltonian is determined by λ and ω separately, versus only the combination λ D .Figure 8 addresses this issue by replotting the data ofFigs. 7(a,b) as a function of λ D for two values of theinverse temperature. For λ D (cid:38) .
8, the data collapsewell, whereas at small λ D S cdw can vary by as much asa factor of two even though λ D is identical. It is likelythat this sensitivity to the individual values of λ and ω is associated with proximity to the QCP. / t FIG. 10. Heat map of the ground state values of S cdw in the( λ, ω ) plane. We compare the semi-metal to CDW transition withincreasing λ D for the π -flux phase and honeycomblattices in Fig. 9. These data are at lower temperaturesthan those of Fig. 8, so that the ground state values of S cdw have been reached for the system sizes shown. Ground State in the ( λ, ω ) Plane Figure 10 provides another perspective on thedependence of the CDW order on λ and ω individually,by giving a heat map of S cdw in the ( λ, ω ) plane at lowtemperature. The bright yellow in upper-left indicatesa strong CDW phase, whereas the dark purple region inlower-right indicates the Dirac semi-metal phase. Thephase boundary is roughly linear, as would be expectedif only the combination λ D = λ / ( ω W ) is relevant. Wenote, however, that this statement is only qualitativelytrue. The more precise line graphs of Fig. 8 indicate thatalong the line λ = (cid:112) λ D, crit W ω ∼ . ω , the separatevalues of λ and ω are relevant. Finite Size Scaling: Finite T Transition
A quantitative determination of the finite temperatureand quantum critical points can be done with finitesize scaling (FSS). Figure 11 gives both raw and scaleddata for S cdw for different lattice sizes L = 4 , , , λ = 2 . ω = 1 . t as a function of β . Unscaleddata are in panel (a): S cdw is small and L -independentat small β (high T ) where c ( r ) is short ranged. Onthe other hand, S cdw is proportional to N = L at S C D W (a) L=4L=6L=8L=103.0 3.5 4.0 4.5 5.00.00.51.0 S C D W / L / (b) -10 -5 0 5 ( c ) L S C D W / L / (c) FIG. 11. (a) The CDW structure factor S cdw as a function of β for several lattice sizes. (b) The scaled CDW structurefactor S cdw / L γ/ν as a function of β using Ising criticalexponents γ = 7 / ν = 1, showing a crossing of differentL at β c = 3 . /t . (c) S CDW / L γ/ν versus ( β − β c )L, giving abest data collapse at β c = 3 . /t . Here the parameters are λ = 2 . ω = 1 . t . large β (low T ), reflecting the long-range CDW orderin c ( r ). Panel (b) shows a data crossing for different L occurs when S cdw /L γ/ν is plotted versus β . A universalcrossing is seen at β t ∼ . ± .
02, giving a precisedetermination of critical temperature T c . The 2D Isingcritical exponents γ = 7 / ν = 1 were used in thisanalysis, since the CDW phase transition breaks a similardiscrete symmetry. Panel (c) shows a full data collapsewhen the β axis is also appropriately scaled by L /ν . Thebest collapse occurs at β c = 3 . /t , consistent with theresult from the data crossing.In the region immeditely above the QCP, the DQMCvalues for T c are roughly five times lower than thoseobtained in MFT, and, indeed, the MFT over-estimationof T c can be made arbitrarily large at strong coupling.This reflects both the relatively low dimensionality ( d =2) and the fact that MFT fails to distinguish moment-forming and moment-ordering temperature scales. Quantum Phase Transition
Analysis of the renormalization group invariant Bindercumulant[49], B = 32 (cid:18) − < S >< S cdw > (cid:19) , (10)can be used to locate the quantum critical point precisely.Only lattice sizes L = 4 n where n is an integer can beused, for other L the Dirac points are not one of theallowed k values and finite size effects are much more D B L=4L=8L=12 D , c r i t ( L ) FIG. 12. Main panel: Binder cumulant as a function of EPCstrength λ D for three lattice sizes. Inverse temperature is β = 2 L and ω is fixed at ω = 1 . t . Inset: Extrapolation ofthe crossings for pairs of sizes as a function of 1 /L to get theQCP in the thermodynamic limit. D T C D W / t Semi-metal CDW honeycomb-flux
FIG. 13. Critical temperature T c for CDW phase transition,obtained from DQMC for both π -flux phase square lattice(blue line) and the honeycomb lattice (red line), in a rangeof coupling strength. λ D is varied by changing λ at fixed ω = 1 . t for both models. Quantum critical point isdetermined using Binder cumulant analysis (discussed below).Data for the honeycomb lattice are taken from [37] Error barsare smaller than symbol size for π -flux data. significant. As exhibited in Fig. 12, for L = 4 , B exhibits a set of crossings in a range about λ D ≈ . /L , as shown in the inset of Fig. 12,gives λ D, crit = 0 . ± . Phase Diagram
Location of the finite temperature phase boundary,Fig. 11, and the QCP, Fig. 12, can be combined into thephase diagram of Fig. 13. Results for the π -flux geometry(blue circles) are put in better context by compared withthose of the honeycomb lattice (red triangles). Data wereobtained at fixed ω = 1 . t . In both geometries, phasetransitions into CDW order happen only above a finite λ D, crit . Beyond λ D, crit , T c rises rapidly to its maximalvalue before decaying. For π -flux model, T c reaches amaximum T c, max /t ∼ .
26 at λ D ∼ .
7, whereas for thehoneycomb lattice T c reaches its maximum T c, max /t ∼ .
20 at λ D ∼ .
5. Similarly λ D,crit for π -flux is largerthan that of the honeycomb lattice, as λ D,crit = 0 . .
27 respectively. When measured in terms ofthe relative Fermi velocities v F = 2 t, . t for the π -flux and honeycomb respectively, these values becomevery similar: λ D, crit /v F = 0 .
21 and 0 .
18 for π -flux andhoneycomb; T c, max /v F = 0 .
13 and 0 . VI. CONCLUSIONS
This paper has determined the quantitative phasediagram for Dirac fermions interacting with local phononmodes on the π -flux lattice. A key feature, shared withthe honeycomb geometry, is the presence of a quantumcritical point λ D, crit below which the system remainsa semi-metal down to T = 0. The values of T c and λ D, crit for the two cases, when normalized to the Fermivelocities, agree to within roughly 10%.We have also considered the question of whetherthe properties of the model can be described in termsof the single ratio λ /ω . We find that qualitativelythis is indeed the case, but that, quantititively, thecharge structure factor can depend significantly onthe individual values of EPC and phonon frequency,especially in the vicinity of the QCP. However this morecomplex behavior is masked by the fact that T c rises sorapidly with λ in that region. In investigating this issuewe have studied substantially smaller values of ω thanhave typically been investigated in QMC treatments ofthe Holstein Hamiltonian.Acknowledgments: The work of Y.-X.Z. and R.T.S. wassupported by the grant DESC0014671 funded by theU.S. Department of Energy, Office of Science. H.G. wassupported by NSFC grant No. 11774019. The authorswould like to thank B. Cohen-Stead and W.-T. Chiu foruseful conversations. [1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, “The electronic propertiesof graphene,” Rev. Mod. Phys. , 109–162 (2009).[2] A.K. Geim, “Graphene: Status and prospects,” Science , 1530 (2009).[3] Wonbong Choi, Indranil Lahiri, RaghunandanSeelaboyina, and Yong Soo Kang, “Synthesis ofgraphene and its applications: A review,” Crit. Rev. inSolid State and Mat. Sci. , 52–71 (2010).[4] Konstantin S Novoselov, VI Fal, L Colombo, PR Gellert,MG Schwab, K Kim, et al. , “A roadmap for graphene,”nature , 192–200 (2012).[5] A Brooks Harris, Tom C Lubensky, and Eugene J Mele,“Flux phases in two-dimensional tight-binding models,”Phys. Rev. B , 2631 (1989).[6] Yongfei Jia, Huaiming Guo, Ziyu Chen, Shun-QingShen, and Shiping Feng, “Effect of interactions on two-dimensional dirac fermions,” Phys. Rev. B , 075101(2013).[7] Y. Otsuka and Y. Hatsugai, “Mott transition in the two-dimensional flux phase,” Phys. Rev. B , 073101 (2002).[8] Yuichi Otsuka, Seiji Yunoki, and Sandro Sorella,“Mott transition in the 2d hubbard model with π -flux,” in Proceedings of the International Conference onStrongly Correlated Electron Systems (SCES2013) (2014)p. 013021.[9] Zi-Xiang Li, Yi-Fan Jiang, and Hong Yao, “Fermion-sign-free majarana-quantum-monte-carlo studies ofquantum critical phenomena of dirac fermions in twodimensions,” New J. of Phys. , 085003 (2015).[10] Francesco Parisen Toldin, Martin Hohenadler, Fakher FAssaad, and Igor F Herbut, “Fermionic quantumcriticality in honeycomb and π -flux hubbard models:Finite-size scaling of renormalization-group-invariantobservables from quantum monte carlo,” Phys. Rev. B , 165108 (2015).[11] Yuichi Otsuka, Seiji Yunoki, and Sandro Sorella,“Universal quantum criticality in the metal-insulatortransition of two-dimensional interacting diracelectrons,” Phys. Rev. X , 011029 (2016).[12] Thomas C. Lang and Andreas M. L¨auchli, “Quantummonte carlo simulation of the chiral heisenberg gross-neveu-yukawa phase transition with a single dirac cone,”Phys. Rev. Lett. , 137602 (2019).[13] H-M Guo, Lei Wang, and RT Scalettar, “Quantum phasetransitions of multispecies dirac fermions,” Phys. Rev. B , 235152 (2018).[14] Huaiming Guo, Ehsan Khatami, Yao Wang, Thomas PDevereaux, Rajiv RP Singh, and Richard T Scalettar,“Unconventional pairing symmetry of interacting diracfermions on a π -flux lattice,” Phys. Rev. B , 155146(2018).[15] Lei Wang, Philippe Corboz, and Matthias Troyer,“Fermionic quantum critical point of spinless fermionson a honeycomb lattice,” New J. of Phys. , 103008(2014).[16] Z.Y. Meng, T.C. Lang, S. Wessel, F.F. Assaad, andA. Muramatsu, “Quantum spin liquid emerging in two-dimensional correlated dirac fermions,” Nature , 847–851 (2010).[17] Y. Otsuka, S. Yunoki, and S. Sorella, “Quantum monte carlo study of the half-filled hubbard model on thehoneycomb lattice,” J. of Phys.: Conf. Series , 012045(2013).[18] Zhichao Zhou, Congjun Wu, and Yu Wang, “Motttransition in the π -flux su (4) hubbard model on a squarelattice,” Phys. Rev. B , 195122 (2018).[19] Chia-Chen Chang and Richard T. Scalettar, “Quantumdisordered phase near the mott transition in thestaggered-flux hubbard model on a square lattice,” Phys.Rev. Lett. , 026404 (2012).[20] T Holstein, “Studies of polaron motion: Part i. themolecular-crystal model,” Annals of Physics , 325 – 342(1959).[21] P.E. Kornilovitch, “Continuous-time quantum montecarlo algorithm for the lattice polaron,” Phys. Rev. Lett. , 5382 (1998).[22] P. E. Kornilovitch, “Ground-state dispersion and densityof states from path-integral monte carlo: Application tothe lattice polaron,” Phys. Rev. B , 3237–3243 (1999).[23] A. S. Alexandrov, “Polaron dynamics and bipolaroncondensation in cuprates,” Phys. Rev. B , 12315–12327(2000).[24] M. Hohenadler, H. G. Evertz, and W. von der Linden,“Quantum monte carlo and variational approaches to theholstein model,” Phys. Rev. B , 024301 (2004).[25] Li-Chung Ku, S. A. Trugman, and J. Bonˇca,“Dimensionality effects on the holstein polaron,” Phys.Rev. B , 174306 (2002).[26] P. E. Spencer, J. H. Samson, P. E. Kornilovitch, andA. S. Alexandrov, “Effect of electron-phonon interactionrange on lattice polaron dynamics: A continuous-timequantum monte carlo study,” Phys. Rev. B , 184310(2005).[27] A. Macridin, G. A. Sawatzky, and Mark Jarrell, “Two-dimensional hubbard-holstein bipolaron,” Phys. Rev. B , 245111 (2004).[28] Aldo H. Romero, David W. Brown, and KatjaLindenberg, “Effects of dimensionality and anisotropyon the holstein polaron,” Phys. Rev. B , 14080–14091(1999).[29] R. T. Scalettar, N. E. Bickers, and D. J. Scalapino,“Competition of pairing and peierls˘charge-density-wave correlations in a two-dimensional electron-phononmodel,” Phys. Rev. B , 197–200 (1989).[30] F Marsiglio, “Pairing and charge-density-wavecorrelations in the holstein model at half-filling,”Physical Review B , 2416 (1990).[31] M. Vekic, R.M. Noack, and S.R. White, “Charge-densitywaves versus superconductivity in the holstein modelwith next-nearest-neighbor hopping,” Phys. Rev. B ,271 (1992).[32] Parhat Niyaz, J. E. Gubernatis, R. T. Scalettar, andC. Y. Fong, “Charge-density-wave-gap formation in thetwo-dimensional holstein model at half-filling,” Phys.Rev. B , 16011–16022 (1993).[33] M. Veki´c and S. R. White, “Gap formation in the densityof states for the holstein model,” Phys. Rev. B , 7643–7650 (1993).[34] JK Freericks, M Jarrell, and DJ Scalapino, “Holstein model in infinite dimensions,” Phys. Rev. B , 6302–6314 (1993).[35] H. Zheng and S. Y. Zhu, “Charge-density-wave andsuperconducting states in the holstein model on a squarelattice,” Phys. Rev. B , 3803–3815 (1997).[36] Eric Jeckelmann, Chunli Zhang, and Steven R. White,“Metal-insulator transition in the one-dimensionalholstein model at half filling,” Phys. Rev. B , 7950–7955 (1999).[37] Y.-X. Zhang, W.-T. Chiu, N. C. Costa, G. G. Batrouni,and R. T. Scalettar, “Charge order in the holstein modelon a honeycomb lattice,” Phys. Rev. Lett. , 077602(2019).[38] Chuang Chen, Xiao Yan Xu, Zi Yang Meng, and MartinHohenadler, “Charge-density-wave transitions of diracfermions coupled to phonons,” Phys. Rev. Lett. ,077601 (2019).[39] Ian Affleck and J Brad Marston, “Large-n limit of theheisenberg-hubbard model: Implications for high-t csuperconductors,” Physical Review B , 3774 (1988).[40] Monika Aidelsburger, Marcos Atala, SylvainNascimbene, Stefan Trotzky, Y-A Chen, and ImmanuelBloch, “Experimental realization of strong effectivemagnetic fields in an optical lattice,” Physical reviewletters , 255301 (2011).[41] Snir Gazit, Mohit Randeria, and Ashvin Vishwanath,“Emergent dirac fermions and broken symmetries inconfined and deconfined phases of z 2 gauge theories,”Nature Physics , 484–490 (2017).[42] Elliott H. Lieb, “Flux phase of the half-filled band,” Phys.Rev. Lett. , 2158–2161 (1994).[43] Gabriel Mazzucchi, Luca Lepori, and AndreaTrombettoni, “Semimetal–superfluid quantum phasetransitions in 2d and 3d lattices with dirac points,” J.of Phys,. B: Atomic, Mol. and Opt. Phys. , 134014(2013).[44] Ya-Jie Wu, Jiang Zhou, and Su-Peng Kou, “Stronglyfluctuating fermionic superfluid in the attractive π -fluxhubbard model,” Phys. Rev. A , 013619 (2014).[45] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar,“Monte carlo calculations of coupled boson-fermionsystems. i,” Phys. Rev. D , 2278–2286 (1981).[46] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh,J. E. Gubernatis, and R. T. Scalettar, “Numerical studyof the two-dimensional hubbard model,” Phys. Rev. B , 506–516 (1989).[47] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White,D. J. Scalapino, and R. L. Sugar, “Sign problem in thenumerical simulation of many-electron systems,” Phys.Rev. B , 9301–9307 (1990).[48] Chunhan Feng, Huaiming Guo, and Richard T Scalettar,“Charge density waves on a half-filled decoratedhoneycomb lattice,” Physical Review B , 205103(2020).[49] Kurt Binder, “Finite size scaling analysis of ising modelblock distribution functions,” Zeitschrift f¨ur Physik BCondensed Matter43