@prefix vivo: .
@prefix edm: .
@prefix ns0: .
@prefix dcterms: .
@prefix skos: .
vivo:departmentOrSchool "Science, Faculty of"@en, "Physics and Astronomy, Department of"@en ;
edm:dataProvider "DSpace"@en ;
ns0:degreeCampus "UBCV"@en ;
dcterms:creator "Adolphs, Clemens"@en ;
dcterms:issued "2016-04-26T02:02:39"@en, "2016"@en ;
vivo:relatedDegree "Doctor of Philosophy - PhD"@en ;
ns0:degreeGrantor "University of British Columbia"@en ;
dcterms:description """In this thesis, we investigated a set of theoretical models frequently used in the field of solid state physics. These models describe coupling of charge carriers to bosonic modes such as phonons or magnons. In particular, the Holstein model describes coupling of charge carriers to dispersionless phonons, whereas the Emery model describes coupling of charge carriers to magnons in hole-doped antiferromagnets.
For the Holstein-like models, we studied how extending the model of the coupling beyond terms that are merely linear in the lattice distortion affects the ground state properties. Using appropriate extensions of the momentum average approximation, we could show that even small nonlinearities have a dramatic effect on the resulting quasi-particle's properties. We further investigated a particular type of nonlinear coupling, the double-well coupling model. After studying the properties of a single quasi-particle, we also showed that this system allows the formation of bound states between two charge carriers and a phonon cloud, the so-called bi-polaron. In contrast to the linear variation of the Holstein model, the resulting bi-polaron can be strongly bound yet lightweight. For the Emery model, we consider an experimentally relevant extension. The original model describes a single layer of CuO₂, relevant for the hole-doped cuprate superconductors. We consider recently synthesized layers of CuO, which can be viewed as two intercalated layers of CuO₂. The resulting system is similar to CuO₂, but different in important aspects. We use a variational method similar in spirit to MA but applicable to magnons instead of phonons to obtain the system's dispersion and compare it to that of the original CuO₂ layer. We observe a discrepancy between these dispersions that cannot be accounted for with a single-band model that is commonly used to model the CuO₂ dispersion. However, it has been a long-standing question whether or not this and other single-band models are appropriate for the description of cuprate physics. With our study of CuO, we demonstrated how a careful experimental analysis of this system can resolve that question."""@en ;
edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/57841?expand=metadata"@en ;
skos:note "Extensions Beyond Standard ModelsbyClemens AdolphsA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2016c© Clemens Adolphs 2016AbstractIn this thesis, we investigated a set of theoretical models frequently used inthe field of solid state physics. These models describe coupling of charge car-riers to bosonic modes such as phonons or magnons. In particular, the Hol-stein model describes coupling of charge carriers to dispersionless phonons,whereas the Emery model describes coupling of charge carriers to magnonsin hole-doped antiferromagnets.For the Holstein-like models, we studied how extending the model of thecoupling beyond terms that are merely linear in the lattice distortion affectsthe ground state properties. Using appropriate extensions of the momentumaverage approximation, we could show that even small nonlinearities havea dramatic effect on the resulting quasi-particle’s properties. We furtherinvestigated a particular type of nonlinear coupling, the double-well couplingmodel. After studying the properties of a single quasi-particle, we alsoshowed that this system allows the formation of bound states between twocharge carriers and a phonon cloud, the so-called bi-polaron. In contrast tothe linear variation of the Holstein model, the resulting bi-polaron can bestrongly bound yet lightweight.For the Emery model, we consider an experimentally relevant extension.The original model describes a single layer of CuO2, relevant for the hole-doped cuprate superconductors. We consider recently synthesized layersof CuO, which can be viewed as two intercalated layers of CuO2. Theresulting system is similar to CuO2, but different in important aspects. Weuse a variational method similar in spirit to MA but applicable to magnonsinstead of phonons to obtain the system’s dispersion and compare it tothat of the original CuO2 layer. We observe a discrepancy between thesedispersions that cannot be accounted for with a single-band model that isiiAbstractcommonly used to model the CuO2 dispersion. However, it has been along-standing question whether or not this and other single-band modelsare appropriate for the description of cuprate physics. With our study ofCuO, we demonstrated how a careful experimental analysis of this systemcan resolve that question.iiiPreface• A version of the work discussed in chapter 2 is published as “C.P.J.Adolphs and M. Berciu, Europhysics Letters 102, 47003 (2013)”. Itrelies on and advances further the techniques introduced previouslyby Professor M. Berciu published in [1] as well as techniques jointlydeveloped in our group by M. Moeller, A. Mukherjee, C.P.J. Adolphs,D.J.J. Marchand and M. Berciu in [2].• The work discussed in chapter 3 is published as “C.P.J. Adolphs andM. Berciu, Physical Review B 89, 035122 (2014). It is based on thesame techniques as cited above.• The work presented in chapter 4 is published as “C.P.J. Adolphs andM. Berciu, Physical Review B 90, 085149 (2014). It is based on thesame previous work as cited above.• A manuscript of the work discussed in 5 has been accepted for publi-cation by Physical Review Letters and is expected to be published bythe end of February 2016. It presents an extension of the work donein [3] and [4].I carried out all the necessary analytical and numerical work for all fourprojects. I wrote the first drafts for all the manuscripts. Prof. M. Berciuassisted in the preparation of the final drafts of these works. The draft forthe final project, currently under review for publication in Physical ReviewLetters, was also assisted by S. Moser and G.A.W. Sawatzky.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Model Hamiltonians . . . . . . . . . . . . . . . . . . . . . . . 21.2 The Holstein model for electron-phonon coupling . . . . . . . 31.3 The three-band Emery model for a cuprate layer . . . . . . . 211.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262 Nonlinear Holstein model . . . . . . . . . . . . . . . . . . . . 292.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.2.1 Summary of model assumptions . . . . . . . . . . . . 332.3 Formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 382.4.1 Atomic limit . . . . . . . . . . . . . . . . . . . . . . . 392.4.2 Negative nonlinear term . . . . . . . . . . . . . . . . 442.4.3 Small carrier concentrations . . . . . . . . . . . . . . 45vTable of Contents2.4.4 Other types of el-ph coupling . . . . . . . . . . . . . . 482.5 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . 483 Quadratic Holstein model . . . . . . . . . . . . . . . . . . . . 493.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.2.1 Summary of model assumptions . . . . . . . . . . . . 563.3 Formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.3.1 The even sector . . . . . . . . . . . . . . . . . . . . . 593.3.2 The odd sector . . . . . . . . . . . . . . . . . . . . . . 623.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.4.1 Atomic limit: t = 0 . . . . . . . . . . . . . . . . . . . 623.4.2 Finite hopping . . . . . . . . . . . . . . . . . . . . . . 673.5 Summary and discussions . . . . . . . . . . . . . . . . . . . . 744 Bipolarons in the quadratic Holstein model . . . . . . . . . 764.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.2.1 Summary of model assumptions . . . . . . . . . . . . 794.3 Formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.3.1 Momentum average approximation . . . . . . . . . . 814.3.2 Exact diagonalization . . . . . . . . . . . . . . . . . . 844.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.5 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . 915 Single-polaron dispersion in the insulating limit of tetrago-nal CuO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 945.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 945.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.2.1 Summary of model assumptions . . . . . . . . . . . . 1015.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . 1045.4.1 The ZRS Bloch state . . . . . . . . . . . . . . . . . . 1065.4.2 Low-energy effective model . . . . . . . . . . . . . . . 107viTable of Contents5.4.3 Comparison to ARPES . . . . . . . . . . . . . . . . . 1136 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166.1 Summary of this work . . . . . . . . . . . . . . . . . . . . . . 1166.2 Further developments . . . . . . . . . . . . . . . . . . . . . . 117Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120AppendicesA Various formal methods . . . . . . . . . . . . . . . . . . . . . . 136B Appendix for Nonlinear Holstein Model . . . . . . . . . . . 137B.1 Free propagator for the quadratic model . . . . . . . . . . . 137B.2 Equations of motion for quartic model . . . . . . . . . . . . . 138C Appendix for Quadratic Holstein Model . . . . . . . . . . . 142C.1 Details for the even-sector . . . . . . . . . . . . . . . . . . . 142C.1.1 Coupling matrices . . . . . . . . . . . . . . . . . . . . 142C.1.2 Manipulation of the EOMs . . . . . . . . . . . . . . . 143C.2 Details for the odd-sector . . . . . . . . . . . . . . . . . . . . 144C.2.1 Equations of Motion . . . . . . . . . . . . . . . . . . . 144C.2.2 Momentum space Green’s functions . . . . . . . . . . 146C.3 Quadratic e-ph coupling with g2 > 0 . . . . . . . . . . . . . . 147D Appendix for T-CuO . . . . . . . . . . . . . . . . . . . . . . . . 150D.1 Exact Diagonalization with Lanczos . . . . . . . . . . . . . . 150viiList of Tables1.1 Relevant parameters for the 3-band Emery model in units ofeV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.1 Some example values of the bipolaron binding energy andeffective mass. . . . . . . . . . . . . . . . . . . . . . . . . . . 89viiiList of Figures1.1 Sketch of the origin of the electron-phonon coupling term.Without an electron present, the nuclear potential Up(x) isgiven by the blue curve I. When an electron is present, thenuclear potential is instead given by the red curve II. . . . . 61.2 Sketch of the first two orders of the Dyson series representedby Feynman diagrams. Straight lines represent the free propa-gator of the charge carrying fermion. Squiggly lines representthe free propagator of the boson (the optical phonon in thiscase). A complete diagram would include labels at the ver-tices denoting the energies and momenta (or other relevantquantum numbers) of the particles. . . . . . . . . . . . . . . . 91.3 Spectral function of the one-dimensional Holstein model asobtained with MA. The effective coupling λ = −g2/2dtΩ isa dimensionless measure for the strength of el-ph interaction.The strong coupling regime occurs at around λ ∼ 1. . . . . . 191.4 Phase diagram of the high-temperature superconducting cupratesfor both electron and hole doping. The undoped compoundsare antiferromagnetic insulators which become superconduct-ing upon either hole or electron doping. Here, we focus on thehole-doped side of the phase diagram. Public domain figurefrom https://en.wikipedia.org/wiki/File:Cuphasediag.png. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22ixList of Figures1.5 (a) Structure of a typical high-temperature superconductingcuprate compound, YBCO. Layers of copper and oxygen arecontained within a host structure of rare-earth and transi-tion metal elements. Figure used under the GFDL, as ob-tained from https://commons.wikimedia.org/wiki/File:Ybco002.svg. (b) The most relevant orbitals for cupratephysics are the copper dx2−y2 orbital (round shapes) and theligand oxygen px and py orbitals (lobes), forming a unit cellwith three orbitals. . . . . . . . . . . . . . . . . . . . . . . . 231.6 E(k) along several cuts of the Brillouin zone for the three-band model in the variational approximation described in[3, 5] with (a) at most 2 magnons and (b) at most 3 magnons.Circles show exact diagonalization results from [6] for a clus-ter of 32 copper and 64 oxygen orbitals. Full lines show theresults from the variational approximation, neglecting spinfluctuations. The overall agreement is excellent and can befurther improved by including local spin fluctuations [5]. . . 272.1 (a) Sketch of a 1D chain of polar molecules; (b) The potentialof the pair with (II) or without (I) an extra charge carrier (fulllines) is approximated by a polynomial (thick dashed lines) . 322.2 GS quasiparticle weight (left panel) and GS average phononnumber (right panel) vs. ζ, in the quadratic (n = 2, lines)and quartic (n = 4, symbols) models, for various values of λand Ω = 0.5t, in one dimension. . . . . . . . . . . . . . . . . 402.3 (left) Zat, and (right) N(at)ph vs. ζ, for g1 =√2 and Ω = 0.5(full lines). Dashed lines show the mean-field estimates, whilethe dot-dashed lines show the results of fitting g˜/Ω˜ to exactlyreproduce the other quantity. See text for more details. . . . 42xList of Figures2.4 Estimate of the bipolaron phase diagram in 1D for Ω/t = 0.5and for different values of ζ, based on second order perturba-tion theory in t. In all four panels, the solid lines show thetransition from S0 (on-site) stable bipolarons to S1 (nearest-neighbor) stable bipolarons, while the dashed lines show theunbinding transition above which bound polarons are unsta-ble. Note that panels (c) and (d) have a significantly rescaledy-axis. See text for more details. . . . . . . . . . . . . . . . . 463.1 Sketch of the crystal structures discussed in this work: (a)1D chain, and (b) 2D plane, consisting of light atoms (filledcircles) intercalated between heavy atoms (empty circles). Inthe absence of carriers, the ionic potential of a light atomis a simple harmonic well. In the presence of a carrier, theionic potential of the light atom hosting it remains an evenfunction of its longitudinal displacement, so the linear e-phcoupling vanishes. In suitable conditions the effective ionicpotential becomes a double well (see text for more details). . 513.2 (a) Overlap between the undoped ground-states with andwithout anharmonic corrections, and (b) the average numberof phonons per site in the undoped system, due to anharmoniccorrections, as a function of θ/Ω. . . . . . . . . . . . . . . . 553.3 Polaron ground-state properties in the atomic limit, for sev-eral values of the g4: a) quasiparticle weight, and b) averagenumber of phonons in the phonon cloud. Other parametersare Ω = 0.5, t = 0. . . . . . . . . . . . . . . . . . . . . . . . . 643.4 Relative error in the ground state energy when computed inthe semiclassical approximation (see text for details). Thecoupling g2 < −Ω/4 is restricted to values for which there isa double-well potential. Other parameters are like in Fig. 3.3. 66xiList of Figures3.5 (color online) Polaron ground-state properties in one dimen-sion for various values of the quartic coupling term g4 as afunction of the quadratic coupling g2: a) quasiparticle weight,and b) average number of phonons in the phonon cloud. Otherparameters are t = 1, Ω = 0.5t. . . . . . . . . . . . . . . . . . 683.6 Polaron ground-state properties in two dimensions for variousvalues of the quartic coupling term g4 as a function of thequadratic coupling g2: a) quasiparticle weight, and b) averagenumber of phonons in the phonon cloud. Other parametersare t = 1, Ω = 0.5t. . . . . . . . . . . . . . . . . . . . . . . . 693.7 A(k, ω) in 1D, for g4 = 0.05, Ω = 0.5 and t = 1, for variousvalues of g2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 703.8 A(k, ω) in 2D, for g4 = 0.05, Ω = 0.5 and t = 1, for variousvalues of g2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 713.9 Real-space diagonal spectral function Aiii(ω) at g4 = 0.05for various values of (negative) g2 in one dimension and forΩ = 0.5. The y-axis has a logarithmic scale. The verticalbars indicate the position of Eeven0 + Ω. . . . . . . . . . . . . . 724.1 Sketch of the crystal structures discussed in this work: (a)1D chain, and (b) 2D plane, consisting of light atoms (filledcircles) intercalated between heavy atoms (empty circles). Inthe absence of carriers, the ionic potential of a light atomis a simple harmonic well. In the presence of a carrier, theionic potential of the light atom hosting it remains an evenfunction of its longitudinal displacement, so the linear e-phcoupling vanishes. In suitable conditions the effective ionicpotential becomes a double well. (Reproduced from the pre-vious chapter) . . . . . . . . . . . . . . . . . . . . . . . . . . 804.2 (a) Bipolaron ground-state energy, and (b) inverse effectivemass for t = 1,Ω = 0.5, g4 = 0.1, computed with ED (solidblack line) and MA (red dots). . . . . . . . . . . . . . . . . . 85xiiList of Figures4.3 Ground-state properties (total energy and inverse mass) ofthe S0 bipolaron and two independent polarons for t = 1,Ω =0.5 and g4 = 0.2, 0.1, and 0.05 for a), b), and c), respectively.For all panels, U = 0. . . . . . . . . . . . . . . . . . . . . . . 874.4 Ground-state properties (total energy and inverse mass) ofthe S0 bipolaron and two independent polarons for t = 1,Ω =2 and g4 = 0.2, 0.1, and 0.05 for a), b), and c), respectively.For all panels, U = 0. . . . . . . . . . . . . . . . . . . . . . . 884.5 Binding energy ∆ and effective-mass ratio mbp/2mp of thebipolarons for Ω = 0.5 and Ω = 2 at different values of g4. . 884.6 Ionic potential (above) and ionic ground-state wavefunction(below) in the single-well and double-well models. Solid linescorrespond to the situation without an additional carrier,dashed lines to the situation with an additional carrier. . . . 904.7 Bipolaron energy (left) and effective mass (right) as a functionof the Hubbard U , for t = 1,Ω = 0.4 and g4 = 0.1. Similarresults are found for other parameters if U is not large enoughto lead to bipolaron dissociation. . . . . . . . . . . . . . . . . 925.1 Structure of a layer of (a) T-CuO, and (b) CuO2. Full cir-cles are Cu, empty squares are O. The Cu 3dx2−y2 orbitalsare drawn at a few sites, with white/dark lobes showing ourchoice for positive/negative signs. The corresponding ligandO 2p orbitals are also indicated on neighboring O sites. TheT-CuO layer can be thought of as two intercalated CuO2 lay-ers sharing common O. The coppers of the two sublatticeshybridize with different O 2p orbitals. Panels (c) and (d)show the two degenerate ground-states of the undoped T-CuO layer. Different colors are used for the Cu spins on thetwo sublattices for better visibility. . . . . . . . . . . . . . . 98xiiiList of Figures5.2 Qp dispersion in units of Jdd for (a) nm = 1, (b) nm = 2,and (c) nm = 3 with full/dashed lines for T-CuO/CuO2. TheBrillouin zone for the magnetic order of Fig. 5.1(c) is shown inred in (b). The shaded area is the smaller BZ for CuO2. Thepoints marked by circles and empty/full squares are equiva-lent in CuO2 but not in T-CuO. (d) Hopping between twoadjacent ZRSs, and (e) between a ZRS (red) and one withx− y symmetry (blue). See text for more details. . . . . . . 1055.3 (a) Unit cell for CuO2, with two Cu spins and four ligand Oorbitals. We use the location i of the down-spin Cu as the ref-erence point. The white/shaded areas indicate our choice forpositive/negative lobs. (b) Magnetic Brillouin zone (shadedregion) vs. full Brillouin zone (unshaded). (c) ZRS between ahole occupying the linear combination of ligand orbitals withx2− y2 symmetry, and the spin of the central Cu. The Blochstate is obtained from its translations on the correspondingmagnetic sublattice. . . . . . . . . . . . . . . . . . . . . . . . 1065.4 Comparison between the one-magnon numerical dispersionand the effective calculations involving both the low-energyand first excited state quasiparticles. The results are in quali-tative agreement, in particular the shift of the minimum alongkx = ky. The dashed lines are a guide to the eye to demon-strate that the peak along Γ − M gets shifted away fromk = (pi/2, pi/2) yet remains there along X −X ′. . . . . . . . 111xivList of Figures5.5 The first two panels show the Brillouin zones correspondingto the two possible magnetic orders of Fig. 5.1(c) and (d).The symbols mark the deep minima (full squares), displacedshallower minima (circles) and very shallow minima (emptysquares). The two figures are related by a C4 rotation. Thethird panel shows the average of these two patterns, where ateach special point the deepest local minimum was selected.This pattern has a restored C4 symmetry and a Brillouinzone (blue line) corresponding to one Cu site/unit cell. Thepattern of deep/shallower minima is like that found experi-mentally. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1125.6 Quasiparticle dispersion for the two orientations of the mag-netic background, along the contour considered in Ref. [20].Top panel shows the results in the hole picture used through-out this work, whereas the low panel shows the same resultsin the electron picture relevant for experiments, with furthertuning of the parameters. Also, we note that the slight dis-placements of the shallower minima have not been observedexperimentally. This may be because of the very broad widthsof the quasiparticle peaks and the loss of spectral weight onone side of these points, which may mask them. In “un-twinned” samples the scattering rate should be lower, whichmay make the observation of these displacements towards Γmore easily visible. Of course, in that case the lack of C4symmetry and existence of really shallow minima should alsobecome visible. . . . . . . . . . . . . . . . . . . . . . . . . . . 115C.1 a) Quasiparticle weight Z, and (b) average number of phononsfor a quadratic model with g2 > 0, g4 = 0 in the atomic limitt = 0, for Ω = 1. . . . . . . . . . . . . . . . . . . . . . . . . . 148xvList of FiguresC.2 Sketch of the lattice potential for i) Holstein, and ii) g2 > 0quadratic models. Full (dashed) lines indicate ionic potentialand ground state wavefunction without (with) an extra chargeon the site. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148xviList of SymbolsH, H Generic HamiltonianMA Momentum average, Momentum average approximationel-ph electron-phononqp quasiparticlet A charge carrier’s hopping amplitudec†i , ci Real-space fermionic creation/annihilation operatorsb†k, bk Momentum-space bosonic creation/annihilation operatorsg Electron-phonon coupling constantgn n-th order electron-phonon coupling constantsΩ Phonon frequencyG(. . . ) Generic symbol for Green’s functionsA(. . . ) Generic symbol for spectral functionsARPES Angular resolved photoemission spectroscopyZ Quasi-particle weightm∗ Effective massΣ(. . . ) Generic symbol for self-energyFn(k, i, ω) Generalized Green’s functionsζ Non-linearity parameter in non-linear Holstein modelEOM Equation of motiong¯0 Momentum-averaged approximate version of the free propa-gatorWn Vectors containing several generalized Green’s functionsαn, βn, γn Matrices in the EOM for the Wnλ Dimensionless coupling-constant of the linear Holstein model,E0(t = 0)/E0(g = 0)xviiList of SymbolsΩat Effective phonon frequency√Ω(Ω + 4g2)U On-site Coulomb repulsionΘ Energy-scale of free-phonon anharmonic termsED Exact diagonalization∆ Binding energy of a bipolaron (Chapter 4) OR Charge-transfer energyAFM AntiferromagneticFM FerromagneticZRS Zhang-Rice singletT-CuO Tetragonal copper oxide CuOCn n-fold rotational symmetryBZ Brillouin zonep†, p Creatio/annihilation operator for holes on oxygen p-orbitalsΓ The (0, 0) point in the Brillouin zoneM The (pi, pi) point in the Brillouin zoneX,X ′ The (0, pi) and (pi, 0) points in the Brillouin zonexviiiAcknowledgementsI want to thank my supervisor, Dr. Mona Berciu, for guiding me through theprocess of my PhD. From my very first inquiry about open PhD positionsto the final drafts of my thesis, she has been extremely helpful at everystep of my academic path. I enjoyed working with her. Always availableyet never micro-managing, she provided the perfect environment for me togrow as a researcher. I appreciate her focus on actual results as opposed toproxy measures such as time spent at the desk. That way I could explorethe beauty of Vancouver and its surrounding mountains while still finishingmy PhD at a good pace and with publications in major journals for each ofmy projects.I also thank the supervisory committee, Dr. Sarah Burke, Dr. MosheRozali and Dr. George Sawatzky for their constructive feedack regardingmy thesis.An important part of my experiene at UBC was the Varsity OutdoorClub and the many friendships forged in the snow or on the rock. They aretoo numerous to list them all here, but I have to especially mention JensVent-Schmidt, Ruanne Lai and Caroline Jung for their support in difficulttimes and for good times in the mountains.I also thank my parents for their support and understanding along theway.Finally, I want to thank Kasia Celler for her love and support and herpatience and understanding during the final months of thesis preparation.xixChapter 1IntroductionSolid states physics is concerned with the physics of many, many particles –on the order of 1023 – interacting with each other. From a purely theoreticalpoint of view, this problem is merely one of applied math: The equation gov-erning most of the phenomena occurring in materials is the non-relativisticSchro¨dinger equation,H |Ψ〉 = E |Ψ〉 (1.1)for the Hamiltonian of electrons and nuclei,H = −Nn∑α=1P2α2Mα−Ne∑j=1p2j2m−Ne∑j=1Nn∑α=1Zαe2|ri −Rα|+Ne∑j 0 is a small convergence factor.Good introductions to Green’s functions are found in most advanced solid-state physics textbooks such as [15]. We note here that we are interested inthe interaction of a single carrier with the phonons; therefore, the vacuumstate |0〉 is understood to contain no other carriers, i.e., ci |0〉 = 0. Weintroduce the spectral functionA(k, ω) = − 1piImG(k, ω) (1.6)which is directly related to a quantity accessible to experimentalists viaAngular Resolved Photoemission Spectroscopy (ARPES)2 By inserting afull set of eigenstates |ψn〉 of H into (1.5), it is a textbook derivation to2For an in-depth discussion of ARPES, see the excellent review in [16]. For a derivationof the ARPES intensity expressed through the spectral function, see the PhD Thesis ofDr. Ebrahimnejad [17].71.2. The Holstein model for electron-phonon couplingshow thatA(k, ω) =∑n∣∣∣〈ψn|c†k|0〉∣∣∣2 δ(ω − En(k)). (1.7)That is, the spectral function is a set of peaks3 located at the Hamiltonian’seigenenergies, with their weights given by the overlap of the correspondingeigenstate with that of a free particle. In particular, the lowest-lying (inenergy) peak tells us where the groundstate of the system is, and its weightis called the quasi-particle weight Z. This quantity can range from 0 to 1and indicates how similar to a free particle’s state the ground state is. In asystem where a carrier interacts with bosonic degrees of freedom, we expectthe free carrier to become dressed with a cloud of bosons and thus acquirea larger effective mass, expressed through a reduced qp weight Z.Let us at this point introduce another quantity that is often mentionedin solid state physics, the self-energy. We begin by noting that for a freeparticle with dispersion E0(k) the momentum-space Green’s function – alsocalled the free propagator – is given byG0(ω,k) =1ω − E0(k) + iηfor small convergence factor η > 0. For an interacting particle, the Green’sfunction will be modified. Via simple algebra, one can, however, alwayswrite it asG(ω,k) =1ω − E0(k)− Σ(ω,k) + iηfor some complex-valued function Σ, the self-energy. The real part of Σrenormalizes the energy of the quasi-particle whereas the imaginary partchanges its lifetime.An advantage of working with the Green’s function G instead of theHamiltonian H is that the Green’s function has a diagrammatic expansionthat allows us to separate out the ”easy”, free part from of the Hamiltonian3The δ-peaks are, in fact, Lorentzians with width given by η81.2. The Holstein model for electron-phonon couplingf+kfyfk+kfyfzfk+kkf{yfkk (1.9)Figure 1.2: Sketch of the first two orders of the Dyson series representedby Feynman diagrams. Straight lines represent the free propagator of thecharge carrying fermion. Squiggly lines represent the free propagator of theboson (the optical phonon in this case). A complete diagram would includelabels at the vertices denoting the energies and momenta (or other relevantquantum numbers) of the particles.from the ”hard”, interacting part. If we write this separation asH = H0 +H1and introduce the operator version of the Green’s function as Gˆ(ω) = [ω −H+ iη]−1, then simple algebra shows thatGˆ(ω) = Gˆ0(ω) + Gˆ(ω)H1Gˆ0(ω) (1.8)which is the so-called Dyson’s identity and Gˆ0(ω) = [ω −H0 + iη]−1 is thefree propagator. The decomposition of H is ideally chosen such that G0 iseasy to compute.Dyson’s identity prescribes an iterative way to compute the Green’s func-tion: By reinserting the identity into itself, an infinite series is obtainedwhere each term is of the form Gˆ0[H1Gˆ0]n. When going to a particular ba-sis, such as momentum space, these terms can be represented as Feynmandiagrams, with vertices representing the interaction H1 and edges represent-ing the free propagator Gˆ0(ω).A good introduction to the use of Feynman diagrams in many-bodyphysics is found in [18]. We show a sketch of them in Fig. 1.2. It representsthe first few iterations of Dyson’s identity in a simplified form. In a com-plete representation, the various propagating lines would be labeled with theenergies and momenta (or other relevant quantum numbers of the particles.91.2. The Holstein model for electron-phonon couplingThe first term in Fig. 1.2 describes a carrier propagating freely. The secondterm describes the creation and subsequent absorption of a phonon. Thefinal two terms describe the creation and absorption of two phonons, oneterm where the phonons get destroyed in the order they were created andone term where they get destroyed in the inverse order. For completeness,we note that if we had defined our vacuum |0〉 as a state with a numberof charge carriers below the Fermi level instead of a state with no carriersat all, the diagrammatic expansion would also include diagrams where thecarrier creates a phonon and that phonon then creates a carrier-hole pairby exciting a carrier from below the Fermi level into a state above it. Sincewe are only concerned with the physics of a single carrier interacting with acloud of bosons, the vacuum is defined to contain no extra carriers and thusthese terms do not arise.While these diagrams are an excellent tool for visualizing many-bodyinteraction processes, the sum over the infinite number of diagrams cannotbe carried out except in the most basic cases. Instead, several differentapproaches exist:• If the interactions in H1 are sufficiently weak, the diagrammatic ex-pansion can be terminated after a reasonably low number of iterations.This works very well in other fields such as quantum electrodynamics,but usually fails in many-body physics due to the stronger interactions.In our specific case, it works if g ≤ Ω and provides the weak-couplinglimit perturbation expansion of the Holstein model.• More common in solid state physics are approaches that sum an infinitesubset of diagrams. These diagrams are chosen in such a way that theinfinite sum can be carried out analytically. An example relevant forthe Holstein model is the Self-Consistent Born Approximation (SCBA)[19] where only those diagrams are included where phonon lines donot overlap (such as the last term shown in Fig. 1.2). While theseapproaches are very popular, a common problem is that they are notcontrolled because they lack a ”small parameter” and the diagramsthat are omitted may have contributions comparable to the ones that101.2. The Holstein model for electron-phonon couplingare included.In short, the first approach sums all diagrams up to a certain order in theinteraction H1 whereas the second approach sums diagrams of all orders,but only a subset of them. In both approaches, the diagrams that do get in-cluded are included with their exact contributions. The momentum averageapproximation provides a unique third approach: It sums all diagrams ofall orders, but it individually approximates each diagram with a simplifiedexpression in a way that allows us to analytically carry out the summation.Let us demonstrate this now. Dyson’s equation for the Holstein modelyieldsG(k, ω) = 〈0|ckGˆ(ω)c†k|0〉= 〈0|ck(Gˆ0(ω) + Gˆ(ω)H1Gˆ0(ω))c†k|0〉 .In the following, we will omit the hats and the argument ω whenever thiscan be done unambiguously.The first part of the expansion of G yields just the non-interactingGreen’s function. Since the vacuum has no phonons and c†k only adds afermion, this gives just 1/(ω − εk + iη) for η → 0+, where εk is the freecarrier dispersion. It depends of course on the particular form Hel. Ford-dimensional nearest-neighbor hopping in a cubic lattice, the dispersion is−2t∑i∈{x,y,z} cos(kia) where k is the momentum and a the lattice constant.For the second part, we insert 1 =∑k |k〉 〈k| and use the fact that G0 isdiagonal in k. This finally yieldsG(k, ω) =1ω − εk + iη ·[1 + 〈0|ckGH1c†k|0〉]. (1.10)To study how H1 acts on c†k |0〉, we transform the momentum to real spaceviac†k =1√N∑ieik·ric†i111.2. The Holstein model for electron-phonon couplingand obtainG(k, ω) =1ω − εk + iη ·[1 +1√N∑ieik·ri 〈0|ckGH1c†i |0〉].Let us focus on just the evaluation of the braket. H1 contains a sum oversites j; because of the term nj , its contribution is zero except for the site iwhere we have created an electron. Since we have no bosons in the system,only the term with b†i gives a contribution. The resulting state then isH1c†i |0〉 = gc†ib†i |0〉 .The resulting matrix element with the Green’s function gets its own name:F1(k, i, ω) := 〈0|ckGˆ(ω)c†ib†i |0〉 (1.11)where the subscript 1 denotes that there is one phonon in the system.We use the Dyson equation a second time to derive an equation of motion(EOM) for F1.F1(k, i, ω) = 〈0|ckGˆ(ω)c†ib†i |0〉= 〈0|ckGˆ(ω)Hˆ1Gˆ0(ω)c†ib†i |0〉 .Note that we dropped the single term Gˆ0(ω) from the Dyson expansion.This is valid because F1 is a matrix element between states with differingnumber of phonons (1 and 0); since G0 conserves the phonon number, thismatrix element is 0. As an intermediate result, we obtainF1(k, i, ω) =∑jG0(j − i, ω − Ω) · 〈0|ckGˆ(ω)Hˆ1c†jb†i |0〉 .where we have switched the basis representation of G0 from momentumspace to position space,G0(j − i, ω) = 〈0|cjGˆ0(ω)c†i |0〉 .121.2. The Holstein model for electron-phonon couplingNow we have to study the action of H1. This time, both the bosonic anni-hilation and creation part have a contribution. Let us treat the destructionpart first. Although cumbersome to write down, it should be clear whathappens: We can only destroy a phonon on sites where there is one, andonly if there is also an electron on that site. This introduces a term δij , sowe obtain 〈0|ckGˆ(ω)c†i |0〉 = G(k, i, ω), where G(k, i, ω) is just the Green’sfunction in mixed basis representation: We add a particle on site i and re-move it with momentum k. The creation part will create another boson atthe site the electron is at.Up to this point, the derivation has been exact. Now we introduce thecrucial approximation: Bosons will be allowed only at a single site i. Toget bosons at another site, one has to first get rid of all the bosons at sitei. This approximation leads to a form for the creation part that is similarto the form for the destruction part, i.e. 〈0|ckGˆ(ω)c†i (b†i )2|0〉 =: F2(k, i, ω).Putting everything together gives usF1 = gG0(0, ω − Ω) · [G(k, i, ω) + F2(k, i, ω)]where G0 is given in real space. This is the Fourier transform of G0(k, ω−Ω),and the dc-value is just the integral over the transformed function, i.e.G0(0, ω − Ω) = 1N∑kG0(k, ω − Ω) =: g¯0(ω − Ω).In other words, the exact form of the free propagator is replaced by itsmomentum-averaged form, giving the method its name. The variationalinterpretation of this approximation is that the full Hilbert space is replacedby a restricted Hilbert space containing only states of the nature|i, j, n〉 = c†i(b†j)n |0〉 .This restriction can be systematically relaxed by allowing states with 1 cloudand one extra phonon, which gives rise to the first systematic improvementof MA, called MA(1). In turn, the basic form of MA is called MA(0). Any131.2. The Holstein model for electron-phonon couplingnumber of improved versions MA(n) can be formulated, where n gives thenumber of additional phonons allowed outside the “main” cloud4. Obviouslyin the limit n→∞ the full Hilbert space is recovered and the approximationconverges to the exact result. However, it turns out that for the Holsteinmodel, adequate accuracy for quasi-particle properties is obtained alreadyfor MA(0). For now, we return to the evaluation of the equations of motion.We see that the equation of motion for F1 contains G(= F0) and F2. Wecan see that an equation of motion for F2 would contain F1 and F3, sinceeach additional application of H1 will either remove or add a phonon. Thus,let us now derive the general expression for Fn. It isFn(k, i, ω) = 〈0|ckGˆ(ω)c†i (b†i )n|0〉= 〈0|ckGˆ(ω)H1Gˆ0(ω)c†i (b†i )n|0〉=∑jG0(j − i, ω − nΩ) 〈0|Gˆ(ω)Hˆ1c†j(b†i )n|0〉The destruction part ofH1 will give us again a δij . Next, note that [bi , f(b†i )] =∂/∂b†if(b†i ) for any function f . Hence [bi , (b†i )n] = n(b†i )n−1. For the creationpart, we make again use of the approximation that additional bosons canbe added only to sites where bosons are already present. It is then straight-forward to arrive at the recursion relation for Fn,F0(k, i, ω) = 〈0|ckGˆ(ω)c†i |0〉Fn(k, i, ω) = gg¯0(ω − nΩ) · [n · Fn−1(k, i, ω) + Fn+1(k, i, ω)] . (1.12)By introducing a cut-off in n, we can find an explicit solution for theserelations. The index n denotes the number of phonons we add to the system.The quantity Fn(k, i, ω) then is the amplitude for putting n phonons and4We remark that MA(n) as formulated in, e.g., [13] are not exactly equivalent tovariational models with states c†i b†a1 . . . b†an(b†j)m; instead, they are defined by how manyorders in the Dyson expansion retain the full form of the free propagator: MA(0) keepsthe free propagator only in the 0-th order and replaces it with its momentum averagedversion beginning at the first order, whereas MA(1) would treat the free propagators inboth the 0-th order and the 1st order exactly and only replace it with the momentumaveraged version beginning at the second order.141.2. The Holstein model for electron-phonon couplingone electron at site i into the system and later having a state where theelectron has momentum k and there are no phonons left in the system. Wecan therefore expect that Fn → 0 for n → ∞. For a cut-off N , we setFn = 0 for all n > N . How to determine a good value for the cut-off willbe discussed later. Let us for now just see how we can solve the recursionrelations for Fn. With the cut-off, we haveFN (k, i, ω) = gg¯0(ω −NΩ) ·N · FN−1(k, i, ω)since FN+1 = 0. This means that FN (k, i, ω) is proportional to FN−1(k, i, ω)and we introduce AN = gg¯0(ω −NΩ) ·N to writeFN (k, i, ω) = AN · FN−1(k, i, ω).With FN now being given in terms of FN−1 only, we can continue. Let usomit the arguments (k, i, ω) from the Fn as they do not change anyway. Letus also introduce the shorthand g¯0(N) for g¯0(ω −NΩ).FN−1 = gg¯0(N − 1) · [(N − 1)FN−2 + FN ]= gg¯0(N − 1) · [(N − 1)FN−2 +ANFN−1] .This equation contains only FN−1 and FN−2. It can be rearranged to giveFN−1 = AN−1FN−2. We see that this continues all the way down to F1,which then will be given as F1 = A1G. Let us derive a general expressionfor the An. For AN , we have the explicit expression from above. The otherAn are defined viaFn = g · g¯0(n) · [nFn−1 + Fn+1]= g · g¯0(n) · [nFn−1 +An+1Fn]151.2. The Holstein model for electron-phonon couplingThis can be rearranged to(1− gg¯0(n)An+1) · Fn = gg¯0(n)nFn−1⇔ Fn = n · gg¯0(n)1− gg¯0(n)An+1 · Fn−1.Thus, by definition,An =n · gg¯0(n)1− gg¯0(n)An+1 . (1.13)We are interested in A1, because via F1 = A1G, we can obtain an explicitexpression for G. A1 is given in terms of a continued fraction.A1 =gg¯0(1)1− gg¯0(1)2 · gg¯0(2)1− gg¯0(2)3 · gg¯0(3)1− . . .We will come back to the evaluation of this expression in an instant, but firstlet us see how knowledge of A1 helps us determining the Green’s function.Inserting what we have derived so far into (1.10) yieldsG(k, ω) =1ω − εk + iη ·[1 +g√N∑ieik·riA1G(k, iω)]=1ω − εk + iη · [1 + gA1G(k, ω)]We can rearrange this to obtainG(k, ω) =1ω − εk + iη − gA1(ω) . (1.14)A closer look reveals that now we not only have an explicit expression forG(k, ω) in terms of A1, we also see that gA1 is the self-energy of the system.Since A1 involves only the momentum-averaged expression g0(ω), the self-energy is k-independent.There are several ways to solve the continued fraction. Let us define Ak1161.2. The Holstein model for electron-phonon couplingas the value for A1 that is obtained when choosing N = k as the cut-off,i.e., all An with n > N are zero. Then one can compute Ak1 for increasingvalues of k until the change Ak+11 −Ak1 is below a certain threshold. This isthe straight-forward approach. The disadvantage is that it involves a lot ofcomputation, and to make the step from k to k + 1 one has to repeat thecomplete computation as the result from the computation of Ak1 cannot beused to arrive at Ak+11 . A better method is the modified Lentz method asdescribed in [20]. The advantage of that method is that we do not have toset a cut-off N in advance. Instead, the n-th convergent can be computedusing the n − 1-th convergent. To use this method, we have to write ourcontinued fraction in the canonical form,x = b0 +a1b1 +a2b2 +a3b3 +.. ..In our case, we haveb0 = 0 a1 = gg¯0(1)bn = 1 an = −g2 · n · g¯0(n− 1) · g¯0(n).The calculation of the bn and an does not require any recursion, so it is trivialto write a function that returns for a given n the appropriate coefficient anor bn. The modified Lentz method is given in Alg. 1. The value of 10−30can be any small value and is needed to avoid division by zero. The valueof ε denotes the desired accuracy of our computation.Figure 1.3 shows some typical spectral functions of the one-dimensionalHolstein model as obtained via MA. It is customary to introduce the di-mensionless coupling constant λ = g2/2dtΩ. As we’ll show later, this is theratio of the model’s ground state energies in the strong and the weak inter-action limits. In the first panel, λ = 0, we just recover the free electron’sdispersion −2t cos(ka), with lattice constant a set to 1. At weak coupling,171.2. The Holstein model for electron-phonon couplingAlgorithm 1: Modified Lentz Methodbeginx := b0if b0 = 0 thenx := 10−30endC := xD := 0j := 0∆ := 0while |∆− 1| > ε doj := j + 1b := b(j, ω)a := a(j, ω)D := b+ a ·Dif D = 0 thenD := 10−30endC := b+ a/Cif C = 0 thenC := 10−30endD := 1/D∆ := C ·Dx = x ·∆endend181.2. The Holstein model for electron-phonon coupling(a) λ = 0 (b) λ = 0.2,Ω = 0.5(c) λ = 1,Ω = −0.5 (d) λ = 10,Ω = 1Figure 1.3: Spectral function of the one-dimensional Holstein model as ob-tained with MA. The effective coupling λ = −g2/2dtΩ is a dimensionlessmeasure for the strength of el-ph interaction. The strong coupling regimeoccurs at around λ ∼ 1.191.2. The Holstein model for electron-phonon couplingλ = 0.2, we observe two things: The dispersion of the lowest lying state stillhas roughly a cosine shape −2t∗ cos(ka) but with t∗ < t. This is because thecarrier becomes dressed with phonons, which increases its effective mass. Inaddition, the spectral function also shows a continuum of energies. Theseare the “polaron+1-phonon” states: For some momentum k, the system’sstate can either be a polaron with energy E0(k) or it could be a polaronwith energy E0(k− q) and, far away, a phonon with momentum q and en-ergy Ω. The energy of this state would be E0(k − q) + Ω. Since q can beany momentum, at momentum k the polaron + phonon states span all theenergies from min(E0)+Ω to max(E0)+Ω. That is, the continuum starts atexactly Ω above the polaronic ground state. Close inspection of the spectralfunction as obtained from MA, however, shows that the continuum insteadoriginates at −2t+ Ω, i.e., precisely Ω above the free-electron ground state.This discrepancy is easiest explained in the variational picture of MA. Re-call that MA(0) is equivalent to restricting the Hilbert space to states whereall phonons are located on the same site. Therefore, within this Hilbertspace we cannot describe states that have a phonon cloud in one locationand a single phonon in another location. Higher orders of MA obtain moreaccurate estimates of the start of the continuum by allowing states wherethere are phonons away from the cloud. A more detailed discussion is foundelsewhere [13]. The main observation is that MA(n) predicts the continuumto start an energy Ω above the polaronic ground state of MA(n−1). Furtherincreasing the coupling λ leads to a further reduction of the polaronic band-width and the appearance of additional continua separated by Ω. Finally,in the very strong interacting limit, the spectrum consists of flat bands withspacing Ω. In fact, the atomic limit t→ 0 can be solved exactly. We do notshow the derivation here because we will solve a more general case – butwith the same underlying method – in Chapter 2. To summarize, the mainproperties observed for the Holstein model are:• Interaction with phonons leads to a bound state, the polaron.• The dressing of a free carrier with a cloud of phonons increases itseffective mass201.3. The three-band Emery model for a cuprate layer• The larger the coupling, the larger the effective mass.• MA works best for ground-state properties.For more in-depth discussion and some of the finer points of MA, such as itsdifference to Dynamical Mean Field Theory (DMFT), we refer the reader tothe literature [1, 12, 13, 21].These results set the stage for our extensions to the standard (linear)Holstein model. In our derivation and motivation of the electron-phononcoupling term, there was no a-priori reason for assuming that the two ionicpotentials I and II in Fig. 1.1 only differ in the equilibrium position butnot in the phonon frequency. Furthermore, the harmonic approximation forpotentials I and II is only valid for small lattice distortions, whereas largecoupling, as we have seen, predicts large lattice distortions. Hence, the linearmodel’s predictions at strong coupling invalidate its assumptions and a morethorough approach is required. In Chapter 2, we study the effect of includingterms up to fourth order in the lattice coordinates. We will show how evensmall non-linear terms have a big qualitative and quantitative effect on thepolaron properties. Next, we study a particularly interesting special caseof a non-linear electron-phonon model where the lattice potential containsonly even terms in the phonon coordinates. For certain parameters, thisleads to a model where the carrier’s presence turns the local single-welllattice potential into a double-well. This is studied in the single-carrier casein Chapter 3 and for the two-carrier case in Chapter 4. We find that thepeculiar nature of the el-ph interaction leads to quasiparticles and boundpairs of quasiparticles (so-called bipolarons) that have properties that aredistinctly different from those of linear models.1.3 The three-band Emery model for a cupratelayerOne of the longest standing puzzles in solid state physics is the microscopicorigin of high-temperature superconductivity in the cuprates [22]. We show211.3. The three-band Emery model for a cuprate layerFigure 1.4: Phase diagram of the high-temperature superconductingcuprates for both electron and hole doping. The undoped compounds are an-tiferromagnetic insulators which become superconducting upon either holeor electron doping. Here, we focus on the hole-doped side of the phase dia-gram. Public domain figure from https://en.wikipedia.org/wiki/File:Cuphasediag.png.221.3. The three-band Emery model for a cuprate layerCopper(a) YBCO crystal structureOCu(b) Unit cell of the 3-band Emery modelFigure 1.5: (a) Structure of a typical high-temperature superconductingcuprate compound, YBCO. Layers of copper and oxygen are containedwithin a host structure of rare-earth and transition metal elements. Fig-ure used under the GFDL, as obtained from https://commons.wikimedia.org/wiki/File:Ybco002.svg. (b) The most relevant orbitals for cupratephysics are the copper dx2−y2 orbital (round shapes) and the ligand oxygenpx and py orbitals (lobes), forming a unit cell with three orbitals.231.3. The three-band Emery model for a cuprate layertheir crystal structure and typical phase diagram in Fig. 1.4. The par-ent compounds are antiferromagnetic insulators. Upon doping electrons orholes into then, they become superconducting. Here, we focus on the case ofdoping a single hole into the insulating parent compound. While supercon-ductivity only arises at finite doping, the “simple” case of having a singlecarrier doped into the system is already very complex, as we shall see below.With the goal of arriving at a model Hamiltonian that is simple enough,but not too simple, we will discuss the progression of proposed model Hamil-tonians. Experiments have shown that the important physics happen in two-dimensional CuO2 layers of copper and oxygen. In the un-doped case, therelevant orbitals close to the Fermi energy are the copper’s in-plane dx2−y2orbital and the oxygen’s px and py orbitals. The cuprates are so-calledcharge-transfer insulators. This means that inserting an electron (remov-ing a hole) takes place on the copper sites whereas removing an electron(inserting a hole) takes place on the oxygen sites [23]. A popular modelHamiltonian, then, is the so-called p-d or Emery model [24]. We give herethe so-called 3-band version, where as a further approximation we note thatthe copper orbital has the largest overlap with those oxygen orbitals that arealigned with it. A unit cell for this model thus contains three sites and threeorbitals as shown in Fig. 1.5(b). Its Hamiltonian is then given by (comparealso its extensive discussion in [6])H3B = Tpd+Tpp+∆pd∑nl+,σ+Upp∑nl+,↑nl+,↓+Udd∑nl↑nl↓ (1.15)withTpd = tpd∑(−p†l+,σ + p†l−,σ)dlσ + h.c. (1.16)Tpp = tpp∑sδp†l++δpl+,σ − t′pp∑(p†l−,σ + p†l+3,σ)pl+,σ (1.17)nl,σ = d†l,σdl,σ (1.18)nl+,σ = p†l+,σpl+,σ (1.19)In this Hamiltonian, Tpd describes hopping between oxygen and copper,241.3. The three-band Emery model for a cuprate layerTable 1.1: Relevant parameters for the 3-band Emery model in units of eV.tpd tpp t′pp ∆pd Upp1.3 0.65 0.38 3.6 4Tpp describes oxygen-oxygen nearest (tpp) and next-nearest (t′pp) neighborhopping. The next-nearest neighbor hopping is, in fact, mediated by thecopper 4s orbitals. Other hoppings, such as y-direction hopping betweentwo px orbitals, is neglected. sδ is the sign of the matrix elements betweenvarious oxygen p orbitals. It depends on our arbitrary choice of the orbitals’phases. We make our choice such that sδ = 1 for hopping along the maindiagonal (upper-right and lower-left) and −1 for the other directions. Thecharge-transfer energy, i.e., the energy cost of moving a hole from a copperto an oxygen site, is given by ∆pd. The on-site Coulomb repulsions are givenby Udd and Upp, respectively. Numerical values for these parameters havebeen obtained from ab-initio calculations in [25] and are presented in Table1.1.Despite being a relatively simple Hamiltonian, the model in (1.15) al-ready contains a variety of approximations, such as neglecting longer-rangeCoulomb interaction and longer range hopping. Despite this, it is still re-garded as too complex to be solved directly due to its many degrees offreedom. Hence, a simpler model Hamiltonian is required. Zhang and Rice,in their seminal work, argued that although hole insertion occurs on the oxy-gen sites, a particular linear combination of the four ligand oxygen orbitalsaround a copper orbital forms a singlet with that orbital, the so-called ZhangRice singlet [26]. By discarding all other linear combinations of the oxygenorbitals, one arrives at an effective single-band model where the Zhang Ricesinglet (ZRS) combines the hole’s charge degree of freedom with the copperion’s spin degree of freedom into a singlet object. This object then hops inthe antiferromagnetic (AFM) background of the copper sites and gives riseto the so-called tJ-model, where t refers to the single-band hopping and J tothe exchange interaction of the AFM background. As it turns out, the tJ-model itself does not adequately describe the single-particle dispersion of the251.4. Outlinehole-doped cuprates. To obtain agreement with experiments, longer rangehopping is added to the model to arrive at the t-t′-t′′-J model. With thisextended model, quantitative agreement with the experimental dispersion isreached for a suitable choice of the parameters.While these models have received considerable attention from the physicscommunity [26–34], a clear and convincing picture of superconductivity inthe cuprates has still not emerged. There is some evidence that one-bandmodels fail to capture the essential physics of the cuprates. However, we stillneed to simplify the Hamiltonian (1.15). This is achieved in two steps: First,we apply the Udd →∞ limit to the three-band Emery model instead of thesingle-band Hubbard model. The resulting model contains a variety of inter-esting higher-order effective terms which will be discussed in Chapter 5. Ithas been studied numerically via exact diagonalization in Bayo Lau’s disser-tation [6]. A further significant approximation was made by Hadi Ebrahim-nejad [3, 5]; it neglects spin-fluctuations between the copper holes. Whilethese fluctuations are important in the t − J model, he showed that theyplay only a very little role in the Emery model. This approximations allowsfor a variational approach much like MA in the phonon case, where the vari-ational space restricts the number of magnons present in the system. Goodqualitative convergence is obtained for even just one magnon, and excellentquantitative convergence is obtained for two magnons, as demonstrated inFig. 1.6.Curiously, despite being based on vastly differing physics, these single-and many-band models agree on the single-particle dispersion. It is there-fore not quite clear which of them more appropriately describes the cupratephysics. In Chapter5, we consider an extension to the cuprates based on thenewly synthesized tetragonal copper oxide T-CuO. We show how careful ex-periments with this material can settle this longstanding modeling question.1.4 OutlineThe remainder of this thesis is organized as follows. Chapter 2 discussesin detail the single-polaron physics of the nonlinear Holstein model. We261.4. Outline-21-20-19E(k)/J dd(pi,0) (0,0) (pi,pi) (0,pi) (pi,0) k-22-21-20E(k)/J dd(a) nm = 2,r(b) nm = 3,rFigure 1.6: E(k) along several cuts of the Brillouin zone for the three-bandmodel in the variational approximation described in [3, 5] with (a) at most2 magnons and (b) at most 3 magnons. Circles show exact diagonalizationresults from [6] for a cluster of 32 copper and 64 oxygen orbitals. Fulllines show the results from the variational approximation, neglecting spinfluctuations. The overall agreement is excellent and can be further improvedby including local spin fluctuations [5].271.4. Outlinefind that at moderate to strong linear coupling, the addition of even smallquadratic terms drastically alters the polaron properties. Chapter 3 intro-duces a special variation of the nonlinear Holstein model, the double-wellelectron-phonon coupling model. In this model, the presence of a charge car-rier switches the local lattice potential from a single- to a double-well. Westudy in detail the single-polaron properties of this model and find that theycannot be described with an effective linear Holstein model. We continuethe study of this model in Chapter 4, where we investigate the bi-polaron,a bound state made of two charge carriers and a cloud of phonons. Ourresults here indicate that realistic parameter regimes exist for which the bi-polaron is strongly bound yet lightweight, in contrast to the results of otherel-ph coupling models where the bi-polaron’s mass increases drastically thestronger it is bound. Finally, in Chapter 5 we study a system where therelevant bosons are magnons instead of phonons; for a layer of tetragonalcopper oxide (T-CuO) we study the dispersion of a single hole doped into theoxygen orbitals. Here, our finding is that the small inter-sublattice hoppingsufficiently changes the symmetry properties of the resulting quasiparticle (aspin-polaron) from those of CuO2 that experimental studies of T-CuO willbe able to settle longstanding questions regarding the single-particle physicsof the underdoped cuprates.28Chapter 2Nonlinear Holstein modelIn this chapter, we introduce a generalized Holstein model with nonlinearelectron-phonon couplings. The model is then studied with the momentum-average approximation. We find that the inclusion of even small non-linearterms can have drastic effects on the properties of the system. This chapteralso provides further review of the momentum average approximation.2.1 IntroductionCoupling of carriers to phonons and the properties of the resulting quasipar-ticles, the polarons, are important for many materials, e.g. organic semicon-ductors [35, 36]. cuprates [37–42], manganites [43], two-gap superconductorslike MgB2 [44–47], etc. In some cases the effective electron-phonon (el-ph)coupling λ is known quite accurately. For others, like the cuprates, esti-mates range from very small (λ ∼ 0.3) to very large (λ ∼ 10) [48, 49]. Onepossible explanation for this is that, especially for stronger couplings wheresimple perturbational expressions are no longer valid, properly fitting theexperimental data to theoretical models can be quite involved [50].Here we consider another possible explanation, namely that at strongel-ph coupling, simple theoretical models may not be valid anymore.As mentioned in the introduction to this thesis, all widely-used models[11, 51] assume at the outset that the displacements xi of the atoms outof equilibrium are small enough to justify expanding the electron-lattice in-teractions to linear order in xi. These linear models generically predict theformation of small polarons or bipolarons at strong coupling, with the car-rier(s) surrounded by a robust phonon cloud. As a result, lattice distortions〈xi〉 are considerable near the carrier(s). Hence, the linear models are based292.1. Introductionon assumptions which are in direct opposition to their predictions.In this chapter we investigate this issue in the single polaron limit, rele-vant for the study of weakly doped materials like very underdoped cuprates[52, 53] and organic semiconductors [35, 36], and for cold atoms/moleculessimulators [54–57]. We study the ground-state (GS) of a single polaron in ageneralized Holstein model including el-ph coupling up to quartic order inxi to test the importance of the higher order terms. We find that for stronglinear coupling even very small quadratic terms drastically change the prop-erties of the polaron. Moreover, we show that these effects go beyond amere renormalization of the parameters of the linear Holstein model. As aresult, attempts to find effective parameters appropriate for a linear modelby using its predictions to fit the properties of real systems are doomed tofailure, as different values will be obtained from fitting different properties.This offers another possible explanation for the wide range of estimates ofthe el-ph coupling in some materials. More importantly, it means that wemust seriously reconsider how to characterize such interactions when theyare strong. Furthermore, this calls for similar investigations of the validityof these linear models at finite carrier concentrations, since it is reasonableto expect that they also fail in the strong coupling limit.To the best of our knowledge, we present here the first systematic, non-perturbative study of the importance of higher-order el-ph coupling termson single polaron properties. We note that in previous work going beyondlinear models, purely quadratic (no linear term) but weak el-ph couplingwas discussed for organic metals using perturbation theory [58, 59], whilelinear and quadratic el-ph coupling was studied in the context of high-TCsuperconductivity in Ref. [60]. A semi-classical study of some non-linearcoupling potentials was carried out in Ref. [61].This chapter is organized as follows: in Section 2.2 we introduce andmotivate the nonlinear Holstein model in more detail than we did in theintroduction. In Section 2.3 we show the modifications to the momentumaverage approximation necessary to solve the model. The results are dis-cussed in Section 2.4. An outlook and concluding remarks are given inSection 2.5.302.2. Model2.2 ModelWe start this section with stating the model Hamiltonian for the generalizednonlinear Holstein model introduced in [62] and the introduction.H = Hel +Hph +Hel−ph, (2.1)The remainder of this section will derive, motivate and explain these terms.The Holstein Hamiltonian models a charge carrier in a molecular crystallike the 1D example sketched in Fig. 2.1(a). A charge carrier introduced insuch a crystal hops between “molecules”, as described byHel =∑kkc†kck, with k = −2td∑α=1cos(kα) (2.2)for nearest-neighbor hopping on a d-dimensional simple cubic lattice.Fig. 2.1(b) illustrates how the lattice part is handled. In the absence ofa carrier, the potential has some form (curve I) which is approximated as aparabola and leads to Hph = Ω∑i b†ibi . This describes harmonic oscillationsof each “molecule” about its equilibrium distance R. If a carrier is present,the potential has some other form (curve II). The difference between I andII leads to Hel−ph. Its details are material specific; here we propose twomodels and choose a generic form based on them.The first model assumes that the carrier occupies an orbital of the ionwith opposite charge. The attraction between them is then some constant,whereas the Coulomb repulsion between the carrier and the ion of like chargeisU(xi) =U0ni1− xiR= U0ni∞∑n=0(xiR)nwhere ni = 1 (ni = 0) if the carrier is (is not) present and U0 > 0 is thecharacteristic energy. Using the expression of the position operator in termsof the phonon creation and annihilation operators, xi =√h¯2µΩ(bi+b†i ) whereµ is the reduced mass of the molecule, and truncating the series at n = 4312.2. ModelpU (x)−+i-1 i+1i(a) (b)RIIIRFigure 2.1: (a) Sketch of a 1D chain of polar molecules; (b) The potentialof the pair with (II) or without (I) an extra charge carrier (full lines) isapproximated by a polynomial (thick dashed lines)leads to:Hel−ph =4∑n=1H(n)el−ph =4∑n=1gn∑ic†ici(bi + b†i )n, (2.3)where gn = g1ζn−1 with ζ = A/R and A =√h¯/(2µΩ) the zero-pointamplitude of the harmonic oscillator.The second model assumes that the carrier is an electron (hole) thatoccupies an anti-bonding (bonding) orbital of the molecule; all bonding or-bitals are initially full since the parent crystal is an insulator. In both casesthe energy increases by an overlap integral which decreases exponentiallywith the distance:U(xi) ∼ nie−R−xiaB (2.4)where aB is the Bohr radius. A Taylor expansion to fourth order in xi leadsagain to Eq. (2.3) but nowgng1= 2n−1ζn−1n!for ζ =g2g1=A2aB, (2.5)where again A =√h¯/(2µΩ).We define the following special cases:Linear model The case where only g1 6= 0. This is the standard Holsteinmodel.322.2. ModelQuadratic model The case where only g1 6= 0 and g2 6= 0.Quartic model The case where all gn 6= 0 for n ≤ 4.The case with only g4 = 0 is not considered because it is unstable: Apolynomial with degree 3 will always diverge to −∞ for one of the directionsx→ ±∞.As discussed in the introduction, the linear Holstein model is charac-terized by two dimensionless parameters: the effective coupling strengthλ = g21/(2dtΩ), where d is the dimension of the lattice, and the adiabaticityratio Ω/(4dt). As long as the latter is not very small, the former controlsthe phenomenology, with the crossover to small polaron physics occurringfor λ ∼ 1 [63]. For ease of comparison, we continue to use these parameterswhen characterizing the higher order models. For the quadratic model, thenew energy scale g2 results in a third dimensionless parameter ζ = g2/g1.For the quartic model there are two more parameters gn/g1, n = 3, 4. Bothscale like ζn−1 but with different prefactors. We use gn/g1 = ζn−1 like inthe first model since for the second model the prefactors are less than 1,making these terms smaller and thus less important.For specificity, from now we assume ζ > 0 (ζ < 0 is briefly discussedat the end of this chapter, and extensively in later chapters). As we showbelow, in this case we find that while quadratic terms are important whenthe linear coupling is large, addition of the n = 3, 4 terms only leads tosmall quantitative changes and can be ignored. This justifies a posterioriwhy we do not include anharmonic corrections in Hph and/or higher orderterms with n > 4 in the electron-phonon coupling.2.2.1 Summary of model assumptionsRecalling the caveat for model Hamiltonians, we summarize here the as-sumptions that went into the model Hamiltonian and give their justification.First, quartic terms in the free-phonon part of the Hamiltonian wereneglected. This will be justified a posteriori by the results, which will showthat quartic terms in the el-ph coupling have only a small effect on thepolaron properties.332.3. FormalismSecond, the coefficients gn follow a power-law dependence on the linearcoupling, gn ∼ ζn−1g1. This is motivated by observing this relation for twodifferent physical origins of electron-phonon coupling.Finally, we assume ζ > 0. This is not so much a model assumptions asit is a focus on one region of parameter space. The case of ζ < 0 will beexplored in the next chapter.2.3 FormalismWe now describe in detail the MA solution for the quadratic model. Thecalculations for the quartic model are analogous but much more tedious.We want to find the single particle Green’s functionG(k, ω) = 〈0|ckGˆ(ω)c†k|0〉 (2.6)where Gˆ(ω) = [ω − H + iη]−1 is the resolvent for this Hamiltonian, withη → 0 a small positive number and |0〉 the vacuum state. From this wecan extract all the polaron’s GS properties [12]. We rewrite the quadraticHamiltonian as H = H0 +H1, whereH0 = Hel +Hph + g2∑ic†ici(2b†ibi + 1)(2.7)H1 =∑ic†ici[g1(b†i + bi ) + g2(b†2i + b2i )]. (2.8)The equation of motion (EOM) for the propagator is obtained recursivelyfrom Dyson’s identity,Gˆ(ω) = Gˆ0(ω) + Gˆ(ω)H1Gˆ0(ω) (2.9)where Gˆ0(ω) = [ω −H0 + iη]−1 is the resolvent for H0. Using it in G(k, ω)yields the EOMG(k, ω) = G0(k, ω)[1 +2∑n=1∑ieik·ri√NgnFn(k, ω; i)](2.10)342.3. Formalismwhere Fn(k, ω; i) = 〈0|ckGˆ(ω)c†i (b†i )n|0〉, which we call a generalized propa-gator.So far, everything has been exact and very similar to the treatment ofthe linear Holstein model in Chapter1, with the exception that our H0 isslightly more complicated, since in addition to free electron propagation italso contains a potential scattering term off the phonons. Following thebasic steps of the momentum average approximation (MA) [12, 13] we letg˜0(ω) =1N∑k1ω − k + iη (2.11)denote the momentum averaged propagator of a free particle with dispersionk. Now consider the propagator associated with H0 in real space:G0(i− j, ω;n) = 1n!〈0|cibni Gˆ0(ω)(b†i )nc†i |0〉 . (2.12)This is the on-site real-space propagator of an otherwise free carrier beingscattered by the on-site potential 2g2n while there are n phonons in thesystem, as described by (2.7). We show in B.1 that for i = j, we haveG0(i− j = 0, ω;n) =[1g˜0(ω − nΩ− g2) − 2g2n]−1. (2.13)The momentum average approximation now consists of settingG0(i− j, ω;n) ≈ δijG0(i− j, ω;n) = g¯0(ω;n) for any n ≥ 1. (2.14)with g¯0(ω;n) = g˜0(ω − nΩ − g2). In words, the real-space free propagatorappearing in the higher-order iterations of Dyson’s identity is approximatedas being purely on-site. This is justified because the polaron GS energylies below the free particle spectrum, and for such energies the free-particlepropagator decreases exponentially with |i− j|. Thus, MA keeps the largestcontribution and ignores the exponentially smaller ones. This becomes exactin the strong-coupling limit t → 0. Since phonon creation and annihilationis contained within H1, the effect of this approximation is that new phonons352.3. Formalismcan only be created at a site where there already are phonons. At this levelof approximation, MA can be viewed as a variational method where onlystates of the form ∑ijc†i(b†j)n |0〉 (2.15)are considered, i.e., states with at most one single-site phonon cloud.By eliminating propagators that take the carrier away from the phononcloud, MA allows us to obtain a simplified hierarchy of EOM involving onlythe generalized Green’s functions Fn. For any n ≥ 1, they read:Fn(k, ω; i) = g¯0(ω;n) ·[n(n− 1)g2Fn−2(k, ω; i)+ ng1Fn−1(k, ω; i) + g1Fn+1(k, ω; i) + g2Fn+2(k, ω; i)].Since the arguments of all Fn propagators are the same, we suppress themin the following for simplicity. The important achievement is that now eachgeneralized propagator is linked via its EOM to only a finite set of generalizedpropagators.In the original MA approach for the linear Holstein model, the EOMcoupled the generalized propagator of order n to those of orders n ± 1.Here, now the connection is to propagators of orders n ± 1, n ± 2. As aconsequence, the EOM cannot be reduced to a simple continued fractionanymore. Instead, we will obtain a continued matrix fraction: Following thetechnique introduced in Ref. [2], we reduce these EOMs to a simple recursiverelation. This is achieved by introducing a vector Wn = (F2n−1, F2n). TheEOM for Wn are γnWn = αnWn−1 + βnWn+1, where the αn, βn and γnare 2× 2 matrices whose coefficients are read off of the EOM, namelyαn|11 = (2n− 1)(2n− 2)g2g¯0(ω; 2n− 1), (2.16)αn|12 = (2n− 1)g1g¯0(ω; 2n− 1), (2.17)αn|21 = 0, (2.18)αn|22 = 2n(2n− 1)g2g¯0(ω; 2n), (2.19)362.3. Formalismwhileβn =(g2g¯0(ω; 2n− 1) 0g1g¯0(ω; 2n) g2g¯0(ω; 2n)), (2.20)γn =(1 −g1g¯0(ω; 2n− 1)−2ng1g¯0(ω; 2n) 1). (2.21)This simple recursive relation for Wn has the solution Wn = AnWn−1 forany n ≥ 1, where An are 2×2 matrices obtained from the infinite continuedfractionAn = [γn − βnAn+1]−1 αn. (2.22)In practice, we start with AN = 0 for a sufficiently large cutoff N , chosen sothat the results are insensitive to further increases in it (N ∼ 100 is usuallysufficient).We find A1 =(0 a120 a22), where a12 and a22 are obtained after usingEq. (2.22) N − 1 times. As a result, F1 = a12F0, F2 = a22F0, whereG(k, ω) =∑ieik·ri√NF0(k, ω; i). (2.23)Using these in Eq. (2.10) leads to a solution of the formG(k, ω) =1ω − k − Σ(ω) + iη , (2.24)with the MA self-energy for the quadratic model:Σ(ω) = g1a12(ω) + g2a22(ω). (2.25)The reason why the self-energy is local at this level of MA is the simplic-ity of this Hamiltonian, whose vertices are momentum independent; thisissue is discussed at length for the linear Holstein model in Ref. [13]. Atthis point, we emphasize that MA is not equivalent – let alone inferior – toDynamical Mean Field Theory. The latter always leads to a momentum-independent self-energy, whereas MA – even at its simplest level – achieves372.4. Results and discussionmomentum dependence for more complicated models, and becomes momen-tum dependent for even the simple Holstein model when higher orders ofthe approximation are used.The quartic model is solved analogously. The main difference is that herethe EOM for Fn involves 9 consecutive terms, from Fn−4 to Fn+4. Thesecan also be rewritten as simple recurrence relations γnWn = αnWn−1 +βnWn+1, but now αn, βn and γn are 4× 4 matrices. Their expressions aretoo long to be listed here, but they can be found in Appendix B.2.2.4 Results and discussionTo gauge the relevance of the higher-order el-ph coupling terms we plot inFig. 2.2 the evolution with ζ of a polaron property that can be directlymeasured, namely the quasiparticle weight Z = m/m∗ where m,m∗ arethe carrier and the polaron mass, respectively. We also show the averagephonon number Nph. The results are for a one-dimensional chain. Resultsin higher dimensions are qualitatively similar to these 1D results for smallλ, and become quantitatively similar to them in the interesting regime oflarge λ where all of them converge towards those of the atomic limit t = 0.First, we note that the ζ = 0 intercepts trace the predictions of the linearmodel: with increased coupling λ, Z decreases while Nph increases as thepolaron acquires a robust phonon cloud [63, 12]. From these intercepts, weestimate that the linear model predicts the crossover to the small polaronregime to occur around λ ∼ 1.5 for this adiabaticity ratio and dimension.The quadratic model, whose predictions are indicated by lines, shows avery strong dependence of ζ for strong linear coupling λ ≥ 1.5: here both Zand Nph vary by about an order of magnitude as ζ increases from 0 to 0.1.For higher ζ, Z and Nph have a slight turnaround towards smaller/largervalues, for reasons explained below, but are still consistent with a largepolaron. These results indicate that the quadratic term can completelychange the behavior of the polaron in the limit of medium and large λ. Forexample, in the quadratic model at λ = 1.5 and ζ ∼ 0.1 the polaron is lightand with a small phonon cloud, in total disagreement with the linear model382.4. Results and discussionprediction of a heavy small polaron at this λ.Of course, this raises the question of how large ζ is. The answer ismaterial specific, but as an extreme case, let H2 be the unit of the molecularcrystal. This case is described by model two, so ζ ∼ 2A/aB, where aB ≈0.5A˚ while A ≈ 0.1A˚ if we use Ω ≈ 0.5eV appropriate for a H2 molecule[64]. This leads to a very large ζ ∼ 0.4. Other atoms are heavier but phononfrequencies are usually much smaller than 0.5eV , so it is not clear whetherA ∼ 1/√µΩ changes much. The Bohr radius (or distance R between atoms,for model 1) is usually larger than 0.5A˚ but not by a lot, maybe up to afactor 5 for R; thus we expect smaller ζ in real materials but the change islikely not by orders of magnitude. Fig. 2.2 shows that values as small asζ ∼ 0.05 already lead to significant quantitative changes in m∗.Inclusion of cubic and quartic terms (the symbols show the results ofthe quartic model) further changes Z and Nph, but these changes are muchsmaller for all ζ, of up to ∼ 10% when compared to the quadratic modelvalues, as opposed to order of magnitude changes between the quadraticand the linear models. Thus, these terms are much less relevant and canbe ignored without losing much accuracy. As discussed, their small effectexplains why we do not consider terms with even higher order n, nor n = 4anharmonic terms in the phonon Hamiltonian.2.4.1 Atomic limitTo understand the effects of the quadratic term at large λ, we study itin the atomic limit t = 0 (λ = ∞) where the carrier remains at one siteand interacts only with the phonons of that site. Focusing on this site, itsquadratic Hamiltonian H(2)at = Ωb†b+∑2n=1 gn(b†+b)n is well-studied in thefield of quantum optics, where it describes so-called squeezed coherent states[65]. The extra charge changes the origin and spring constant of the originalharmonic oscillator which means that the Hamiltonian is easily diagonalizedby changing to new bosonic operators γ† = ub† + vb+ w, where u, v and w392.4. Results and discussion0 0.2 0.4 0.6 0.8 1ζ0.00.20.40.60.8Z0 0.2 0.4 0.6 0.8 1ζ0246Nphλ = 0.5, n = 2λ = 1.0, n = 2λ = 1.5, n = 2λ = 2.0, n = 2λ = 0.5, n = 4λ = 1.0, n = 4λ = 1.5, n = 4λ = 2.0, n = 4Figure 2.2: GS quasiparticle weight (left panel) and GS average phononnumber (right panel) vs. ζ, in the quadratic (n = 2, lines) and quartic(n = 4, symbols) models, for various values of λ and Ω = 0.5t, in onedimension.402.4. Results and discussionare such that H(2)at = Ωatγ†γ + E(at)GS . We findΩat =√Ω(Ω + 4g2) (2.26)u =√(Ω + 2g2 + Ωat) /(2Ωat) (2.27)w = g1√Ω/Ω3at (2.28)v = sgn(g2)√(Ω + 2g2 − Ωat) /(2Ωat). (2.29)From these, we obtainEatGS = −g21ΩΩ2at+12(Ωat − Ω) (2.30)N(at)ph =12[Ω + 2g2Ωat− 1]+g21(Ω + 4g2)2(2.31)Zat =1uexp[−w2(1− vu)]. (2.32)The latter result requires the expansion of the squeezed coherent states inthe number state basis [66].Figure 2.3 shows Zat and N(at)ph vs. ζ (thick lines), which agree well withthe corresponding λ = 2 results of Fig. 2.2. In particular, for ζ → 0 we findΩat = Ω + 2g1ζ +O(ζ2) (2.33)N(at)ph =g21Ω2[1− 8g1Ωζ +O(ζ2)](2.34)explaining their linear increase/decreases for small ζ.The slight turnaround of the Z and Nph curves at larger values of ζis also observed in the atomic limit of the quadratic model. The reasonis that the first term in N(at)ph increases whereas the second term decreaseswith ζ. As discussed above, for small ζ the second term dominates andthe overall number of phonons decreases. For large ζ, however, the secondterm vanishes whereas the first term diverges as√g2 =√ζg1. Hence, as ζincreases Natph has a minimum, and then starts to increase with ζ. Basically,here the g2(b†2+b2) coupling dominates over the linear coupling g1(b†+b) and412.4. Results and discussion0 0.2 0.4 0.6 0.8 1ζ0.00.20.40.60.81.0ZatExactMFTMatching Nph0 0.2 0.4 0.6 0.8 1ζ02468Nph(at)ExactMFTMatching ZFigure 2.3: (left) Zat, and (right) N(at)ph vs. ζ, for g1 =√2 and Ω = 0.5 (fulllines). Dashed lines show the mean-field estimates, while the dot-dashedlines show the results of fitting g˜/Ω˜ to exactly reproduce the other quantity.See text for more details.422.4. Results and discussionchanges the trend. Physically, at small ζ the dominating effect is a narrowingof the harmonic potential, making phonons more costly to generate. At evenlarger ζ, the effect of shifting the origin to a new position dominates.This leads us to pose the question whether these exact results of thequadratic atomic model can be fit well by an effective linear model H(1)at =Ω˜b†b+ g˜(b†+b), for some appropriate choice of the effective parameters Ω˜, g˜.One way to achieve this is with a mean-field ansatzb†2 ≈ 2〈b†〉b† − 〈b†〉2, (2.35)with 〈b†〉 the GS expectation value of b†. The self-consistency condition〈b†〉 = −g1 + 2g2〈b†〉Ω + 2g2(2.36)leads to the mean-field estimatesΩ˜MF = Ω + 2g2 g˜MF = g1 − 2g1g2Ω + 4g2. (2.37)Thus, for small ζ = g2/g1, Ω˜MF increases whereas g˜MF decreases with in-creasing ζ so the effective coupling λ˜ = g˜2/(2dtΩ˜) decreases with ζ. This isconsistent with the observed move away from the small polaron limit withincreasing ζ. Quantitatively, however, these mean-field results (dashed linesin Fig. 2.3) are not very accurate for small ζ, and fail to capture even qual-itatively the correct behavior when ζ 1, since here N (at)ph → ∞ whileN(MF )ph = g˜2MF/Ω˜2MF → 0.In fact, there is no choice for effective linear parameters g˜ and Ω˜ thatreproduces the results of the quadratic model. This is because in the linearmodel, both Z˜ and N˜ph are functions of g˜/Ω˜ only. Fig. 2.3 shows that if onechooses this ratio so that N(at)ph = N˜ph, then Z˜ (dot-dashed line in the leftpanel) disagrees with Zat, and vice versa. Even more significant is the factthat even if one could find a way to choose g˜, Ω˜ so that the overall agreementis satisfactory for all GS properties, the linear model’s prediction for higherenergy features would still be completely wrong. For example, it would432.4. Results and discussionpredict the polaron+one-phonon continuum to occur at EGS + Ω˜ insteadof the proper EGS + Ω threshold. Since in the atomic limit the predictionsof the quadratic model cannot be reproduced with a renormalized linearmodel, we conclude that this must hold true at finite hopping t as well, atleast for large λ where the quadratic terms are important.So far we discussed moderate values of the adiabaticity ratio Ω/t = 0.5,as well as the anti-adiabatic (atomic) limit. MA predicts similar results inthe adiabatic limit Ω/t → 0 for large λ, where it remains accurate, but isunsuitable to study small and moderate couplings [13]. We expect that herethe quadratic coupling is essential even for small couplings λ → 0, becausethe term 2g2∑i b†ibi ensures that phonons are gapped even though Ω = 0.2.4.2 Negative nonlinear termSo far we also only discussed the case ζ > 0. The behavior of models withζ < 0 can be glimpsed at from the exact results in the atomic limit. Forsmall negative ζ, the results listed above show that the average phononnumber N(at)ph increases with |ζ| while the qp weight Zat decreases fast, i.e.the polaron moves more strongly into the small polaron limit. This is inagreement with the MA predictions for the quadratic model (not shown).Here, however, we must limit ourselves to values |ζ| < Ω/(4g1) so that Ωatremains a real quantity (a similar threshold is found for the full quadraticmodel. Note that the value of this threshold decreases with increasing λ).For values of |ζ| above this threshold the quadratic model becomes unsta-ble. This, of course, is unphysical. In reality, here one is forced to includehigher order (anharmonic) terms in the phonon Hamiltonian Hph since theyguarantee the stability of the lattice if the quadratic terms fail to do so.Such anharmonic terms may have little to no effect in the absence of thecarrier, but clearly become important in its presence, in this limit. Theycan be treated with the same MA formalism we used here. Their effects,as well as a full analysis of all possible signs of the non-linearities and theresulting polaron physics will be presented in the next chapter. For now, itis obvious that in the case ζ < 0, higher order terms in el-ph coupling also442.4. Results and discussionplay a key role in determining the polaron properties unless λ is very small,and therefore cannot be ignored.2.4.3 Small carrier concentrationsThe results presented so far clearly demonstrate the importance of non-linearel-ph coupling terms if the linear coupling λ is moderate or large, throughtheir significant effects on the properties of a single Holstein polaron.A reasonable follow-up question is whether such dramatic effects are lim-ited to the single polaron limit or are expected to extend to finite carrierconcentrations. While the limit of large carrier concentrations cannot betackled with our approach, here we present strong evidence that quadraticterms are likely to be equally important at small but finite carrier concen-trations.Of course, for finite carrier concentrations one needs to supplement theHamiltonian with a term describing carrier-carrier interactions. The sim-plest such term is an on-site Hubbard repulsion HU = U∑i ni↑ni↓, andgives rise to the Hubbard-Holstein Hamiltonian. The linear version of thisHamiltonian has been studied extensively by a variety of numerical meth-ods [63]. In particular, for low carrier concentrations and focusing on thesmall polaron/bipolaron limit, the phase diagram has been shown to con-sist of three regions: (i) for large λ and small U , the deformation energyfavours the formation of on-site bipolarons, also known as the S0 bipolarons;(ii) increasing U eventually makes having two carriers at the same site tooexpensive, and the S0 bipolarons evolve into weakly-bound S1 bipolarons,where the two carriers sit on neighboring sites. The binding is now providedby virtual hopping processes which allow each carrier to interact with thecloud of its neighbor. However, at smaller λ and larger U this binding mech-anism is insufficient to stabilize the S1 bipolaron, and instead one finds (iii)a ground state consisting of unbound polarons.This phase diagram has been found numerically in 1D [67] and 2D [68]for the linear Hubbard-Holstein model. Some results in 3D have also becomeavailable very recently [69]. In 1D and 2D, the separation lines between the452.4. Results and discussion0 1 2 3 4 5λ05U/t0 1 2 3 4 5λ0102030U/t0 1 2 3 4 5λ0102030U/t0 1 2 3 4 5λ05U/tS0S1unbound(a) ζ = 0 (b) ζ = 0.01(d) ζ = 0.1(c) ζ = 0.05unboundunboundunboundS1S1S1S0S0 S0Figure 2.4: Estimate of the bipolaron phase diagram in 1D for Ω/t = 0.5and for different values of ζ, based on second order perturbation theoryin t. In all four panels, the solid lines show the transition from S0 (on-site) stable bipolarons to S1 (nearest-neighbor) stable bipolarons, while thedashed lines show the unbinding transition above which bound polarons areunstable. Note that panels (c) and (d) have a significantly rescaled y-axis.See text for more details.462.4. Results and discussionvarious phases are found to be close to those estimated using second orderperturbation theory in the hopping t, starting from the atomic limit [67, 68].This is expected since for large linear coupling λ, the results always convergetoward those predicted by the atomic limit.Since the quadratic Hamiltonian can be diagonalized exactly in theatomic limit, we use second order perturbation theory to estimate the lo-cation of the separation lines for various values of ζ > 0. The results areshown in Fig. 2.4. Panel (a) shows the rough phase diagram for ζ = 0, inagreement with the asymptotic estimates shown in Refs. [67, 68] (note thatthe definition of the effective coupling used in those works differs by vari-ous factors from our definition for λ). Panels (b)-(d) show a very significantchange with increasing ζ. Even the presence of an extremely small quadraticterm ζ = 0.01 moves the two lines to considerably lower U values, as shownin panel (b), while for ζ = 0.05 and 0.1, the bipolarons are stable only in avery narrow region with small values of U (note that the vertical axes arerescaled for panels (c) and (d)).The dramatic change with increasing ζ in the location of these asymp-totic estimates for the various bipolaron transitions/crossovers strongly sug-gests that non-linear el-ph coupling terms remain just as important in thelimit of small carrier concentrations as they have been shown to be in thesingle polaron limit. In particular, these results suggest that the presenceof non-linear el-ph coupling terms leads to a significant suppression of thephonon-mediated interaction between carriers, so that the addition of a smallrepulsion U suffices to break the bipolarons into unbound polarons (whoseproperties are also strongly affected by the non-linear terms, as alreadyshown).Since publication of the work reported in this chapter, research carriedout by Steve Johnston has used quantum Monte Carlo methods to studythe nonlinear Holstein model at finite carrier concentration [70, 71]. Theyfind that in a two-dimensional Hubbard model, Charge-Density-Wave corre-lations are dramatically suppressed, even for small non-linear terms. Theirfindings confirm our earlier remarks that non-linear terms remain importantat finite carrier concentration.472.5. Concluding remarks2.4.4 Other types of el-ph couplingThe Holstein model is the simplest example of a g(q) model, i.e. a modelwhere the electron-phonon interaction depends only on the momentum ofthe phonon. Physically, such models appear when the coupling to the latticemanifests itself through a modulation of the on-site energy of the carrier.The Fro¨hlich model is another famous example of g(q) coupling. Modelsof this type are found to have qualitatively similar behavior, with smallpolarons forming when the effective coupling increases. These small polaronsalways have robust clouds, with significant distortions of the lattice in theirvicinity. We therefore expect that non-linear terms become important forall such models at sufficiently large linear coupling.2.5 Concluding remarksTo summarize, we have shown that non-linear terms in the el-ph couplingmust be included in a Holstein model if the linear coupling is large enoughto predict small polaron formation, and that doing so may very significantlychange the results. We also argued that these changes cannot be accountedfor by a linear Holstein model with renormalized parameters. These resultsshow that we have to (re)consider carefully how we model interactions withphonons (more generally, with any bosons) in materials where such inter-actions are expected to be strong, at least for models where this couplingmodulates the on-site energy of the carriers.48Chapter 3Quadratic Holstein modelIn this chapter, we show that in crystals where light ions are symmetricallyintercalated between heavy ions, the electron-phonon coupling for carrierslocated at the light sites cannot be described by a Holstein model. Weintroduce the double-well electron-phonon coupling model to describe themost interesting parameter regime in such systems, and study it in the singlecarrier limit using the momentum average approximation. For sufficientlystrong coupling, a small polaron with a robust phonon cloud appears at lowenergies. While some of its properties are similar to those of a Holsteinpolaron, we highlight some crucial differences. These prove that the physicsof the double-well electron-phonon coupling model cannot be reproducedwith a linear Holstein model.3.1 IntroductionWe recall that in most standard electron-phonon coupling models, it is as-sumed that the electron-phonon coupling is linear in the lattice displace-ments. This is a natural assumption because if the displacements are small,the linear term is the most important contribution. However, the coefficientof the linear term of the Holstein-like local electron phonon coupling mayvanish due to symmetries of the crystal. In such cases, the most importantcontribution is the quadratic term.Here we introduce, motivate and study in detail a Hamiltonian describingsuch quadratic electron-phonon (e-ph) coupling relevant for many commoncrystal structures, consisting of intercalated sublattices of heavy and lightatoms. We focus on the single carrier limit and the parameter regime wherethe carrier dynamically changes the effective lattice potential from a single-493.2. Modelwell to a double-well; hence, we call this the double-well e-ph coupling. Weuse the momentum-average approximation[1, 12] to compute the propertiesof the resulting polaron with high accuracy. We find that although thepolaron shares some similarities with the Holstein polaron, it also differsin important aspects. Indeed, we show that the physics of the double-welle-ph coupling model cannot be described by a renormalized linear Holsteinmodel.To the best of our knowledge, this is the first systematic, non-perturbativestudy of such a quadratic model. In the previous chapter and in [62] westudied the effect of quadratic (and higher) corrections added to a linearterm. Weak, purely quadratic coupling was studied using perturbation the-ory in Refs. [58, 59]. Other works considered complicated non-linear latticepotentials and couplings but treated the oscillators classically,[72, 73, 61]or discussed anharmonic lattice potentials but for purely linear coupling[74, 75]. Away from the single-carrier limit, the Holstein-Hubbard modelin infinite dimensions was shown to have parameter regions where the ef-fective lattice potential has a double-well shape;[76–78] this was then usedto explain ferroelectricity in some rare-earth oxides. See [79] and referencestherein. However, the effect of a double-well e-ph coupling on the propertiesof a single polaron were not explored in a fully quantum-mechanical modelon a low-dimensional lattice.This work is organized as follows: in Section 3.2 we introduce the Hamil-tonian, motivate its use for relevant systems, and discuss all approximationsmade in deriving it. In Section 3.3 we review the theoretical means bywhich we study our Hamiltonian. In Section 3.4 we present our results, andin Section 3.5 we give our concluding discussion and an outlook for futurework.3.2 ModelThe crystal structures of interest are illustrated in Fig. 3.1(a) for 1D, andFig. 3.1(b) for 2D cases. The 3D crystal would have a perovskite structurebut we do not discuss it explicitly because, as we show below, dimensionality503.2. Modelionic potentialno extra carrier with extra carrierhopping(a) (b)Figure 3.1: Sketch of the crystal structures discussed in this work: (a) 1Dchain, and (b) 2D plane, consisting of light atoms (filled circles) interca-lated between heavy atoms (empty circles). In the absence of carriers, theionic potential of a light atom is a simple harmonic well. In the presenceof a carrier, the ionic potential of the light atom hosting it remains an evenfunction of its longitudinal displacement, so the linear e-ph coupling van-ishes. In suitable conditions the effective ionic potential becomes a doublewell (see text for more details).plays no strong role in determining the polaron properties.The undoped compound is an insulator made of light atoms, shown asfilled circles, intercalated between heavy ones, shown as empty circles. Tozeroth order, the vibrations of the heavy atoms can be ignored while thoseof the light atoms are described by independent harmonic oscillators Hph =Ω∑i b†ibi , where bi annihilates a phonon at the ith light atom. (We set themass of the light ions M = 1, and also h¯ = 1). In reality there is weakcoupling between these oscillators giving rise to a dispersive optical phononbranch. However, the dispersion can be ignored if its bandwidth is smallcompared to all other energy scales. We do so in the following.Consider now the addition of a carrier. If it occupies orbitals centeredon the heavy atoms, its coupling to the oscillations of the light atoms isdescribed by breathing-mode coupling models [80, 81, 14, 82]. Here weare instead interested in the case where the carrier is located on the lightatoms. Such is the situation for a CuO2 plane as shown in Fig. 3.1(b),since the parent compound is a charge-transfer insulator[23] so that upondoping, the holes reside on the light O sites (of course, there are additional513.2. Modelcomplications due to the magnetic order of the Cu spins; we ignore thesedegrees of freedom in the following). The carrier moves through nearest-neighbor hopping between light atoms: Tˆ = −t∑〈i,j〉 (c†icj + h.c.), whereci is the carrier annihilation operator at light atom i.Given the symmetric equilibrium location of the light ion hosting thecarrier between two heavy ions, it is clear that the on-site e-ph couplingcannot be linear in the displacement δxi of that light ion: the sign of thedisplacement cannot matter. Thus, e-ph coupling in such a material is notdescribed by a Holstein model. This assertion is supported by detailed mod-elling. For simplicity, we assume that the interactions with the neighboringheavy atoms are dominant (longer-range interactions can be easily includedbut lead to no qualitative changes). There are, then, two distinct contribu-tions to the e-ph coupling:Electrostatic coupling: The carrier changes the total charge of the lightion it resides on. If the distance between adjacent light and heavy ions isd, and if U(x) is their additional Coulomb interaction due to the carrier,then the potential increases by U(d + δxi) + U(d − δxi). This is an evenfunction and thus has no linear (or any odd) terms in δxi. The coefficientof the quadratic term (δxi)2 can be either positive or negative, dependingon the charge of the carrier (electron or hole).Hybridization: Even though charge transport is assumed to take placein a light atom band, there is always some hybridization tlh allowing thecarrier to hop onto an adjacent heavy ion. If ∆ is the corresponding energyincrease, assumed to be large, then the carrier can lower its on-site energyby −t2lh/∆ through virtual hopping to a nearby heavy ion and back. Thehopping tlh depends on the distance between ions; for small displacementstlh(δx) ≈ tlh(1 + αδx) where α is some material-specific constant. Becausethe light ion is centered between two heavy ions, such contributions addto−t2lh∆[(1 + αδx)2 + (1− αδx)2] = −2t2lh∆ [1 + α2(δx)2] . The potential isagain even in δx. In this case, the coefficient of the quadratic term is alwaysnegative.523.2. ModelGiven that δxi ∼ bi + b†i , it follows that the largest (quadratic) contri-bution to the e-ph coupling for such a crystal has the general form:H(2)e-ph = g2∑ic†ici(bi + b†i)2where all prefactors have been absorbed into the energy scale g2, and thesum is over all light ions. From the analysis above we know that g2 mayhave either sign.Physically, H(2)e-ph shows that the presence of a carrier modifies the cur-vature of its ion’s lattice potential, and thus changes the phonon frequencyat that site from Ω to Ωat =√Ω2 + 4Ωg2. If g2 > 0 then Ωat > Ω, makingphonon creation more costly. As we show in Appendix C.3, this leads to arather uninteresting large polaron with very weakly renormalized properties.This is why in the following we focus on the case with g2 < 0.We note at this point that what we have shown is that the linear termof a Holstein-like model vanishes in these symmetric cases. The same is notnecessarily true for more complicated models of electron-phonon coupling.As a simple example, if we include coupling of the carrier on one site tophonons of the neighboring light atoms, the coupling will not be symmetricin the neighboring sites’ displacements and thus retains a linear term. Suchlinear terms arise whether this coupling is electrostatic in nature or due tohybridization (see discussion above). These terms, however, are not includedin the original, linear, Holstein model as it is intended for cases where theon-site terms are dominant. In this case, then, it is valid for us to also notinclude them in the non-linear Holstein model. These longer-range termsare likely to be smaller than the on-site terms we keep given that all theseinteractions decrease quickly with distance. In particular, we do not expectthem to change the double-well physics we discuss.For sufficiently negative g2, Ωat vanishes or becomes imaginary, i.e. thelattice is unstable. This is unphysical; in reality the bare ionic potentialcontains higher order terms that stabilize the lattice. This means that forg2 < 0 we must include anharmonic (quartic) terms in the phonon Hamil-533.2. Modeltonian and, for consistency, also in the e-ph coupling, so thatHph = Ω∑ib†ibi + Θ∑i(b†i + bi )4H(4)el-ph =∑n∈{2,4}gn∑ic†ici (b†i + bi )n,where Θ is the scale of the anharmonic corrections. In physical situationsΘ Ω and 0 < g4 |g2|, or the Taylor expansions would not be sensiblestarting points.The anharmonic terms in Hph make the total Hamiltonian unwieldy, be-cause the phonon vacuum |0〉 is no longer the undoped ground-state, andthe new undoped ground state |0˜〉 has no simple analytical expression. Inorder to be able to proceed with an analytical approximation, we argue thatthese terms can be absorbed into the e-ph coupling; this is a key approxi-mation of the model. The reasoning is as follows: At those lattice sites thatdo not have a carrier, the quartic terms have little effect if θ Ω. Thisstatement is verified by exact diagonalization of Hph. Results are shownin Fig. 3.2 where we plot the overlap O = |〈0|0˜〉|2 (per site) between theundoped ground-states with and without anharmonic corrections, as well asthe average number of phonons at a site of the undoped lattice. Even forunphysically large values Θ/Ω ∼ 1, the overlap O remains close to 1 whileNph 1, showing that the undoped ground-state has not changed signifi-cantly in the presence of anharmonic corrections. From now we ignore thesecorrections at sites without an additional carrier.However, for sites that have a carrier present, we cannot ignore theanharmonic term: As discussed, it is crucial for stabilizing the lattice. Sincethis term is similar to the quartic term in the e-ph coupling, they can both begrouped together, resulting in the approximate Hamiltonian for our crystal:H = Tˆ + Ω∑ib†ibi + g2∑ic†ici(b†i + bi)2+ (g4 + Θ)∑ic†ici(b†i + bi)4(3.1)543.2. Model0 1 2Θ / Ω0.00.20.40.60.81.0O0 1 2Θ / Ω0.00.10.20.30.40.5Nph(a) (b)Figure 3.2: (a) Overlap between the undoped ground-states with and with-out anharmonic corrections, and (b) the average number of phonons per sitein the undoped system, due to anharmonic corrections, as a function of θ/Ω.553.2. Modelwith an effective quartic e-ph coupling term g4 + Θ, which from now on wewill simply call g4. This is the Hamiltonian that we investigate in this work.Before proceeding, let us review what we are neglecting when we discardthe anharmonic corrections at the unoccupied sites. Besides ignoring thechange in the undoped ground state from |0〉 to |0˜〉 (which is a reasonableapproximation if θ/Ω 1, as discussed), we also assume that only the e-phcoupling can change the number of phonons in the system, whereas in thefull model the phonon number is also changed by anharmonic corrections.This latter approximation is valid if the timescale for anharmonic phononprocesses τ4 ∼ 1/Θ is much longer than the characteristic polaron timescaleτp ∼ m∗/(mt), where m∗ is the effective polaron mass.Let us briefly summarize the basic properties of the lattice potential,which equals Ve(δx) = Ω2(δx)2/2 for sites without an extra carrier, andVc(δx) = Ω2at(δx)2/2 + 4Ω2g4(δx)4 for sites with one carrier. If g2 > −Ω/4,the first term describes a harmonic well with frequency Ωat and Vc(δx) de-scribes a single well centered at δx = 0. If g2 < −Ω/4, however, Ωat becomespurely imaginary. In this case, Vc(δx) becomes a double-well potential with alocal maximum at δx = 0. The two wells are centered at ±xeq = ±√−Ω−4g216Ωg4.For δx ≈ ±xeq we obtain Vc(δx) ≈ V (xeq) − Ω2at(δx ∓ xeq)2, which locallydescribes a harmonic well of frequency Ω2eff = −2Ω2at. Interestingly, this isindependent of g4, whose only role is to control the location and depth ofthe two wells (they are further apart and deeper for smaller g4).3.2.1 Summary of model assumptionsAs in the previous chapter, we will summarize here the assumptions thatwent into the model.First, quartic terms of the free-phonon part of the Hamiltonian got ab-sorbed into the el-ph coupling. We justify this by explicitly studying theeffect of quartic terms on the phonons and showing that they do not signifi-cantly change the free-phonon behaviour, especially in the physical situationof Θ Ω.Second, it is assumed that the timescale of these anharmonic processes563.3. Formalismis much slower than the polaronic timescale. We justify this as follows.Assuming an overall small anharmonic contribution Θ, a light polaron willcertainly move on a faster timescale than 1/Θ. A heavy polaron, on theother hand, implies that the carrier cannot travel far from the phonon cloud;thus, anharmonic effects predominantly occur when the carrier is present,and again the effect of free-phonon anharmonic terms can be neglected.Finally, we assume that models exist where the quadratic coupling termcan be sufficiently negative. Since the hybridization term is always negative,we show here a back-of-the-envelope calculation for the Coulomb term. Forour example, we assume a chain of hydrogen ions intercalated between someunspecified heavier ions with charge +1, and that the charge carrier is anelectron. Assuming that longer range interactions are screened, the Coulombpotential experienced by the extra charge carrier is then given byU(δx) =−e24pi0[1a+ δx+1a− δx].Here we set a ≈ 1A˚, which is slightly larger than the length of a hydrogenmolecule. Expanding the potential to second order in δx provides us withan estimate for g2. Assuming a phonon frequency of Ω = 0.5eV , just as inthe previous chapter, gives a zero-point amplitude of A ∼ 0.1A˚ and finallyleads to g2 ≈ −0.14eV , so the condition of g2 > −Ω/4 is satisfied. Again,parameters will be slightly different in actual materials, but not by ordersof magnitude. It is thus very plausible that the double-well electron-phononcoupling is realized in some materials.3.3 FormalismWe want to find the single particle Green’s functionG(k, ω) = 〈0|ckGˆ(ω)c†k|0〉,where Gˆ(ω) = [ω − H + iη]−1 is the resolvent of Hamiltonian (3.1). Fromthis, we can obtain all the polaron’s ground state properties as well as itsdispersion [12].Grouping terms in the Hamiltonian according to how they affect the573.3. Formalismphonon number, we rewriteH = H0 +Hp +H2 +H4 (3.2)withH0 = Tˆ + Ω∑ib†ibi + g2 + 3g4 (3.3)Hp =∑inib†ibi(2g2 + 6g4 + 6g4b†ibi) (3.4)H2 =∑ini[(g2 + 6g4)(b†,2i + b2i ) + 4g4(b†,3i bi + b†ib3i )](3.5)H4 = g4∑ini(b†,4i + b4i). (3.6)We note that H0 and Hp conserve phonon number whereas H2 and H4change it by ±2 and ±4, respectively. The constant g2 + 3g4 in H0 isabsorbed into ω in the following derivations, but plots of the spectral weightwill show actual energies.One important property of this Hamiltonian is that it preserves thephonon number parity on each site: because its terms only change the num-ber of phonons by multiples of two, any eigenstate is a sum of basis stateshaving only even (or only odd) number of phonons. The Hilbert space canthus be divided into an even and an odd (phonon number) sector, whichcan be diagonalized separately. We emphasize that this symmetry is dif-ferent from the parity symmetry under a global lattice inversion r → −r.The latter has been studied extensively for the linear Holstein model [83],where it was shown that polaron states with total momentum K = 0, pi havewell defined (spatial) parity. The phonon number parity, on the other hand,corresponds to a unitary transformation b†i → −b†i , i.e., a local inversion ofthe phonon coordinates. The number parity symmetry also correlates withthe local spatial parity of the ions, since the spatial parity operator for sitei can be written as Pˆi = exp(ipib†ibi).583.3. Formalism3.3.1 The even sectorWe compute the Green’s function via the same continued matrix fractionsmethod[2] previously used by us to compute the Green’s function of a gen-eralized Holstein model with linear and higher-order terms[62] within theframework of the momentum average (MA) approximation. This approxi-mation was shown to be highly accurate for models with Holstein coupling[1, 12]. The reasons for this (such as obeying exact sum rules) can be verifiedto hold for this model, too. To be specific, here we implement the MA(2)flavor which allows us to also locate the continuum lying above the polaronband. Recall from the introduction that this version of MA means that theexact form of the free propagator is used for those orders of the Dyson serieswhere there are at most 2 phonons in the system [13].We begin our derivation by dividing the Hamiltonian into H = H0 +H1 with H1 = Hp + H2 + H4. Using Dyson’s identity Gˆ(ω) = Gˆ0(ω) +Gˆ(ω)H1Gˆ0(ω), where Gˆ0(ω) = [ω −H0 + iη]−1, we obtainG(k, ω) = G0(k, ω)[1+∑ieikRi√N(g2 + 6g4)F1(k, ω; i, i) + g4F2(k, ω; i, i)](3.7)with Fn(k, ω; i, j) = 〈0|ckGˆ(ω)ci(b†i )2n−2(b†j)2|0〉 being the generalized prop-agator for a system with 2n phonons in total, 2n− 2 of them on site i withthe other two on site j. The difference between MA(2) and the original MA,which we also call MA(0), is that for F1 we also use its exact equation ofmotion (EOM),F1(k, ω; i, j) = G(k, ω; j)G0(j − i, ω − 2Ω)(2g2 + 12g4)+ F1(k, ω; j, j)G0(j − i, ω − 2Ω)(4g2 + 36g4)+ 8g4F2(k, ω; j, j)G0(j − i, ω − 2Ω)+∑lG0(l − i, ω − 2Ω) [F2(k, ω; l, j)(g2 + 6g4) + F3(k, ω; l, j)g4] . (3.8)593.3. Formalismwhich is obtained by applying Dyson’s identity again, and introducingG(k, ω; j) = 〈0|ckGˆ(ω)c†j |0〉 , G0(j − i, ω) = 〈0|cjGˆ0(ω)c†i |0〉 .The equations of motion for the Fn propagators with n ≥ 2 are approximatedby replacing the free propagator G0(j−i, ω−2nΩ)→ δi,j g¯0(ω−2nΩ), whereg¯0(ω) =1N∑kG0(k, ω) is the momentum averaged free propagator. At lowenergies this is a good approximation because G0(j − i, ω − 2nΩ) decaysexponentially with the distance |j − i| if ω − 2nΩ < −2dt in d dimensions.This is also justified by the variational meaning of the MA approximationsalready discussed. Essentially, MA(2) assumes that all phonons in the cloudare at the same site but also allows for a pair of phonons to be created at asite away from the cloud.The resulting EOMs are different depending on whether i = j or i 6= j.If we define F=n (k, ω; i) = Fn(k, ω; i, i) and F6=n (k, ω; i, j) = Fn(k, ω; i, j) fori 6= j, we obtainF=n (k, ω; i) =g¯0(ω − 2nΩ)[F=n−2(2n)4¯g4 + F=n−1((g2 + 6g4)(2n)2¯ + 4g4(2n)3¯)+(4ng2 + 12ng4 + 24n2g4)F=n + (g2 + 6g4 + 8ng4)F=n+1 + g4F=n+2]. (3.9)F 6=n (k, ω; i, j) = g¯0(ω − 2nΩ)[g4(2n− 2)4¯F 6=n−2+((g2 + 6g4)(2n− 2)2¯ + (2n− 2)3¯ · 4g4)F 6=n−1+[2(2n− 2)g2 + 12(n− 1)g4 + 6(2n− 2)2g4]F 6=n +[g2 + 6g4 + 4(2n− 2)g4]F 6=n+1 + g4F 6=n+2](3.10)where we use the notation xn¯ = x!/(x−n)!. We also omitted the argumentsfrom the Fn appearing on the right hand sides, as they remain unchanged.These EOMs connect generalized Green’s functions Fn with Fn±1 andFn±2. We reduce this to a first order recurrence relation[62] by introducing603.3. Formalismvectors W=n = (F=2n, F=2n+1) and analogously for W6=n . Below, we write Wnwithout the index = or 6= for results that apply to both W=n and W 6=n . Byinserting the EOMs into the definition of Wn, we obtain a matrix EOM forthe Wn,γnWn = αnWn−1 + βnWn+1. (3.11)The coefficients of these matrices are read off from the EOM for the Fn.They are listed in appendix C.1.1.Using the fact that limn→∞An = 0 we can show[62] thatWn = AnWn−1, with An = [γn − βnAn+1]−1 αn.By introducing a suitably large cut-off N where we set AN+1 = 0, we cancompute all An with n ≤ N as continued matrix fractions. Knowledge ofA1 allows us to express F2 and F3 in terms of F1 and F0 = G. Following aseries of steps presented in appendix C.1.2, we obtain a closed equation forF1 in terms of G, which we then finally use to compute G. The end resultof these manipulations is the self energyΣ(ω) =(g2 + 6g4 +A=1 |12g4)g˜0(ω)a=01− g˜0(ω)(a=1 − a6=)+ g4A=1 |11.with g˜0(ω) = g¯0(ω− 2Ω−a 6=) and the other coefficients defined in appendixC.1.2. The independence of the self-energy on momentum is the consequenceof the local form of the coupling and of the non-dispersive phonons, similarto the MA results for the Holstein model [13]. Momentum-dependence wouldbe acquired in a higher flavor of MA, but is likely to be weak. Finally, theGreen’s function is:G(k, ω) =1ω − k − Σ(ω) + iη . (3.12)One can now use the matrices An to generate the generalized propagatorsFn, which allow one to reconstruct the entire polaron wavefunction (withinthis variational space) [84]. For the quantities of interest here, however, thesingle-particle Green’s function suffices.613.4. Results3.3.2 The odd sectorHere we calculate the Green’s function for a state that already has a phononin the system. Since the phonon number can only change by 2 or 4, thissingle phonon can never be moved to another site, so it is natural to computethe Green’s function in real space. The most general such real space Green’sfunction is:Gijl(ω) = 〈0|bl cjGˆ(ω)c†ib†l |0〉 .Applying the Dyson identity leads to the EOMGijl(ω) = G0(j − i, ω − Ω)+∑i′G0(i′ − i, ω − Ω) 〈0|bl cjGˆ(ω)H1c†i′b†l |0〉 .We then split the sum over all lattice sites into a term i′ = l where theelectron is on the same site as the extra phonon, and a sum over all theother sites. The subsequent steps are very similar to those for the even-sector Green’s function. We summarize them in Appendix C.2, where wealso discuss how various propagators that enforce translational symmetry –i.e. propagators defined in momentum space – can be obtained from thesereal-space Green’s functions.The end result for the real-space Green’s functions is Gijl(ω) = G0(j −i, ω˜) + G0(l − i, ω˜)G0(j − l, ω˜)(a=o − a6=o )[1 − g¯0(ω˜)(a=o − a6=o )]−1 where ω˜ =ω − a 6=o − Ω. The coefficients a=o and a 6=o are listed in appendix C.2.3.4 Results3.4.1 Atomic limit: t = 0We begin our analysis with the atomic limit since it is a good starting pointfor understanding the properties of the small polaron, which is the moreinteresting regime. However, we note an important distinction between theHolstein model and our double-well model. In the former, the atomic limit isthe infinite-coupling limit. In the latter, g4 sets an additional energy scale.623.4. ResultsThus, the atomic limit is not the same as the strong coupling limit; thelatter also requires that g4/|g2| be small.Before doing any computations, we can describe some general featuresof the spectrum. As already discussed, the phonon component of the wave-functions has either even or odd phonon number parity. Since this is due tothe spatial symmetry in the local ionic displacement, in any eigenstate theion is equally likely to be found in either well. As usual, the ground statehas even symmetry since it has no nodes in its wavefunction. Subsequenteigenstates always have one more node than the preceding eigenstate, sostates with even and odd parity alternate. The exception is the limit ofinfinite well separation, g4/|g2| → 0+, where the 2nth and 2n + 1st eigen-states become degenerate. The system can then spontaneously break parityto have the ion definitely located in the left or in the right well, like in aferroelectric. For a finite g4 this is not possible in the single carrier limit,but it can be achieved at finite carrier concentration through spontaneoussymmetry breaking.As discussed, our results are obtained with MA. In the atomic limit MAis exact[1] because for t = 0 the free propagator is diagonal in real-space sothe terms ignored by MA vanish. Thus, MA results must be identical here tothose obtained by other exact means. To check our implementation of MA,we used exact diagonalization (ED) with up to a few thousand phonons;this suffices for an accurate computation of the first few eigenstates. EDand MA results agree, as required.Figure 3.3 shows the ground-state quasiparticle weight Z (the overlapbetween the polaron ground-state and the non-interacting carrier ground-state), and the ground-state average number of phonons in the cloud, Nph, asa function of g2 < 0, for various values of g4. Z has an interesting behavior.At g2 = 0 it is slightly below 1 because of the quartic terms. As |g2| isincreased, Z first rises towards a value close to 1 and then sharply drops.This turnaround is caused by the terms that involve both g2 and g4, i.e.(2g2+6g4)∑i nib†ibi fromHp and (g2+6g4)∑i ni(b†,2i +b2i ) fromH2. Startingfrom g2 = 0 and making it increasingly more negative will at first decreasethese coefficients, thereby renormalizing the ground state less. For even633.4. Results0 1 2|g2|0.00.20.40.60.81.0Zg4 = 0.1g4 = 0.05g4 = 0.02g4 = 0.010 1 2|g2|0510152025Nph(a) (b)Figure 3.3: Polaron ground-state properties in the atomic limit, for severalvalues of the g4: a) quasiparticle weight, and b) average number of phononsin the phonon cloud. Other parameters are Ω = 0.5, t = 0.643.4. Resultsmore negative g2, however, Z decreases sharply as the absolute value of thesecoefficients increases; this is paralleled by a strong increase in Nph. Based onthis argument, the peak in Z should occur for −6g4 < g2 < −3g4, which isindeed the case. The strong-coupling limit of a small polaron (correspondingto small Z, large Nph values) is therefore reached either by increasing |g2|or by lowering g4.While this allows us to conclude that in the atomic limit the crossoverinto the small polaron regime occurs at g23g4 ≈ −1.5, it also illustrates thedifficulty in defining an effective coupling for this model. For the Hol-stein model, the dimensionless effective coupling λ is the ratio between theground-state energies in the atomic limit and in the free electron limit; thecrossover to the small polaron regime occurs at λ ∼ 1. For the double-wellmodel the introduction of an effective coupling is not as straightforward,because the atomic limit has vastly different properties depending on theratio g2/g4, so comparing the energy in this limit to that of a free electronis not sufficient. (Moreover, there is no analytic expression for the groundstate energy of the double well potential). For these reasons, we continue touse the bare coupling parameters g2 and g4 to characterize our model.For strong coupling, we can accurately estimate the ground state energyby using the barrier depth and effective harmonic frequency of the double-well potential, E0,sc = Vc(xeq) + Ωeff/2. Fig. 3.4 shows the relative errorof this estimate, which indeed decreases as parameters move deeper intothe small polaron regime. Since here the tunnelling between the two wellsalso becomes increasingly smaller, one may think that we can describe thisregime accurately by assuming that the carrier becomes localized in one ofthe wells (thus breaking parity), i.e. that we can approximate the full latticepotential as being a single harmonic well centered at either xeq or −xeq. Ofcourse, the latter situation can be modelled with a linear Holstein model.It turns out that this is not the case. In the standard Holstein model,the charge carrier cannot change the curvature of the lattice potential andthus cannot account for the difference between Ω and Ωeff. To accountfor the change in the curvature of the well, one would have to consider atleast a Holstein model with both linear and quadratic e-ph coupling terms.653.4. Results0 1 2|g2|10-310-210-11001011021- E0,sc/E0g4 = 0.1g4 = 0.05g4 = 0.02g4 = 0.01Figure 3.4: Relative error in the ground state energy when computed in thesemiclassical approximation (see text for details). The coupling g2 < −Ω/4is restricted to values for which there is a double-well potential. Otherparameters are like in Fig. 3.3.663.4. ResultsAlthough it is possible to find effective parameters g1,eff, g2,eff and Ωeff sothat the resulting lattice potential in the presence of the carrier has the samelocation and curvature as one of the wells of the double-well potential, thecorresponding quasi-particle weight Zeff severely underestimates Z. Thisis because the single well approximation severely overestimates the latticepotential at x = 0, thereby reducing the overlap between the ground state ofthe shifted well and that of the original well. We conclude that the double-well coupling cannot be accurately described by a (renormalized) Holsteincoupling even in this simplest limit.3.4.2 Finite hoppingWe focus on results from the even sector because it describes states accessibleby injecting the carrier in the undoped ground-state. The odd sector isaccessed only if the carrier is injected into an excited state with an oddnumber of phonons present in the undoped system; we briefly discuss thiscase at the end of the section.We begin by plotting the ground-state values of Z and Nph, for 1D and2D lattices, in Figs. 3.5 and 3.6 respectively. Since the MA self-energy islocal, the effective polaron mass m∗ = m/Z, where m is the free carrier mass;we therefore do not plot m∗ separately. Apart from t = 1, the parametersare like in Fig. 3.3. Note that the kinks in the Nph curves for g4 = 0.02 arenot physical; they arise from numerical difficulties in resolving the preciselocation of the ground state peak when Z → 0.Qualitatively, the polaron properties show the same dependence on g2 asin the atomic limit, but the shape and location of the turnarounds is slightlymodified: As one would expect, the presence of finite hopping counteractsthe formation of a robust polaron cloud and increases the quasi-particleweight Z for any given g2 and g4 when compared to the atomic limit.The results in one and two dimensions are strikingly similar. The 2DZ is slightly larger than the 1D Z, and Nph in 2D is slightly lower than in1D. This is expected because in higher dimensions, the polaron formationenergy is competing against a larger carrier kinetic energy. These results673.4. Results0 1 2|g2|0.00.20.40.60.81.0Zg4 = 0.2g4 = 0.1g4 = 0.05g4 = 0.020 1 2|g2|05101520Nph(a) (b)Figure 3.5: (color online) Polaron ground-state properties in one dimen-sion for various values of the quartic coupling term g4 as a function of thequadratic coupling g2: a) quasiparticle weight, and b) average number ofphonons in the phonon cloud. Other parameters are t = 1, Ω = 0.5t.683.4. Results0 1 2|g2|0.00.20.40.60.81.0Zg4 = 0.2g4 = 0.1g4 = 0.05g4 = 0.020 1 2|g2|05101520Nph(a) (b)Figure 3.6: Polaron ground-state properties in two dimensions for variousvalues of the quartic coupling term g4 as a function of the quadratic couplingg2: a) quasiparticle weight, and b) average number of phonons in the phononcloud. Other parameters are t = 1, Ω = 0.5t.suggest that dimensionality is not playing a key role for the double-wellmodel, similar to the situation for the Holstein model. This is why we didnot consider 3D systems explicitly.We now move on to discuss the evolution of the spectral weight A(k, ω) =− 1pi ImG(k, ω) with increasing |g2|, at a fixed value of g4. This is shown inFig. 3.7 for 1D, and in Fig. 3.8 for 2D. Because the evolution is again qual-itatively similar in the two cases, we analyze in more detail the 1D results.Here, at small quadratic coupling g2 = −0.5, we observe the appearanceof a polaron band below a continuum of states. This continuum begins atE0 + 2Ω, and consists of excited states comprising the polaron plus twophonons far away from it. (In our MA(2) approximation, the continuumactually begins at EMA(0)0 + 2Ω, not at EMA(2)0 + 2Ω, for reasons detailed inRef. [13]).693.4. Results0 π/2πk−4−202ω(a) g2 = −0.50 π/2πk−4−202ω(b) g2 = −10 π/2πk−4−202ω(c) g2 = −1.50 π/2πk−4−202ω(d) g2 = −2Figure 3.7: A(k, ω) in 1D, for g4 = 0.05, Ω = 0.5 and t = 1, for variousvalues of g2.703.4. Results(0,0) (π,π)−8−6−4−20246ω(π,0) (0,0)(a) g2 = −0.5(0,0) (π,π)−8−6−4−20246ω(π,0) (0,0)(b) g2 = −1(0,0) (π,π)−8−6−4−20246ω(π,0) (0,0)(c) g2 = −1.5Figure 3.8: A(k, ω) in 2D, for g4 = 0.05, Ω = 0.5 and t = 1, for variousvalues of g2.713.4. Results10-410-210010-410-210010-410-210010-410-2100A(ω)10-410-2100-10 -5 0 5ω10-410-2100|g2|=1.50|g2|=1.25|g2|=1.00|g2|=0.75|g2|=0.50|g2|=0.25Figure 3.9: Real-space diagonal spectral function Aiii(ω) at g4 = 0.05 forvarious values of (negative) g2 in one dimension and for Ω = 0.5. The y-axishas a logarithmic scale. The vertical bars indicate the position of Eeven0 +Ω.Note that due to the parity-preserving nature of the Hamiltonian there isno analog of the polaron+one-phonon continuum starting at E0 + Ω, whichis observed in all linear coupling models. Trying to mimic the results ofthe double-well coupling with a linear model will, therefore, lead to a wrongassignment for the value of Ω.At small |g2|, the polaron band flattens out just below the polaron+two-phonon continuum. With increasing |g2|, its bandwidth decreases as thepolaron becomes heavier, and additional bound states appear below thecontinuum. This is similar to the evolution of the spectrum of a Holsteinpolaron when moving towards stronger effective coupling [12]. However, asalready discussed, this does not mean that the two Hamiltonians can bemapped onto one another.723.4. ResultsFor completeness, let us also discuss some of the features of the oddsector. In particular, we focus on the local Green’s function Giii(ω), whichcan be written asGiii(ω) = g¯0(ω˜) +g¯0(ω˜)2(a= − a6=)1− g¯0(ω˜)(a= − a6=)with ω˜ = ω−Ω−a 6=. One can verify that a 6= equals the MA(0) self-energy forthe even sector, up to a shift by Ω of its frequency. The equation for Giii(ω)then shows that the odd sector spectral function comprises two parts: (i)the first term is just the momentum-averaged spectral function of the even-sector, shifted in energy by Ω due to the presence of the extra phonon. Onecan think of these as states where the even-sector polaron does not interactwith the extra phonon. This contribution therefore has weight starting fromE0 + Ω; (ii) the second part describes interactions between the polaron andthe extra phonon. An interesting question is whether these can lead to abound state, i.e. to a new polaron with odd numbers of phonons in its cloud.This question is answered in Fig. 3.9 where we plotAiii(ω) = − 1pi ImGiii(ω)for different values of |g2| and g4 = 0.05, Ω = 0.5, t = 1, in one dimension.The vertical bars indicate the position of E0 + Ω, where indeed a contin-uum begins, as expected from the previous discussion. At sufficiently strongcoupling |g2| we find a discrete bound state below that continuum, showingthat the polaron can bind the extra phonon. In fact, it is more proper to saythat the extra phonon (which is localized somewhere on the lattice) bindsthe polaron to itself and therefore localizes it. One can think of this as anexample of “self-trapping”, except here there is an external trapping agentin the form of the extra phonon.One might wonder whether this localized bound state in the odd sectorcould ever be at an energy below the polaron ground-state energy E0 ofthe even sector, i.e. become the true ground-state. This is not the case;as explained above, in the atomic limit the ionic states alternate betweeneven and odd symmetry. Introducing a finite hopping allows the polaron tofurther lower its energy by delocalizing, but this is only possible in the evensector. Thus, we always expect the even-sector polaron to have an energy733.5. Summary and discussionsbelow that of this localized state.As stated before, the two subspaces with even and odd phonon numberare never mixed, at least at zero temperature. At finite temperature, theextra charge is inserted not into the phonon vacuum but into a mixed statecontaining a number of thermally excited phonons. We therefore expectthe resulting spectral function to show features of both the even and oddsectors. To be more precise, some spectral weight should be shifted from theeven-sector spectral weight to the odd-sector spectral weight as T increasesand there is a higher probability to find one or more thermal phonons in theundoped state. We plan to study the temperature depend properties of thisdouble-well coupling elsewhere.3.5 Summary and discussionsHere we introduced and motivated a model for purely quadratic e-ph cou-pling, relevant for certain types of intercalated lattices, wherein the carrierdynamically changes the on-site lattice potential from a single well into adouble well potential. All the approximations made in deriving this modelwere analyzed. In particular, we argued that ignoring the anharmonic latticeterms at the sites not hosting the carrier should be a good approximation.However, a more in-depth numerical analysis might be needed to furthervalidate this assumption.We used the momentum average approximation to obtain the model’sground state properties and its spectral function in the single polaron limit,in one and two dimensions. We found that for sufficiently strong quadraticcoupling a small polaron forms. Although the polaron behaves somewhatsimilarly to the polaron of the linear Holstein model, the double-well modelcannot be mapped onto an effective linear model: apart from the differencein the location of the continuum in the even sector, the double-well modelalso has an odd sector that should be visible at finite T , and which is en-tirely absent in the Holstein model. This is due to the double-well potentialmodel’s invariance to local inversions of the ionic coordinate; this symmetryis not found in the Holstein model. The polaron in this odd sector is also743.5. Summary and discussionsqualitatively different from the Holstein polaron, in that it is localized nearthe additional phonon present in the system when the carrier is injected. Ofcourse, if the assumption of an Einstein mode is relaxed, then the phonon ac-quires a finite speed and this polaron would become delocalized, as expectedfor a system invariant to translations. However, this would still be quali-tatively different than a regular polaronic solution because this polaron’sdispersion would be primarily controlled by the phonon bandwidth, not thecarrier hopping.Our results suggest that researchers interpreting their measurementsfrom, e.g., angular-resolved photoemission spectroscopy, must carefully con-sider the nature of their system’s e-ph coupling: if they assume linear cou-pling where the lattice symmetry calls for a quadratic one, the parametersextracted from fitting to such models will have wrong values.While we have laid here the basis for a thorough investigation of theproperties of the double-well e-ph coupling model, much work remains tobe done. We believe that adjusting already existing numerical schemes suchas diagrammatic Monte Carlo to this model is straightforward and lookforward to a comparison of numerically exact results with our MA results.In addition, there are certain ranges of parameters for which MA is not well-suited, such as the adiabatic limit Ω→ 0 at weak coupling, or systems withfinite carrier densities. We anticipate that these regimes will be exploredwith a range of numerical and analytical tools, especially the finite carrierregime which should be relevant for modelling ferroelectric materials.We plan to extend our study of the double-well e-ph coupling beyondthe single-polaron limit. We deem especially interesting the parameter rangewhere the lattice potential remains a single well if only one carrier is present,but changes into a double well when a second charge is added. In this case,we anticipate the appearance of a strongly bound bipolaron while the singlepolarons are relatively light. Such states are not possible in the Holsteinmodel.Finally, extending our MA treatment to finite temperature should yieldinteresting insights into the interplay between the two symmetry sectorsrevealed by the spectral weight.75Chapter 4Bipolarons in the quadraticHolstein modelIn this chapter, we use the Momentum Average approximation (MA) tostudy the ground-state properties of strongly bound bipolarons in the double-well electron-phonon (el-ph) coupling model, whose single polaron solutionwas discussed in the previous chapter. We show that this model predicts theexistence of strongly bound yet lightweight bipolarons in some regions of theparameter space. This provides a novel mechanism for the appearance ofsuch bipolarons, in addition to long-range el-ph coupling and special latticegeometries.4.1 IntroductionAs pointed out in the previous chapters, when a charge carrier becomesdressed by a cloud of phonons, the quasi-particle that forms – the polaron– may have quite different properties from the free particle, such as a largereffective mass and renormalized interactions with other particles. One par-ticularly interesting effect of the latter is the formation of bipolarons, wherean effective attraction mediated by exchange of phonons binds the carrierstogether. If the binding is strong enough, the two phonon clouds mergeinto one, resulting in a so-called S0 bipolaron. Weaker binding, where eachpolaron maintains its cloud and the binding is mediated by virtual visitsto the other carrier’s cloud, is also possible and results in a S1 bipolaron[67, 68].The existence of bipolarons is interesting for many reasons. For in-stance, it has been suggested that Bose-Einstein condensation of bipolarons764.2. Modelmight be responsible for superconductivity in some high-Tc materials. Foran overview, see [85]. For this to occur, the bipolaron must be stronglybound so it can survive up to high temperatures. However, such strongbinding generally requires strong electron-phonon coupling. In most simplemodels of el-ph coupling such as the Holstein model,[11] this also resultsin a large effective mass of the bipolaron[67, 68] which severely reduces itsmobility and makes it likely to become localized by even small amounts ofdisorder.For this reason, much of the theoretical work on bipolarons is focused onfinding models and parameter regimes for which the bipolaron is stronglybound yet relatively light. So far, successful mechanism are based eitheron longer-range electron-phonon interactions[86–88, 69] or on special latticegeometries such as one-dimensional ladders or triangular lattices [89].Here we show that the recently proposed (short-range) double-well el-phcoupling model[90] also predicts the existence of strongly bound bipolaronswith relatively low effective mass in certain regions of the parameter space,thus revealing another possible mechanism for their appearance. Our studyuses the Momentum Average (MA) approximation,[1, 12, 62, 90] which wevalidate with exact diagonalization in an enlarged variational space. Sincein the single-particle case the dimensionality of the underlying lattice hadlittle qualitative impact, we focus here on the one-dimensional case.This work is organized as follows. In Section 4.2 we introduce the Hamil-tonian for the double-well model and in Section 4.3 we discuss the methodswe use to solve it. In Section 4.4 we present results for the bipolaron bindingenergy and effective mass, and in Section 4.5 we summarize our conclusionsand an outlook for future work.4.2 ModelThe double-well el-ph coupling model was introduced in the previous chap-ter, where its single polaron was studied. For ease of reference, we repeatsome of its motivation and introduction here.The model is relevant for crystals whose structure is such that a sublat-774.2. Modeltice of light ions is symmetrically intercalated with one of much heavier ions;the latter are assumed to be immobile. Moreover, charge transport occurson the sublattice of the light ions. An example is the one-dimensional inter-calated chain shown in Fig. 4.1(a). Another example is a two-dimensionalCuO layer, sketched in Fig. 4.1(b), where the doping holes move on the lightoxygen ions placed in between the heavy copper ions. In such structures,because in equilibrium each light ions is symmetrically placed between twoimmobile heavy ions, the potential felt by a carrier located on a light ionmust be an even function of that ion’s longitudinal displacement from equi-librium, i.e. the first derivative of the local potential must vanish. As aresult, the linear electron-phonon coupling is zero by symmetry, and oneneeds to consider the quadratic coupling. This is what the double-well el-phcoupling model does.Starting from the single-polaron Hamiltonian describing double-well el-ph coupling, introduced in Ref. [90], we add the appropriate terms for themany-electron problem to obtainH = Tˆ + Ω∑ib†ibi + U∑inˆi↑nˆi↓+ g2∑iσc†iσciσ(b†i + bi)2+∑ig(ni)4(b†i + bi)4. (4.1)Here, ciσ and bi are annihilation operators for a spin-σ carrier at site iand a phonon at site i, respectively. Tˆ describes hopping of free carrierson the sublattice of light ions in an intercalated lattice like that sketchedin Fig. 4.1. For simplicity, we consider nearest-neighbor hopping only,Tˆ = −t∑〈i,j〉,σ c†iσcjσ+h.c., although our method can also treat longer-rangefinite hopping [2]. The next two terms describe a single branch of disper-sionless optical phonons with energy Ω, and the Hubbard on-site Coulombrepulsion with strength U . The last two terms describe the el-ph couplingin the double-well model. As mentioned, in lattices like that sketched inFig. 4.1, the coupling depends only on even powers of the light-ion displace-ment δxˆi ∝ b†i +bi (the heavy ions are assumed to be immobile). As a result,784.2. Modelthe lowest order el-ph coupling is the quadratic term whose characteristic en-ergy g2 can have either sign, depending on modeling details. As discussed atlength in the previous chapter, the interesting physics occurs when g2 < 0so that the el-ph coupling “softens” the lattice potential. For sufficientlynegative g2 this renders the lattice locally unstable in the harmonic approx-imation and requires the inclusion of quartic terms in the lattice potential.For consistency, one should then also include quartic terms in the el-phcoupling. Under reasonable assumptions the quartic lattice terms can becombined with the quartic el-ph coupling term on sites hosting a carrierand ignored on all other sites. Because the resulting quartic term containscontributions from both the lattice potential and from the el-ph interac-tion, it should not be assumed to be linear in the carrier number, unlike thequadratic term which arises purely from el-ph coupling. Instead, we use thegeneral formg(ni)4 = g4 ·0, if ni = 01, if ni = 1α, if ni = 2where ni =∑σ c†iσciσ is the number of carriers on site i, and α is a constantbetween 1 and 2. Setting α = 2 assumes that quartic lattice effects arenegligible compared to the quartic el-ph terms, whereas α = 1 is the oppositeextreme. For the remainder of this article we set α = 1, so that g(1)4 = g(2)4 =g4. This case leads to stronger coupling, since a lower g4 results in deeperwells that are further apart,[90] and thus represents the parameter regime weare interested in. Physically, this describes the situation where the quarticlattice terms are much larger than the quartic el-ph coupling; however theyare still negligible compared to the quadratic lattice terms and therefore canbe ignored at sites without a carrier.4.2.1 Summary of model assumptionsIn addition to the assumptions that went into the quadratic Holstein modelfrom the previous chapter, our assumptions are as follows:794.3. Formalismionic potentialno extra carrier with extra carrierhopping(a) (b)Figure 4.1: Sketch of the crystal structures discussed in this work: (a) 1Dchain, and (b) 2D plane, consisting of light atoms (filled circles) interca-lated between heavy atoms (empty circles). In the absence of carriers, theionic potential of a light atom is a simple harmonic well. In the presenceof a carrier, the ionic potential of the light atom hosting it remains an evenfunction of its longitudinal displacement, so the linear e-ph coupling van-ishes. In suitable conditions the effective ionic potential becomes a doublewell. (Reproduced from the previous chapter)First, the Coulomb interaction between carriers is on-site only. This isa valid approximation in the presence of reasonable screening, since in theintercalated lattices we consider, the charge-carrying ions are separated bya large distance, with another ion in between.Second, we carry out our calculations for the strongest coupling casefor the quartic coupling term, which is assumed to contain contributionsfrom both the coupling and the free-phonon part. Using a different case willchange the results only quantitatively, not qualitatively.4.3 FormalismWe compute the bipolaron binding energy and effective mass using the mo-mentum average (MA) approximation [1, 12, 62, 90]. Since we are interestedin strongly-bound bipolarons which have a large probability of having bothcarriers on the same site, the version of MA used here is the variationalapproximation that discounts states where the two carriers occupy differ-ent sites. This results in an analytic expression of the two-particle Green’sfunction which is used to efficiently explore the whole parameter space. The804.3. Formalismaccuracy of this flavor of MA is verified by performing exact diagonalizationin a much larger variational subspace (details are provided below). In theregime of interest the agreement is very favorable, showing that the effortrequired to perform the analytical calculation for a flavor of MA describinga bigger variational space is not warranted.4.3.1 Momentum average approximationWe define states with both carriers at the same site, |i〉 = c†i↑c†i↓ |0〉, andstates of given total momentum k with both carriers at the same site,|k〉 = 1√N∑ieik·ri |i〉 .The bipolaron dispersion Ebp(k) is obtained from the lowest energy pole ofthe two-particle Green’s functionG(k, ω) = 〈k|[ω −H+ iη]−1|k〉 ,where η → 0+ is a small convergence factor. The effective bipolaron massis 1/mbp = ∂2Ebp/∂k2|k=0. Throughout this work we set h¯ = 1, a = 1.We split the Hamiltonian into H = H0 +H1 with H0 = Tˆ + Ω∑i b†ibidescribing the free system and H1 containing the interaction terms. Weapply Dyson’s identity Gˆ(ω) = Gˆ0(ω) + Gˆ(ω)H1Gˆ0(ω) whereGˆ0(ω) = [ω −H0 + iη]−1is the resolvent of H0 and we also defineG0(k, ω) = 〈k|Gˆ0(ω)|k〉 = 1N∑q1ω + iη − (k− q)− (q)as a non-interacting two-particle propagator, where (k) is the free carrierdispersion. N → ∞ is the number of light-ion sites of the lattice. In 1D,G0(k, ω) equals the momentum-averaged single-particle free propagator inone dimension for an effective hopping integral 2t cos(k/2), for which an814.3. Formalismanalytic expression is known [91]. In higher dimensions, such propagatorscan be calculated as discussed in Ref. [92].As mentioned, in a variational sense the MA used here amounts to ne-glecting all states where the carriers are not on the same site. This ap-proximation is justified for the description of the strongly bound on-site(S0) bipolaron, which is expected to have most of its weight in the sec-tor where both carriers are on the same site. Another way to look atthis is that the bipolaron ground-state energy in the strongly-bound casemust be well below the non-interacting two-particle continuum, and the freetwo-particle propagator will have vanishingly small off-diagonal matrix el-ements at such energies. Ignoring them, the equation of motion becomesG(k, ω) ≈ G0(k, ω) + 〈k|Gˆ(ω)H1|k〉G0(k, ω), and thus:G(k, ω) =G0(k, ω)(1 +∑ieikRi√N[(g2 + 6g4)F1(k, ω, i)+ g4F2(k, ω, i) + UF0(k, ω, i)]).where Fn(k, ω, i) = 〈k|Gˆ(ω)b†,2ni |i〉 is a generalized two-particle propagator.Equations of motion for the Fn propagators are obtained in the same way,and read:Fn(k, ω, i) = g¯0(ω − 2nΩ)[g4(2n)4¯Fn−2(k, ω, i)+((2g2 + 6g4)(2n)2¯ + 4g4(2n)3¯)Fn−1(k, ω, i)+ (8ng2 + 12ng4 + 24n2g4 + U)Fn(k, ω, i)+ (2g2 + 6g4 + 8ng4)Fn+1(k, ω, i) + g4Fn+2(k, ω, i)]. (4.2)where we use the shorthand notation xn¯ = x!/(x− n)! and have introduced824.3. Formalismthe momentum-averaged free two-carrier propagator,g¯0(ω) := 〈i|Gˆ0(ω)|i〉 = 1N∑kG0(k, ω)=1N2∑k,q1ω − (k− q)− (q) + iη .In 1D, g¯0(ω) equals the diagonal element of the free propagator for a particlein 2D, which can be expressed in terms of elliptical functions and calculatedefficiently [91]. Similar considerations hold in higher dimensions [92].The equations of motions are then solved following the procedure de-scribed at length in Refs. [62, 90]. For consistency, we sketch the main stepshere. First, we introduce vectors Wn = (F2n−1, F2n)T for n ≥ 0 (the ar-guments k, ω, i of the propagators are not written explicitly from now on).Note that with this definition, W0 = (F−1, F0), yet F−1 is not properlydefined. However, the final result has no dependence on F−1, as we showbelow. The equations of motion are then rewritten in terms of Wn to readγnWn = αnWn−1 + βnWn+1. The matrix elements of the 2 × 2 matricesαn, βn, γn, are easily read off Eq. (4.2).Defining An = [γn − βnAn+1]−1αn, the physical solution of these recur-rence equations is Wn = AnWn−1. Introducing a sufficiently large cut-offNc where WNc = 0, we can then compute A1 and have W1 = A1W0, i.e.,(F1F2)= A1(F−1F0)=(a11 a12a21 a22)(F−1F0).One can easily check that a11 = a21 = 0. Thus, we obtain F1 = a12F0 andF2 = a22F0. Substituting these results back into the EOM for G we obtainG(k, ω) = G0(k, ω)[1+∑ieikRi√N((2g2 + 6g4)a12 + g4a22 + U)F0(k, ω, i)].Since, by definition, G(k, ω) =∑ieikRi√NF0(k, ω, i), and given that a12, a22834.3. Formalismare functions of ω only, we find:G(k, ω) =1G−10 (k, ω)− (2g2 + 6g4)a12 − g4a22 − U. (4.3)Note that the coefficients a12 and a22 depend on all parameters of themodel, including U . As a result, the position of the lowest pole of Eq. (4.3)is not simply linear in U , although this is a good approximation for thestrongly bound bipolaron.We emphasize that this MA expression becomes exact in two limitingcases. First, in the atomic limit t→ 0 the free propagator has no off-diagonalterms and thus no error is introduced by dropping them from the equationsof motion. Second, without el-ph interactions (gn = 0) the Hamiltonianreduces to the Hubbard model which is exactly solvable in the two-particlecase [93]. In both cases MA gives the exact solution.4.3.2 Exact diagonalizationThe results obtained via MA as outlined above are checked against exactdiagonalization results in a bigger variational subspace designed to describewell the strongly bound S0 bipolaron. Hence, we only consider states whereall the phonons are located on the same lattice site and at least one of thetwo electrons is close to this cloud. The basis states are of the form|k, n, δ1, δ2〉 =∑ieikri√Nb†,ni c†i+δ1,↑c†i+δ2,↓ |0〉with the constraint that either δ1 or δ2 is below a certain cut-off. In addition,a global cut-off Nc is imposed on n+ δ1 + δ2. The ground state within thevariational space is then computed using standard eigenvalue techniques.The main difference between these ED and MA results is that MA dis-cards contributions from configurations where the carriers are at differentlattice sites. Comparing the two therefore allows us to gauge the impor-tance of such terms, and to decide whether the speed gained from using theanalytical MA expressions counterbalances the loss of accuracy.844.4. Results-1.5-1-0.50g2-15-10-5E bp-1.5-1.0-0.50.0g210-210-11002m/mbpEDMA(a) (b)Figure 4.2: (a) Bipolaron ground-state energy, and (b) inverse effective massfor t = 1,Ω = 0.5, g4 = 0.1, computed with ED (solid black line) and MA(red dots).4.4 ResultsFrom now on we focus on the one-dimensional (1D) case, since our resultsfrom the previous chapter suggests that going to higher dimensions leads toqualitatively similar results.Before discussing the MA results, we first compare them to those ob-tained from ED in the larger variational subspace discussed above. A typicalcomparison (for t = 1,Ω = 0.5 and g4 = 0.1) is shown in Fig. 4.2. The leftpanel shows the ground state energy and the right panel shows the inverseeffective mass of the bipolaron. In the regime where the bipolaron is stronglybound, i.e. where its energy decreases fast and its effective mass increasessharply as |g2| increases, we find excellent agreement for the energy. Themasses also agree reasonably well, but MA systematically overestimates thebipolaron mass. This is a direct result of the more restrictive nature of theMA approximation: By discarding configurations where the carriers occupy854.4. Resultsdifferent sites, the mobility of the bipolaron is underestimated and thus theeffective mass is overestimated. Nonetheless, this error is not very large, andonly means that the bipolarons in the double well model are even lighterthan calculated by MA. Due to similarly good agreement in all cases we ver-ified, for the remainder of this work we only discuss results obtained withthe more efficient MA method.We emphasize that our approximation for computing the Green’s func-tion is only valid in the regime of strong binding and does not describecorrectly the physics at weak coupling. Since neither MA nor ED, as im-plemented here, allow for the formation of two phonon clouds, neither de-scribes the formation of a weakly bound S1 bipolaron (where polarons formon neighboring sites and interact with each other’s clouds via virtual hop-pings), nor the dissociation into two polarons as the coupling is furtherdecreased [67, 68]. Accuracy in these parameter regimes can be improvedby applying more sophisticated – yet much more tedious – versions of MAor ED for suitably expanded variational spaces. Here, however, we want tofocus on the strong-coupling regime, where our results are accurate.We show the ground-state properties of the bipolaron compared to thoseof two single polarons in Figs. 4.3, 4.4 for two different values of Ω. In allthose panels, we have set U = 0 for simplicity; the role of finite U will bediscussed at the end of this section.The ground state energy of the bipolaron behaves qualitatively similarfor all values of Ω and g4 in that it shows a kink at some g2 where the slopebecomes steeper. This signifies the onset of the strong-coupling regime wherethe bipolaron energy is well below the energy of two independent polarons,consistent with a strongly bound bipolaron. At weaker coupling the resultsare not accurate since – as explained above – our version of MA cannotdescribe the dissociation of the bipolaron.Consider now the behavior of the effective masses. For all parametersconsidered here, we see that the single polaron mass mp starts out slightlyabove the free electron mass m, then decreases until it is almost as light asthe free electron, before increasing again. This turnaround in the polaronmass is due to partial cancellation effects of the quadratic and quartic el-864.4. Results-7-6-5-42EpEbp-15-10-5-30-20-10-1.5-1.0-0.50.0g20.00.20.40.60.81.0m/mp2m/mbp-1.5-1.0-0.50.0g2-1.5-1.0-0.50.0g2(a) g4= 0.2 (c) g4= 0.05(b) g4= 0.1Figure 4.3: Ground-state properties (total energy and inverse mass) of theS0 bipolaron and two independent polarons for t = 1,Ω = 0.5 and g4 = 0.2,0.1, and 0.05 for a), b), and c), respectively. For all panels, U = 0.ph coupling terms, as discussed in Ref. [90]. We observe that when thebipolaron is already quite strongly bound, the single polaron can still bevery light. Empirically, we find that in the strong-coupling regime mp ∼m exp(−γ∆p/Ω) where ∆p = −2t−Ep is the single-polaron binding energyand γ is a small numerical prefactor. This behavior is also found in theHolstein model[11] in the strong-coupling limit, where γ = 1,∆p = −g2/Ω.The prefactor γ can be much smaller in the double well model because ofthe nature of the ionic potential. This was explained in detail in Ref. [90],and will be discussed in the context of bipolarons later in this section.The bipolaron effective mass fluctuates around the value of 2m in theweak coupling regime. As explained above, here our method does not de-scribe two independent polarons, but two independent free electrons whoseeffective mass should just be 2m. However, the two-particle spectral func-tion in this case does not have a low-energy quasi-particle peak. Instead, ithas a continuum spanning the allowed two-particle continuum.874.4. Results-8-6-42EpEbp-20-15-10-5-30-20-10-1.5-1.0-0.50.0g20.00.51.0m/mp2m/mbp-1.5-1.0-0.50.0g2-1.5-1.0-0.50.0g2a) g4 = 0.2 b) g4 = 0.1 c) g4 = 0.05Figure 4.4: Ground-state properties (total energy and inverse mass) of theS0 bipolaron and two independent polarons for t = 1,Ω = 2 and g4 = 0.2,0.1, and 0.05 for a), b), and c), respectively. For all panels, U = 0.05101520∆ / Ωg4 = 0.2g4 = 0.1g4 = 0.05-1.5-1.0-0.5g2100101102mbp / 2mp0510-1.5-1-0.50g2100101102Ω=0.5 Ω=2Figure 4.5: Binding energy ∆ and effective-mass ratio mbp/2mp of the bipo-larons for Ω = 0.5 and Ω = 2 at different values of g4.884.4. ResultsTable 4.1: Some example values of the bipolaron binding energy and effectivemass.Ω g4 |g2| ∆/t mbp/2m0.5 0.1 0.9 1.25 8.30.2 1.3 1.48 4.42 0.1 1.3 3.11 5.90.2 1.5 1.03 1.8These issues disappear at stronger coupling where a strongly bound bipo-laron forms and MA becomes accurate. The figures show that here the bipo-laron quickly gains mass with increased coupling strength |g2|, and that thisincrease is stronger the smaller g4 is. Note that a smaller g4 actually meansstronger coupling, because the wells are deeper and further apart [90].The same data is displayed in a different way in Fig. 4.5, where weshow the magnitude of the bipolaron binding energy ∆ = 2Ep − Ebp andthe ratio of bipolaron to single-polaron masses, mbp/2mp. The stronglybound bipolaron regime (where the results are accurate) is reached whenthese quantities vary fast with g2. In particular, the results for mbp/2mpshow that here the bipolaron mass increases much more quickly than thepolaron mass. This is not surprising for models like this, where the phononsmodulate the on-site energy of the carrier. At strong coupling the resultscan be understood starting from the atomic limit t = 0, treating hopping asa perturbation. Since both carriers must hop in order for the bipolaron tomove, one expects that mbp/m ∝ (mp/m)2; indeed, we find this relation tobe valid for a wide range of parameters for our model.Although in this regime the bipolaron quickly gains mass, there are pa-rameter ranges where its mass is still rather light while the bipolaron isstrongly bound. Examples of such parameters are given in Table 4.1. Wenote that qualifiers such as ”strongly bound“ and ”light“ are subjective. Inour case, we take the bipolaron as strongly bound when the binding energy∆/Ω > 1 and the ratio mbp/2m < 10− 20, consistent with other references[89, 86].894.4. Resultssingle well double-wellpotentialionic wfFigure 4.6: Ionic potential (above) and ionic ground-state wavefunction (be-low) in the single-well and double-well models. Solid lines correspond to thesituation without an additional carrier, dashed lines to the situation withan additional carrier.Light but strongly bound bipolarons were previously found for long-range el-ph coupling [86]. The explanation is that in such models, carriersinduce a spatially extended lattice deformation, not one that is located inthe immediate vicinity of the carrier as is the case at strong coupling inlocal el-ph coupling models. Because of their extended nature, the overlapbetween clouds displaced by one lattice site (which controls the effectivehopping) remains rather large, meaning that the polarons and bipolaronsremain rather light in such models.Even though it is due to a local el-ph coupling, the mechanism resultingin light bipolarons in our model is qualitatively similar, as illustrated inFig. 4.6. In the linear Holstein model the effect of an additional carrier addedto a lattice site is to shift the equilibrium position of the ionic potential.The ionic wavefunctions corresponding to an empty and an occupied sitetherefore have only small overlap, which strongly reduces the effective carrierhopping. In the double-well model, in contrast, the ionic wavefunction forthe doubly-occupied site has appreciable overlap with the ionic wavefunctionfor an empty or a singly-occupied site and thus does not reduce the effectivehopping as much.We conclude with a brief discussion of the effects of a finite, repulsiveU . For a very strongly coupled S0 bipolaron, most of the weight is in states904.5. Conclusions and outlookwith both carriers on the same site. In this regime, the binding energydecreases (nearly) linearly with U , ∆bp(U) ≈ ∆bp(U = 0) − U . However,increasing U increases the energy cost of the S0 state and thus encourageshybridization with off-site states, which results in an overall smaller effectivemass. We show results for the bipolaron energy and effective mass as afunction of U in Fig. 4.7. We stay within the regime U < ∆bp where thebipolaron remains strongly bound. As predicted, the energy of the bipolaronincreases linearly with U , which in turn means that the binding energy ∆bpdecreases linearly with U . The effective mass also decreases (approximately)linearly with U . This can be demonstrated for the strong coupling limit viasecond order perturbation theory in the hopping. Following along the linesin Refs. [67, 68], the effective hopping of the S0 bipolaron is of the formm−1bp ∝ teff ∼−t2e−γ∆/Ω2Ep − Ufor some constant γ. We see that the mass itself decreases linearly with U ,with a steeper slope the larger the effective mass at U = 0.In essence, provided that it is not large enough to break the bonding, afinite U does not change the overall picture and merely tunes the balancebetween the bipolaron binding energy and its effective mass.4.5 Conclusions and outlookIn conclusion, we have investigated the bipolaron ground-state properties inthe dilute limit of the double-well el-ph coupling model at strong coupling.We have demonstrated that due to the particular nature of the carrier-induced ionic potential, the double-well bipolaron can be strongly boundwhile remaining light compared to the bipolaron in the Hubbard-Holsteinmodel. This suggests a new route to stabilizing such bipolarons, in additionto previously discussed mechanisms based on long-range el-ph coupling orspecial lattice geometries. We expect that a combination of these mecha-nisms will lead to even lighter bipolarons.In this work, we have used and validated a simple extension of the Mo-914.5. Conclusions and outlook0 1 2 3U-12-10-8-6-4E bp0 1 2 3U050100150200mbp / 2mg2 = -1.00g2 = -1.25Figure 4.7: Bipolaron energy (left) and effective mass (right) as a functionof the Hubbard U , for t = 1,Ω = 0.4 and g4 = 0.1. Similar results arefound for other parameters if U is not large enough to lead to bipolarondissociation.924.5. Conclusions and outlookmentum Average approximation to the two-carrier case. While this general-ization is appropriate to describe a strongly bound S0 bipolaron, it cannotdescribe the off-site (S1) bipolaron that forms at larger Hubbard repulsionU , or the unbinding of the bipolaron at even larger U . A more sophisticatedversion of MA, currently under development, will give us insight into thefull phase diagram of the double-well model.93Chapter 5Single-polaron dispersion inthe insulating limit oftetragonal CuOWe argue that tetragonal CuO (T-CuO) has the potential to finally set-tle one long-standing modelling issue for cuprate physics. We compare theone-hole quasiparticle (qp) dispersion of T-CuO to that of cuprates, in theframework of the strongly-correlated (Udd → ∞) limit of the three-bandEmery model. Unlike in CuO2, magnetic frustration in T-CuO breaks theC4 rotational symmetry and leads to strong deviations from the Zhang-Ricesinglet picture in parts of the reciprocal space. Our results are consistentwith angle-resolved photoemission spectroscopy data but in sharp contra-diction to those of a one-band model previously suggested for them. Thesedifferences identify T-CuO as an ideal material to test a variety of scenariosproposed for explaining cuprate phenomenology.5.1 IntroductionUnderstanding the high-temperature superconductivity in cuprates [22] isone of the biggest challenges in condensed matter physics. These layeredmaterials contain two-dimensional (2D) CuO2 layers which exhibit antifer-romagnetic (AFM) order in the undoped limit, and host the superconductingCooper pairs upon doping. Consequently, it is widely believed that under-standing the behaviour of a doped CuO2 layer is the key to understand theunusual properties of these materials.945.1. IntroductionThe first step in this quest is to understand the nature and dynamics ofthe quasiparticle (qp) that forms when one hole is doped into a CuO2 layer.Despite huge efforts on the theory side, this issue is not yet fully settled.The most relevant orbitals for the physics of the CuO2 layer are theCu 3dx2−y2 and the O 2p ligand orbitals, and the appropriate model isthe three-band Emery Hamiltonian [24]. Zhang and Rice argued, however,that the resulting quasiparticle is a Zhang-Rice singlet5 (ZRS) hopping onthe Cu sublattice and is thus well described by the much simpler one-bandt-J or Hubbard Hamiltonians [26–30]. A lot of effort focusing on these(relatively) simpler one-band models was thus to follow. In the absenceof exact solutions or accurate approximations for these strongly-correlated2D Hamiltonians, progress was made through numerical studies of finite-size clusters. These showed that the qp dispersion is strongly influencedby the quantum fluctuations of the AFM background [94], and that longer-range hopping is needed to obtain a dispersion in quantitative agreementwith what was measured experimentally [95–97, 16]. Similar results are ob-tained within Cluster Dynamical Mean Field Theory [31]. The longer-rangehoppings required to achieve this agreement are similar to those calculatedtheoretically [32–34]. This was taken as proof that these extended one-bandmodels provide the correct description, and the focus shifted to studyingthem at finite hole concentrations. Surprisingly, a lot of this work discardsthe longer-range hopping terms despite their proven relevance. While muchsuch work was done in the past two decades, the lack even of consensus thatthese one-band models support robust, high-temperature superconductivityraises strong doubts about how appropriate they are to describe the hole-doped side of the phase diagram. It is much more likely that the one-bandmodels (with properly adjusted longer-range hopping) describe correctly thephysics of the electron-doped side, because in this case the full O bands areinert spectators.There are two reasons why the one-band models might fail to find thedesired physics at finite doping: (i) they may describe the single qp cor-5We give a proper definition of the ZRS state later in this chapter when discussing thelow-energy effective model. A sketch is given in Fig. 5.3.955.1. Introductionrectly yet fail to appropriately model the effective interactions between qps,responsible for pairing. This was shown to occur when different degreesof freedom on different sublattices are mapped onto an effective single bandmodel [98, 99]. Because in cuprates the doped holes reside on oxygen whereasthe magnons reside on copper sites [100], a one-band model may thus fail tomimic the correct interaction between them; (ii) they may predict the cor-rect qp dispersion for the wrong reasons, by describing different physics evenat the single hole level. Support for the latter view comes from our recentwork on the Udd → ∞ limit of the three-band Emery model; the resultingHamiltonian has spins at the Cu sites and the doped holes move on the Osublattice [100]. In stark contrast to one-band models where spin fluctua-tions are key to obtaining the correct qp dispersion, here it is recovered evenin the absence of spin fluctuations [3]. This qualitative difference shows thatalthough the quasiparticles of these models have similar dispersion, this isdriven by different physics [5].To decisively settle the question of whether these one- and three-bandmodels are equivalent, one must compare them for a material similar enoughto CuO2 that it should be described by similar Hamiltonians, however onefor which they give different predictions so that (at least) one of them canbe falsified experimentally. In this Letter we show that layers of tetragonalCuO (T-CuO) are precisely such a material, whose careful investigation canfinally resolve these fundamental modeling issues.Thin films of several unit cells of T-CuO were recently grown epitaxiallyon a SrTiO3 substrate [101, 102]. They can be thought of as a stack ofweakly-interacting CuO layers, whose structure is depicted in Fig. 5.1(a)and consists of two intercalated CuO2 lattices (sharing the same O). A CuO2layer is sketched in Fig. 5.1(b). Because Cu 3dx2−y2 orbitals only hybridizewith their ligand O 2p orbitals, shown in the same color in Fig. 5.1, the twoCuO2 sublattices would be effectively decoupled if pp hopping between thetwo sets of O 2p orbitals was absent.In this case, a hole doped into one sublattice would evolve just like ina regular CuO2 layer, and the same (now doubly-degenerate) qp dispersionwould be predicted by both one- and three-band models, as already dis-965.2. Modelcussed.However, the CuO2 sublattices are coupled by pp hopping between thetwo sets of O 2p orbitals, which lifts this degeneracy. The resulting qp disper-sion was measured by angle resolved photoemission spectroscopy (ARPES)[4]. It was found to be overall quite similar to that of CuO2 and seemed tobe well described by a small cluster study of a t-t′-t′′-J model. As we shownext, this conclusion is opposite to the one we find for the Udd → ∞ limitof the three-band model. We predict a qualitatively different dispersion forT-CuO and CuO2, but these differences are masked in magnetically twinnedsamples.We give a detailed presentation of the model in section 5.2. In section5.3 we describe the method we use to solve it. Results and discussions aregiven in section 5.4.5.2 ModelWe study the Udd →∞ limit of the Emery model, with spins at the Cu sitesand a single doped hole on the O sublattice. This limit is justified becauseUdd is much larger than the other energy scales [25]. For a CuO2 layer, thecorresponding Hamiltonian, see Fig. 5.1(b), is [100]:Hˆ = Tˆpp + Tˆswap + HˆJpd + HˆJdd . (5.1)For simplicity of notation, we provide here the explicit expressions for theterms in the CuO2 Hamiltonian, assuming that only the ligand O 2p orbitalsare included. Generalization to including both sets of O 2p orbitals, and alsoto T-CuO, is straightforward. With the sign of positive/negative lobes aspictured in Fig. 5.1 and using p†i,σ as the creation operator for a hole in theligand O 2p orbital located at i, we have:Tpp = tpp∑i∈O,δ,σrδp†i,σpi+δ,σ − t′pp∑i∈O,σp†i,σ(pi−,σ + pi+,σ). (5.2)975.2. Model(a) (b)tt’ppJdd pptsw Jpdtpp~(c) (d)Figure 5.1: Structure of a layer of (a) T-CuO, and (b) CuO2. Full circles areCu, empty squares are O. The Cu 3dx2−y2 orbitals are drawn at a few sites,with white/dark lobes showing our choice for positive/negative signs. Thecorresponding ligand O 2p orbitals are also indicated on neighboring O sites.The T-CuO layer can be thought of as two intercalated CuO2 layers sharingcommon O. The coppers of the two sublattices hybridize with different O 2porbitals. Panels (c) and (d) show the two degenerate ground-states of theundoped T-CuO layer. Different colors are used for the Cu spins on the twosublattices for better visibility.985.2. ModelThe lattice constant is set to a = 1. The vectors δ = ±(0.5, 0.5),±(0.5,−0.5)are the distances between any O and its four nn O sites, and rδ = ±1 setsthe sign of each nn pp hopping integral in accordance with the overlap ofthe 2p orbitals involved. Next nn hopping is included only between O 2porbitals pointing toward a common bridging Cu, separated by ε = (1, 0) or(0, 1); hybridization with the 4s orbital of the bridging Cu further booststhe value of this hopping integral.Tswap describes Cu-mediated effective hopping accompanied by a spin-swap. Specifically, the hole at a Cu site adjacent to the doped hole hops toone of its other neighbor O sites, followed by the doped hole falling into thevacated Cu orbital. Because the original doped hole replaces the Cu hole,their spins are swapped. ThusTswap = −tsw∑i∈Cu,u6=u′,σ,σ′su−u′p†i+u,σpi+u′,σ′ |iσ′〉〈iσ|, (5.3)where u,u′ = (±0.5, 0), (0,±0.5) are the distances between a Cu and itsfour nn O sites. It shows the change of the Cu spin located at Ri from σ toσ′ as the doped hole changes its spin from σ′ to σ while moving to anotherO. The sign sη = ±1 is due to the overlaps of the orbitals involved in theprocess, and the overall minus is because of the interchange in the order ofholes6.The origin ofHˆJpd = Jpd∑i,uSi · si+u (5.4)is similar, except that now the Cu hole hops onto the O that is hosting thedoped hole, followed by one of the two holes returning to the Cu. Unlikefor Tˆswap, charge is not moved in this process; instead, it gives rise to AFMexchange between the spin si+u of the doped hole and that of its neighborCu, Si.6We have not specified a canonical order of fermions in the many-body states becausewe only deal with a single explicit charge carrier. However, whatever order for the basisstates we choose, the swap process will exchange an electron on a copper site and anelectron on an oxygen sign, leading to the minus sign.995.2. ModelFinally,HˆJdd = Jdd∑〈i,j〉′Si · Sj (5.5)is the usual AFM coupling between neighbor Cu spins. The sum runs over allnearest-neighbor copper pairs except the one that has the bridging oxygenoccupied by the doped hole. The energy scale Jdd ∼ 150 meV is taken asthe unit of energy, in terms of which tpp = 4.1, t′pp = 0.6tpp, tsw = 3.0 andJpd = 2.8 [25]. Note that the oxygen Hubbard repulsion Upp is not includedin Eq. (5.1) because we consider only the case of a single doped hole.In CuO2 the important O 2p orbitals are the ligand orbitals, but it isstraightforward to generalize the model to also include the in-plane non-ligand orbitals [3]. Because they do not hybridize with the Cu 3dx2−y2orbitals, their addition does not change Tˆswap, HˆJpd or HˆJdd , all of which arisefrom such hybridization. Only Tˆpp must be supplemented accordingly. Bysymmetry, nn hopping between two non-ligand orbitals is the same tpp as forligand orbitals, with signs dictated by the lobes’ overlap. Hopping betweennn ligand and non-ligand orbitals, which we call Tˆmix and is shown by theblack arrow in Fig. 5.1(a), has magnitude t˜pp/tpp = (tpp,σ − tpp,pi)/(tpp,σ +tpp,pi) = 0.6 because tpp,σ = 4tpp,pi [103]. Thus, adding the non-ligand orbitalsdoes not introduce new energy scales. For CuO2, their inclusion has a minoreffect on the qp dispersion [3].The Hamiltonian for T-CuO is a straightforward generalization of Eq.(5.1). The pp hopping is unchanged and described by the same Tˆpp + Tˆmixdiscussed above. Because of the two intercalated Cu sublattices, there aretwo sets of terms Tˆswap, HˆJpd and HˆJdd which couple Cu spins on eachsublattice to each other and to the doped hole – if the latter occupies a2p orbital that has ligand character for that Cu sublattice. We use thesame parameters for T-CuO like for CuO2 (the results remain qualitativelysimilar if the parameters are varied within reasonable ranges) and focus onthe effect of Tˆmix, which moves the hole between the two sets of 2p orbitalsand changes to which Cu sublattice it is coupled.For completeness, we note that in Ref. [4] a small nn exchange J˜dd ∼ 0.04was also included between Cu spins on the different sublattices. We have1005.2. Modelconsidered this term, as well as a weak J˜pd coupling between the doped holeoccupying a non-ligand orbital and its neighbor Cu spins (FM exchange isfavored by Hund’s coupling when the Cu hole hops into the O orbital or-thogonal to that hosting the doped hole). None of these terms were found tolead to qualitative changes, and because of their relatively small magnitude,their quantitative effects are minor. As a result, we ignore them from nowon.5.2.1 Summary of model assumptionsThere are quite a few assumptions going into the derivation of this model.For ease of reference, we summarize them here.The first assumption concern the nature of the hopping terms includedin the Hamiltonian. Longer-range hoppings could be considered, but areneglected because in the case of CuO2 they were shown to have negligibleeffect [17].The next assumption concerns the degree to which the Udd → ∞ limitis treated perturbatively. The resulting terms give rise to magnetic interac-tions and are carried out to leading order. Higher-order terms only serve toslightly renormalize already present terms [6].The biggest assumption is in treating the magnetic interaction betweenholes on the copper ions with an Ising term instead of a Heisenberg term.In effect, this means neglecting the copper’s spin fluctuations; the only waycopper spins can get flipped is through their interaction with the oxygen hole.In essence, this assumption is similar to the assumption made for the double-well electron-phonon coupling model where quartic free-phonon terms wereabsorbed into the quartic el-ph interaction term. The justification, however,is a different one. The timescale of spin fluctuations would be 1/Jdd, which islower than the other timescales in the system. As such, it is not immediatelyclear that neglecting these fluctuations is valid. However, the results in [3]clearly justify this approximation a posteriori.Additional assumptions concern terms arising from the interaction be-tween the two intercalated layers. Here, we only explicitly include nearest-1015.3. Methodneighbor hopping between them. Certain higher-order magnetic terms werefound to have little effect on the ground state properties and are thereforeneglected.5.3 MethodAs discussed in Ref. [3], we extract the qp dispersion Eqp(k) from the one-hole propagator computed variationally in a restricted Hilbert space thatallows up to nm magnons to be created by the doped hole through Tˆswapand HˆJpd processes, assuming that it was injected in a Ne´el-like background.Of course, in reality there are spin fluctuations in the AFM background, butbecause their energy scale Jdd is small, they occur so slowly as to have littleeffect on the qp dispersion: the hole creates and moves its magnon cloud ona timescale faster than that controlling the spin fluctuations, so the lattercan be ignored [3, 5]. If the T-CuO energy scales are similar, and given theweak coupling between the two Cu sublattices, this approximation shouldremain valid.In undoped T-CuO, each Cu sublattice has AFM order due to its HˆJddterm. As a result, any weak coupling J˜dd between the two Cu sublatticesis fully frustrated: a spin on one sublattice interacts with equal numbersof up and down spins from the other sublattice. Nevertheless, order bydisorder selects one of the two degenerate states depicted in Fig. 5.1(c), (d)as the ground-state of the undoped system [104–106]. Because they haveFM chains running along either the x = y or x = −y diagonals, they arerelated by a C4 rotation so it suffices to study one case. Unlike in CuO2, forT-CuO in either of these states we expect that the quasiparticle dispersionEqp(k) is not invariant to C4 rotations, only to C2 ones.Let us now discuss our variational method in a bit more detail. For anygiven momentum k, our corresponding variational Hilbert space contains allstates with at most nm magnons and with a restriction on the maximumallowed magnon-hole distances, as discussed next.Zero-magnon states are the eight Bloch states of momentum k, one foreach of the eight oxygen orbitals in the magnetic unit cell. States with1025.3. Methodmagnons are Bloch states of total momentum k for given hole and magnonsconfigurations with fixed relative distances between holes and magnons. Forexample, the one-magnon states have the form|k, α, δ〉 = 1√N∑ieik·rip†iα,↓S+i+δ |Ne´el〉where i runs over the unit cells, iα denotes oxygen orbital α = 1, . . . , 8within unit cell i, and S+i+δ is the magnon creation operator for a down-spin copper of the magnetic unit cell located a distance δ away from unitcell i. The two-magnon and three-magnon states are constructed similarlyby adding magnons in the units cells located at distances δ2 and δ3 apart,respectively.Our restriction is that any |δ| ≤ Nmax. For the 1- and 2-magnon states,the variational approach together with the Lanczos method [107] can eas-ily handle distances of up to Nmax = 10, although the low-energy (quasi-particle) results are converged even for Nmax = 2. This is not surprisingbecause we are concerned with the physics of the polaronic bound state,where the magnons are bound close to the hole. For the 3-magnon states,we set Nmax = 1, i.e. we only include 3-magnon configurations where the3 magnons are located either in the same unit cell as the hole, or in a unitcell directly adjacent to it. These are the configurations with the highestweight in the quasiparticle cloud but, as the results show, they do not lead tosignificant changes. Less likely 3-magnon configurations, with the magnonsspread further apart, can therefore be safely ignored at these energies.We note that this approximation is valid for calculating low-energy prop-erties. It would fail at describing higher energy features such as the correctlocations for the polaron+one-magnon continuum, since its states necessarilyhave a magnon far away from the polaron, so the corresponding configura-tions need to be included to capture such higher energy features.This approximation is conceptually similar to the momentum averageapproximation: The creation of bosons (phonons or magnons) lowers theenergy available to the charge carrier such that it becomes a good approxi-mation to restrict the presence of bosons to a final (small) number of lattice1035.4. Results and discussionsites. A difference is that while a lattice site in the Holstein model can har-bour an arbitrary number of phonons, the number of magnons is limited toone per lattice site. The result is, depending on the formulation, a finite setof coupled equations of motion for the propagators or a sparse Hamiltonianof manageable size.The resulting Hamiltonian matrix in our basis is sparse, because for eachconfiguration there are only a few allowed hoppings and spin flips. Thus,it readily lends itself to the Lanczos method, which provides us with thespectral function A(k, ω) relying only on computationally cheap matrix-vector products. The lowest-lying peak in the spectral function then givesus the qp dispersion.5.4 Results and discussionFigures 5.2(a)-(c) show Eqp(k) for the magnetic order of Fig. 5.1(c) ob-tained with the variational method for nm = 1, 2, 3, respectively, inside theBrillouin zone (BZ) displayed in Fig. 5.2(d). Full/dashed lines are for T-CuO/CuO2.In CuO2, at the points marked by circles and squares there are equiva-lent nearly isotropic minima [97, 16]. With increasing nm, the bandwidthnarrows and the dispersion flattens below the polaron+one magnon contin-uum (both these effects are due to standard polaronic physics discussed inRef. [3]) but the shape is unchanged. The results are nearly converged atnm = 3 for CuO2, with a bandwidth of ∼ 2Jdd in agreement with availableexact diagonalization results and with experimental data [3]).In T-CuO, we verified that for Tˆmix = 0 the same (but now doubly-degenerate) dispersion is obtained. When Tˆmix is turned on, this degeneracyis lifted. Only the low-energy eigenstate is shown in Fig. 5.2. Again,results are essentially converged for nm = 3. As expected, the dispersionloses its invariance to C4 rotations because the qp is now evolving in amagnetic background that lacks this symmetry. The dispersion in the kx =−ky quadrants again displays deep, isotropic minima around ±(pi2 ,−pi2 ) (fullsquares) and is thus rather similar to that in CuO2. The difference, however,1045.4. Results and discussion-22-20-18-16-14E qp(k)-24-22-20E qp(k)(0,0) (pi,pi) (pi,0) (0,0) (0,pi) (pi,0) k-24-22-20E qp(k)(a) nm = 1(b) nm = 2(c) nm = 3k yk x(d)(f)ππ12345(e)12341243Figure 5.2: Qp dispersion in units of Jdd for (a) nm = 1, (b) nm = 2, and(c) nm = 3 with full/dashed lines for T-CuO/CuO2. The Brillouin zone forthe magnetic order of Fig. 5.1(c) is shown in red in (b). The shaded areais the smaller BZ for CuO2. The points marked by circles and empty/fullsquares are equivalent in CuO2 but not in T-CuO. (d) Hopping between twoadjacent ZRSs, and (e) between a ZRS (red) and one with x− y symmetry(blue). See text for more details.1055.4. Results and discussionkxkypipi(c)i i−(b)(a)iFigure 5.3: (a) Unit cell for CuO2, with two Cu spins and four ligand Oorbitals. We use the location i of the down-spin Cu as the reference point.The white/shaded areas indicate our choice for positive/negative lobs. (b)Magnetic Brillouin zone (shaded region) vs. full Brillouin zone (unshaded).(c) ZRS between a hole occupying the linear combination of ligand orbitalswith x2 − y2 symmetry, and the spin of the central Cu. The Bloch state isobtained from its translations on the corresponding magnetic sublattice.is significant in the kx = ky quadrants near the ±(pi2 , pi2 ) points (circles).Not only are energies here higher than at the ±(pi2 ,−pi2 ) points, but theseminima are shifted toward the Γ point. Note also that the BZ corners (emptysquares) continue to mark local minima, but now lying at high energies justbelow the polaron+one magnon continuum.5.4.1 The ZRS Bloch stateBefore we can continue our discussion of the results, we have to describe inmore detail the nature of the Zhang-Rice singlet state. For the CuO2 lattice,taking into consideration the AFM order of the Cu spins, we can choose theunit cell as shown in Fig. 5.3(a), with the corresponding magnetic Brillouinzone shown in Fig. 5.3(b).To define a ZRS Bloch state, we first introduce:p†x2−y2,i,σ =12[p†i+x2,σ + p†i+y2,σ− p†i−x2,σ − p†i−y2,σ](5.6)1065.4. Results and discussionwhich describes the doped hole occupying a linear combination of ligandorbitals with x2 − y2 symmetry, centered on the Cu located at i. The ZRSis obtained when a hole occupying such a state is locked in a singlet withthe Cu spin, therefore a ZRS Bloch state can be defined as7:|d,k〉 =∑i∈Cu↓eik·Ri√Np†x2−y2,i,↑ − p†x2−y2,i,↓S+i√2|Ne´el〉. (5.7)Here, k is any momentum in the magnetic Brillouin zone and the sum isonly over sites in the spin-down Cu sublattice (since a spin-up doped holecan form a ZRS only with these spins).Of course, one can also define Bloch states associated with singlets thathave other symmetries for the linear combination of O orbitals. Of all thesestates, in CuO2 the ZRS Bloch state is found to have the largest overlapwith the low-energy quasiparticle wavefunction.Its first excited state, on the other hand, is found to have the largestoverlap with Bloch states based on the singlet with x − y symmetry, i.e.the singlet obtained using p†x−y,i,σ =12[p†i+x2,σ + p†i+y2,σ+ p†i−x2,σ + p†i−y2,σ]instead of p†x2−y2,i,σ in Eq. (5.7). We call this state |p,k〉.5.4.2 Low-energy effective modelWe now prove that the unusual dispersion for T-CuO involves physics be-yond the Zhang-Rice singlet. As such, it cannot be described by one-bandmodels obtained through a projection onto these states.We start by estimating the effect of Tˆmix on the two CuO2-like degenerateeigenstates that appear in its absence, and whose energy E0(k) is shown bythe dashed lines in Fig. 5.2. Especially near the (±pi2 ,±pi2 ) points, the CuO2quasiparticle indeed has a large overlap with a ZRS Bloch state [5], andthe hole occupies the x2 − y2 linear combination of O 2p ligand orbitals7On occasion, reviewers have asked whether this expression misses an S−i operator,corresponding to the S+i in the second term within the sum. We emphasize that theexpression is correct as it stands and is not lacking a S−i term. It correctly describes astate that is a superposition of intact AFM order and one-magnon AFM order.1075.4. Results and discussionsketched for two nn sites in Fig. 5.2(e). For T-CuO, these two Bloch statescombine into one Bloch state |d,k〉 with a momentum k in the bigger, T-CuO Brillouin zone. If we use it as an approximation for the low-energyeigenstate, then the T-CuO dispersion becomes Eqp(k) ≈ E0(k) + δE(k),where δE(k) = 〈d,k|Tˆmix|d,k〉 is calculated in the following.Previously, we have viewed T-CuO as two copies of a CuO2 lattice, byadding two more Cu ions and four more oxygen orbitals to the unit cellwhile keeping the same lattice vectors. An alternative view is to use thesame unit cell as for CuO2, i.e., two copper ions and four oxygen orbitals,but have one of the lattice vectors shortened by a factor 2, depending on therelative magnetic ordering of the sublattices. If the copper spins are alignedferromagnetically along the main diagonal, the proper lattice vectors area1 = a(1/2, 1/2)T and a2 = a(−1, 1)T .While both views are valid, the latter view has the smallest possible unitcell and thus does not exhibit folding of the true Brillouin zone.For a very simple effective theory, we use perturbation theory to in-vestigate the influence of the inter-sublattice hopping Tˆmix. Because thisoperator moves the hole between the two sets of O orbitals but cannot movemagnons between the different Cu sublattices, only the magnon-free partof the quasi-particle wavefunction will contribute to matrix elements of thisoperator. Consider a two-dimensional Hilbert space containing the d- andp-wave ZRS states of momentum k. Without Tˆmix, we haveH =(E0 00 E1)where E0(k) and E1(k) are the CuO2 dispersions of the groundstate andthe first excited state.We now apply Tˆmix to each of the four oxygen basis states from whichthe ZRS-type states are built. Let|1,k〉 :=∑i∈Cu↓eik·Ri√Np†i+x2,↑ |Ne´el〉 .1085.4. Results and discussionand |2,k〉, |3,k〉 and |4,k〉 defined analogously for the other oxygen orbitalsenumerated counter-clockwise around the same copper ion. Note that herethe sum is over all down-spins in the T-CuO lattice. Taking care of theproper phases for the hoppings, we haveTˆmix |1,k〉 = −2t˜pp cos(k · a1) |1,k〉 − t˜pp[e−ik·a1 + e−ik·(a1−a2)]|3,k〉Tˆmix |2,k〉 = −2t˜pp cos(k · a1) |2,k〉 − t˜pp[e−ik·a1 + e−ik·(a1+a2)]|4,k〉Tˆmix |3,k〉 = −2t˜pp cos(k · a1) |3,k〉 − t˜pp[eik·a1 + eik·(a1−a2)]|1,k〉Tˆmix |4,k〉 = −2t˜pp cos(k · a1) |4,k〉 − t˜pp[eik·a1 + eik·(a1+a2)]|2,k〉From this, we can now compute 〈d/p,k|Tˆmix|d/p,k〉.〈d,k|Tˆmix|d,k〉 = t˜pp cos(k · a1) [cos(k · a2)− 1]〈p,k|Tˆmix|p,k〉 = t˜pp cos(k · a1) [cos(k · a2)− 3]〈p,k|Tˆmix|d,k〉 = it˜pp sin(k · a1) [cos(k · a2) + 1]Strictly speaking, these matrix elements should be weighted by the appro-priate quasiparticle weights Zk corresponding to projection of the true qpeigenstate onto these non-interacting Bloch states. However, these weightsare known to be rather featureless near the (pi/2, pi/2) points that are mostrelevant in this discussion, so we ignore them in the following. This willaffect results quantitatively, but not qualitatively.If we ignore the p-wave states and only take into consideration the d-waveZRS state, then the change in energy due to Tˆmix becomesδE(k) = −t˜pp cos kx + ky2[1− cos(kx − ky)] . (5.8)Because δE(kx = −ky) = −2t˜pp sin2 kx and δE(ky = kx∓pi) = −2t˜pp sin |kx|,the minima at ±(pi2 ,−pi2 ) (full squares) are pushed to lower energies. Sim-ilarly, the minima at the corners of the BZ (empty squares) are pushed tohigh energies, which agrees with the results of Fig. 5.2.However, things are different near the ±(pi2 , pi2 ) points (circles). Because1095.4. Results and discussionδE(kx = ky) = 0, here the dispersion should remain unchanged instead ofthe minima moving toward the Γ point. Moreover, we find that the overlapof the T-CuO quasiparticle wavefunction with the ZRS Bloch state |d,k〉vanishes at k = ±(pi2 , pi2 ). These facts clearly prove that the changes nearthe ±(pi2 , pi2 ) points cannot be due to Zhang-Rice singlet physics.Indeed, Tˆmix hopping between x2−y2 linear combinations centred at nnCu sites is suppressed, see Fig. 5.2(e): eg., a hole at site 1 of the lower Cu(red) hops into p†1 + p†3 of the upper Cu (blue), which is orthogonal to itsx2− y2 linear combination. Instead, here hopping between adjacent x2− y2and x−y combinations is enhanced, see Fig. 5.2(f). The shift of the ±(pi2 , pi2 )minima toward Γ is due to a large mixing of the singlet with x−y symmetryinto the quasiparticle eigenstate, which thus loses its ZRS nature. Note thatexperiments like Refs. [108–110], which are sensitive only to the local singletcharacter, cannot distinguish between a ZRS singlet and one of such mixedsymmetry.Since the d-wave ZRS state alone cannot explain the dispersion of T-CuO, we now turn our focus to the Hilbert space spanned by the d-waveand p-wave ZRS states. It turns out that considering Tˆmix in this two-dimensional Hilbert space is already sufficient to explain qualitatively thequasi-particle physics of T-CuO.Let us consider the special case of the kx = ky = k line, where k · a2 = 0and k · a1 = ka. We then haveTmix(k) =(0 −2it˜pp sin(ka)2it˜pp sin(ka) −2t˜pp cos(ka).)The full numerical results showed that the minimum along this line getsshifted from k = pi/2 closer to the Γ-point, and this readily follows fromthe simple form here: The off-diagonal elements provide mixing of the d-and p-states, and the energy −2t˜pp cos(ka) of the p-state then moves theminimum to a lower k.On the other hand, for kx = pi − ky = k we find that Tdd = Tpp = 0and Tpd = it˜pp [1 + cos(pi − 2k)]. In this case, the Hamiltonian is symmetric1105.4. Results and discussion(0,0) (1,1) (1.5, 0.5) (0.5, -0.5) (0,0) (0,1) (1,0)(kx, ky) / pi/a-25-20-15E qp / JddfulleffectiveFigure 5.4: Comparison between the one-magnon numerical dispersion andthe effective calculations involving both the low-energy and first excitedstate quasiparticles. The results are in qualitative agreement, in particularthe shift of the minimum along kx = ky. The dashed lines are a guide tothe eye to demonstrate that the peak along Γ −M gets shifted away fromk = (pi/2, pi/2) yet remains there along X −X ′.around k = pi/2, and thus the minimum gets only shifted down in energywhile remaining at k = pi/2. Other directions can be considered similarly.In Fig. 5.4 we compare the full numerical results to the effective low-energy results from the 2 × 2 Hamiltonian. The qualitative agreement isstriking, in particular the shift of the minimum from (pi2 ,pi2 ) towards the Γpoint emerges as predicted. Of course, quantitative disagreements resultfrom the crudeness of our approximation: we neglected the quasiparticleweights which would narrow the bandwidth for this effective low-energyresult, we ignored other higher-energy eigenstates as well as the fact that theCuO2 low-energy and first excited quasiparticles are not pure x2−y2- and x−y ZRS-like Bloch states, respectively. Nonetheless, the results satisfactorilyshow that the changes near (pi2 ,pi2 ) arise as a result of mixing between stateswith d- and p-symmetry. This mixing, in turn, results from a desire to gainkinetic energy from hopping between the sublattices, which is not possiblein the pure ZRS subspace: for this relative arrangement of the two Ne´el1115.4. Results and discussionkxkypikykxkxkyFigure 5.5: The first two panels show the Brillouin zones corresponding tothe two possible magnetic orders of Fig. 5.1(c) and (d). The symbols markthe deep minima (full squares), displaced shallower minima (circles) andvery shallow minima (empty squares). The two figures are related by a C4rotation. The third panel shows the average of these two patterns, whereat each special point the deepest local minimum was selected. This patternhas a restored C4 symmetry and a Brillouin zone (blue line) correspondingto one Cu site/unit cell. The pattern of deep/shallower minima is like thatfound experimentally.1125.4. Results and discussionCu sublattices, Tˆmix cannot hop a ZRS with momentum kx = ky betweenneighbor Cu sites, see Fig. 2(d) of the main text. Higher-energy physics, ofnon-ZRS origin, then becomes relevant.We checked that adding terms like J˜dd and J˜pd has no qualitative effect:the dispersion remains like in Fig. 5.2. This is expected because their ma-trix elements are small and/or featureless near (±pi2 ,±pi2 ). We are thereforeconfident that our prediction is robust.5.4.3 Comparison to ARPESAngular resolved photoelectron spectroscopy (ARPES) finds the T-CuO qpdispersion to obey C4 symmetry and to have a large Brillouin zone, corre-sponding to a unit cell containing one Cu and one O atom [4]. Both featuresare very surprising for the long-range magnetic orders of Figs. 5.1(c), (d),each of which break the C4 symmetry. Moreover, any AFM-type order hasat least two magnetically non-equivalent Cu atoms and thus can have a BZlike in Fig. 5.2(d) or smaller, not larger. Our results become consistentwith the ARPES data if we assume the presence of magnetic domains inboth ground-states, so that their average is measured experimentally. In-deed, as shown below, averaging the band-structure of Fig. 5.2(d) with itscounterpart rotated by 90o leads to an apparent doubling of the Brillouinzone and a new pattern of minima with two different energies, in agreementwith those found experimentally.The result of averaging over domains with both possible orientationsis demonstrated in Fig. 5.5. Panel (a) shows the T-CuO Brillouin zones(red rectangles) from Fig. (d) of the main text, which corresponds to themagnetic order of Fig. 5.1(c). The symbols mark the deep minima (fullsquares), displaced shallower minima (circles) and very shallow minima(empty squares). Panel (b) is obtained by a C4 rotation and correspondsto the magnetic order of Fig. 5.1(d). Their average is shown in panel (c),where the symbols now mark the lowest-energy local minimum. The result-ing pattern agrees with that measured experimentally for T-CuO. The largeBrillouin zones (blue squares) corresponding to a unit cell with one Cu per1135.4. Results and discussionbasis emerges naturally, as do the patterns of deep and shallower minima.Fig. 5.6 shows the quasiparticle dispersion along the same contour dis-cussed in Fig. 1 of Ref. [4], for both magnetic orientations. Note that ourresults are in hole language, so to compare with their ARPES data the en-ergies should be reversed, Eqp → −Eqp. For ease of comparison, the lowerpanel shows the same results in this electron picture, with the symmetrypoints also labelled like in Ref. [4].We predict that a dispersion like in Fig. 5.2 should appear in the ARPESof “magnetically untwinned” T-CuO films in the insulating limit. This isvery different and therefore easily distinguishable from the one-band modelprediction [4]. The observation of this pattern, with shallower displacedminima in two of the quadrants, will provide a clear proof of low-energyphysics beyond the ZRS, and of the superiority of three-band models tomodel such materials. If T-CuO films can be doped, for example by gating,this new pattern of minima will open extraordinary opportunities to testmany ideas relating the shape of the Fermi surface, location of “hot spots”and possibility of nesting, to much of the cuprate phenomenology, includingthe symmetry of the pairing and of the superconducting gap, formation ofstripes, appearance and relevance of various other ordered phases, etc.While the feasibility of such experiments remains to be determined, animportant lesson from this study is that low-energy physics of non-ZRSnature can arise in such materials in suitable circumstances/symmetries.The presence of disorder, of other nearby quasiparticles, of stripes, charge-density wave or other ordered phases may have a similar effect in CuO2layers.1145.4. Results and discussion(0,0) (-pi,0) (-2pi,0) (-pi,−pi) (0,0) k-24-23-22-21-20-19E qp / JddOrientation 1c)Orientation 1d)Γ B X’ B’ M A’ X/M’ A Γ k-21-20-19-18-17-16E qp / JddOrientation 1c)Orientation 1d)Figure 5.6: Quasiparticle dispersion for the two orientations of the magneticbackground, along the contour considered in Ref. [20]. Top panel shows theresults in the hole picture used throughout this work, whereas the low panelshows the same results in the electron picture relevant for experiments, withfurther tuning of the parameters. Also, we note that the slight displacementsof the shallower minima have not been observed experimentally. This maybe because of the very broad widths of the quasiparticle peaks and the lossof spectral weight on one side of these points, which may mask them. In“untwinned” samples the scattering rate should be lower, which may makethe observation of these displacements towards Γ more easily visible. Ofcourse, in that case the lack of C4 symmetry and existence of really shallowminima should also become visible.115Chapter 6Conclusions6.1 Summary of this workIn this thesis, we have studied a variety of simple extensions to commonlyused model Hamiltonians. Throughout, we have seen that even seeminglyminor alterations can lead to drastic quantitative and, more importantly,qualitative changes in the system’s behaviour. Our results highlight thatlongstanding and commonly held assumptions about certain systems andmodels can be wrong when they are based on extrapolations of the simplestandard Hamiltonian.In Chapter 2, we first introduced higher-order non-linear terms into thestandard Holstein model’s electron-phonon coupling. We demonstrate thatthe most dramatic qualitative changes are due to the quadratic term andthus a-posterio justify that we focus mainly on this case. By accessingthe polaron’s ground-state properties via an extension to the momentum-average approximation, it is seen that even a small quadratic term leads toa complete change in the polaron’s characteristics, making it much lighter.We argue that these results invalidate a large portion of what is commonlyassumed regarding polaronic behaviour: The linear model, by necessity, as-sumes that lattice deformations are small. In the strong coupling regime,it then predicts that lattice deformations are large. The quadratic model,on the other hand, suffers from no such contradiction: By assuming mod-erate lattice deformations and thus keeping terms up to second order, thequadratic model predicts only small to moderate lattice deformations andthus remains internally consistent.In Chapter 3, we consider a different type of non-linear Holstein model,one where the linear term is non-existent due to symmetry. We call this1166.2. Further developmentsthe double-well electron-phonon coupling model. For this model, we showthat its results cannot be reproduced by an effective linear model. We thenhighlighted similarities and differences between this model’s polaron andthat of the linear model.The double-well model’s bipolaron is studied in Chapter 4, where weshow that the particular nature of the double-well electron-phonon couplingcan lead to bipolarons that are strongly bound yet lightweight. Bound quasi-particles with this property were previously believed to be restricted to muchmore complicated models, either to those with peculiar lattice geometriesor much more complicated electron-phonon interactions. In studying thisproblem, we have also introduced a first step to extending the momentum-average approximation to systems with more than one particle.Finally, in Chapter 5 we study a novel type of copper-oxide layer, the re-cently grown tetragonal copper oxide T-CuO. This model can be thought ofas an extension of the thoroughly studied CuO2 layer, extended with a sec-ond intercalated copy weakly coupled via inter-sublattice hopping. We useda variational method similar to MA and solved it via exact diagonalizationto study the spectral function and dispersion of the resulting quasi-particle.We discover that the dispersion is quite significantly changed from that ofCuO2 in a way that cannot be explained by a one-band Zhang-Rice singlettype of state. Instead, the inclusion of the small inter-band coupling requiresus to use a multi-band description. The importance of this result comes fromthe fact that a long-standing theoretical issue regarding the cuprates is pre-cisely whether the correct model for the cuprates is a single-band model(where spin-fluctuations play an important part) or a multi-band model(where spin-fluctuations are relatively unimportant). We demonstrate howour results establish that T-CuO can be an experimental test case for thisimportant modelling question.6.2 Further developmentsOur calculations for the non-linear Holstein model in Chapter 2 were per-formed for the single-polaron case. A follow-up question to our results is1176.2. Further developmentswhether quadratic terms also have a strong impact at finite density. Thisrequires different techniques. A former member of our group has pursuedthis avenue for the two-dimensional Hubbard-Holstein model [70, 71] usingdeterminant Quantum Monte Carlo methods. Their results confirm thatsmall non-linear terms remain important, in that they strongly influencethe effective electron-lattice coupling and charge-density-wave correlations.Similarly, the double-well electron-phonon coupling model should bestudied at finite temperature and density. An open question would be if thethermoynamic limit shows spontaneous symmetry breaking, which wouldmanifest itself as (anti-)ferroelectricity. Furthermore, with the results ofChapter 4 in mind, a study of the superconducting properties would revealif the strongly-bound yet lightweight bipolarons of the double-well modelcan condensate and give rise to polaronic superconductivity.The immediate follow-up to our study of T-CuO is clear: Experimental-ists should aim to grow untwinned samples of T-CuO and carefully studyits dispersion, to answer whether the results of our three-band model or ofa one-band tJ-model are those realized in the material. Another interest-ing question is whether there are other effects present in CuO2 that cangive rise to similar non-ZRS physics. Such effects would likely be basedon symmetries that are incompatible with a d-wave object, similar to howthe inter-sublattice hopping leads to destructive interference for the d-waveZRS.Another avenue of exploration is to follow up on the premise that it isthree-band models, not one-band models, that correctly capture the inter-action between quasi-particles. Currently ongoing work in our group stud-ies the effective interaction between two holes in a CuO2 layer in the sameUdd →∞ limit of the three-band Emery model as we used in our study of T-CuO. The question is whether this model, in contrast to the tJ-model, showsan attractive effective interaction between two holes mediated by magnons.Apart from the particular model Hamiltonians we have studied in thisthesis, our general approach is flexible and can be applied to other models. Itproves particularly useful when studying a small number of carriers couplingto gapped bosonic modes. The underlying idea of both MA and the various1186.2. Further developmentsvariational approaches is that creation of these bosonic modes costs energy.In the MA view, this means that the carrier loses energy when it createsa boson. At energies below the free carrier’s ground state, its real-spacepropagator decays exponentially, justifying its replacement by a simplifiedversion. In the variational view, subspaces with different boson numbersare well separated, justifying a cut-off at relatively small boson number.In addition, a combination of the MA and variational arguments justifiesvariational approximations where the carrier is forced to stay close to thebosonic modes it excites.119Bibliography[1] Mona Berciu. Green’s function of a dressed particle. Phys. Rev.Lett., 97:036402, Jul 2006. doi: 10.1103/PhysRevLett.97.036402. URLhttp://link.aps.org/doi/10.1103/PhysRevLett.97.036402.[2] M Mo¨ller, A Mukherjee, CPJ Adolphs, DJJ Marchand, and M Berciu.Efficient computation of lattice green functions for models with longerrange hopping. Journal of Physics A: Mathematical and Theoretical,45(11):115206–115214, 2012.[3] Hadi Ebrahimnejad, George A Sawatzky, and Mona Berciu. The dy-namics of a doped hole in a cuprate is not controlled by spin fluctua-tions. Nature Physics, 2014.[4] S. Moser, L. Moreschini, H.-Y. Yang, D. Innocenti, F. Fuchs, N. H.Hansen, Y. J. Chang, K. S. Kim, A. L. Walter, A. Bostwick, E. Roten-berg, F. Mila, and M. Grioni. Angle-resolved photoemission spec-troscopy of tetragonal cuo: Evidence for intralayer coupling betweencupratelike sublattices. Phys. Rev. Lett., 113:187001, Oct 2014. doi:10.1103/PhysRevLett.113.187001. URL http://link.aps.org/doi/10.1103/PhysRevLett.113.187001.[5] H Ebrahimnejad, G A Sawatzky, and M Berciu. Differences betweenthe insulating limit quasiparticles of one-band and three-band cupratemodels. Journal of Physics: Condensed Matter, 28(10):105603, 2016.[6] Bayo Lau. Modeling Polarons in Transition-Metal Oxides. PhD thesis,The University Of British Columbia (Vancouver), 2011.120Bibliography[7] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev.,136:B864–B871, Nov 1964. doi: 10.1103/PhysRev.136.B864. URLhttp://link.aps.org/doi/10.1103/PhysRev.136.B864.[8] W. Kohn and L. J. Sham. Self-consistent equations including exchangeand correlation effects. Phys. Rev., 140:A1133–A1138, Nov 1965. doi:10.1103/PhysRev.140.A1133. URL http://link.aps.org/doi/10.1103/PhysRev.140.A1133.[9] M. Born and R. Oppenheimer. Zur quantentheorie der molekeln. An-nalen der Physik, 389(20):457–484, 1927. ISSN 1521-3889. doi: 10.1002/andp.19273892002. URL http://dx.doi.org/10.1002/andp.19273892002.[10] Neil W Ashcroft, N David Mermin, and Sergio Rodriguez. Solid-statephysics. American Journal of Physics, 46(1):116–117, 1978.[11] T Holstein. Studies of polaron motion: Part i. the molecular-crystalmodel. Annals of Physics, 8(3):325 – 342, 1959. ISSN 0003-4916. doi:http://dx.doi.org/10.1016/0003-4916(59)90002-8. URL http://www.sciencedirect.com/science/article/pii/0003491659900028.[12] Glen L Goodvin, Mona Berciu, and George A Sawatzky. Greens func-tion of the holstein polaron. Physical Review B, 74(24):245104, 2006.[13] Mona Berciu and Glen L Goodvin. Systematic improvement of themomentum average approximation for the greens function of a holsteinpolaron. Physical Review B, 76(16):165109, 2007.[14] Glen L. Goodvin and Mona Berciu. Momentum average approxi-mation for models with electron-phonon coupling dependent on thephonon momentum. Phys. Rev. B, 78:235120, Dec 2008. doi:10.1103/PhysRevB.78.235120. URL http://link.aps.org/doi/10.1103/PhysRevB.78.235120.[15] Gerald D Mahan. Many-particle physics. Springer Science & BusinessMedia, 2013.121Bibliography[16] Andrea Damascelli, Zahid Hussain, and Zhi-Xun Shen. Angle-resolvedphotoemission studies of the cuprate superconductors. Rev. Mod.Phys., 75:473–541, Apr 2003. doi: 10.1103/RevModPhys.75.473. URLhttp://link.aps.org/doi/10.1103/RevModPhys.75.473.[17] Hadi Ebrahimnejad. Variational studies of dressed quasiparticles prop-erties and their interactions with external potentials. PhD thesis, Uni-versity of British Columbia, Dec 2014. URL https://open.library.ubc.ca/cIRcle/collections/24/items/1.0166096.[18] Richard D Mattuck. A guide to Feynman diagrams in the many-bodyproblem. Courier Corporation, 2012.[19] C. Slezak, A. Macridin, G. A. Sawatzky, M. Jarrell, and T. A. Maier.Spectral properties of holstein and breathing polarons. Phys. Rev.B, 73:205122, May 2006. doi: 10.1103/PhysRevB.73.205122. URLhttp://link.aps.org/doi/10.1103/PhysRevB.73.205122.[20] William H Press. Numerical recipes 3rd edition: The art of scientificcomputing. Cambridge university press, 2007.[21] Mona Berciu. Berciu replies:. Phys. Rev. Lett., 98:209702, May 2007.doi: 10.1103/PhysRevLett.98.209702. URL http://link.aps.org/doi/10.1103/PhysRevLett.98.209702.[22] J.G. Bednorz and K.A. Mller. Possible hight c superconductivity inthe balacuo system. Zeitschrift fr Physik B Condensed Matter, 64(2):189–193, 1986. ISSN 0722-3277. doi: 10.1007/BF01303701. URLhttp://dx.doi.org/10.1007/BF01303701.[23] J. Zaanen, G. A. Sawatzky, and J. W. Allen. Band gaps and electronicstructure of transition-metal compounds. Phys. Rev. Lett., 55:418–421,Jul 1985. doi: 10.1103/PhysRevLett.55.418. URL http://link.aps.org/doi/10.1103/PhysRevLett.55.418.[24] V. J. Emery. Theory of high-tc superconductivity in oxides. Phys. Rev.122BibliographyLett., 58:2794–2797, Jun 1987. doi: 10.1103/PhysRevLett.58.2794.URL http://link.aps.org/doi/10.1103/PhysRevLett.58.2794.[25] Masao Ogata and Hidetoshi Fukuyama. The t j model for the oxidehigh- t c superconductors. Reports on Progress in Physics, 71(3):036501, 2008. URL http://stacks.iop.org/0034-4885/71/i=3/a=036501.[26] F. C. Zhang and T. M. Rice. Effective hamiltonian for the su-perconducting cu oxides. Phys. Rev. B, 37:3759–3761, Mar 1988.doi: 10.1103/PhysRevB.37.3759. URL http://link.aps.org/doi/10.1103/PhysRevB.37.3759.[27] H. Eskes and G. A. Sawatzky. Tendency towards local spin compen-sation of holes in the high-Tc copper compounds. Phys. Rev. Lett.,61:1415–1418, Sep 1988. doi: 10.1103/PhysRevLett.61.1415. URLhttp://link.aps.org/doi/10.1103/PhysRevLett.61.1415.[28] K A Chao, J Spalek, and A M Oles. Kinetic exchange interaction in anarrow s-band. Journal of Physics C: Solid State Physics, 10(10):L271,1977. URL http://stacks.iop.org/0022-3719/10/i=10/a=002.[29] K. A. Chao, J. Spa lek, and A. M. Oles´. Canonical perturbation ex-pansion of the hubbard model. Phys. Rev. B, 18:3453–3464, Oct 1978.doi: 10.1103/PhysRevB.18.3453. URL http://link.aps.org/doi/10.1103/PhysRevB.18.3453.[30] Masao Ogata and Hiroyuki Shiba. Holes in two-dimensional cuo2 sys-tem pairing mechanism emerging from small-cluster studies. Journalof the Physical Society of Japan, 57(9):3074–3088, 1988. doi: 10.1143/JPSJ.57.3074. URL http://dx.doi.org/10.1143/JPSJ.57.3074.[31] A. I. Lichtenstein and M. I. Katsnelson. Antiferromagnetism andd -wave superconductivity in cuprates: A cluster dynamical mean-field theory. Phys. Rev. B, 62:R9283–R9286, Oct 2000. doi:10.1103/PhysRevB.62.R9283. URL http://link.aps.org/doi/10.1103/PhysRevB.62.R9283.123Bibliography[32] E. Pavarini, I. Dasgupta, T. Saha-Dasgupta, O. Jepsen, and O. K.Andersen. Band-structure trend in hole-doped cuprates and cor-relation with tcmax. Phys. Rev. Lett., 87:047003, Jul 2001. doi:10.1103/PhysRevLett.87.047003. URL http://link.aps.org/doi/10.1103/PhysRevLett.87.047003.[33] O.K. Andersen, A.I. Liechtenstein, O. Jepsen, and F. Paulsen.{LDA} energy bands, low-energy hamiltonians, t, t, t (k), and j.Journal of Physics and Chemistry of Solids, 56(12):1573 – 1591,1995. ISSN 0022-3697. doi: http://dx.doi.org/10.1016/0022-3697(95)00269-3. URL http://www.sciencedirect.com/science/article/pii/0022369795002693. Proceedings of the Conference on Spectro-scopies in Novel Superconductors.[34] H. Eskes and G. A. Sawatzky. Single-, triple-, or multiple-bandhubbard models. Phys. Rev. B, 44:9656–9666, Nov 1991. doi:10.1103/PhysRevB.44.9656. URL http://link.aps.org/doi/10.1103/PhysRevB.44.9656.[35] Hiroyuki Matsui, Andrei S. Mishchenko, and Tatsuo Hasegawa. Distri-bution of localized states from fine analysis of electron spin resonancespectra in organic transistors. Phys. Rev. Lett., 104:056602, Feb 2010.doi: 10.1103/PhysRevLett.104.056602. URL http://link.aps.org/doi/10.1103/PhysRevLett.104.056602.[36] S. Ciuchi and S. Fratini. Band dispersion and electronic lifetimesin crystalline organic semiconductors. Phys. Rev. Lett., 106:166403,Apr 2011. doi: 10.1103/PhysRevLett.106.166403. URL http://link.aps.org/doi/10.1103/PhysRevLett.106.166403.[37] A Lanzara, PV Bogdanov, XJ Zhou, SA Kellar, DL Feng, ED Lu,T Yoshida, H Eisaki, Atsushi Fujimori, K Kishio, et al. Evidencefor ubiquitous strong electron–phonon coupling in high-temperaturesuperconductors. Nature, 412(6846):510–514, 2001.124Bibliography[38] KM Shen, F Ronning, DH Lu, WS Lee, NJC Ingle, W Meevasana,F Baumberger, A Damascelli, NP Armitage, LL Miller, et al. Missingquasiparticles and the chemical potential puzzle in the doping evolu-tion of the cuprate superconductors. Physical review letters, 93(26):267002, 2004.[39] Dmitri Reznik, L Pintschovius, M Ito, S Iikubo, M Sato, H Goka,M Fujita, K Yamada, GD Gu, and JM Tranquada. Electron–phononcoupling reflecting dynamic charge inhomogeneity in copper oxide su-perconductors. Nature, 440(7088):1170–1173, 2006.[40] Jinho Lee, K Fujita, K McElroy, JA Slezak, M Wang, Y Aiura,H Bando, M Ishikado, T Masui, J-X Zhu, et al. Interplay of electron–lattice interactions and superconductivity in bi2sr2cacu2o8+ δ.Nature, 442(7102):546–550, 2006.[41] C Gadermaier, AS Alexandrov, VV Kabanov, P Kusar, T Mertelj,X Yao, C Manzoni, D Brida, G Cerullo, and D Mihailovic. Electron-phonon coupling in high-temperature cuprate superconductors deter-mined from electron relaxation rates. Physical review letters, 105(25):257001, 2010.[42] O Gunnarsson and O Ro¨sch. Interplay between electron–phonon andcoulomb interactions in cuprates. Journal of Physics: Condensed Mat-ter, 20(4):043201, 2008.[43] Norman Mannella, Wanli L Yang, Xing Jiang Zhou, Hong Zheng,John F Mitchell, Jan Zaanen, Thomas P Devereaux, Naoto Na-gaosa, Zahid Hussain, and Z-X Shen. Nodal quasiparticle in pseu-dogapped colossal magnetoresistive manganites. Nature, 438(7067):474–478, 2005.[44] T Yildirim, O Gu¨lseren, JW Lynn, CM Brown, TJ Udovic, Q Huang,N Rogado, KA Regan, MA Hayward, JS Slusky, et al. Giant anhar-monicity and nonlinear electron-phonon coupling in mgb 2: a com-125Bibliographybined first-principles calculation and neutron scattering study. Physi-cal review letters, 87(3):037001, 2001.[45] Amy Y Liu, II Mazin, and Jens Kortus. Beyond eliashberg super-conductivity in mgb 2: anharmonicity, two-phonon scattering, andmultiple gaps. Physical Review Letters, 87(8):087005, 2001.[46] II Mazin and VP Antropov. Electronic structure, electron–phononcoupling, and multiband effects in mgb 2. Physica C: Superconductiv-ity, 385(1):49–65, 2003.[47] Hyoung Joon Choi, David Roundy, Hong Sun, Marvin L Cohen, andSteven G Louie. The origin of the anomalous superconducting prop-erties of mgb2. Nature, 418(6899):758–760, 2002.[48] A Kaminski, M Randeria, JC Campuzano, MR Norman, H Fretwell,J Mesot, T Sato, T Takahashi, and K Kadowaki. Renormalization ofspectral line shape and dispersion below t c in bi 2 sr 2 cacu 2 o 8+δ. Physical Review Letters, 86(6):1070, 2001.[49] TK Kim, AA Kordyuk, SV Borisenko, A Koitzsch, M Knupfer,H Berger, and J Fink. Doping dependence of the mass enhancementin (p b, b i) 2 s r 2 c a c u 2 o 8 at the antinodal point in the super-conducting and normal states. Physical review letters, 91(16):167002,2003.[50] C. N. Veenstra, G. L. Goodvin, M. Berciu, and A. Damascelli. Spec-tral function tour of electron-phonon coupling outside the migdallimit. Phys. Rev. B, 84:085126, Aug 2011. doi: 10.1103/PhysRevB.84.085126. URL http://link.aps.org/doi/10.1103/PhysRevB.84.085126.[51] H. Frhlich. Electrons in lattice fields. Advances in Physics, 3(11):325–361, 1954. doi: 10.1080/00018735400101213. URL http://dx.doi.org/10.1080/00018735400101213.126Bibliography[52] AS Mishchenko and Naoto Nagaosa. Electron-phonon coupling anda polaron in the t- j model: From the weak to the strong couplingregime. Physical review letters, 93(3):036402, 2004.[53] AS Mishchenko, N Nagaosa, Z-X Shen, G De Filippis, V Cataudella,TP Devereaux, Christian Bernhard, Kyung Wan Kim, and J Zaanen.Charge dynamics of doped holes in high t c cuprate superconductors: aclue from optical conductivity. Physical review letters, 100(16):166401,2008.[54] Felipe Herrera and Roman V Krems. Tunable holstein model withcold polar molecules. Physical Review A, 84(5):051401, 2011.[55] JP Hague and C MacCormick. Quantum simulation of electron–phonon interactions in strongly deformable materials. New Journalof Physics, 14(3):033019, 2012.[56] Vladimir M Stojanovic´, Tao Shi, C Bruder, and J Ignacio Cirac. Quan-tum simulation of small-polaron formation with trapped ions. Physicalreview letters, 109(25):250501, 2012.[57] Felipe Herrera, Kirk W. Madison, Roman V. Krems, and MonaBerciu. Investigating polaron transitions with polar molecules. Phys.Rev. Lett., 110:223002, May 2013. doi: 10.1103/PhysRevLett.110.223002. URL http://link.aps.org/doi/10.1103/PhysRevLett.110.223002.[58] O Entin-Wohlman, H Gutfreund, and M Weger. Mass and frequencyrenormalization for quadratic electron-phonon coupling. Solid StateCommunications, 46(1):1–5, 1983.[59] O Entin-Wohlman, H Gutfreund, and M Weger. Band narrowing andsusceptibility enhancement by a quadratic electron-phonon interac-tion. Journal of Physics C: Solid State Physics, 18(3):L61, 1985.[60] Vincent H Crespi and Marvin L Cohen. Anharmonic phonons and127Bibliographyhigh-temperature superconductivity. Physical Review B, 48(1):398,1993.[61] VM Kenkre. What do polarons owe to their harmonic origins? PhysicaD: Nonlinear Phenomena, 113(2):233–241, 1998.[62] C. P. J. Adolphs and M. Berciu. Going beyond the linear approxi-mation in describing electron-phonon coupling: Relevance for the hol-stein model. EPL (Europhysics Letters), 102(4):47003, 2013. URLhttp://stacks.iop.org/0295-5075/102/i=4/a=47003.[63] H Fehske and SA Trugman. Numerical solution of the holstein polaronproblem. In Polarons in Advanced Materials, pages 393–461. Springer,2007.[64] G Herzberg. Dissociation energy and ionization potential of molecularhydrogen. Physical review letters, 23(19):1081, 1969.[65] Christopher Gerry and Peter Knight. Introductory quantum optics.Cambridge university press, 2005.[66] David Stoler. Equivalence classes of minimum uncertainty packets.Physical Review D, 1(12):3217, 1970.[67] J Bonca, T Katrasnik, and SA Trugman. Mobile bipolaron. Physicalreview letters, 84(14):3153, 2000.[68] Alexandru Macridin, GA Sawatzky, and Mark Jarrell. Two-dimensional hubbard-holstein bipolaron. Physical Review B, 69(24):245111, 2004.[69] AR Davenport, JP Hague, and PE Kornilovitch. Mobile small bipo-larons on a three-dimensional cubic lattice. Physical Review B, 86(3):035106, 2012.[70] Shaozhi Li and S. Johnston. The effects of non-linear electron-phononinteractions on superconductivity and charge-density-wave correla-tions. Europhysics Letters, 109:27007, 2015.128Bibliography[71] Shaozhi Li, EA Nowadnick, and S Johnston. Quasiparticle proper-ties of the nonlinear holstein model at finite doping and temperature.Physical Review B, 92(6):064301, 2015.[72] Y. Zolotaryuk, P. L. Christiansen, and J. Juul Rasmussen. Polarondynamics in a two-dimensional anharmonic holstein model. Phys. Rev.B, 58:14305–14319, Dec 1998. doi: 10.1103/PhysRevB.58.14305. URLhttp://link.aps.org/doi/10.1103/PhysRevB.58.14305.[73] S Raghavan and V M Kenkre. Quantum mechanical bound rotator asa generalized harmonic oscillator. Journal of Physics: Condensed Mat-ter, 6(47):10297, 1994. URL http://stacks.iop.org/0953-8984/6/i=47/a=013.[74] N. K. Voulgarakis and G. P. Tsironis. Stationary and dynamicalproperties of polarons in the anharmonic holstein model. Phys. Rev.B, 63:014302, Dec 2000. doi: 10.1103/PhysRevB.63.014302. URLhttp://link.aps.org/doi/10.1103/PhysRevB.63.014302.[75] P. Maniadis, G. Kalosakas, K. Ø. Rasmussen, and A. R. Bishop. Po-laron normal modes in the peyrard-bishop-holstein model. Phys. Rev.B, 68:174304, Nov 2003. doi: 10.1103/PhysRevB.68.174304. URLhttp://link.aps.org/doi/10.1103/PhysRevB.68.174304.[76] J. K. Freericks, M. Jarrell, and D. J. Scalapino. Holstein model ininfinite dimensions. Phys. Rev. B, 48:6302–6314, Sep 1993. doi: 10.1103/PhysRevB.48.6302. URL http://link.aps.org/doi/10.1103/PhysRevB.48.6302.[77] D. Meyer, A. C. Hewson, and R. Bulla. Gap formation and soft phononmode in the holstein model. Phys. Rev. Lett., 89:196401, Oct 2002.doi: 10.1103/PhysRevLett.89.196401. URL http://link.aps.org/doi/10.1103/PhysRevLett.89.196401.[78] Keisuke Mitsumoto and Yoshiaki no. Phonon softening and double-well potential formation due to electronphonon interaction in heavy-fermion systems. Physica C: Superconductivity, 426431, Part 1:330129Bibliography– 334, 2005. ISSN 0921-4534. doi: http://dx.doi.org/10.1016/j.physc.2004.12.021. URL http://www.sciencedirect.com/science/article/pii/S092145340500273X. Proceedings of the 17th Interna-tional Symposium on Superconductivity (ISS 2004)Advances in Super-conductivity {XVIIProceedings} of the 17th International Symposiumon Superconductivity (ISS 2004).[79] A. Bussmann-Holder, A. Simon, and H. Bu¨ttner. Possibility of a com-mon origin to ferroelectricity and superconductivity in oxides. Phys.Rev. B, 39:207–214, Jan 1989. doi: 10.1103/PhysRevB.39.207. URLhttp://link.aps.org/doi/10.1103/PhysRevB.39.207.[80] O. Ro¨sch and O. Gunnarsson. Electron-phonon interaction in thet-j model. Phys. Rev. Lett., 92:146403, Apr 2004. doi: 10.1103/PhysRevLett.92.146403. URL http://link.aps.org/doi/10.1103/PhysRevLett.92.146403.[81] Bayo Lau, Mona Berciu, and George A. Sawatzky. Single-polaronproperties of the one-dimensional breathing-mode hamiltonian. Phys.Rev. B, 76:174305, Nov 2007. doi: 10.1103/PhysRevB.76.174305. URLhttp://link.aps.org/doi/10.1103/PhysRevB.76.174305.[82] S. Johnston, F. Vernay, B. Moritz, Z.-X. Shen, N. Nagaosa, J. Zaanen,and T. P. Devereaux. Systematic study of electron-phonon couplingto oxygen modes across the cuprates. Phys. Rev. B, 82:064513, Aug2010. doi: 10.1103/PhysRevB.82.064513. URL http://link.aps.org/doi/10.1103/PhysRevB.82.064513.[83] O. S. Bariˇsic´. Calculation of excited polaron states in the holsteinmodel. Phys. Rev. B, 69:064302, Feb 2004. doi: 10.1103/PhysRevB.69.064302. URL http://link.aps.org/doi/10.1103/PhysRevB.69.064302.[84] Mona Berciu. Momentum average approximations for the holsteinpolaron. Canadian Journal of Physics, 86(4):523–527, 2008.130Bibliography[85] A.S. Alexandrov. Superconducting polarons and bipolarons. InA.S. Alexandrov, editor, Polarons in Advanced Materials, vol-ume 103 of Springer Series in Materials Science, pages 257–310. Springer Netherlands, 2007. ISBN 978-1-4020-6347-3. doi:10.1007/978-1-4020-6348-0 7. URL http://dx.doi.org/10.1007/978-1-4020-6348-0_7.[86] J. Boncˇa and S. A. Trugman. Bipolarons in the extended hol-stein hubbard model. Phys. Rev. B, 64:094507, Aug 2001. doi:10.1103/PhysRevB.64.094507. URL http://link.aps.org/doi/10.1103/PhysRevB.64.094507.[87] J. P. Hague, P. E. Kornilovitch, J. H. Samson, and A. S. Alexan-drov. Superlight small bipolarons in the presence of a strong coulombrepulsion. Phys. Rev. Lett., 98:037002, Jan 2007. doi: 10.1103/PhysRevLett.98.037002. URL http://link.aps.org/doi/10.1103/PhysRevLett.98.037002.[88] A. S. Alexandrov, J. H. Samson, and G. Sica. Superlight smallbipolarons from realistic long-range coulomb and fro¨hlich interac-tions. Phys. Rev. B, 85:104520, Mar 2012. doi: 10.1103/PhysRevB.85.104520. URL http://link.aps.org/doi/10.1103/PhysRevB.85.104520.[89] J. P. Hague and P. E. Kornilovitch. Light and stable triplet bipo-larons on square and triangular lattices. Phys. Rev. B, 82:094301, Sep2010. doi: 10.1103/PhysRevB.82.094301. URL http://link.aps.org/doi/10.1103/PhysRevB.82.094301.[90] Clemens P. J. Adolphs and Mona Berciu. Single-polaron properties fordouble-well electron-phonon coupling. Phys. Rev. B, 89:035122, Jan2014. doi: 10.1103/PhysRevB.89.035122. URL http://link.aps.org/doi/10.1103/PhysRevB.89.035122.[91] Eleftherios N Economou. Green’s functions in quantum physics.Green’s Functions in Quantum Physics:, Springer Series in Solid-State131BibliographySciences, Volume 7. ISBN 978-3-540-28838-1. Springer-Verlag BerlinHeidelberg, 2006, 1, 2006.[92] Mona Berciu. Few-particle green’s functions for strongly correlatedsystems on infinite lattices. Phys. Rev. Lett., 107:246403, Dec 2011.doi: 10.1103/PhysRevLett.107.246403. URL http://link.aps.org/doi/10.1103/PhysRevLett.107.246403.[93] G. A. Sawatzky. Quasiatomic auger spectra in narrow-band metals.Phys. Rev. Lett., 39:504–507, Aug 1977. doi: 10.1103/PhysRevLett.39.504. URL http://link.aps.org/doi/10.1103/PhysRevLett.39.504.[94] S. A. Trugman. Interaction of holes in a hubbard antiferromagnet andhigh-temperature superconductivity. Phys. Rev. B, 37:1597–1603, Feb1988. doi: 10.1103/PhysRevB.37.1597. URL http://link.aps.org/doi/10.1103/PhysRevB.37.1597.[95] P. W. Leung and R. J. Gooding. Dynamical properties of the single-hole t - J model on a 32-site square lattice. Phys. Rev. B, 52:R15711–R15714, Dec 1995. doi: 10.1103/PhysRevB.52.R15711. URL http://link.aps.org/doi/10.1103/PhysRevB.52.R15711.[96] P. W. Leung, B. O. Wells, and R. J. Gooding. Comparison of 32-site exact-diagonalization results and arpes spectral functions for theantiferromagnetic insulator sr2cuo2cl2. Phys. Rev. B, 56:6320–6326,Sep 1997. doi: 10.1103/PhysRevB.56.6320. URL http://link.aps.org/doi/10.1103/PhysRevB.56.6320.[97] B. O. Wells, Z. X. Shen, A. Matsuura, D. M. King, M. A. Kastner,M. Greven, and R. J. Birgeneau. e versus k relations and many bodyeffects in the model insulating copper oxide sr2cuo2cl2. Phys. Rev.Lett., 74:964–967, Feb 1995. doi: 10.1103/PhysRevLett.74.964. URLhttp://link.aps.org/doi/10.1103/PhysRevLett.74.964.[98] Mirko Mo¨ller, George A. Sawatzky, and Mona Berciu. Magnon-mediated interactions between fermions depend strongly on the lat-132Bibliographytice structure. Phys. Rev. Lett., 108:216403, May 2012. doi: 10.1103/PhysRevLett.108.216403. URL http://link.aps.org/doi/10.1103/PhysRevLett.108.216403.[99] Mirko Mo¨ller, George A. Sawatzky, and Mona Berciu. Role of thelattice structure in determining the magnon-mediated interactionsbetween charge carriers doped into a magnetically ordered back-ground. Phys. Rev. B, 86:075128, Aug 2012. doi: 10.1103/PhysRevB.86.075128. URL http://link.aps.org/doi/10.1103/PhysRevB.86.075128.[100] Bayo Lau, Mona Berciu, and George A. Sawatzky. High-spin polaronin lightly doped cuo2 planes. Phys. Rev. Lett., 106:036401, Jan 2011.doi: 10.1103/PhysRevLett.106.036401. URL http://link.aps.org/doi/10.1103/PhysRevLett.106.036401.[101] Wolter Siemons, Gertjan Koster, Dave H. A. Blank, Robert H. Ham-mond, Theodore H. Geballe, and Malcolm R. Beasley. Tetragonalcuo: End member of the 3d transition metal monoxides. Phys. Rev.B, 79:195122, May 2009. doi: 10.1103/PhysRevB.79.195122. URLhttp://link.aps.org/doi/10.1103/PhysRevB.79.195122.[102] D. Samal, Haiyan Tan, Y. Takamura, W. Siemons, Jo Verbeeck,G. Van Tendeloo, E. Arenholz, C. A. Jenkins, G. Rijnders, and GertjanKoster. Direct structural and spectroscopic investigation of ultrathinfilms of tetragonal cuo: Six-fold coordinated copper. EPL (Euro-physics Letters), 105(1):17003, 2014. URL http://stacks.iop.org/0295-5075/105/i=1/a=17003.[103] W.A. Harrison. Elementary Electronic Structure. Elementary Elec-tronic Structure. World Scientific, 1999. ISBN 9789810238957. URLhttps://books.google.ca/books?id=eAfi3ajSTcUC.[104] Christopher L. Henley. Ordering due to disorder in a frustrated vec-tor antiferromagnet. Phys. Rev. Lett., 62:2056–2059, Apr 1989. doi:133Bibliography10.1103/PhysRevLett.62.2056. URL http://link.aps.org/doi/10.1103/PhysRevLett.62.2056.[105] P. Chandra, P. Coleman, and A. I. Larkin. Ising transition in frus-trated heisenberg models. Phys. Rev. Lett., 64:88–91, Jan 1990.doi: 10.1103/PhysRevLett.64.88. URL http://link.aps.org/doi/10.1103/PhysRevLett.64.88.[106] Ce´dric Weber, Luca Capriotti, Gre´goire Misguich, Federico Becca,Maged Elhajal, and Fre´de´ric Mila. Ising transition driven by frus-tration in a 2d classical model with continuous symmetry. Phys.Rev. Lett., 91:177202, Oct 2003. doi: 10.1103/PhysRevLett.91.177202. URL http://link.aps.org/doi/10.1103/PhysRevLett.91.177202.[107] Cornelius Lanczos. An iteration method for the solution of the eigen-value problem of linear differential and integral operators. UnitedStates Governm. Press Office, 1950.[108] L. H. Tjeng, B. Sinkovic, N. B. Brookes, J. B. Goedkoop, R. Hesper,E. Pellegrin, F. M. F. de Groot, S. Altieri, S. L. Hulbert, E. Shekel, andG. A. Sawatzky. Spin -resolved photoemission on Anti -ferromagnets:Direct observation of zhang-rice singlets in cuo. Phys. Rev. Lett.,78:1126–1129, Feb 1997. doi: 10.1103/PhysRevLett.78.1126. URLhttp://link.aps.org/doi/10.1103/PhysRevLett.78.1126.[109] N. B. Brookes, G. Ghiringhelli, O. Tjernberg, L. H. Tjeng, T. Mi-zokawa, T. W. Li, and A. A. Menovsky. Detection of zhang-ricesinglets using spin-polarized photoemission. Phys. Rev. Lett., 87:237003, Nov 2001. doi: 10.1103/PhysRevLett.87.237003. URL http://link.aps.org/doi/10.1103/PhysRevLett.87.237003.[110] N. B. Brookes, G. Ghiringhelli, A.-M. Charvet, A. Fujimori,T. Kakeshita, H. Eisaki, S. Uchida, and T. Mizokawa. Stability ofthe zhang-rice singlet with doping in lanthanum strontium copper ox-ide across the superconducting dome and above. Phys. Rev. Lett.,134115:027002, Jul 2015. doi: 10.1103/PhysRevLett.115.027002. URLhttp://link.aps.org/doi/10.1103/PhysRevLett.115.027002.[111] Erik Koch. The lanczos method. In Eva Pavarini, Erik Koch, DieterVollhardt, and Alexander Lichtenstein, editors, The LDA+DMFT ap-proach to strongly correlated materials, volume 1, 2011. URL http://www.cond-mat.de/events/correl11/manuscript/Koch.pdf.135Appendix AVarious formal methods136Appendix BAppendix for NonlinearHolstein ModelB.1 Free propagator for the quadratic modelThe Hamiltonian whose diagonal (in real space) Green’s function we needto compute isH0 = Hel +Hph + g2∑ic†ici(2b†ibi + 1).We look at the Green’s functionG0(j − i, ω − nΩ) = 〈0|cjbni Gˆ0(ω)(b†i )nc†i |0〉 ,which describes free propagation of an electron from the side of the phononcloud to another site j. We split the Hamiltonian H0 asH00 = Hel +Hph + g2, H01 = 2g2∑ic†icib†ibi137B.2. Equations of motion for quartic modeland use Dyson’s identity to obtainG0(j − i, ω − nΩ) = G00(j − i, ω − nΩ− g2)+∑l〈0|bni cjGˆ0(ω)H1c†l (b†i )n|0〉 ·G00(l − i, ω − nΩ− g2)= G00(j − i, ω − nΩ− g2)+ 2g2G0(j − i, ω;n)G00(i− i, ω − nΩ− g2)=G00(j − i, ω − nΩ− g2)1− 2g2g˜0(ω − nΩ− g2) .G00(j − i, ω) = 1N∑keik·(Rj−Ri)ω−k+iη is the free propagator for the electron, andg˜0(ω) = G00(i − i, ω) is a short-hand notation. The only dependence ofG0 on the distance of sites j − i is in the free electron propagator G00.For energies below the free tight-binding model’s ground state energy, thisquantity decays exponentially with distance j − i, which justifies the MAapproximation [1].G0(j − i, ω) ≈ δijG0(j − i, ω).With g¯0(ω − nΩ) = G0(i− i, ω − nΩ), we obtaing¯0(ω, n) =1[g˜0(ω − nΩ− g2)]−1 − 2g2n. (B.1)We can obtain g¯0(ω;n) for any lattice model whose single-electron propaga-tor we can compute. This includes the tight-binding model in any dimensionwith next-neighbor and finite-range hopping [2].B.2 Equations of motion for quartic modelIn the quartic model, the el-ph coupling term is given by∑4n=1 gn(b+ b†)n.The equation of motion for the generalized Green’s functions then dependson the expansion of these terms.The first two terms have already been treated in the main text; they138B.2. Equations of motion for quartic modelreproduce the quadratic model. For the third and fourth term, we have toexpand (b + b†)3 and (b + b†)4, respectively. This is a tedious and boringexercise in correctly applying the bosonic commutation relations to groupall terms together. Once this is done, we see that the third-order term willlink a propagator Fn to propagators Fn±3 as well as Fn±1. The latter comesfrom terms such as b†b2 and bb†b from the expansion, which changes theoverall phonon number by −1. Expanding the fourth-order bracket, then,shows that this term couples Fn to Fn±4, Fn±2 and Fn itself.Hence, instead of having vectors Wn of size 2 and matrices αn, βn, γnof size 2 × 2, we now have vectors of size 4 and matrices of size 4 × 4,with Wn = (F4n−3, F4n−2, F4n−1, F4n). Inserting the expansions for thevarious coupling terms into Dyson’s identity then allows us again to collectthe generalized propagators of the appropriate order to obtain an EOM inmatrix form, where we give the matrix elements below, using the followingnotational conventions:• xn¯ = x(x− 1)(x− 2) . . . (x− n+ 1).• gn0 = g0(ω, n), withg0(ω, n) =[1g¯0(ω − nΩ− g2 − 3g4) − n× (2g2 + 6g4 + 6ng4)]−1and g¯0 being the bare free electron propagator.139B.2. Equations of motion for quartic modelWith this, we can then write the non-zero matrix elements as follows:α11 = g4n−30 · (4n− 3)4¯ · g4α12 = g4n−30 · (4n− 3)3¯ · g3α13 = g4n−30 · (4n− 3)2¯ · (g2 + (16n− 14)g4)α14 = g4n−30 · (4n− 3) · (g1 + (12n− 9)g3)α22 = g4n−20 · (4n− 2)4¯ · g4α23 = g4n−20 · (4n− 2)3¯ · g3α24 = g4n−20 · (4n− 2)2¯ · (g2 + (16n− 10)g4)α33 = g4n−10 · (4n− 1)4¯ · g4α34 = g4n−10 · (4n− 1)3¯ · g3α44 = g4n0 · (4n)4¯ · g4β11 = g4n−30 · g4β21 = g4n−20 · g3β22 = g4n−20 · g4β31 = g4n−10 · (g2 + (16n− 4)g4)β32 = g4n−10 · g3β33 = g4n−10 · g4β41 = g4n0 · (g1 + (12n+ 3)g3)β41 = g4n0 · (g2 + (16n+ 6)g4)β41 = g4n0 · g3β41 = g4n0 · g4140B.2. Equations of motion for quartic modelγ11 = 1γ12 = −g4n−30 · (g1 + (12n− 6)g3)γ13 = −g4n−30 · (g2 + (16n− 6)g4)γ14 = −g4n−30 · g3γ21 = −g4n−20 · (4n− 2) · (g1 + (12n− 6)g3)γ22 = 1γ23 = −g4n−20 · (g1 + (12n− 3)g3)γ24 = −g4n−20 · (g2 + (16n− 2)g4)γ31 = −g4n−10 · (4n− 1)2¯ · (g2 + (16n− 4)g4)γ32 = −g4n−10 · (4n− 1) · (g1 + (12n− 3)g3)γ33 = 1γ34 = −g4n−10 · (g1 + 12ng3)γ41 = −g4n0 · (4n)3¯ · g3γ42 = −g4n0 · (4n)2¯ · (g2 + (16n− 2)g4)γ43 = −g4n0 · (4n) · (g1 + 12ng3)γ44 = 1Then, the recusion defined by γnWn = αnWn−1 + βnWn+1 is solved inthe same manner as outlined in the main text.141Appendix CAppendix for QuadraticHolstein ModelC.1 Details for the even-sectorC.1.1 Coupling matricesThe matrices appearing in Eq. (3.11) are:γ=n |11 = 1− g¯0(ω − 4nΩ)(8ng2 + 24ng4 + 96n2g4)γ=n |12 = −g¯0(ω − 4nΩ)(g2 + 6g4 + 16ng4)γ=n |21 = −g¯0(ω − (4n+ 2)Ω)((g2 + 6g4)(4n+ 2)2¯+4g4(4n+ 2)3¯)γ=n |22 = 1− g¯0(ω − (4n+ 2)Ω)((8n+ 4)g2+(24n+ 12)g4 + 24(2n+ 1)2g4)α=n |11 = g¯0(ω − 4nΩ)(g4(4n)4¯)α=n |12 = g¯0(ω − 4nΩ)((g2 + 6g4)(4n)2¯ + 4g4(4n)3¯)α=n |21 = 0α=n |22 = g¯0(ω − (4n+ 2)Ω(g4(4n+ 2)4¯)142C.1. Details for the even-sectorβ=n |11 = g¯0(ω − 4nΩ)g4β=n |12 = 0β=n |21 = g¯0(ω − (4n+ 2)Ω)(g2 + 6g4 + (16n+ 8)g4)β=n |22 = g¯0(ω − (4n+ 2)Ω)g4The matrices for 6= sector are the same if we substitute n→ n− 1/2 every-where except in the argument of g¯0(ω).C.1.2 Manipulation of the EOMsWe can rewrite the EOM of F1 by inserting the matrices A=1 and A6=1 andcollecting terms. This results inF1(ij) = G0(j − i, ω − 2Ω) [a=0 G(j) + a=1 F=1 (j)]+∑l 6=jG0(l − i, ω − 2Ω)a 6=F 6=1 (lj). (C.1)where we omit the arguments k and ω for shorter notation. We give ex-pressions for the various coefficients below. For now, we rewrite the EOMasF1(ij) = G0(j − i, ω − 2Ω)×[a=0 G(j) + (a=1 − a 6=1 )F=1 (j)]+∑lG0(l − i, ω − 2Ω)a6=F1(lj). (C.2)Defining G0(ω)ij := G0(j − i, ω), we can write this as a matrix product:∑l[δil − a 6=G0(ω − 2Ω)il]F1(lj) =G0(ω − 2Ω)ij[a=0 G(j) + (a=1 − a6=1 )F1(jj)].143C.2. Details for the odd-sectorWe multiply this from the left with G−10 (ω − 2Ω) and obtain∑l[G−10 (ω − 2Ω)rl − a6=δrl]F1(lj) = δrj[a 6=0 G(j) + (a=1 − a 6=)F1(jj)].Next, we use the fact that G−10 (ω−2Ω)rl = δrl(ω−2Ω)−Hˆrl, so subtractinga 6=δrl from this just shifts its frequency to obtain G−10 (ω− 2Ω− a 6=)rl. As aresult:F1(ij) = G0(ω − 2Ω − a 6=)ij ×[a=0 G(j) + (a=1 − a6=1 )F1(jj)].Since in the EOM for G we only require F1(jj), we solve for that diagonalelement and obtainF1(jj) =g¯0(ω − 2Ω− a 6=)a=0 G(j)1− g¯0(ω − 2Ω− a 6=)(a=1 − a6=1 ).The coefficients are obtained by just inserting the appropriate matrices Aninto the EOM and collecting terms:a=0 = 2g2 + 12g4 + (g2 + 14g4)A=1 |11 + g4A=1 |21a=1 = 4g2 + 36g4 + (g2 + 14g4)A=2 |12 + g4A=2 |22a6= = (g2 + 6g4)A6=1 |12 + g4A 6=1 |22Finally, F1(jj) are used in Eq. (3.7) to obtain G(k, ω).C.2 Details for the odd-sectorC.2.1 Equations of MotionStarting from the EOM for Gijl(ω), we let H1 act on the states in thosesums, to find for the diagonal state:H1c†l b†l |0〉 = (2g2 + 12g4)c†l b†l |0〉 + (g2 + 10g4)c†l b†,3l |0〉 + g4c†l b†,5l |0〉144C.2. Details for the odd-sectorwhile for the off-diagonal ones:H1c†i′b†l |0〉 = (2g2 + 6g4)c†i′b†l |0〉+ (g2 + 6g4)c†i′b†,2l b†l |0〉+ g4c†i′b†,4l b†l |0〉 .We now define the generalized Green functions as:Fn(k, i, j, ω) = 〈k|Gˆ(ω)cib†,2ni bj |0〉so we always have the extra phonon at site j. The equation of motion forG then becomes: Gijl(ω) = G0(j − i, ω − Ω) +[(2g2 + 12g4)F=0 (l) + (g2 +10g4)F=1 (l)+g4F=2 (l)]Gill+∑i′ 6=l[(2g2 +6g4)F6=0 (i′, l)+(g2 +6g4)F6=1 (i′, l)+g4F6=2 (i′, l)]G0(i′− i, ω−Ω). Again, we start by separating the cases F=n andF 6=n . The resulting equations of motion for F=n are like those of the even-sector F=n with n → n + 1/2, while those for F 6=n are like those of theeven-sector F 6=n with n→ n+ 1.In the spirit of MA(2), only the EOM for G, which already has onephonon present, is kept exact, while in the EOMs for all the Fn with n ≥ 1we approximate G0(i − j, ω) → δij g¯0(ω). We introduce matrices Wn =(F2n−1, F2n). Again we obtain an equation like Eq. (3.11), where now:γ=11 = 1− g¯0(ω − (4n− 1)Ω)((4n− 1)(2g2 + 6g4+ 6g4(4n− 1))γ=12 = −g¯0(ω − (4n− 1)Ω) (g2 + 6g4 + 4g4(4n− 1))γ=21 = −g¯0(ω − (4n+ 1)Ω)((4n+ 1)2¯(g2 + 6g4)+ (4n+ 1)3¯ · 4g4)γ=22 = 1− g¯0(ω − (4n+ 1)Ω)(4n+ 1)× (2g2 + 6g4 + 6g4(4n+ 1))145C.2. Details for the odd-sectorα=11 = g¯0(ω − (4n− 1)Ω)(4n− 1)4¯g4α=12 = g¯0(ω − (4n− 1)Ω)×((4n− 1)2¯(g2 + 6g4) + (4n− 1)3¯ · 4g4)α=21 = 0α=22 = g¯0(ω − (4n+ 1)Ω)(4n+ 1)4¯g4β=11 = g¯0(ω − (4n− 1)Ω)g4β=12 = 0β=21 = g¯0(ω − (4n+ 1)Ω) (g2 + 6g4 + 4g4(4n+ 1))β=22 = g¯0(ω − (4n+ 1)Ω)g4The matrices for W 6=n are obtained from these by replacing n → n − 1/4everywhere except in the argument of g¯0. The remaining steps are in closeanalogy to those for obtaining the even-sector Green’s function and notreproduced here.The coefficients occurring in the final results for the odd-sector Green’sfunction area=o = 2g2 + 12g4 + (g2 + 10g4)A=1 |1,2 + g4A=1 |2,2a6=o = (g2 + 6g4)A6=1 |1,2 + g4A6=1 |2,2.C.2.2 Momentum space Green’s functionsRather than having the phonon present at a lattice site l, we can constructan electron-phonon state of total momentum K as|K,n〉 =∑ieiKRi/√Nc†ib†i+n |0〉where n is the relative electron-phonon distance. It is easy to show that〈K,m|Gˆ(ω)|K,n〉 = Gi,i+n−m,i+n(ω) exp(iKa(n−m)) where a is the latticeconstant. In particular, the odd-polaron propagator n = m = 0 is just the146C.3. Quadratic e-ph coupling with g2 > 0completely local real space propagator Giii(ω). In other words, the odd-sector polaron shows no dispersion at all.Another Green’s function of interest is given by〈k′, q′|Gˆ(ω)|k, q〉 = 〈0|ck′bq′Gˆ(ω)b†qc†k|0〉where we insert an electron of momentum k into a system where the phononhas momentum q. Conservation of total momentum demands that k + q =k′ + q′. It is again easy to show that the resulting propagator is〈k′, q′|Gˆ(ω)|k, q〉 =δkk′δqq′G0(k, ω˜) +1NG0(k′, ω˜)G0(k, ω) · a=o − a 6=01− g¯0(ω˜)(a=o − a 6=0 ).Since the latter term vanishes in the thermodynamic limit N →∞, we areleft with just the even-sector polaron propagator. This is to be expected:In an infinite system, an electron does not scatter off a single impurity. Ifinstead we assume a finite but low density np of phonons, the prefactor 1/Nin the scattering term is replaced with np.This brief analysis shows that the interesting physics of the odd phononnumber sector are best observed in real space.C.3 Quadratic e-ph coupling with g2 > 0Fig. C.1 shows that for g2 > 0, g4 = 0, the e-ph coupling has an extremelyweak effect even in the atomic limit t = 0, since the quasiparticle weight Zremains very close to 1 while the average number of phonons is very small.An explanation for this behaviour is sketched in Fig. C.2: in the linearHolstein model, the carrier displaces the harmonic lattice potential of itssite, as sketched in the left panel. The overlap between the ground statewavefunctions of the original and the displaced potentials is then the overlapbetween the tails of two Gaussians with different centers, which decreasesexponentially with increasing displacement. Indeed, in the atomic limit for147C.3. Quadratic e-ph coupling with g2 > 00 1 2g20.00.20.40.60.81.0Z0 1 2g20.00.10.20.30.40.5Nph(a) (b)Figure C.1: a) Quasiparticle weight Z, and (b) average number of phononsfor a quadratic model with g2 > 0, g4 = 0 in the atomic limit t = 0, forΩ = 1.i) ii) ionic potentialw/o extra chargewith extra chargeground state wave functionsw/o extra chargewith extra chargeFigure C.2: Sketch of the lattice potential for i) Holstein, and ii) g2 > 0quadratic models. Full (dashed) lines indicate ionic potential and groundstate wavefunction without (with) an extra charge on the site.148C.3. Quadratic e-ph coupling with g2 > 0the linear Holstein model Z ∼ exp[−(g/Ω)2]. In the purely quadratic modelwith positive g2, however, the electron merely changes the shape of the wellby increasing Ω to Ωat. The overlap between the ground states of the originaland modified potential is that of two Gaussians with the same center butdifferent widths. We can calculate this overlap analytically to findZ =√1−(Ω− ΩatΩ + Ω2at)2(C.3)For Ω = 1.0, even for g2 = 100Ω we still have Z ≈ 0.42. We conclude that apositive, purely quadratic electron-phonon coupling has negligible effect onthe dynamics of a charge carrier. In particular, no crossover into the smallpolaron regime occurs for positive g2 for any reasonable coupling strength.Finite t results (not shown) fully support this conclusion.149Appendix DAppendix for T-CuOD.1 Exact Diagonalization with LanczosThe Lanczos algorithm is an efficient method for computing the ground-stateand spectral function of a Hamiltonian, provided that the Hamiltonian issparse, i.e., has a low number of matrix elements [107].As an iterative method, the Hamiltonian H enters the Lanczos methodonly in the computation of matrix-vector products H |ψ〉. For many relevantmodels, the number of non-zero matrix elements per row is a small constantthat does not grow with the system size. Therefore, the computational costof the matrix vector product grows only linearly with the size of the Hilbertspace. In contrast, the diagonalization of a dense matrix has a computationalcost that grows as the third power of the Hilbert space [20]. Here we givejust a very brief introduction to the Lanczos method and refer the reader tothe extensive literature that exists on this algorithm.At its core, the Lanczos method is a variational approach. Let H denotethe matrix representation of some Hamiltonian H. Then, H is diagonalizedin a specially constructed subspace: We start with an arbitrary, randomlygenerated starting vector q0. Then, we setqi+1 =Hqi − αiqi − βi−1qi−1βiwith αi = 〈qi|H|qi〉 , βi = ||Hqi − αiq− βi−1qi−1||2.One can show that in the basis of the subspace spanned by the qi, the150D.1. Exact Diagonalization with LanczosHamiltonian has matrix representation T where T is tridiagonal,T =α0 β1β1 α1 β2. . .. . .. . .βm−1 alpham−1 βmβm αm. (D.1)The advantage of the Lanczos method, then, is that the matrix T can bediagonalized very efficiently with algorithms specialized for tridiagonal ma-trices [20], and that the extremal eigenvalues of T converge rapidly towardsthose of H, with m much smaller than the size of the full Hilbert space.The Lanczos method can also be used to compute the spectral functionby using as the starting vector not a randomly generated vector but insteadthe vector for which the spectral function is to be computed. Details of thismethod can be found in the literature. See for example the lecture notesfrom [111].151"@en ;
edm:hasType "Thesis/Dissertation"@en ;
vivo:dateIssued "2016-09"@en ;
edm:isShownAt "10.14288/1.0300230"@en ;
dcterms:language "eng"@en ;
ns0:degreeDiscipline "Physics"@en ;
edm:provider "Vancouver : University of British Columbia Library"@en ;
dcterms:publisher "University of British Columbia"@en ;
dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@* ;
ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@* ;
ns0:scholarLevel "Graduate"@en ;
dcterms:title "Extensions beyond standard models"@en ;
dcterms:type "Text"@en ;
ns0:identifierURI "http://hdl.handle.net/2429/57841"@en .