Variational studies of dressedquasiparticles’ properties and theirinteractions with external potentialsbyS. Hadi Ebrahimnejad RahbariB.Sc., Tabriz University, 2004M.Sc., Sharif University of Technology, 2007M.Sc., The University of British Columbia, 2009A 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)December 2014c© S. Hadi Ebrahimnejad Rahbari 2014AbstractIn this thesis, we investigated the spectral properties of polaronic quasiparti-cles resulting from the coupling of a charge carrier to the bosonic excitationsof ordered environments. Holstein polarons, describing an electron locallycoupled to dispersionless phonons, and spin polarons in hole-doped antifer-romagnets are considered. The Green’s function of the Holstein polaron ina lattice with various kinds of impurities is calculated using the momentumaverage approximation. The main finding is that the scenario where themere effect of the coupling to phonons is to enhance the electron’s effectivemass is incomplete as the phonons are found to also renormalize the impuritypotential the polaron interacts with. This formalism is applied first to thecase of a single impurity and the range of parameters where the polaron’sground state is localized is identified. The lifetime of the polaron due toscattering from weak but extended disorder is then studied and it is shownthat the renormalization of the disorder potential leads to deviation of thestrong coupling results from Fermi’s golden rule’s predictions.Motivated by the hole-doped cuprate superconductors, the motion ofa hole in a two-dimensional Ising antiferromagnet and its binding to anattractive impurity is studied next, based on a variational scheme that allowsfor configurations with a certain maximum number of spin flips. The roleof the magnetic sublattices in determining the symmetry of the resultingbound states is discussed. Next, a more realistic model describing a hole ina CuO2 layer which retains the O explicitly, is considered. By neglecting thefluctuations of the Cu spins and using a variational principle similar to thatof the previous chapter, a semi-analytic solution for the Green’s function ofthe hole in an infinite 2D lattice is constructed. The resulting quasiparticledispersion shows the proper ground state and other features observed inexperiments. The lack of importance of the background spin fluctuationsis justified based on the hole’s ability to move on the O sublattice withoutdisturbing the Cu spins. Finally, the model is generalized to gauge theimportance of other relevant O orbitals.iiPreface• A version of the work discussed in chapter 2 is published as “H.Ebrahimnejad and M. Berciu, Physical Review B 85, 165117 (2012)”.It is a development of a previous work by Prof. M. Berciu which waspublished in Ref. [1].• The work discussed in chapter 3 is published as “H. Ebrahimnejad andM. Berciu, Physical Review B 86, 205109 (2012)”.• The work discussed in chapter 4 is published as “H. Ebrahimnejad andM. Berciu, Physical Review B 88, 104410 (2013)”. The variationalmethod used here was suggested before by M. Berciu in Ref. [2].• The work discussed in chapter 5 is published as “H. Ebrahimnejad, G.A. Sawatzky and M. Berciu, Nature Physics 10, 951-955 (2014)”. Amore detailed manuscript is under review for publication in anotherjournal. The model we used here was introduced in Ref. [3].I carried out analytical and numerical work for all four projects and wrotethe first draft of the manuscript for the first three, which were then revisedby M. Berciu before submission. The preparation of the draft for the paperpublished in Nature Physics as well as calculations for the three-magnoncase are done by M. Berciu.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Hamiltonian for lattice polarons . . . . . . . . . . . . . . . . 31.1.1 Polarons in systems with other forms of electron-bosoncoupling . . . . . . . . . . . . . . . . . . . . . . . . . 61.2 Studying spectral properties using Green’s function formalism 71.3 Green’s function of the Holstein model and its variationalcalculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.4 Overview of thesis and further reading . . . . . . . . . . . . 92 Trapping of three-dimensional Holstein polarons by variousimpurities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Momentum average approximation for inhomogeneous sys-tems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.3 Polaron near a single impurity . . . . . . . . . . . . . . . . . 212.3.1 Impurity changing the on-site energy . . . . . . . . . 222.3.2 Impurity changing the electron-phonon coupling . . . 292.3.3 Isotope impurity . . . . . . . . . . . . . . . . . . . . . 312.4 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 35ivTable of Contents3 A perturbational study of the lifetime of a Holstein polaronin the presence of weak disorder . . . . . . . . . . . . . . . . 383.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.2 The model and its solution . . . . . . . . . . . . . . . . . . . 393.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.3.1 Weak electron-phonon coupling . . . . . . . . . . . . 453.3.2 Strong electron-phonon coupling . . . . . . . . . . . . 473.4 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 524 Binding carriers to a non-magnetic impurity in a two-dimensionalsquare Ising anti-ferromagnet . . . . . . . . . . . . . . . . . . 554.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.2 The Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.3 Propagation of the hole in the clean system . . . . . . . . . . 594.4 The effect of the impurity . . . . . . . . . . . . . . . . . . . . 664.4.1 Quasiparticle and impurity on the same sublattice . . 674.4.2 Quasiparticle and impurity on different sublattices . . 714.5 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 765 Accurate variational solution for a hole doped in a CuO2layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.2 Three-band model for a cuprate layer and its t-J analogue . 805.2.1 Doping the CuO2 layer with a hole . . . . . . . . . . 845.3 Green’s functions . . . . . . . . . . . . . . . . . . . . . . . . 875.3.1 One-magnon approximation . . . . . . . . . . . . . . 885.3.2 Two-magnon approximation . . . . . . . . . . . . . . 935.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 945.4.1 Dispersion relations . . . . . . . . . . . . . . . . . . . 945.4.2 Effect of various terms in the Hamiltonian . . . . . . 975.5 Quasiparticle weight . . . . . . . . . . . . . . . . . . . . . . . 995.5.1 Comparison with ARPES . . . . . . . . . . . . . . . . 995.6 Spin polaron vs the Zhang-Rice singlet . . . . . . . . . . . . 1025.7 From three-band to five-band model: effect of other in-plane2p orbitals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1055.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1096 Summary and development . . . . . . . . . . . . . . . . . . . 1106.1 Summary of this work . . . . . . . . . . . . . . . . . . . . . . 1106.2 Outlook and further developments . . . . . . . . . . . . . . . 112vTable of ContentsBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114AppendicesA MA0 solution for the clean system . . . . . . . . . . . . . . . 121B IMA1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123C MA0 solution for F (n)k,si(ω) and W nm(ω) . . . . . . . . . . . . . 126D Disorder self-energy in Average T-matrix approximation 128E The equations of motion for G0,R(ω) . . . . . . . . . . . . . . 129F Derivation of ARPES intensity . . . . . . . . . . . . . . . . . 131G Solving for Gα,β(k, ω) in the five-band model . . . . . . . . . 135viList of Figures2.1 Diagrammatic expansion for (a) the disorder GF Gdij(ω) (boldred line); (b) the polaron GF in a clean system G0ij(ω˜) (doublethin line), and (c) the “instantaneous” approximation for thepolaron GF in a disordered system, Gdij(ω˜) (double bold redline). The thin black lines depict free electron propagators,the wriggly lines correspond to phonons, and scattering onthe disorder potential is depicted the dashed lines ending withcircles. See text for more details. . . . . . . . . . . . . . . . . 162.2 (color online) (a) Diagrammatic expansion for the full MAsolution Gij(ω) (thick dashed blue line) in terms of the cleanpolaron Green’s function (double thin line) and scattering onthe renormalized disorder potential ǫ∗l (ω), depicted by dashedlines ended with squares. (b) Diagrammatic expansion ofǫ∗l (ω). For more details, see text. . . . . . . . . . . . . . . . . 192.3 LDOS at the impurity site ρ(0, ω) in units of 1/t, vs. theenergy ω/t, for (a) MA0 with lc = 0 for U = 1.8, 1.9, 1.95 and2.0. The dashed line shows the DOS for the clean system,times 100; (b) Bound polaron peak in MA0 at U/t = 2, andcutoffs in the renormalized potential lc = 0, 1, 2. Panels (c)and (d) are the same as (a) and (b), respectively, but usingMA1. Panel (e) shows the DMC results from Ref. [4], forthe same parameters as (a) and (b). Other parameters areΩ = 2t, λ = 0.8, η/t = 10−3. . . . . . . . . . . . . . . . . . . . 242.4 Phase diagram separating the regime where the polaron ismobile (below the line) and trapped (above the line). Theeffective coupling is λ = g2/(6tΩ), and the threshold trappingpotential Uc is shown for several values of Ω/t. The MAresults (filled symbols) compare well with the DMC results ofRef. [4] (empty symbols). . . . . . . . . . . . . . . . . . . . . 25viiList of Figures2.5 (color online) Effective value of the impurity attraction U∗/U ,extracted from the scaling E∗B/t∗ = f(U∗/t∗), for λ = 0.5, 1.5and Ω = 3t. In order to have the best fit to the data, weplot each curve starting from slightly larger U/t than thecorresponding Uc/t given in Fig. (2.4). For more details, seetext. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6 (color online) Real part of the additional on-site MA0 disorderpotential v0(0, ω) over a wide energy range, for U = 2t,Ω =2t, λ = 0.8 and two values of η. . . . . . . . . . . . . . . . . . 282.7 (color online) Phase diagram separating the regime where thepolaron is delocalized (below the line) and trapped (abovethe line), as a function of the difference between the impu-rity and the bulk electron-phonon coupling, gd − g. Symbolsshow MA0 results, while the dashed lines correspond to theinstantaneous approximation. . . . . . . . . . . . . . . . . . . 302.8 For large λ, in the clean system (dashed red line) the firstpolaron band is separated by an energy gap from the nextfeatures in the spectrum (here, the band associated with thesecond bound state). For gd < g, an “anti-bound” impuritystate is pushed inside this gap (black full line). Parametersare Ω = t = 1, λ = 1.2 and η = 10−3. . . . . . . . . . . . . . . 322.9 (color online) Critical effective coupling λ∗ above which an im-purity state may appear for a sufficiently light isotope. Belowthis line, polarons cannot be bound near isotopes. Symbolsshows MA0 results. The dashed line is the analytic low-boundfor λ∗ discussed in the text. . . . . . . . . . . . . . . . . . . . 332.10 Phase diagram separating the regime where the polaron isdelocalized (below the line) and trapped (above the line), asa function of the difference between the Ωd − Ω, on a loga-rithmic scale. Symbols show MA0 results for Ω = 4t, 8t. AsΩd → ∞, these critical lines converge towards their corre-sponding λ∗ (dashed lines), below which polaron states arealways delocalized. . . . . . . . . . . . . . . . . . . . . . . . . 352.11 LDOS at the impurity site near an isotope with Ωd = Ω + 7t(full line). Two discrete states, one above and one below thebulk polaron band, are seen. For comparison, the DOS inthe clean system (multiplied by 10) is shown as a dashed line.Parameters are Ω = 4t, λ = 2.5 and η = 10−2. . . . . . . . . . 36viiiList of Figures3.1 (a) Inverse polaron lifetime 1/τk vs its peak energy Ek, and(b) average density of states ρ(ω), for a weak electron-phononcoupling and two values of the disorder ∆. The solid anddashed lines are the corresponding Fermi golden rule (FGR)and ATA results, respectively (see text for more details).Other parameters are Ω = t, η/t = 10−2 in (a) and η/t =5 × 10−3 in (b). . . . . . . . . . . . . . . . . . . . . . . . . . . 443.2 Inverse lifetime of the polaron of momentum k = (π/8, 0, 0)and λ = 0.5 vs the strength of disorder, ∆/t (full squares).The inverse lifetime for a free electron with the same momen-tum is shown by triangles. Empty squares show 1/(t∗τk) vs∆/t∗ for the polaron (for this λ, t∗ = 0.881t). Other param-eters are Ω = t, η/t = 10−2. . . . . . . . . . . . . . . . . . . . 463.3 (Color online). Top panels: 1/(tτk) vs Ek/t for three levelsof disorder: ∆/t = 0.05, 0.1 and 0.15. The symbols shows theMA result for a strong coupling λ = 1.2, while the full anddashed lines show Fermi’s golden rule and the ATA predic-tions, respectively. Bottom panels: The average DOS ρ(ω)for that ∆ (full line) and the DOS in the clean system, ρH(ω)(dashed line) vs ω. Parameters are Ω = t, η = 10−3t. . . . . . 473.4 (Color online) The same average DOS vs ω displayed in thelower panels of Fig. 3.3, but now shown for the entire polaronband. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.5 (i) FGR approximation for the disorder-averaged Green’s func-tion of a carrier (double thin line), in terms of that of the freecarrier (thin line) and uncorrelated disorder (dashed line) inthe absence of electron-phonon coupling; (ii) FGR approx-imation for the disorder-averaged GF of a polaron (doublethick line), in terms of that of the clean polaron (thick line);(iii) clean polaron GF in terms of free carrier (thin lines) andphonon (curly lines) propagators; (iv) the first few terms inthe MA approximation for the disorder averaged GF of a po-laron. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.6 (Color online) (a) Inverse lifetime vs disorder, and (b) inverselifetime vs energy Ek, as disorder is turned on, for two mo-menta k1 = (2π/9, 0, 0) and k2 = (π/6, 0, 0), for a polaronwith λ = 1.2,Ω = t, η = 10−3t. The asterisk in panel (b)marks the clean polaron GS energy in the clean system, EGS ,for these parameters. See text for more details. . . . . . . . . 52ixList of Figures3.7 (color online) Same data as shown in Fig. 3.6(a) but withrescaled axes; 1/(t∗τk) vs ∆/t∗ for momentum k2 = (2π/9, 0, 0)(empty squares) is compared with the free electron lifetime(triangles) at the same momentum, for λ = 1.2 where t∗ =0.148t. Other parameters are Ω = t and η/t = 10−3. . . . . . 534.1 Square lattice of Nee´l-ordered Ising spins. Removing one ofthe spin-up particles creates a hole and results in a configu-ration with sz = −12 . . . . . . . . . . . . . . . . . . . . . . . 594.2 Effective first- and second-nearest-neighbor hoppings of thehole (the blue circle) achieved with loops involving only upto three magnons. The latter is realized when the hole startsa second loop before removing the last magnon it createdduring the first loop. The magnons are shown by red circles.The properly oriented spins are not shown. . . . . . . . . . . 614.3 The choice of coordinate systems for the lattice with impurity.The impurity, shown as purple square, is at the origin of thexy axes that span the original lattice with unit vectors x,y.The XY axes are rotated by 45 ◦ and span the sublattice(black dots) on which the quasiparticle propagates via theelementary vectors y± x. . . . . . . . . . . . . . . . . . . . . 634.4 The total density of states (top panels), effective hoppingst1(ω), t2(ω) (middle panels) and on-site energy ε(ω) (bottompanels) in the clean system for two different values of t/J¯ .The effective parameters are relatively constant within theenergy band, explaining why the DOS has the generic formexpected for a bare particle with nearest-neighbor hopping ona square lattice. Here J¯ = 1, η = 10−3 and t = 3 (left panels)and t = 6 (right panels), respectively. . . . . . . . . . . . . . . 664.5 (Top): the energy dispersion and the quasiparticle weight forthe lowest polaron band. The ground state is at k = (0, 0),while k = (pi2 , pi2 ) is a saddle point. The top-right quadrantof the full Brillouin zone (BZ) and that of the correspond-ing magnetic BZ (dashed green) is shown at the inset. Theband folding due to AFM order leads to symmetric dispersionalong (0, 0)− (π, π). The quasiparticle weight has only minorvariations which can be related to the fact that the polaronbandwidth (≈ 0.3) is considerably smaller than the energy ofa typical spin defect (few J¯) in the magnon cloud. t = 3J¯ = 3and η = 10−3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67xList of Figures4.6 (Top) LDOS at the impurity site for various values of U . Thedashed line is the DOS in the clean system, times 4. At finiteU , a single bound state splits from the continuum and itsbinding energy increases with U . Curves are shifted verticallyto help visibility. (bottom) The effective on-site energy at theimpurity site is essentially equal to U . Parameters are t = 6,J¯ = 1 and η = 10−3. . . . . . . . . . . . . . . . . . . . . . . . 684.7 Relative amplitude of the wave functions corresponding tothree of the bound states shown in Fig. 4.6, at various dis-tances from the impurity site. Lines are exponential fits.States with bigger binding energies have shorter decay lengths. 704.8 (Top) LDOS ρ(r = x;ω) with curves shifted vertically to helpvisibility; (Bottom) Ueff(r = x;ω) at the quasiparticle’s sub-lattice site located nearest to the impurity. Up to three boundstates split from the continuum upon increasing U . The pres-ence of the impurity at r = 0 induces a finite effective on-siteattraction at r = x, whose value is significantly smaller thanU (Bottom). Parameters are t = 6, J¯ = 1, and η = 10−3. . . . 714.9 Relative amplitude of the upper and lower bound states forU = −2J¯ , at various distances from the impurity site. Linesare exponential fits. . . . . . . . . . . . . . . . . . . . . . . . 724.10 Binding energy of the s-wave bound state at U = −J¯ vs. t/J¯ ,when the quasiparticle and the impurity are on the same sub-lattice (top panel) and different sublattices (bottom panel).The smaller binding energy at strong hopping is due to thereduction in the quasiparticle’s effective mass, which makesit harder to trap. When the quasiparticle and the impurityare on different sublattices, the enhancement of Ueff at smallt dominates over the effective mass decrease, explaining thegrowth of the binding energy here. . . . . . . . . . . . . . . . 744.11 The gap between the px and either of the s or d states, for afixed U . Its enhancement as a function of t/J¯ reflects the ro-tational kinetic energy gain of the quasiparticle as it becomeslighter with increasing t. . . . . . . . . . . . . . . . . . . . . . 75xiList of Figures5.1 (a): Cu, shown in red, is surrounded by six O ligands inoctahedral geometry. The static electric field of these lig-ands partially removes the degeneracy of Cu 3d orbitals intoeg and t2g manifolds. This degenerate configuration distortsto remove the remaining degeneracies and ends up having3dx2−y2 as its highest partially-occupied orbital, (b). Thevalues of splittings for La2CuO4 are given. In the case ofSr2CuO2Cl2, splitting within eg is such strong that 3d3z2−r2orbital is pushed below the entire t2g manifold. . . . . . . . . 815.2 (a) The ground state of charge-transfer insulator (left) andits lowest excitation (right) ∆pd higher in energy. ↑ and ↓represent electrons and © stands for a hole (electron repre-sentation). (b) In a Mott insulator, however, charge excita-tion with the lowest energy corresponds to an electron movingfrom one Cu to another at an energy cost of Udd. . . . . . . . 825.3 One set of processes giving rise to the superexchange cou-pling, Jdd, between two neighbouring Cu holes. Arrows rep-resent the holes and empty O levels corresponds to them beingdoubly-occupied with electrons (hole notation). The energiesof intermediate configurations are given. . . . . . . . . . . . . 835.4 Energy diagram of the undoped CuO2 layer, illustrating acharge-transfer insulator. The filled bands are shaded, abovethem lies the Fermi energy in the gap. Since ∆pd < Udd,removing electrons created holes on O’s rather than on Cu’sand this moves the Fermi energy into the O 2p-band. . . . . . 845.5 (a): The CuO2 layer with the spin-up hole occupying one ofthe O 2p orbitals. nn (δ) and next nn hops (ǫ) between O’sthat are considered in our model are shown. For the hole, nnhopping is negative along y = x direction where 2p orbitalshave their lobes with similar phases pointing towards eachother and positive along y = −x. (b): Cu and O ions thatmake the unit cell are shown in green. (c): Processes thatgive rise to nn (c1) and next nn (c2) spin-swap. (d): The fullBrillouin zone of the CuO2 lattice (the outer square) whichencloses the magnetic Brillouin zone (shaded). We study thequasiparticles with momenta lying along the dashed red line,starting from (0, π) and ending at (π, 0). . . . . . . . . . . . . 86xiiList of Figures5.6 Numbering various one-magnon Green’s function VMi at M= 1, 3: the one with the hole lying to the far left is i = 0 andthe index i increases in the counterclockwise direction. Thecurrent configuration corresponds to V 21 . . . . . . . . . . . . . 895.7 Numbering various two-magnon GFs at M = 1, 3. The con-figuration shown corresponds to W 23 . . . . . . . . . . . . . . . 925.8 The spectral weight at nm = 2 approximation, Apx1 (k, ω) incolor scale as a function of momentum and energy (in unitsof Jdd ≈ 0.1 − 0.15 eV), showing several coherent bands. Forthe rest of this chapter, only the lowest energy band will bestudied. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.9 (a): Energy dispersion (in units of Jdd ≈ 0.1−0.15 eV) of thehole without allowing for the creation of magnons, nm = 0.The ground state is still at k = (pi2 , pi2 ). (b): Comparisonbetween one- and two-magnon approximation results for thehole’s energy dispersion along various cuts in the full BZ. Thedots are relevant exact-diagonalization calculations for a finitesize sample of the same model, shifted in energy in order toease the comparison. . . . . . . . . . . . . . . . . . . . . . . . 975.10 The lowest quasiparticle dispersion (nm = 2 and in units ofJdd ≈ 0.1 − 0.15 eV) when various terms in the Hamiltonian(5.3) are on and off. Jpd has almost no effect on the dispersion(a), whereas without ts (c) the model fails to properly predictthe expected ground state k = (pi2 , pi2 ). For comparable valuesof ts and tpp, the correct shape of energy dispersion re-appears(d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.11 top: quasiparticle weights extracted from the sublattice GFApx1 (k, ω) in one- and two-magnon approximation. The cor-responding energy dispersion (for the two-magnon approxi-mation, in units of Jdd ≈ 0.1 − 0.15 eV) is presented in thelower panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.12 The quasiparticle weights in nm = 2 approximation extractedfrom spectral weights that correspond to projections onto or-bitals with different polarizations, 2px and 2py. . . . . . . . . 1015.13 Comparison between one- and two-magnon approximation re-sults for ARPES quasiparticle weight. This shows disagree-ment with experiments for momenta along (0, 0)− (π, π). Seetext for further discussions. . . . . . . . . . . . . . . . . . . . 102xiiiList of Figures5.14 Inside the dashed circle are the copper and oxygens involvedin Zhang-Rice singlet. Labelling of one-magnon GFs are suchthat once the copper and hole’s spin are flipped, the currentconfiguration becomes V 21 . The sign of each singlet contribu-tion to ZRS is given. . . . . . . . . . . . . . . . . . . . . . . . 1035.15 top: overlap between the states in the lowest polaron band|φ〉 whose energy dispersion (in units of Jdd ≈ 0.1−0.15 eV) isshown at the lower panel, and the corresponding Zhang-Ricesinglet states. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045.16 The CuO2 layer with two in-plane 2p orbitals per O. The nnhopping between new (blue) and original (green) orbitals, t˜pp,and the next nn hopping between the new orbitals, t˜′pp, areshown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1065.17 Energy dispersion (in units of Jdd ≈ 0.1 − 0.15 eV) andARPES quasiparticle weight of the five-band model that in-cludes two 2p orbitals per O. Whereas the energy dispersionsare qualitatively similar to those given by the three-bandmodel, the quasiparticle weight is now less sensitive to thestructure of magnon cloud. . . . . . . . . . . . . . . . . . . . 1075.18 Constant-energy map of ARPES spectral weight for the five-band model. It clearly shows the asymmetry of the quasipar-ticle weight between the inside and outside of the magneticBZ, which is a robust feature for the generic values of param-eters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108F.1 The effect of finite K/λ (lowest order correction) on the ARPESquasiparticle weight of the five-band model. The differencebetween the quasiparticle dispersions is invisible on the sameenergy scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . 133xivAcknowledgementsI would like to use this opportunity to thank my adviser, Prof. Mona Berciu,for all her efforts in guiding me through my PhD years, a task that she didwith patience, enthusiasm and good manners. It was through working withher that I learned how asking the right questions is an important trait forany serious researcher, a quality she is also quite good at it. She gave me theright problems which eventually resulted in successful projects that madeparts of my thesis. I am grateful for her always making herself available byall means, making it possible for me to better focus my attention and effortstowards accomplishing my research goals in a timely manner.I would also like to acknowledge the members of my PhD committeefor attending scheduled meetings, reading the thesis and their suggestionsfor its improvement. Special thanks goes to Prof. George Sawatzky for hiscollaboration and insightful remarks which turned the project discussed inchapter five into a more impactful piece of work.Last but not least, I thank my parents and my brother Habib, who wasmy inspiration to pursue physics, for their support in all stages of my life.Finally, I gratefully thank my wife, Parisa, for all her support, caring andlove which has been a real source of encouragement to me and will certainlybe so for years to come.xvChapter 1IntroductionUnderstanding various macroscopic physical properties exhibited by solidstate materials generally requires studying their microscopic structure, thedynamics of their constituent electrons and their interaction with the restof the solid. The simplest description of electrons inside a solid (apart fromthe particle in a box model) is given by the band theory in which mobileelectrons are thought of as interacting with a perfectly periodic lattice ofimmobile ions1. In this approximation, the ε(k) ∼ k2 energy dispersion ofan electron in free space is replaced by a series of energy bands, εn(k), thatshow the periodicity of the corresponding reciprocal lattice [5]. If one canignore electron-electron correlations, the many-electron states are then con-structed by populating these single-particle energy bands according to Fermistatistics which forbids multiple occupancy of any quantum state. The bandtheory successfully explains why some solids are metallic whereas others areelectrical insulators: insulators are those that have their highest occupiedenergy band completely filled, as electrons in completely filled bands do notcarry a net electrical current and exciting them to empty bands costs a finiteamount energy. Materials with partially filled valence bands are metallic,since a net current results when electrons are excited into empty states ofthe same band with an infinitesimal energy cost.With many electrons to fill the energy bands, one has to account for theCoulomb repulsion between them beyond mean-field approximations. Atlow temperatures and for weak electron-electron interactions, it results inscattering of electrons into low-lying unoccupied states near the Fermi level.The Fermi liquid theory predicts a one-to-one correspondence between thek-states of the interacting system and those of the non-interacting system[6]. This analogy introduces the notion of Landau quasiparticles as the ele-mentary excitations in systems of interacting fermions. These quasiparticleshave a different effective mass m∗ 6= mel and, due to phase space restrictions,they are much more weakly interacting than the original electrons the closerthey get to the Fermi surface where they become long-lived and accord-1The potential includes a mean-field (Hartree-Fock) contribution from the other elec-trons in the solid.1Chapter 1. Introductioningly well-defined. This gives a single particle description of the interactingsystem where the quasiparticles replace the bare electrons.In the strongly correlated systems, however, the motion of individualelectrons depends on the state of the other electrons. This is generally thecase for systems where the interaction between electrons is strong comparedto their kinetic energy that enable electrons to avoid each other. The sup-pression of kinetic energy in such systems gives rise to a variety of fascinat-ing ordered states such as those observed in magnetic materials, fractionalquantum Hall and other interacting low-dimensional systems and variousunconventional superconductors. The quasiparticles in these systems haveeither exotic properties or their mere existence is still debated.Coming back to the single electron limit, the assumption that the ionsare at rest violates the Heisenberg uncertainty principle by neglecting theirzero-point motion. Furthermore, even though ions are much heavier thanelectrons, the attractive interaction among the two causes the lattice to dis-tort. The excess density of positive charges in the distorted part of thelattice acts like a local potential well for the electron. It was even proposedthat in a strongly distorted lattice an electron traps itself into the poten-tial well it creates [7]. In a clean lattice, however, translational invarianceprevents the electron from being localized at any particular site. Therefore,it remains mobile but its effective mass becomes larger due to the latticedistortions it drags behind.In quantum mechanical terms, where the lattice distortions are describedas phonons, a many-body picture emerges where the electron plus the ac-companying lattice distortion is viewed as a quasiparticle comprised of theelectron surrounded by a fluctuating cloud of phonons. This is known asa lattice polaron. They are a useful way of thinking about charge carriers(electrons or holes) in ionic crystals and organic semiconductors [8].In many-electron systems, the lattice distortion created by an electroncan attract another electron and result in an effective attraction betweenthe two electrons which, in suitable conditions, can bind them together intoa bipolaron state. If the electrons are much faster than the lattice ions,however, the distortion created by the first electron persists even after itmoves away. When the second electron arrives in this region, the first one isalready gone. Therefore, there will not be a strong repulsive force betweenthem and the effective attraction prevails even if it is rather weak. This mayresult in bound states in the momentum space, well-known for their role inestablishing conventional (BCS) superconductivity [9].21.1. Hamiltonian for lattice polarons1.1 Hamiltonian for lattice polaronsThe generic Hamiltonian for a system of an electron in a lattice isH = − ~22m∇2r +∑jVel−a(r−Rj) + Hlat, (1.1)where r is the position of the electron, Rj are those of the lattice sites andHlat is the lattice Hamiltonian which, for harmonic lattice deformations, canbe cast into the usual phonon Hamiltonian (see below). Vel−a(r−Rj) is theinteraction between the electron (with the bare mass m) and the latticeion at Rj . For small lattice deformations, uj = Rj −R0j ≪ a, where a isthe lattice constant and R0j are the equilibrium position of lattice sites, thesecond term can be approximated to its first-order expansion in ujH = − ~22m∇2r +∑jVel−a(r−R0j)︸ ︷︷ ︸Hb−∑juj .∇rVel−a(r−R0j) + Hlat. (1.2)This linear approximation is sufficient for most situations. However, it isexpected to fail when the coupling strength between the electron and thelattice is too strong [10] or the physical system has some kind of symmetrythat forbids the linear dependence of energy on lattice distortions [11]. Thefirst two terms which describe the electron in a periodic potential are theonly ones considered in the band theory: Hbψk(r) = ε(k)ψk(r), where theeigenstates are Bloch waves. The spin is irrelevant for a single electron (po-laron) and its index is therefore dropped. In the tight-binding approximationand using the language of second quantization, Hb is given as∑k ε(k)c†kck,where for a simple cubic lattice in d dimensions, ε(k) = ε0−2t∑di=1 cos(ki).t is the hopping amplitude and it describes the possibility of tunnelling be-tween overlapping atomic states and ε0 is the on-site energy. k is known asthe crystal momentum and it is unique only inside the first Brillouin zone(BZ), −π/a < ki ≤ π/a.The next term describes the coupling to lattice distortions and its secondquantized form depends on the specific from of the potential Vel−a. However,it can be written asVel−ph =∑i,jg(|R0i −R0j |)c†i ci(b†j + bj) (1.3)since it couples the lattice distortion at R0j , uj ∝ (b†j + bj), to the electrondensity at r = R0i and is a function of r − R0j . In momentum space, the31.1. Hamiltonian for lattice polaronscorresponding electron-phonon couplings strength depends on the phononmomentum q onlyH =∑kε(k)c†kck +1√N∑k,qg(q)c†k+qck(b†−q + bq) +∑qω(q)b†qbq, (1.4)where N → ∞ is the number of lattice sites and the total momentum isassumed to be conserved. bq is the phonon annihilation operator. It is as-sumed that there is only a single branch of phonons with energy dispersionω(q), but generalizations to include multiple phonon modes is straightfor-ward [12]. In this thesis we only focus on dispersionless optical phonons,ω(q) = Ω, which in the long wavelength limit, correspond to out of phasemotion of the neighbouring lattice sites. They are the relevant vibrationalmodes of ionic crystals, and are called optical because they can be excitedby light and result in oscillating dipole moments. With further assump-tion of constant electron-phonon coupling, g(q) = g, and nearest-neighbourhopping only, this reduces to the Holstein model [13] and in real space itcorresponds to coupling of an electron to lattice distortions at the same siteHH = −t∑<i,j>(c†i cj + H.c.) + g∑ic†i ci(b†i + bi) + Ω∑ib†i bi. (1.5)This model was initially proposed to study molecular crystals with trans-verse vibrational degrees of freedom. Although it looks quite simple, thereis no exact solution valid for the whole range of parameters. It has beenstudied using a variety of numerical and approximate analytical methods,including the Momentum-Average approximation (MA) [14] which we willuse in this thesis. It shows the generic features associated with the polaronphysics, such as the dependence of polaron properties on the structure ofaccompanying boson cloud and the crossover from a large to small polaronas the strength of electron-phonon coupling increases. When the couplingis large, the electron strongly distorts the lattice at the same site and thepolaron becomes very heavy. This is a small polaron, whose phonon cloud islimited to the close vicinity of the electron. Weak electron-lattice coupling,on the other hand, only causes slight lattice distortions that are extendedover several lattice sites. The polaron is accordingly much lighter in thiscase.Note that the previous derivation, although typical of that found in text-books, ignores a second important type of el-ph coupling: phonons can alsocouple to the kinetic energy of the electron by modifying its hopping ampli-tude. The latter depends exponentially on the distance between sites, which41.1. Hamiltonian for lattice polaronsvaries as a result of lattice vibrations: a → a + ui+1 − ui (one dimension).Such effects give rise to new coupling functions and will be encountered ina different context within this work (see chapter 4).The Holstein model has exact solutions for the asymptotic limits of in-finitely weak or strong coupling. For g = 0, there are no phonons and theGreen’s function is that of the free electron in the lattice,G(k, ω) = 1ω − ε(k) + iη (1.6)where ε(k) is the tight-binding energy dispersion. The other limit is the socalled impurity limit, t = 0, where the solution is given by the Lang-Firsovformula [15],G(ω) = e−g2/Ω2∞∑n=01n! (gΩ)2n 1ω − nΩ + g2/Ω + iη . (1.7)With t = 0 the polaron is localized at one site and has infinite mass; theGreen’s function is therefore independent of the momentum. The phononnumber obeys a Poisson distribution, |ψ(n)|2 = e−g2/Ω2 1n!(gΩ)2n, and there-fore the structure of phonon cloud is that of a coherent state, b|ψ〉 = gΩ |ψ〉.The spectrum is quite different for these two limits, going from a contin-uous band of energies ε(k) into a set of discrete states ε(n) = nΩ−g2/Ω, andit is of great interest to understand how the ground state and other excitedstates evolve as both g and t become finite. Perturbation theory can beused to extend the applicability range of the asymptotic results, but it failsfor the intermediate parameter regimes. Various numerical schemes, such asthe exact diagonalization (ED) [16–18] and quantum Monte Carlo (QMC)method [19–21], have been used to study Holstein polarons in this case. TheMA approximation, which we will use in this thesis, is a semi-analytic ap-proach that works reasonably well for the entire parameter range. It has avariational interpretation and its accuracy can be systematically improved.Other examples of couplings include the Fro¨hlich [22], g(q) ∝ 1/q, thedeformation coupling, g(q) ∝ q and the breathing mode Hamiltonian forwhich g(q) ∝ sin(q/2) in one dimension [23, 24].All this discussion assumes that the typical size of the polaron cloud iscomparable with the lattice constant a. A different limit that has receivedmuch attention is that of the large polaron [25], whose cloud extends overa length scale exceeding a by orders of magnitude, so that the electron’sdynamics can be treated in the continuous approximation where its kinetic51.1. Hamiltonian for lattice polaronsenergy has again quadratic dependence on the momentum (but with a renor-malized mass). While this is an interesting topic on its own, here we onlyconsider lattice models and small polarons.1.1.1 Polarons in systems with other forms ofelectron-boson couplingIn terms of the bare electron parameters, the electron-phonon coupling con-sidered in Eq. (1.3) can be interpreted as the one modifying the on-siteenergy of electron since it is proportional to c†i ci. This causes the wavefunc-tion of atomic orbitals occupied by the electron to change accordingly. Asa result, the hopping amplitudes which are related to the overlap betweenthese orbitals are also modified due to coupling to phonons. In the linearapproximation, this gives rise to an additional electron-phonon coupling ofthe following formH′el−ph =∑<i,j>,ltij,lc†i cj(b†l + bl) + H.c. . (1.8)When expressed in the momentum representation, it corresponds to a cou-pling strength that depends on the momenta of both electron and phonon,g(k,q). Therefore, the coupling to phonons in this case modifies the kineticenergy of the electron. This type of coupling is encountered, for example,in the SSH model of conducting polymers [26, 27].Although we introduced polarons in the context of electron-phonon cou-pling, one has to note that polaron physics is more general than this andis relevant to any situation where a charge carrier interacts with the col-lective excitations of an ordered environment. In cuprate superconductors,for example, doping holes move in the CuO2 layer in which Cu spins areanti-ferromagnetically ordered. This motion disturbs the order of Cu spinsand creates (bosonic) magnetic excitations, known as magnons, that dressthe hole in the same way as phonons did with an electron in a distortedlattice. This is known as a spin or magnetic polaron and, based on differentmodels and variational methods inspired by MA (see below), we will studythem in chapters 4 and 5 of this thesis. We will see that this is a case ofboson-modulated hopping.61.2. Studying spectral properties using Green’s function formalism1.2 Studying spectral properties using Green’sfunction formalismGreen’s functions are convenient mathematical tools for studying variousstatic or dynamical properties of quantum systems, such as energy eigensys-tems and response functions. They are particularly useful for many-bodyproblems where there are generally no exact solutions and perturbation orsome other approximation is needed. There are different types of Green’sfunctions, but the one we exclusively use in this thesis is the one-particleretarded Green’s function defined as followsG(k, τ) = −i〈0|T [ck(τ)c†k(0)]|0〉= −iΘ(τ)〈0|ck(τ)c†k(0)|0〉 (1.9)where ck(τ) = eiHτ cke−iHτ is the fermion annihilation operator in Heisen-berg picture and |0〉 is the vacuum state. T is the time-ordering operatorand Θ(τ) is the Heaviside function, reflecting the fact that ck|0〉 = 0. Theusefulness of this quantity becomes apparent in its Lehmann representation,G(k, ω) =∫ ∞−∞dτeiωτG(k, τ) = limη→0+∑α|〈ψα|c†k|0〉|2ω − Eα + iη, (1.10)where α are quantum numbers that identify a complete set of one-particleeigenstates {|ψα〉} of H with energies {Eα}: H|ψα〉 = Eα|ψα〉. In a closedsystem where the total momentum is conserved, [H,P] = 0 and |ψα〉 istherefore a momentum eigenstate with the same k, otherwise 〈ψα|c†k|0〉 = 0.Eq. (1.10) shows that the poles of the Green’s function give the energyspectrum of the Hamiltonian and the associated residues contain informationabout the energy eigenfunctions. In many-body problems, the latter is alsoknown as the quasiparticle weight, Zk,α = |〈ψα|c†k|0〉|2, and it correspondsto the overlap between eigenstates of H and those of the non-interactingproblem. It hence contains information about the degree of renormalizationof states of the system due to many-body interactions. It can also be used toextract various information about the quasiparticles such as their effectivemass and the structure of the accompanying boson cloud [14]. In addition,the imaginary part of the Green’s function is the spectral functionA(k, ω) = − 1π ImG(k, ω) =∑αZk,αδ(ω −Eα) (1.11)and, in angle-resolved photo-emission experiments, it is related to the cur-rent of outcoming photo-electrons with momentum k [28]. The total density71.3. Green’s function of the Holstein model and its variational calculationof states is given by integrating A(k, ω) over entire momenta in the first BZρ(ω) =∫BZdkA(k, ω). (1.12)When there is disorder in the lattice, the crystal momentum k is not con-served and one has to work in the real-space basis,Gij(ω) = 〈0|ciGˆ(ω)c†j |0〉 (1.13)where Gˆ(ω) = 1/(ω −H + iη) is the resolvent of Hamiltonian, H. It givesthe local density of state (LDOS)ρ(Ri, ω) = −1π ImGii(ω) =∑α|ψα(Ri)|2δ(ω − Eα), (1.14)which can be studied to find the spatial profile of eigenstates ψα(Ri). Wewill use this quantity in order to identify localized states of polarons nearvarious types of impurities.1.3 Green’s function of the Holstein model andits variational calculationVariational methods will be extensively used throughout this thesis as a wayof doing nonperturbative calculations. The general idea is to exclude statesthat are unlikely to be realized by the system, due to energetic reasons forexample, and only keep a particular subspace of the Hilbert space that in-volves the more relevant states. This is generally accomplished by separatingthe Hamiltonian into an easy part which can be trivially diagonalized, andthe rest of the Hamiltonian which usually involves many-body interactions,H = H0 + H1, and then using the Dyson’s identity:Gˆ(ω) = Gˆ0(ω) + Gˆ(ω)H1Gˆ0(ω), (1.15)where Gˆ0(ω) = 1/(ω −H0 + iη) is the resolvent associated with the easypart, H0. Acting by H1 upon any state that is within the subspace of interestwill eventually link to states that lie outside it. Discarding such states resultsin a simplified set of equations between various Green’s function, which canbe solved efficiently to find the desired Green’s functions. One advantage ofthis variational approach is its controlled accuracy: one can always enlargethe variational space until an acceptable level of accuracy is achieved.81.4. Overview of thesis and further readingThe Green’s function of the Holstein model has been evaluated using avariational scheme based on this general approach, known as the momentum-average approximation (MA) [14], where the spatial extent of the phononcloud can be varied to improve the accuracy. MA has been successful inaccurately calculating various properties of the Holstein polaron in cleansystems. In the next two chapters, we will extend it to systems with variouskinds of disorder.1.4 Overview of thesis and further readingSo far, we have assumed that the lattice in which the electron moves is cleanof any sort of disorder and its only deviation from complete periodicity isdue to lattice distortions. In reality, however, disorder is introduced intolattice structures during sample growth and preparation and it hinders themotion of electrons even in the absence of lattice distortions. When dis-order is weak compared to kinetic energy of electrons, the main effect ofdisorder is to scatter electrons between different Bloch states and this limitsthe lifetime of these states. For strong disorder, however, it was proposedin seminal work by Anderson [29] that, in three spatial dimensions, elec-trons with certain energies become exponentially localized in finite regionsof the lattice as a result of destructive interference between multiply scat-tered electron waves from various impurity sites. As such localized electronslack itinerant properties, the solid turns into an electrical insulator whenlocalization happens for states near Fermi energy [30]. While it takes a fi-nite amount of disorder to localize all electronic states in three dimensions,this is achieved for any amount of disorder in one and two dimensions 2.Once again, this picture applies to bare electrons which move freely in thelattice in the absence of disorder. In systems of strongly interacting elec-trons, there are other insulating states which arise even in the absence ofany disorder. They are known as Mott insulators in which the electrons areconfined into atomic-like states due to strong repulsion with other electrons.Therefore, the charge fluctuations become strongly suppressed and magneticorder develops at low energies.Compared to bare electrons, polarons are expected to be more suscepti-ble to external disorder because of their enhanced effective mass. In otherwords, the same amount of disorder is effectively stronger when measuredagainst the polaron’s reduced kinetic energy. Although this is in principletrue, it implicitly assumes that the heavier polaron interacts with the same2In two dimensions, states are only logarithmically localized when disorder is week [31].91.4. Overview of thesis and further readingdisorder potential as the bare electron does.In chapter 2, based on the extension of the momentum-average approxi-mation to inhomogeneous systems, we will challenge this picture by closelyexamining the interplay between disorder and coupling to phonons for thecase of Holstein polaron [32]. As the electron-phonon coupling becomesstronger, the simple picture where the effect of coupling to phonons can becast into the enhancement of the effective mass starts to fail and modifica-tions in the form of renormalization of the disorder potential by phononsbecome important. These modifications reflect the fact that, due to thebinding between the electron and its phonon cloud, the electron’s abilityto interact with existing disorder is affected in a way that depends on thestructure of the cloud surrounding it. This reveals previously unexplored as-pects of interaction of polaronic quasiparticles with external potentials andmakes predictions that defy naive expectations [33]. This formalism will beapplied to study bound states of a three-dimensional Holstein polaron nearvarious kinds of single impurities. The accuracy of MA approximation isconfirmed by comparing our results to that of exact numerical treatments,when available. For a detailed review of the historical development of latticepolarons and various techniques developed for studying them, including theMA approximation, the reader may refer to Ref. [34] and [35].We continue on this track in chapter 3, where the effect of weak butextended disorder on the lifetime of Holstein polarons is studied. Differentbehaviours are recovered in the weak and strong coupling regimes and therole of disorder renormalization and the validity of results at strong-couplinglimit will be discussed.Next, we study the spin polaron description of holes doped into antifer-romagnetic insulators. The motivation is to study the spectral propertiesof a single charge carrier in hole-doped copper-oxide superconductors. Thisis an intrinsically multi-band problem where there are several Cu and Oorbitals involved, but in chapter 4 we use an effective one-band description(based on the Zhang-Rice singlet picture [36]) in order to study the trappingof polarons of this model near an attractive impurity. Renormalization ofimpurity potential due to coupling to magnons is shown to be relevant tothis problem as well and the role of magnetic sublattices in this regard willbe discussed.In the last chapter, a more realistic model for a hole in the CuO2 layeris considered. It has been originally proposed and solved using exact diago-nalization techniques in Ref. [3]. Here, we propose a semi-analytic solutionwhich is based on a variational principle similar to that used to generate theMA approximation for lattice polarons. Our approach will prove to be very101.4. Overview of thesis and further readingfast and efficient numerically and, unlike the ED method, it does not sufferfrom finite-size effects. We recover the hole’s energy dispersion in agreementwith ARPES measurements of the parent compound and determine the roleof various terms in the Hamiltonian in giving rise to its generic features.The efficiency of our method enables us to include orbitals that are usu-ally neglected in other studies of cuprates and gauge their effects on thequasiparticle’s spectral properties.11Chapter 2Trapping ofthree-dimensional Holsteinpolarons by variousimpurities2.1 IntroductionThe challenge to understand the effects of disorder on the behavior of par-ticles strongly coupled to bosons from their environment is commonly en-countered in correlated electron systems. For example, high-Tc cupratesare doped anti-ferromagnetic insulators in which, beside strong couplingto magnons, ARPES measurements have also provided evidence for strongelectron-phonon coupling [37]. At the same time, charge carriers moving inthe CuO2 layers are subject to random disorder potentials from the adja-cent dopant layers. Substituting only a few percent of the Cu atoms withimpurities suppresses superconductivity by localizing the low energy elec-tronic states [38]. Inhomogeneities in the superconducting gap measured inhigh resolution tunnelling experiments have been attributed to atomic scaledisorder in the phonon energy and the electron-phonon coupling strength inthese materials [39].The result of the interplay between disorder and coupling to such phononsdepends on their relative strengths. Disorder that is considered weak forfree electrons can be strong enough to localize a polaron, that is the dressedquasiparticle which consists of the electron together with its cloud of phonons,because of its heavy effective mass. On the other hand, whereas in the weakdisorder regime electron-phonon coupling hinders the motion of electrons,such coupling actually facilitates the electron mobility in the strongly dis-ordered regime where the Anderson localization prevails [40].Certain aspects of this problem have been studied with various approxi-mations, most of which rely on sophisticated non-perturbative methods [41]122.2. Momentum average approximation for inhomogeneous systemssuch as the statistical dynamic mean field theory (DMFT) [42, 43], or dy-namical coherent potential approximation (DCPA) [44, 45].In this chapter we focus on Holstein polarons which are quasiparticlescomprising a charge carrier plus a cloud of Einstein phonons describing thelattice distortions in the vicinity of the charge carrier [13]. We considerthe simplest case of a single impurity and study the formation of polaronicbound states near that impurity. This problem has been studied using var-ious numerical techniques, such as diagrammatic Monte Carlo (DMC) [4]and continuous time quantum Monte Carlo techniques [46]. Here, we usea generalization of the Momentum Average (MA) approximation to inho-mogeneous systems [1] in order to study this problem. The method can, inprinciple, be used to study all possible types of disorder for various kindsof electron-phonon coupling. Here, we stick to the Holstein-type coupling,but consider different types of disorder due to single impurities, namely avariation in the on-site energy, in the electron-phonon coupling and/or inthe phonon energy, whereas the previous studies only considered impuritiesof the first type. The accuracy of our method for this case is demonstratedby comparison with DMC results.The chapter is organized as follows: In Sec. 2.2 we describe our methodand discuss its meaning (full details are provided in the Appendix A andB.). Sec. 2.3 presents our results for the three types of impurities, and Sec.2.4 contains the summary.2.2 Momentum average approximation forinhomogeneous systemsWe start with deriving the MA approximation for a lattice that has disorderin the form of random on-site energies, inhomogeneities in the strength ofelectron-phonon coupling and phonon frequencies. A simpler case (withonly on-site disorder) has been briefly discussed in Ref. [1]. The HolsteinHamiltonian for such a lattice isH = Hd + Vˆel−ph = H0 + Vˆd + Vˆel−ph, (2.1)where Hd describes the non-interacting part of the Hamiltonian, includingH0 = −t∑〈i,j〉(c†i cj + H.c.) +∑iΩib†i bi132.2. Momentum average approximation for inhomogeneous systemswhich contains the kinetic energy of the particle and energy of the bosonmodes (~ = 1), andVˆd =∑iǫic†i cidescribing the on-site disorder. The interaction partVˆel−ph =∑igic†i ci(b†i + bi)describes the inhomogeneous Holstein-like coupling [13] between the particleand the bosons. Here, i indexes sites of a lattice which can be in anydimension, and of finite or infinite extent. The operators ci and bi describe,respectively, particle and boson annihilation from the corresponding stateassociated with lattice site i. The spin of the particle is ignored because itis irrelevant in this case, but generalizations are straightforward [47]. Forsimplicity, we assume nearest-neighbour hopping; this approximation canalso be trivially relaxed. Depending on the model of interest, the on-siteenergies ǫi, electron-phonon couplings gi and phonon frequencies Ωi can beassumed to be random variables.Our goal is to calculate the zero temperature single particle retardedGreen’s function (GF) in real spaceGij(ω) = 〈0|ciGˆ(ω)c†j |0〉 =∑α〈0|ci|α〉〈α|c†j |0〉ω − Eα + iη, (2.2)where |0〉 is the vacuum state where there are no phonons or charge carrierin the lattice, Gˆ(ω) = (ω−H+ iη)−1 is the resolvent of H with η → 0+, andEα, |α〉 are single polaron eigenenergies and eigenstates of the Hamiltonian,H|α〉 = Eα|α〉. Note that α refers to a set of quantum numbers that togetheruniquely label the eigenstates of H. Knowledge of this GF allows us tofind the spectrum from the poles, and the local density of states (LDOS)which is measured in scanning tunnelling microscopy (STM) experiments,ρ(i, ω) = − 1pi ImGii(ω).To calculate this quantity, we repeatedly use Dyson’s identity: Gˆ(ω) =Gˆd(ω)+Gˆ(ω)Vˆel−phGˆd(ω), where Gˆd(ω) = (ω−Hd+iη)−1. The first equationof motion generated this way reads asGij(ω) = Gdij(ω) +∑lglF(1)il (ω)Gdlj(ω), (2.3)whereGdij(ω) = 〈0|ciGˆd(ω)c†j |0〉 (2.4)142.2. Momentum average approximation for inhomogeneous systemsare, in principle, known quantities and we have introduced the generalizedGFsF (n)ij (ω) = 〈0|ciGˆ(ω)c†jb†nj |0〉.Note that F (0)ij (ω) = Gij(ω). Next, we generate equations of motion forthese generalized GFs. Using the Dyson identity, for any n ≥ 1, results inF (n)ij (ω) =∑l 6=jglGdlj(ω − nΩj)〈0|ciGˆ(ω)c†l b†l b†nj |0〉+gjGdjj(ω − nΩj)[F (n+1)ij (ω) + nF(n−1)ij (ω)]. (2.5)This equation relates F (n) not only to other GFs of a similar type, but alsointroduces new propagators with phonons at two different sites. Equations ofmotions can be written for these new GFs, linking them to yet more generalGFs, and so on and so forth. The resulting hierarchy of coupled equationsdescribes the problem exactly, but is unmanageable. Approximations areneeded to simplify it and find a closed-form solution.The main idea behind the MA approximations is to simplify this set ofequations by dropping exponentially small contributions in each equation ofmotion. At the simplest level – the so-called MA0 approximation – we ignorethe first term in Eq. (2.5) for any n ≥ 1. This is reasonable at low energies,ω ∼ EGS , where ω−nΩj is well below the energy spectrum of Hd and, there-fore, Gdlj(ω − nΩj) is guaranteed to decrease exponentially with increasingdistance |l− j|. Therefore, we keep only the largest propagator correspond-ing to l = j, and ignore exponentially smaller l 6= j contributions. Althoughthis is the simplest possible such approximation, it is already accurate atlow energies, as shown in the results section. It can also be systematicallyimproved, as discussed below. First, however, we complete this MA0-levelsolution.The simplified equation of motion now reads asF (n)ij (ω) = gjGdjj(ω − nΩj)[F (n+1)ij (ω) + nF(n−1)ij (ω)]. (2.6)On physical grounds, we expect that F (n)ij (ω) must vanish for sufficientlylarge n, because its Fourier transform corresponds to the probability am-plitude that a particle injected in the system will generate n phonons at alater time (or annihilate all n phonons that were initially present), and thisbecomes unlikely as n increases. Therefore, this recursive relation admits asolution asF (n)ij (ω) = An(j, ω)F(n−1)ij (ω), (2.7)152.2. Momentum average approximation for inhomogeneous systemsFigure 2.1: Diagrammatic expansion for (a) the disorder GF Gdij(ω) (boldred line); (b) the polaron GF in a clean system G0ij(ω˜) (double thin line),and (c) the “instantaneous” approximation for the polaron GF in a disor-dered system, Gdij(ω˜) (double bold red line). The thin black lines depict freeelectron propagators, the wriggly lines correspond to phonons, and scatter-ing on the disorder potential is depicted the dashed lines ending with circles.See text for more details.where using Eq. (2.7) in (2.6) gives An(j, ω) as the following continuedfractionAn(j, ω) =ngjGdjj(ω − nΩj)1 − gjGdjj(ω − nΩj)An+1(j, ω). (2.8)This can be efficiently evaluated starting from a sufficiently large cutoff Ncsuch that ANc+1(j, ω) = 0. Generally speaking, Nc must be much larger thanthe average number of phonons expected at site j; in practice, the cutoffis increased until convergence is reached to within the desired accuracy.Substituting F (1)ij (ω) = A1(j, ω)Gij(ω) in Eq. (2.3) leads to a closed system162.2. Momentum average approximation for inhomogeneous systemsof linear equations for the original GF:Gij(ω) = Gdij(ω) +∑lGil(ω)glA1(l, ω)Gdlj(ω). (2.9)This equation has a similar structure to the equation linking the disorderGF, Gdij(ω), to that of the free particle G0ij(ω) = 〈0|ci(ω − H0 + iη)−1c†j |0〉(in the absence of coupling to phonons), which is depicted diagrammaticallyin Fig. 2.1(a), and which reads asGdij(ω) = G0ij(ω) +∑kGdik(ω)ǫkG0kj(ω). (2.10)This analogy shows that coupling to phonons renormalizes the on-site dis-order as ǫl → ǫl + glA1(l, ω). Note that A1(l, ω) depends not only on thelocal phonon frequency Ωl, coupling gl and ǫl, but also on all bare on-siteenergies {ǫk} through the disorder GF Gdll. This is the simplest example ofthe emergence of a renormalized potential for this problem, that is madevery transparent within the MA approximation.While Eq. (2.9) can be solved directly for a finite-size system, we canimprove its efficiency and reveal a different physical interpretation by explic-itly removing the “average” contribution due to the electron-phonon inter-actions. Let g and Ω be the average values of the electron-phonon couplingand phonon energies. We assume that the on-site energy average is zero,ǫ¯i = 0, (a finite value results in a trivial shift of all energies). Then, letAn(ω) =ngg0(ω − nΩ)1 − gg0(ω − nΩ)An+1(ω)(2.11)be the continued fractions corresponding to these average parameters, wherewe use the short-hand notationg0(ω) = G0ii(ω) =1N∑k1ω − ε(k) + iηfor diagonal elements of the free propagator in the absence of disorder. Itis given by the momentum average of the free propagator over the Brillouinzone, where ε(k) is the free particle energy dispersion: ε(k) = −2t∑3i=1 cos ki.The “average” renormalization of the on-site energy is now recognized torepresent the corresponding MA0 self-energy for a “clean” system, i.e. ahomogeneous system with average coupling and phonon frequency,ΣMA0(ω) = gA1(ω),172.2. Momentum average approximation for inhomogeneous systemssee, for instance, Eqs. (11) and (12) of Ref. [48].Introducing the effective disorder potentialv0(l, ω) = glA1(l, ω) − ΣMA0(ω), (2.12)Eq. (2.9) can be rewritten asGij(ω) = Gdij(ω˜) +∑lGil(ω)v0(l, ω)Gdlj(ω˜), (2.13)where ω˜ = ω − ΣMA0(ω). This energy renormalization, ω → ω˜, reflects thefact that processes describing the formation of polaron in the “clean” systemhave been explicitly summed.Note that solving Eq. (2.13) is numerically more efficient compared toEq. (2.9), especially when disorder is confined to a finite region in thesystem, away from which v0(l, ω) vanishes rapidly. As explained below, Eq.(2.12) reveals a different interpretation of the interplay between disorder andelectron-phonon coupling.Consider first the meaning of Gdij(ω˜), which would be the solution if wecould ignore v0(j, ω). In the absence of on-site disorder this term equalsG0ij(ω˜), that is the expected solution for a polaron in the clean system,depicted diagrammatically in Fig. 2.1(b) (of course, here the exact self-energy is approximated to ΣMA0(ω)). Comparing Figs. 2.1(a) and 2.1(b),it follows that Gdij(ω˜) is the sum of the diagrams shown in Fig. 2.1(c).At first sight, this seems to be the full answer for this problem, since thesediagrams sum the contributions of all the processes in which the polaronscatters once, twice, etc., on the disorder potential. This is certainly theanswer obtained in the limit of an “instantaneous” approximation valid whenΩ → ∞, i.e. when the ions are very light and respond instantaneouslyto the motion of electrons. In this case, one can perform a Lang-Firsovtransformation and after an additional averaging over phonons, one arrivesat the following approximative effective Hamiltonian [46, 49]Hinst = −t∗∑〈i,j〉(c†i cj + H.c.) +∑i(ǫi −g2Ω )c†i ci, (2.14)where t∗ = t exp (−g2/Ω2) is the renormalized polaron hopping, and −g2/Ωis the polaron formation energy (for simplicity, here we assume that thephonon energies and electron-phonon coupling are homogeneous, gi → g,Ωi →Ω). The GF of this Hamiltonian is also given by Fig. 2.1(c), if the polaron182.2. Momentum average approximation for inhomogeneous systemspropagator is approximated byG0ij(ω˜) →1N∑keik·(Ri−Rj)ω − ε∗(k) + g2Ω + iη,where ε∗(k) uses the renormalized kinetic energy. Of course, using the fullexpression for G0ij(ω˜) is preferable since the self-energy ΣMA0(ω) describesmuch more accurately the overall energy shift and effective mass renormal-ization than those asymptotic expressions, and it also includes the quasipar-ticle weight.That Gdij(ω˜) cannot be the full answer becomes obvious when we rewritethe clean polaron propagators in terms of free particle and phonon lines, i.e.,we substitute the expansion corresponding to Fig. 2.1(b) in 2.1(c). Doingso reveals that, within this “instantaneous” approximation, scattering ofthe electron on the disorder potential is allowed only when no phononsare present, see zoom-in in Fig. 2.1(c). However, for moderate and largeelectron-phonon coupling, one expects the probability to find no phononsin the system to be exponentially small, therefore the processes summed inFig. 2.1(c) have very low probabilities.Figure 2.2: (color online) (a) Diagrammatic expansion for the full MA so-lution Gij(ω) (thick dashed blue line) in terms of the clean polaron Green’sfunction (double thin line) and scattering on the renormalized disorder po-tential ǫ∗l (ω), depicted by dashed lines ended with squares. (b) Diagram-matic expansion of ǫ∗l (ω). For more details, see text.What is missing in Fig. 2.1(c) are diagrams describing the scattering ofthe electron on the disorder potential in the presence of the phonons from192.2. Momentum average approximation for inhomogeneous systemsthe polaron cloud. Their contribution is included through the renormalizedpotential v0(l, ω) in the second term of Eq. (2.13). Indeed, the full MAsolution shows that the polaron scatters not on the bare disorder ǫl, but onthe renormalized disorder potentialǫ∗l (ω) = ǫl + v0(l, ω), (2.15)as depicted in Fig. 2.2(a). The diagrammatic expansion of the additionalterm v0(l, ω), shown in Fig. 2.2(b), verifies that it indeed describes theeffective scattering in the presence of various number of phonons.Taken together, the diagrams summed in Fig. 2.2 represent all possi-ble contributions to the polaron propagator in the disordered system. TheMA0 approximation consists of discarding exponentially small contributionsfrom each of these diagram, as already discussed. MA0 also has an exactvariational meaning, namely of assuming that the polaron cloud can havephonons only on a single site, in direct analogy with the clean case [50, 51].It is quite remarkable that all diagrams corresponding to this variationalapproximation can still be summed analytically in closed form, even in thepresence of disorder.As it is the case for the clean system, MA can be systematically im-proved by keeping more terms in Eq. (2.5). In particular, at the MA1level, we treat the equation for the F (1) functions exactly, and discard ex-ponentially small off-diagonal propagators only for n ≥ 2. The logic hereis that the propagators appearing in F (1) have the highest energy, thereforethe slowest exponential decay. The MA1 equations can also be solved inclosed form. The details are presented in the Appendix B. The final solu-tion looks identical to Eq. (2.13), except the renormalized energy is nowω˜ = ω − ΣMA1(ω) while the renormalized potential v0(l, ω) is replaced bya more complicated, yet more accurate expression v1(l, ω). The meaning ofall these quantities, however, is the same.To summarize, MA reveals that the role of electron-phonon coupling istwo-fold. On one hand, it renormalizes the quasiparticle properties due topolaron formation, just like in the clean system (as revealed by the explicitappearance of the “clean” system self-energy). However, this coupling alsorenormalizes the disorder potential experienced by the particle, ǫl → ǫ∗l (ω).As we show later for various types of disorder, this renormalization is non-trivial in that it has strong retardation effects, and has significant conse-quences. “Instantaneous” approximations completely ignore this renormal-ization, and therefore miss essential physics. In contrast, MA gives a closedanalytical expression for all quantities of interest, in a formulation whose202.3. Polaron near a single impuritymeaning is very transparent and whose accuracy can be systematically im-proved.2.3 Polaron near a single impurityWe now apply the general formalism described above to the simplest type of“disorder”, namely an otherwise clean three-dimensional simple cubic latticewith a single impurity. Impurities which modulate the on-site energy, thestrength of the electron-phonon coupling and/or the phonon frequency areseparately considered. We investigate under what conditions such impuritiescan trap the polaron.As a reference case, we first review briefly the solution in the absenceof electron-phonon coupling. In this case, the impurity can only modulatethe on-site potential, ǫi = −Uδi0, and the Hamiltonian reduces to the Wolff-Clogston model [52, 53]:Hd = −t∑<i,j>(c†i cj + H.c.) − Uc†0c0 = H0 + Vˆd. (2.16)For this form of Vˆd, Eq. (2.10) reads asGdij(ω) = G0ij(ω) − UGdi0(ω)G00j(ω), (2.17)which is trivially solved to giveGdij(ω) = G0ij(ω) − UG0i0(ω)G00j(ω)1 + UG000(ω). (2.18)Because of translational and time reversal symmetry, G0ij(ω) = G0i−j,0(ω) =G0j−i,0(ω), etc. The LDOS at the impurity site i = 0 is then found to beρd(0, ω) = −1π ImGd00(ω) =ρ0(0, ω)|1 + UG000(ω)|2, (2.19)where ρ0(0, ω) = − 1pi ImG000(ω) is the LDOS in the clean system (same astotal DOS, because of the translational invariance). As a result, a boundstate below the continuum, signalled by a delta-function peak in ρ(0, ω),occurs if and only if the denominator of Eq. (2.19) vanishes. For a 3Dsimple cubic lattice with an attractive impurity, this means that an impuritybound state appears if there is an energy E < −6t such that ReG000(E) =−1/U (below the continuum, the imaginary part of G000(E) vanishes). This212.3. Polaron near a single impurityequation can be solved graphically to find that a bound state appears forany U ≥ Uc = − 1ReG000(−6t) ∼ 3.96t.In the presence of electron-phonon coupling the equations are more com-plicated, but the idea is the same: we calculate the LDOS at the impuritysite and compare it to the DOS of the clean system. If the former has apeak below the threshold of the latter, then a bound state exists at thatenergy. We then vary U to find the threshold value above which a boundstate is guaranteed. More details about the impurity state, such as its local-ization length, statistics for the phonon cloud, etc., can be extracted fromthe LDOS at sites in the neighbourhood of the impurity. Here we focus onidentifying when bound impurity states are stable.2.3.1 Impurity changing the on-site energyThe Hamiltonian describing this case isH = Tˆ + Ω∑ib†ibi + g∑ic†i ci(bi + b†i ) − Uc†0c0, (2.20)where Tˆ is the particle’s tight-binding kinetic energy. We are interested inattractive impurities, U > 0, when an impurity state can be bound near theimpurity site. To find the LDOS at the impurity site, we need to solve Eq.(2.13) to find G00(ω). Note that now Gdij(ω) is known, being given by Eq.(2.18). The free-particle propagators G0ij(ω) = 1N∑keik·(Ri−Rj)ω−ε(k)+iη , where forthe simple cubic lattice ε(k) = −2t∑3i=1 cos ki, can be calculated by doingthe integrals over the Brillouin zone. A more efficient approach, which weuse, is discussed in Ref. [54], and also in chapters 4 and 5.Since the disorder GF approaches to the bulk GF away from the impurity,Gdll(ω) → G0ll(ω) when |l| → ∞, the renormalized impurity potential decaysfast at these sites as A1(l, ω) → A1(ω). Physically, this is because theimpurity has less and less influence at sites far from where it is located.Mathematically, this follows from Eq. (2.18) and the fact that G0l0(ω − nΩ)decreases exponentially with the distance between site l and the origin, atenergies below the free particle continuum which are of interest here. As aresult, in Eq. (2.13) we only need to sum over sites l close to the impurity,and this small system of equations can be solved very efficiently.The appropriate value for this cutoff varies depending on the variousparameters, but generically it decreases as the energies of interest becomelower. Physically, this can be understood as follows. First, let us explainwhy is the renormalized disorder potential non-local, even though the bare222.3. Polaron near a single impurityimpurity potential is local. The answer is provided by the diagrams whichcontribute to it, see Fig. 2.2(b). Consider, for simplicity, the MA0 approxi-mation. In this case, all phonon lines appearing in these diagrams start andend at the site l for which v0(l, ω) is being calculated – this is the site wherethe phonon cloud is located. However, the electron is found with variousprobabilities away from the polaron cloud, so it can scatter on the impurityif this is located within the “radius” of the polaron, where the electron ismost likely to be found. In other words, the range of the renormalized po-tential is controlled by the polaron size. From medium to large couplings,as the polaron becomes smaller, the renormalized potential becomes morelocal. At small couplings, though, the distance between the electron andits phonon cloud can be appreciable, and the range of the effective disorderpotential increases. All this is clearly manifest in the expansion of v0(l, ω)to the lowest order in electron-phonon coupling g,v0(l, ω) ≈ g2[Gdl0(ω−Ω)−G0l0(ω −Ω)] =−Ug21 + Ug0(ω)G0l0(ω−Ω)G00l(ω−Ω),which shows direct proportionality to the ability of particle to visit theimpurity in the presence of a phonon, G0l0(ω − Ω).In Fig. 2.3(a) we plot ρ(0, ω) over a wider energy range, for severalvalues of U , within MA0 approximation. The dashed line shows the DOSof the clean system, multiplied by 100 for visibility. For U/t = 1.8, 1.9, theimpurity attraction is not sufficient to bind a state below the continuum,although the LDOS is pushed towards the lower band-edge. For U = 1.95t,there is a peak just below the continuum. Because of the finite value ofη, the two features are not completely separated and the continuum onsetlooks like a “shoulder”, however lowering η allows us to clearly separate thetwo features. Finally, for U = 2t the bound state peak is clearly belowthe continuum, so in this case threshold value of impurity attraction for abound state to appear is Uc ≈ 1.95t. This Uc value coincides with the DMCresult, although our polaron energies are slightly higher than the exact DMCvalues shown in Fig. 2.3(e), as expected for a variational approximation.Increasing the variational space by going to MA1 level significantly improvesthe agreement with the DMC since all features move toward lower energies,see Fig. 2.3(c). The critical value of Uc ≈ 1.95t for which an impurity stateappears below the continuum is only slightly affected, in spite of the overallshift of the spectral weights.The dependence of the bound state energy on the cutoff is shown in Figs.2.3(b) and (d) for MA0 and MA1, respectively, when U = 2t, Ω = 2t andthe effective coupling λ = g26tΩ = 0.8. At these energies, keeping only the232.3. Polaron near a single impurity-7.53 -7.52ω/t0102030ρ(0,ω)lc = 0lc = 1lc = 2-7.5 -7.4 -7.3ω/t0246ρ(0,ω)U = 2.00tU = 1.95tU = 1.90tU = 1.80tU = 0, x100(b)(a)-7.554 -7.548 -7.542ω/t0102030ρ(0,ω)lc = 0lc = 1lc = 20246ρ(0,ω)U = 2.00tU = 1.95tU = 1.90tU = 1.80tU = 0, x100-7.6 -7.5 -7.4 -7.3ω/t0246ρ(0,ω)(c) (d)(e)Figure 2.3: LDOS at the impurity site ρ(0, ω) in units of 1/t, vs. the energyω/t, for (a) MA0 with lc = 0 for U = 1.8, 1.9, 1.95 and 2.0. The dashed lineshows the DOS for the clean system, times 100; (b) Bound polaron peakin MA0 at U/t = 2, and cutoffs in the renormalized potential lc = 0, 1, 2.Panels (c) and (d) are the same as (a) and (b), respectively, but using MA1.Panel (e) shows the DMC results from Ref. [4], for the same parameters as(a) and (b). Other parameters are Ω = 2t, λ = 0.8, η/t = 10−3. 242.3. Polaron near a single impurity0 1 2 3 4Uc/t0123456λΩ = .5t (MA1)Ω = 2t (MA0)Ω = 4t (MA0)Ω = 8t (MA0)Ω = .5t (DMC)Ω = 2t (DMC)Ω = 4t (DMC)Ω = 8t (DMC)Figure 2.4: Phase diagram separating the regime where the polaron is mobile(below the line) and trapped (above the line). The effective coupling isλ = g2/(6tΩ), and the threshold trapping potential Uc is shown for severalvalues of Ω/t. The MA results (filled symbols) compare well with the DMCresults of Ref. [4] (empty symbols).local part in MA0, i.e. setting v0(l, ω) → δl,0v0(0, ω), is already a very goodapproximation. Including the correction from the 6 nearest neighbor sites(lc = 1) lowers the energy somewhat, but the contribution from the secondnearest neighbor sites (lc = 2) is no longer visible on this scale in either theMA0 or MA1 level.By repeating this process for different values of the parameters we cantrace Uc in the parameter space. This is shown in Fig. 2.4, for Ω/t =0.5, 2, 4, 8 and various effective couplings. The MA results (filled symbols)are in good quantitative agreement with the DMC results (empty symbols)for larger values of phonon frequency, Ω ≥ t. Shown here are MA0 resultsfor cutoff lc = 0. Using MA1 and/or increasing the cutoff changes the valuesof Uc by less than 1% everywhere we checked. For the smaller frequenciessuch as Ω = 0.5t, MA is known to become quantitatively less accurate252.3. Polaron near a single impurityat intermediary couplings [14, 50], and indeed, here we see a discrepancywith the DMC data even for the MA1 results. To improve the quantitativeagreement here, one should use a 2- or 3-site MA variational approximationfor the phonon cloud, as discussed in Ref. [55].As expected, Uc approaches the threshold value for a bare particle forλ→ 0 which is roughly 3.96t. Uc decreases for stronger effective couplings,but the lines never intersect the λ-axis: Uc = 0 is impossible, since thepolaron cannot be trapped (localized) in a clean system, although it canhave an extremely heavy effective mass when λ is strong. This is also thereason why Uc decreases with increasing λ: a heavier quasiparticle is easierto trap by the impurity [4]. However, we claim that this is not the fullstory, and that the renormalization of the trapping potential also plays anon-trivial role, as detailed below.0 2 4 6 8 10U/t1.01.21.41.61.82.0U*/ Uλ = 0.5λ = 1.5Figure 2.5: (color online) Effective value of the impurity attraction U∗/U ,extracted from the scaling E∗B/t∗ = f(U∗/t∗), for λ = 0.5, 1.5 and Ω = 3t.In order to have the best fit to the data, we plot each curve starting fromslightly larger U/t than the corresponding Uc/t given in Fig. (2.4). Formore details, see text.Consider Hamiltonian (2.16), which describes the impurity problem in262.3. Polaron near a single impuritythe absence of electron-phonon coupling. The binding energy of the impu-rity state (once formed) is a monotonic function of the only dimensionlessparameter of this problem, U/t: EBt = f(Ut)for any Ut ≥ Uct ≈ 3.96. In thepresence of electron-phonon coupling, if one views the polaron as a quasi-particle with an effective hopping t∗ which scatters on the same impurity Uas the bare particle (instantaneous approximation), then the polaron bind-ing energy should be E∗Bt∗ = f(Ut∗)for any Ut∗ ≥ Uct∗ ≈ 3.96. In particular,this predicts Uc/t = 3.96t∗/t to decrease with increasing λ, in qualitativeagreement with Fig. 2.4. This hypothesis can be tested. The function f(x)is easy to calculate numerically; we can also extract the binding energy E∗bfor the trapped states by comparing their energy to the ground state energyof the polaron in the clean system, and the effective hopping t∗ ∼ 1/m∗ isdirectly linked to the effective polaron mass m∗ in the clean system [14].What we find is that such a scaling is not obeyed. Instead, one needsto also rescale the impurity potential, i.e. use E∗Bt∗ = f(U∗t∗)where U∗ 6= U .Of course, this scaling assumes that the scattering potential is local. As wediscussed, this is not true although it is a good approximation for mediumand large couplings.In Fig. 2.5 we show the renormalized value U∗/U extracted this way,as a function of U/t above the corresponding threshold values Uc/t, for amedium and a large effective coupling λ = 0.5, 1.5 and Ω = 3t. Qualitativelysimilar curves are found for other parameters. We see that U∗ → U onlywhen U →∞ and becomes the dominant energy scale (hence anything elseis a small perturbation). For fixed values of U and Ω, U∗ is found to increasewith increasing coupling λ – this is also expected, since the renormalizationis directly caused by the electron-phonon coupling, see Fig. 2.2(b). Thisrenormalization is a direct illustration of the general result of Eq. (2.15):the electron-phonon coupling changes not only the properties of the polaron(its effective mass), but also the effective disorder potential it experiences.However, it is very wrong to expect that the potential renormalizationcan always be described by a simple rescaling by some overall value. Indeed,Eq. (2.15) shows that the renormalized potential is expected to be a func-tion of energy, because of retardation effects. This function is not roughlyconstant, instead it has significant and very non-trivial dependence on ω,as illustrated by plotting v0(0, ω) over a wide energy range, in Fig. 2.6.Similar curves are found for other values of the parameters. Over a narrowrange of energies around ω ∼ −7.5t where the bound state forms for theseparameters (see Fig. 2.3), v0(0, ω) varies slowly and can be approximatedas an overall negative constant. This explains why here we can approximate272.3. Polaron near a single impurity-7 -6 -5 -4 -3 -2ω/t-80-60-40-20020406080Re[v0(0,ω)]η = 0.05η = 0.01Figure 2.6: (color online) Real part of the additional on-site MA0 disorderpotential v0(0, ω) over a wide energy range, for U = 2t,Ω = 2t, λ = 0.8 andtwo values of η.ǫ∗0(ω) = −U + v0(0, ω) ≈ −U∗, with U∗ > U , as discussed above. At higherenergies, however, v0(0, ω) goes through singularities and changes sign fromnegative to positive and back. Although at first sight these singularities aresurprising, they should be expected based on Eq. (2.12). The self-energyof the Holstein polaron is known to have such singularities, especially atmedium and higher couplings where an additional second bound state formsand the continuum above shows strong resonances spaced by Ω [14]. Inparticular, as λ → ∞ and the spectrum evolves towards the discrete Lang-Firsov limit En = − g2Ω + nΩ, the self-energy has a singularity at the top ofeach corresponding band. The renormalized potential of Eq. (2.12) is thedifference between two such curves, displaced from each other. It is thereforenot surprising that it has such nontrivial behavior.Physically, such strong retardation effects are not surprising either, sincethe additional potential v0(i, ω) describes the scattering of the electron inthe presence of the phonon cloud. The structure of the phonon cloud varieswith energy; for instance, one expects quite different clouds at low energies282.3. Polaron near a single impuritywithin the polaron band vs at higher energies, where there is a continuumof incoherent states with finite lifetimes. This suggests that the diagramsof Fig. 2.2(b) that contribute most to the series change with energy, andso does the total result. As a final note, we mention that at these higherenergies, MA1 should be used. It is well-known that MA0 fails to describeproperly the location of the polaron+one phonon continuum, since it doesnot include the corresponding variational states [51]. This problem is fixedat the MA1 and higher levels [50].To summarize, the MA approximation is found to agree well with re-sults from DMC in describing the trapping of Holstein polaron near a singleattractive impurity. Although it is not quantitatively as accurate as DMC,besides efficiency its main advantage is that the analytic equations that de-scribe MA allow us to understand the relevant physics. In particular, weshowed that coupling to bosons renormalizes the disorder in a very non-trivial way.2.3.2 Impurity changing the electron-phonon couplingWe now assume that the impurity does not change the on-site energy, butinstead it modifies the value of electron-phonon coupling at its siteH = Tˆ + Ω∑jb†jbj +∑jgjc†jcj(bj + b†j), (2.21)where gj = g+(gd−g)δj,0. Thus, gd and g are the electron-phonon couplingsat the impurity site and in the bulk of the system, respectively. Since Vˆd = 0in this case, the non-interacting part of the Hamiltonian is Hd = H0. Thus,Gdjj(ω − nΩ) → G0jj(ω − nΩ) ≡ g0(ω − nΩ) in the continued fractions, Eq.(2.8), whose dependance on the site index j is now through the coupling gjonly. As a result, the effective disorder potential v0(j, ω) vanishes everywhereexcept at the impurity site, j = 0:v0(j, ω) ≡ ∆(ω)δj,0, (2.22)where ∆(ω) = gdA1(0, ω) − gA1(ω), and A1(0, ω) is like in Eq. (2.11) butwith g → gd. This shows that even though ǫi = 0 in this case, the inho-mogeneity gives rise to an effective potential ǫ∗i (ω) = δi,0∆(ω). This is nowlocal because only when the particle is at the impurity site it can experiencea different coupling. Equation (2.13) can now be solved analytically to giveρ(0, ω) = − 1π Im( g0(ω˜)1 − ∆(ω)g0(ω˜)). (2.23)292.3. Polaron near a single impurity0 1 2 3 4 5 6(gd- g)/t024 λΩ = 3t, MAΩ = 6t, MAΩ = 8t, MAΩ = 3t, IAΩ = 6t, IAΩ = 8t, IAFigure 2.7: (color online) Phase diagram separating the regime where thepolaron is delocalized (below the line) and trapped (above the line), asa function of the difference between the impurity and the bulk electron-phonon coupling, gd− g. Symbols show MA0 results, while the dashed linescorrespond to the instantaneous approximation.We can now find the values of gd− g for an impurity state to emerge belowthe continuum, for given values of g and Ω. The results are shown in Fig. 2.7for gd > g, when the polaron formation energy at the impurity site, −g2d/Ω, islower than the bulk value −g2/Ω, and a bound state may be expected to formeven within the instantaneous approximation. Symbols show MA0 results,while the dashed lines are for the instantaneous approximation (DMC resultsare not available for this case). The two agree quantitatively only in the limitλ → 0. This proves that the additional renormalization included in MA issignificant for this type of impurity, as well. We note that all critical linesintercept the x-axis at a finite value, i.e. for any value of Ω and g = 0,there is a critical finite value gd above which an impurity state forms. Forexample, for Ω = 3t this critical value is gd ≈ 3.9t and it increases withgrowing Ω, as expected.302.3. Polaron near a single impurityUnlike for an impurity which changes the on-site potential, and whichcan bind at most one impurity state, impurities which change the electron-phonon coupling can bind multiple impurity levels. As gd increases and theenergy of the impurity level moves towards lower energies, additional boundstates, spaced by roughly Ω, emerge whenever the distance between the lastone and the bulk polaron band is of order Ω. The origin of this sequence ofbound states is straightforward to understand in the limit gd ≫ g, t, wherethe dominant term in the Hamiltonian isH ≈ gdc†0c0(b†0 + b0) + Ωb†0b0,with c†0c0 ≈ 1 as the weight of the bound state is concentrated at the im-purity site. This Hamiltonian can be exactly diagonalized using the Lang-Firsov transformation [15]. Using new operators B0 = b0 + α, where α is areal number, this Hamiltonian can be written asH = ΩB†0B0 + (gd − αΩ)(B†0 +B0) + α2Ω − 2αgd.Choosing α = gd/Ω diagonalizes this Hamiltonian, resulting in a series ofequally spaced eigenenergies En = nΩ− g2d/Ω. For finite gd, all states whichlie below the bulk polaron continuum become impurity states, and basicallydescribe excited bound states with additional phonons at the impurity site.So far we have considered gd > g, where a ground-state impurity levelcan emerge. It is important to note that discrete peaks can also appear forgd < g, although not at low energies. This happens when λ is sufficientlylarge that there is a gap between the bulk polaron band and higher featuresin the spectrum, such as the polaron+one phonon continuum, or the bandassociated with the second bound state, once it forms [56]. A typical exampleis shown in Fig. 2.8, where in the presence of an impurity with a weakercoupling gd < g (full line), a discrete state appears above the polaron band.Since most of its weight is removed from the bulk polaron band (dashedline), we interpret this as being an “anti-bound” polaron state. A similarstate is also expected to appear for a repulsive on-site U < 0 potential.2.3.3 Isotope impurityThe last case we consider is an isotope impurity. Because of the differentmass of an isotope substitution, Md 6= M , its phonon frequency Ωd ∼M−1/2das well as electron-phonon coupling gd ∝ 1/√MdΩd ∼ M−1/4d are different.Interestingly, both the effective coupling λd = g2d/(6tΩd) = g2/(6tΩ) = λand the polaron formation energy −g2d/Ωd = −g2/Ω show no isotope effect312.3. Polaron near a single impurity-7.6 -7.4 -7.2ω/t00.10.2ρ(0,ω)gd - g = -0.028tgd = gFigure 2.8: For large λ, in the clean system (dashed red line) the firstpolaron band is separated by an energy gap from the next features in thespectrum (here, the band associated with the second bound state). Forgd < g, an “anti-bound” impurity state is pushed inside this gap (black fullline). Parameters are Ω = t = 1, λ = 1.2 and η = 10−3.[57]. As a result, within the instantaneous approximation of Eq. (2.14) onewould predict that the isotope is “invisible” and the polaron spectrum isbasically unaffected by its presence.We consider a single isotope impurity located at the originH = Tˆ +∑jgjc†jcj(bj + b†j) +∑jΩjb†jbj , (2.24)where Ωj = Ω + (Ωd − Ω)δj,0 and gj = g + (gd − g)δj,0 are chosen such thatλd = λ. Just like in the previous section, because there is no on-site disorderǫi = 0, Gdjj(ω) = G0jj(ω) and the effective disorder potential is again local,i.e. it vanishes everywhere other than the impurity site. As a result, theLDOS at the impurity site is given by Eq. (2.23), except that here∆(ω) = Σd(ω) − ΣMA0(ω), (2.25)322.3. Polaron near a single impuritywhere Σd(ω) has the same functional form as ΣMA0(ω), but with g → gdand Ω → Ωd.By investigating this LDOS for different values of parameters we find thatthere exists a threshold value of the effective coupling, λ∗, below which low-energy bound states do not form irrespective of how small or large Md/Mis. In other words, for λ < λ∗ the behavior agrees with the prediction of theinstantaneous approximation.We can estimate λ∗ as follows. Consider the case of a very light isotope,so that Ωd, gd ≫ Ω, g, t. In this limit, Σd(ω) → −g2d/Ωd = −6tλ. The boundstate appears when the LDOS is singular as its denominator vanishes1 − ∆(ω)g0(ω˜) = 0 → ΣMA0(ω) + 1/g0(ω˜) = −6tλ, (2.26)where Eq. (2.25) and Σd(ω) ≈ −6tλ are used.0 2 4 6 8Ω/t12λ∗λ*l = 0.66Figure 2.9: (color online) Critical effective coupling λ∗ above which an im-purity state may appear for a sufficiently light isotope. Below this line,polarons cannot be bound near isotopes. Symbols shows MA0 results. Thedashed line is the analytic low-bound for λ∗ discussed in the text.Consider now the limiting case when a bound state emerges just belowthe bulk polaron ground state, i.e. Eq. (2.26) has a solution at ω ≤ EGS .In the clean system, the polaron ground state energy EGS is the lowest pole332.3. Polaron near a single impurityof G(k = 0, ω) = [ω − ε(k = 0) − ΣMA0(ω)]−1, so it satisfies:EGS = −6t + ΣMA0(EGS).Using this in Eq. (2.26) suggests that a solution can exist if λ > λ∗, whereλ∗ =∣∣∣∣EGS6t∣∣∣∣− 1 − 16tg0(−6t), (2.27)in which g0(−6t) ≈ −1/3.96t for 3D tight-binding model. Since EGS → −6tas λ→ 0, we find that λ∗ → 0.66 in this limit, and that it increases as EGSbecomes more negative, for example with increasing λ. These considerationsare confirmed by the data shown in Fig. 2.9. Here, the symbols show valuesof λ∗ found numerically with MA0, and the dashed line is the lower boundof 0.66, discussed above.For λ > λ∗, bound impurity states can appear near isotopes if gd andΩd = g2d/(6tλ) are sufficiently large. In Fig. 2.10 we show critical lines fortwo cases, Ω = 4t, 8t. The symbols show the MA0 results, which convergetowards their corresponding λ∗ values as Ωd → ∞, as expected. Of course,the largest values considered for Ωd are unphysical; we use them only toillustrate the convergence towards λ∗.The existence of a region of the parameter space where bound polaronstates appear near an isotope is in direct contradiction of the instantaneousapproximation, and again illustrates the importance of the renormalizeddisorder potential ∆(ω) which makes their trapping possible. In this context,it is worth mentioning that there is clear evidence for electronic states boundnear isotope O16 defects in CuO2 planes [58], although the precise nature ofthese states has not been clarified and the measurements are certainly not inthe extremely underdoped regime where our single-polaron results are valid.Interestingly, when such bound states form near an isotope, the spectrumis different than that for the other two types of impurities. As shown inFig. 2.11, bound states now appear simultaneously both below and abovethe bulk polaron continuum, not just below it. This provides a possible“fingerprint” for polarons trapped near isotopes. Finally, we note that evenwhen no low-energy impurity state is observed, it is again possible to havehigher energy bound states inside the gaps opening between various featuresin the bulk polaron spectrum.To summarize, in the presence of isotope defects, polarons in the weaklycoupled regime λ < λ∗ always remain delocalized. Only for λ > λ∗ it is pos-sible to trap polarons near an isotope. This makes the isotope substitutionquite distinct from the other two cases, where bound states exist for any λif the impurity is strong enough.342.4. Summary and conclusions1 10 100 1000(Ωd - Ω)/t1234λΩ = 4tΩ = 8tFigure 2.10: Phase diagram separating the regime where the polaron isdelocalized (below the line) and trapped (above the line), as a function ofthe difference between the Ωd − Ω, on a logarithmic scale. Symbols showMA0 results for Ω = 4t, 8t. As Ωd → ∞, these critical lines converge towardstheir corresponding λ∗ (dashed lines), below which polaron states are alwaysdelocalized.2.4 Summary and conclusionsWe studied the threshold for the emergence of polaron bound states near var-ious types of single impurities, using the momentum-average approximation.Electron-phonon coupling was shown to strongly renormalize the impuritypotential in a nontrivial way which includes strong retardation effects. Thisis a feature that is completely absent in the instantaneous approximation,which is the only other available “simple” description of this problem.We considered the simplest models of impurities that change the strengthof the on-site energy, the local electron-phonon coupling, or are isotopesubstitutions that modify both the coupling and the phonon energy. Wecalculated the polaron binding phase diagrams for each case. The first casehad been considered previously by numerical methods [4, 46], and our resultsare in good quantitative agreement with their predictions. To our knowledge,the other two cases have not been investigated before. We showed that in the352.4. Summary and conclusions-15.6 -15.3 -15ω/t0123ρ(0, ω)Ωd - Ω = 7tΩd - Ω = 0, x10Figure 2.11: LDOS at the impurity site near an isotope with Ωd = Ω + 7t(full line). Two discrete states, one above and one below the bulk polaronband, are seen. For comparison, the DOS in the clean system (multiplied by10) is shown as a dashed line. Parameters are Ω = 4t, λ = 2.5 and η = 10−2.first two cases bound states always exist for a sufficiently strong impurity,however the polaron remains delocalized for the case of isotope substitutionof arbitrary strength if the effective coupling is weaker than a thresholdvalue, λ∗. Differences in the LDOS at the impurity site have also beenfound, such as the possibility to bind multiple states near an impurity thatchanges the coupling, or the unusual fingerprint of discrete states both belowand above the bulk polaron continuum, for an isotope bound state.Of course, a realistic description of an impurity in a real system maycombine several of these inhomogeneities, and even the form of the electron-phonon coupling could be affected. MA gives an efficient yet quite accurateway to deal with such cases, and can be easily generalized to other types ofcouplings where MA has been used succesfully to describe bulk properties.Whereas we expect the single impurity results to remain valid for a sys-tem with multiple impurities if the mean free path is long and the polaroninteracts with one impurity at a time, in the presence of significant dis-order, when multiple scattering processes become important, the polaronscan undergo Anderson localization. This limit has been addressed withina generalized DMFT approach [42], however we believe that our simpler362.4. Summary and conclusionsformulation might provided additional insight and uncover previously unex-plored aspects of Anderson localization for polarons. The opposite limit ofweak (and extended) disorder is discussed in the next chapter.37Chapter 3A perturbational study ofthe lifetime of a Holsteinpolaron in the presence ofweak disorder3.1 IntroductionStudying the behavior of solid state systems under the simultaneous actionof disorder and interactions is a significant challenge in condensed matterphysics. Strong correlations in interacting systems often give rise to sharpquasiparticles. Scattering of such quasiparticles from weak disorder shouldjust limit their lifetime. For strong disorder and no interactions, it is wellunderstood [59] that constructive interference of the backscattered waves canlocalize single particles such that they lose their itinerancy. If interactionsare turned on, there is no consensus about the effect of disorder on thequasiparticles of interacting systems.In chapter 1 I studied the bound state of a Holstein polaron around sin-gle impurities of arbitrary strength. For the case where disorder is extendedall over the lattice, if one views the polaron as a particle with a renormal-ized mass, strong disorder should result in Anderson localization. However,phonon-assisted hopping of carriers between localized states is well estab-lished as a conduction process in lightly doped semiconductors [60]. Thissuggests that, for suitable electron-phonon couplings, polarons may still beitinerant in a disorder potential that would localize particles with the sameeffective mass. A possible explanation for this was offered in the first chap-ter, where the momentum average (MA) approximation was used to showthat the electron-phonon coupling renormalizes the disorder potential in astrongly energy-dependent manner, so that the effective disorder seen bypolarons can be drastically different from the bare disorder.A possible approach for studying this problem, valid for weak disorder,383.2. The model and its solutionis to use perturbation theory and perform the disorder average analytically,similarly to the Born approximation widely employed for charge carriers inthe absence of electron-phonon coupling [61]. Because localization cannotbe described within a perturbational calculation, the polaron eigenstatesremain extended and self-averaging over all disorder realizations is appro-priate.In this chapter, I follow this approach and calculate, using MA andfor weak disorder, the disorder-averaged Green’s function of the Holsteinpolaron and the resulting polaron lifetime and energy shift. Because here wefocus only on the lowest-energy polaron states, which are already accuratelydescribed at the MA0 level, in the following we restrict ourselves to thisflavor and call it MA for simplicity.The chapter is organized as follows. In Sec. 3.2, I present the gen-eralization of MA to include disorder perturbationally. Sec. 3.3 containsthe results and their analysis, and Sec. 3.4 is for conclusion. Some of thecomputational details are given in appendixes.3.2 The model and its solutionWe start again with the Hamiltonian for a Holstein polaron in a lattice withrandom on-site energiesH = Hd + Vˆel−ph = H0 + Vˆd + Vˆel−ph, (3.1)where the non-interacting part of the Hamiltonian, Hd, is divided into Vˆd =∑i ǫic†i ci, describing the on-site disorder potential experienced by the chargecarrier, andH0 = −t∑〈i,j〉(c†i cj + H.c.) + Ω∑ib†ibi (3.2)describing the kinetic energy of charge carrier plus the (optical) phononenergies (~ = 1). The interaction termVˆel−ph = g∑ic†i ci(b†i + bi) (3.3)describes the Holstein coupling between a single charge carrier and phonons.ci and bi are annihilation operators for the charge carrier and phonons,respectively, in a simple cubic lattice whose sites are indexed by i. The spinof the carrier is irrelevant and is therefore ignored.393.2. The model and its solutionThe on-site energies, {ǫi}, are taken from an uncorrelated symmetricrandom distributionP({ǫi}) = ΠiP(ǫi). (3.4)For Anderson disorder, P(ǫi) is the completely random, flat distributionP(ǫi) ={1/(2∆) if − ∆ ≤ ǫi ≤ ∆0 otherwise,(3.5)where ∆ can be thought of as the strength of disorder. For a binary alloy,P(ǫi) = xδ(ǫi − ǫA) + (1 − x)δ(ǫi − ǫB) , with x being the concentration ofA-type atoms and energies are shifted so that xǫA + (1 − x)ǫB = 0.Our aim is to calculate the single polaron Green’s function (GF) of thisproblem and average it analytically over all disorder configurations given byEq. (3.4) and (3.5). From now on, we use an overbar to denote disorder-averaged quantities. The strength of disorder, σ ≡ (ǫ¯2i )1/2, is taken to beweak compared to polaron bandwidth in the clean system, so that it can betreated perturbationally. As a result, the polaronic picture remains valid butits lifetime is expected to become finite due to scattering from the disorderpotential Vˆd.Since Vˆd is weak compared to other terms in the Hamiltonian, we treat itas a perturbation. Dividing the Hamiltonian as H = HH + Vˆd, where HH isthe Hamiltonian of the Holstein polaron in the clean lattice, we use Dyson’sidentity, Gˆ(ω) = GˆH(ω) + Gˆ(ω)VˆdGˆH(ω), to relate Gˆ(ω) of the system withdisorder to that of the clean system, GˆH(ω). To the second order in Vˆd, wefindGˆ(ω) ≈ GˆH(ω) + GˆH(ω)VˆdGˆH(ω)+GˆH(ω)VˆdGˆH(ω)VˆdGˆH(ω). (3.6)Because disorder breaks translational invariance, the eigenstates for anysingle disorder realization are not labelled by the lattice momentum k.However, averaging over all disorder configurations restores the transla-tional invariance and makes k a good quantum number again. As a result,〈0|ckGˆ(ω)c†k′ |0〉 = δk,k′G¯(k, ω) and we only need to calculate the diagonalmatrix elementG¯(k, ω) = GH(k, ω) +∑iǫ¯i〈0|ckGˆH(ω)c†i ciGˆH(ω)c†k|0〉+∑i,jǫiǫj〈0|ckGˆH(ω)c†i ciGˆH(ω)c†jcjGˆH(ω)c†k|0〉. (3.7)403.2. The model and its solutionHere, GH(k, ω) = 〈0|ckGˆH(ω)c†k|0〉 is the polaron GF in the clean system:GH(k, ω) =1ω − ε(k) − ΣMA(ω) + iη(3.8)where ε(k) = −2t∑i cos(ki). For completeness, its MA solution is brieflyreviewed in Appendix A.Since ǫ¯i = 0 for symmetrically distributed disorder, the first-order con-tribution vanishes3. Furthermore, because the disorder is uncorrelated,ǫiǫj = ǫ¯2i δij ≡ σ2δij and the disorder-averaged GF becomesG¯(k, ω) −GH(k, ω) = σ2∑i〈0|ckGˆH(ω)c†i ciGˆH(ω)c†i ciGˆH(ω)c†k|0〉. (3.9)The challenge is to use MA to calculate the matrix elements appearing inthe sum. We do this to the same level of accuracy as GH(k, ω) is evaluated[48].These matrix elements can be broken into products of generalized GFby inserting identity operators, I, between the creation and annihilationoperators. Since the MA flavor we use here is equivalent with assuming thatthe phonon cloud only extends over one site [48], at this level of accuracyit suffices to truncate I ≈ ∑l,n(1/n!)b†nl |0〉〈0|bnl , i.e., to ignore states withphonons at two or more sites (such states can be added systematically inhigher flavors of MA). Doing so results inG¯(k, ω) −GH(k, ω) = σ2∑i,l,s,n,m1n!m!〈0|ckGˆH(ω)c†i b†nl |0〉×〈0|bnl ciGˆH(ω)c†i b†ms |0〉〈0|bms ciGˆH(ω)c†k|0〉.This expression involves two different generalized propagators that we haveto evaluate, namely 〈0|ckGˆH(ω)c†i b†nl |0〉 and 〈0|bnl ciGˆH(ω)c†i b†ms |0〉.First, as detailed in Appendix C, to the level of accuracy of MA0, forn ≥ 1 both propagators vanish unless i = l. Therefore, we only have toevaluate F (n)ki (ω) ≡ 〈0|ckGˆH(ω)c†i bi†n|0〉 for n ≥ 1. Note that F(0)ki (ω) isalready known: F (0)ki (ω) = 〈0|ckGˆH(ω)c†i |0〉 = GH(k, ω) exp(−ik ·Ri)/√N ,where we have used 〈0|ckGˆH(ω)c†k′ |0〉 = δkk′GH(k, ω). Here, N → ∞ is thenumber of sites in the lattice. The details of the calculation for n ≥ 1 arepresented in Appendix C. The final result isF (n)ki (ω) = Γn(ω)F(0)ki (ω), (3.10)3A finite average can be removed trivially by an overall energy shift.413.3. Resultswhere Γn(ω) are easy to calculate products of continued fractions.Next, we calculate W nm(ω) = 〈0|bni ciGˆH(ω)c†i b†mi |0〉. Note that becauseof the invariance to translations in the clean system, this quantity is actuallyindependent of i. The detailed derivation of these functions is also presentedin Appendix C.With these expressions in hand, the disorder averaged GF, to secondorder in σ, becomesG¯(k, ω) = GH(k, ω) + σ2 [GH(k, ω)]2∞∑n,m=0Γn(ω)W nm(ω)Γm(ω)n!m! . (3.11)To the same order, this identifies the sum in Eq. 3.11 as the disorder self-energy in the presence of coupling to phonons,Σdis(ω) = σ2∞∑n,m=0Γn(ω)W nm(ω)Γm(ω)n!m! . (3.12)This is our main result. From a computational point of view, because thefactorials in the denominator grow rapidly with increasing index, the infinitesums can be safely truncated at finite values for n and m. Cutoffs of 20proved sufficient for all cases we examined. Using GH(k, ω) = 1/(ω− ε(k)−ΣMA(ω) + iη) we can finally writeG¯(k, ω) = 1ω − ε(k) − Σtot(ω) + iη, (3.13)where the total self-energy is Σtot(ω) = ΣMA(ω) + Σdis(ω). This implicitsummation gives a more accurate expression for the disorder-averaged GFthan Eq. (3.11), with which it agrees to O(σ2).3.3 ResultsWe are now prepared to study the effect of weak disorder on the polaronlifetime and energy. At this level of perturbation theory, disorder only en-ters through its standard deviation σ. In the following, we assume Andersondisorder of width 2∆, for which σ = ∆/√3. We will use either ∆ or σ tocharacterize the disorder, as convenient, but we emphasize that any othertype of disorder that has the same σ would lead to the same answer withinthis perturbational approximation. To characterize the electron-phonon cou-pling, it is convenient to use the effective coupling λ = g2/(6tΩ).423.3. ResultsOnce the GF is known, the energy broadening of an otherwise coher-ent polaron state of momentum k, which is inversely proportional to itslifetime, is given by the width of the lowest peak in the spectral function,A¯(k, ω) = − 1pi ImG¯(k, ω). This broadening measures the rate at which thepolaron leaves that momentum state due to scattering from the impurity po-tential Vˆd. In a clean system, the polaron states below the phonon-emissionthreshold, EGS + Ω, are infinitely long lived, therefore the low-energy spec-tral weight is a Dirac delta function (in fact, a Lorentzian of width η → 0).Mathematically, this is a consequence of the fact that (in the absence ofdisorder) the polaron self-energy ΣMA(ω) has a vanishing imaginary partfor all energies inside the polaron band.Disorder-induced finite lifetime broadens the delta functions into almostLorentzians. As a reference, we first review the case without electron-phononcoupling, λ = 0, for which the only nonzero term in Eq. (3.12) correspondsto m = n = 0. ThereforeΣdis(ω) = σ2W 00(ω) = σ2g0(ω), (3.14)where g0(ω) is again the momentum-averaged single particle GF,g0(ω) =1N∑k1ω − ε(k) + iη . (3.15)The resulting spectral weight has a peak of width τ−1k centred at energy Ekwhich is given by the pole conditionω − ε(k) − Σdis(ω) + iη = 0, (3.16)where ω = Ek − iτ−1k and η → 0+.Because τ−1k ∼ σ2 is small for weak disorder, we approximate Σdis(ω) ≈Σdis(Ek) +O(σ4). Using this in Eq. (3.16) gives Ek and τ−1k as follows [62]:Ek = ε(k) + ReΣdis(Ek)τ−1k = −ImΣdis(Ek).(3.17)The first expression determines the energy shift compared to the electronenergy in the clean system, ε(k). Since Re[g0(ω)] is negative for ω < 0 andpositive for ω > 0, this implies that the energy band slightly widens up inthe presence of disorder.Using Eq. (3.14), the inverse lifetime becomes τ−1k = −σ2Im[g0(Ek)].However, Im[g0(Ek)] is proportional to the total density of states (DOS) for433.3. Resultsthe clean system,Im[g0(Ek)] = −πN∑k′δ(Ek − ε(k′)) = −πρ0(Ek),and thereforeτ−1k = πσ2ρ0(Ek) = πσ2ρ0(ε(k)) + O(σ4). (3.18)This is simply Fermi’s golden rule. Since the density of states vanishesoutside the energy band of the clean system, ε(k) ∈ [−6t, 6t], this resultpredicts infinite lifetime for all states with |Ek| ≥ 6t, and finite lifetime forall states within the free particle band. We will return to this point later.-6.5 -6.0 -5.5Ek/t0.0000.0020.0040.0061/tτk FGR, ∆=0.2tFGR, ∆=0.4tATA, ∆=0.2tATA, ∆=0.4tMA, ∆=0.2tMA, ∆=0.4t-6.5 -6.0 -5.5ω/t0.000.010.020.030.04ρ(ω)∆ = 0∆ = 0.2t∆ = 0.4t(a) (b)Figure 3.1: (a) Inverse polaron lifetime 1/τk vs its peak energy Ek, and (b)average density of states ρ(ω), for a weak electron-phonon coupling and twovalues of the disorder ∆. The solid and dashed lines are the correspondingFermi golden rule (FGR) and ATA results, respectively (see text for moredetails). Other parameters are Ω = t, η/t = 10−2 in (a) and η/t = 5 × 10−3in (b).443.3. Results3.3.1 Weak electron-phonon couplingThe analysis is performed similarly in the presence of electron-phonon cou-pling, but then using the appropriate total self-energy, Eq. (3.13). We firstconsider weak electron-phonon coupling, λ = 0.5. In Fig. 3.1(a) we plotthe polaron inverse lifetime for states in the lowest polaron band, for twodifferent values of the disorder strength, ∆ = 0.2t and 0.4t (squares andcircles, respectively). These values are extracted from Lorentzian fits to thelowest peak in A¯(k, ω). The broadening η is decreased until Ek and τ−1kconverged to the values presented in Fig. (3.1) that are independent of η.For this small λ, the MA ground state energy of the polaron in the cleansystem is EGS = −6.534t. The weak disorder does not shift the eigenstatessignificantly. In fact, as shown in Fig. 3.1(b), the average density of statesρ(ω) = − 1pi Im∑k G¯(k, ω) is nearly identical to that of the clean system,although the band becomes slightly broader with increasing ∆. The inverselifetime vanishes below the clean system band edge, EGS , and it increases as√Ek − EGS above, which is the expected DOS for the clean system at thebottom of the band. This is very similar to the λ = 0 result, except for therenormalization of the DOS by the electron-phonon interaction. Indeed, ifwe think of the polaron as a simple quasiparticle whose density of states isρ(ω) [renormalized from ρ0(ω) for a free electron], the inverse lifetimes wefind at the bottom of the polaron band are in good agreement with thosepredicted by Fermi’s golden rule (FGR), πσ2ρ(Ek) (see full lines).While the agreement between the two is good near the bottom of theband, it becomes systematically worse at higher energies. To verify that thisamount of disorder is still sufficiently small so that the disagreement is notdue to using perturbational results outside their validity range, we also showaverage T-matrix (ATA) results (dashed lines). ATA is a simple way to treatdisorder beyond the lowest order in perturbation theory, for a system withλ = 0. We briefly discuss it in Appendix D, as well as how we extend it forfinite λ. ATA converges to πσ2ρ(ω) in the limit of small σ, therefore theagreement between FGR and ATA confirms that the contribution of higherorder terms in σ is indeed negligible. The disagreement with MA at higherenergies is, therefore, not an artifact of using perturbation theory.The meaning of this disagreement at higher energies should, however,be treated with some caution. It is well known [51, 63] that this flavor ofMA fails to reproduce the correct polaron+one phonon continuum, whichshould start at EGS + Ω (this problem is fixed by MA1 and higher flavors).One consequence is that MA overestimates the bandwidth of the polaron atweak couplings. Indeed, in Fig. 3.1(a) we see that the polaron band extends453.3. Resultswell past EGS + Ω. In other words, we know that at these higher energiesMA is not accurate enough, so the results shown in Fig. 3.1 should only betrusted close to the bottom of the band, where the agreement is good.0 0.1 0.2 0.3 0.4 0.5∆/t0.0000.0010.0021/tτkpolaron, λ = 0.5free electron, λ = 0polaron, λ = 0.5, axes rescaled by t/t*Figure 3.2: Inverse lifetime of the polaron of momentum k = (π/8, 0, 0) andλ = 0.5 vs the strength of disorder, ∆/t (full squares). The inverse lifetimefor a free electron with the same momentum is shown by triangles. Emptysquares show 1/(t∗τk) vs ∆/t∗ for the polaron (for this λ, t∗ = 0.881t).Other parameters are Ω = t, η/t = 10−2.The monotonic increase of the polaron’s inverse lifetime with increasingdisorder strength is shown in Fig. 3.2, for the lowest polaron state of momen-tum k = (π/8, 0, 0) and λ = 0.5 (full squares). For comparison, also shownis the corresponding lifetime of a bare electron (λ = 0, triangles) with thesame momentum. Both curves show the expected ∝ σ2 increase predictedby Fermi’s golden rule, but the polaron lifetime is somewhat shorter. Themost likely reason for this is the renormalization of the polaron mass byelectron-phonon coupling. Indeed, if instead we plot 1/(t∗τk) vs ∆/t∗ usingthe polaron lifetimes, the results are much closer to those of the free electron,especially for small values of the disorder (empty squares).The conclusion, thus far, is that Fermi’s golden rule agrees well with our463.3. Resultsresults at energies where this flavor of MA can be trusted. In other words, atweak electron-phonon coupling, the effect of disorder can be quantitativelyunderstood if we think of the polaron as a simple particle with a renormalizedmass (or DOS), and use Fermi’s golden rule.-7.75 -7.70Ek/t01×10-42×10-43×10-44×10-41/tτk -7.75 -7.70Ek/t-7.75 -7.70Ek/t-7.75 -7.70ω/t0.100.050.0ρ(ω)-7.75 -7.70ω/t-7.75 -7.70ω/t∆ = 0.05t ∆ = 0.10t ∆ = 0.15tFigure 3.3: (Color online). Top panels: 1/(tτk) vs Ek/t for three levels ofdisorder: ∆/t = 0.05, 0.1 and 0.15. The symbols shows the MA result for astrong coupling λ = 1.2, while the full and dashed lines show Fermi’s goldenrule and the ATA predictions, respectively. Bottom panels: The averageDOS ρ(ω) for that ∆ (full line) and the DOS in the clean system, ρH(ω)(dashed line) vs ω. Parameters are Ω = t, η = 10−3t.3.3.2 Strong electron-phonon couplingWe now check whether this also holds true at strong electron-phonon cou-pling. We use λ = 1.2 for which a robust small polaron appears in the cleansystem. The top panels in Fig. 3.3 show the polaron inverse lifetime vs.its energy for increasing strength of disorder (symbols), as well as Fermi’sgolden rule (full lines) and the ATA (dashed lines) predictions. The latter473.3. Resultstwo are indistinguishable, confirming that these levels of disorder are in-deed perturbationally small. The bottom panels show the average densityof states, ρ(ω) (full lines). For comparison, the polaron DOS in the absenceof disorder, ρH(ω), is also shown (dashed line).Let us first consider the DOS. As in the other cases we see that, byincreasing disorder, the band edge shifts down to lower energies. However,the effect is quantitatively much more significant here compared to the weakcoupling case, Fig. 3.1(b), because the polaron band is much narrower. Thisis seen in Fig. 3.4, where we show the same densities of states but over thewhole polaron band. One surprise is that the entire band moves to lowerenergies with increasing disorder. This is different from what happens for abare particle, where the band broadens symmetrically on both sides. Thedifferent behavior at the upper edge is likely due to the difference in theirspectra. While for a single particle the energy band is the only feature inits spectrum, the spectrum of the polaron is quite complicated, with manyother features, such as a band associated with the second bound state, thepolaron+one-phonon continuum, etc., lying above the lowest polaron band[48]. With increasing disorder all these features should move toward lowerenergies. Level repulsion from these high energy states would explain whythe upper edge of the polaron band moves to lower energies.The inverse lifetime again vanishes for states whose energy lies belowthe band edge of the clean system, Ek < EGS ≈ −7.73t. Because theshift of the disorder-averaged DOS is now significant, this means that, fora quite large energy range at the bottom of the band, the polaron has aninfinite lifetime despite the presence of disorder. We emphasize that thisis qualitatively similar to the result for a bare particle at the bottom of itsband; it is simply more pronounced here. The meaning of this (un-physical)infinite lifetime for these low-energy states is discussed extensively below;briefly, we believe that it signals a failure of the perturbation theory at theseenergies. These low-energy states are most susceptible to localization, so theperturbational calculation and its predictions are not reliable here.For high energy polaron states with Ek > EGS , the lifetime in the pres-ence of disorder becomes finite, as expected. However, here the MA resultsdisagree quantitatively with the FGR and ATA results at all energies. Thelatter two are nearly indistinguishable, suggesting again that these levels ofdisorder are small enough that perturbation theory should be valid. Thedisagreement cannot be blamed on MA, either: at such strong couplingsand correspondingly low energies, MA is extremely accurate for the entirepolaron band [14, 48]. The disagreement is, therefore, physically meaning-ful and its origin can be quite easily traced. If we explicitly separate the483.3. Results-7.75 -7.7 -7.65ω/t0.000.040.080.12ρ(ω)∆ = 0∆ = 0.05t∆ = 0.10t∆ = 0.15tFigure 3.4: (Color online) The same average DOS vs ω displayed in thelower panels of Fig. 3.3, but now shown for the entire polaron band.n = m = 0 contribution in the disorder self-energy, Eq. (3.12) becomesΣdis(ω) = σ2g0(ω − ΣMA(ω)) + σ2∑n+m>0Γn(ω)W nm(ω)Γm(ω)n!m! .If the contribution of terms with n + m > 0 can be ignored, this resultreduces to Fermi’s golden rule, since ρ(ω) = − 1pi Img0(ω − ΣMA(ω)). Thedisagreement between MA and FGR, then, comes from the contribution ofthe terms with n+m > 0. These terms cannot be ignored at strong electron-phonon coupling. Consider, for instance, F (n)ki (ω) = 〈0|ckGˆH(ω)c†i b†ni |0〉,which is proportional to Γn(ω). If we Fourier transform to real times, F (n)ki (τ)is proportional to the amplitude of probability that if an electron is injectedin the system, it is found in the presence of n bosons at a later time τ ,all at the same site. At strong couplings, the electron dresses itself with alarge phonon cloud to become a polaron, so the probability of finding it withmany phonons in its vicinity should be considerable, while the probabilityof finding the electron without any phonons (n = 0) is exponentially small.In the large λ limit, the terms which are expected to contribute most are493.3. Results= +(i)(ii) = +=+ + . . .++ ++ ++ + + + . . .=FGRMAFGR(iii)(iv)Figure 3.5: (i) FGR approximation for the disorder-averaged Green’s func-tion of a carrier (double thin line), in terms of that of the free carrier (thinline) and uncorrelated disorder (dashed line) in the absence of electron-phonon coupling; (ii) FGR approximation for the disorder-averaged GF ofa polaron (double thick line), in terms of that of the clean polaron (thickline); (iii) clean polaron GF in terms of free carrier (thin lines) and phonon(curly lines) propagators; (iv) the first few terms in the MA approximationfor the disorder averaged GF of a polaron.those with n ≈ g2/Ω2, i.e., values close to the average number of phononsin the polaron cloud. In contrast, for small λ the phonon cloud is veryfragile and, in fact, most of the time the electron is alone (resulting in alarge quasiparticle weight). This is why, for weak coupling, keeping only then = m = 0 term in the sum is a good approximation.These considerations are illustrated diagramatically in Fig. 3.5. Panel503.3. Results(i) shows the Born approximation for the disorder averaged GF of a car-rier, which leads to Fermi golden’s rule expression for the lifetime, as al-ready discussed. Panel (ii) shows the same for the disorder averaged GFof the polaron; as discussed above, this is equivalent with keeping only then = m = 0 term in the disorder self-energy, Eq. (3.12). Since each cleanpolaron propagator starts and ends with a free carrier propagator, this ap-proximation means that the electron can scatter from disorder only in theabsence of phonons; this is why this approximation fails at large electron-phonon coupling, where the chance of having no phonons around during animpurity scattering is extremely weak. In contrast, the full MA expressionincludes diagrams such as shown in panel (iv), where phonon and disorderlines cross. One can think of these as leading to an effective renormalizationof the disorder strength, especially since these diagrams are very similar tothose which result in the renormalization of a single impurity potential [1],details discussed in chapter 2.The complete evolution of τ−1k and Ek with disorder ∆ at strong electron-phonon coupling is shown in Fig. 3.6 for two different momenta k1 =(2π/9, 0, 0) and k2 = (π/6, 0, 0), whose free electron energies are ε(k1) ≈−5.5t and ε(k2) ≈ −5.7t, respectively. In Fig. 3.6(a), we trace their in-verse lifetimes as a function of disorder. At small ∆, the energies of boththese states are well above EGS and their scattering rates increase mono-tonically with ∆, as one would expect on general grounds. However, theinverse lifetimes reach a maximum after which they begin to decrease fastand vanish eventually. The value of ∆ where τ−1k vanishes corresponds tothe disorder at which Ek = EGS . For larger disorder, Ek moves below thefree polaron band edge, and its lifetime remains infinite. This is more clearlyshown by Fig. 3.6(b), where the inverse lifetimes are plotted vs. the corre-sponding eigenenergy Ek for these two momenta, as disorder is increased.This confirms that the scattering rates for both polaron states vanish whentheir energy drops below EGS , whose location is marked by the asterisk.The value of ∆ where this happens depends on how far above EGS was theenergy Ek of this polaron, in the limit of ∆ → 0.Figure 3.3 shows that using the FGR expression, i.e. τ−1k = πσ2ρ(Ek),is quantitatively wrong. We can also compare the polaron lifetime, wherefinite, with that of a free particle of renormalized mass, similar to the com-parison in Fig. 3.2. This is shown in Fig. 3.7, where we compare 1/(t∗τk)vs ∆/t∗ for the polaron, with the inverse lifetime of a free electron at thesame momentum. While roughly quadratic dependence is observed for thepolaron at small disorder, the coefficient is quite different from that for thefree electron. For strong disorder, the disagreement is even worse.513.4. Summary and conclusions0.00 0.05 0.10 0.15∆/t01×10-42×10-41/tτkk1k2-7.70-7.72Ek/t01×10-42×10-41/tτkk1k2(a) (b)increasingdisorder∗Figure 3.6: (Color online) (a) Inverse lifetime vs disorder, and (b) inverselifetime vs energy Ek, as disorder is turned on, for two momenta k1 =(2π/9, 0, 0) and k2 = (π/6, 0, 0), for a polaron with λ = 1.2,Ω = t, η =10−3t. The asterisk in panel (b) marks the clean polaron GS energy in theclean system, EGS , for these parameters. See text for more details.This shows that for intermediate and large electron-phonon coupling,where a heavy small polaron forms, its lifetime in the presence of disorderis not described quantitatively by the predictions corresponding to a bareparticle with renormalized mass. The polaron has an internal structurewhich manifests itself as a significant deviation from Fermi’s golden rulepredictions even for weak disorder. The scattering of the electron in thepresence of its phonon cloud is quite different from that of a bare particleof the same effective mass, but without a cloud [1, 32].3.4 Summary and conclusionsUsing MA to deal with the electron-phonon coupling and perturbation the-ory to deal with the weak disorder, we derived an expression for the disorder-averaged GF of the Holstein polaron in a simple cubic lattice with random523.4. Summary and conclusions0 0.1 0.2 0.3 0.4 0.5∆/t*00.0010.0021/t*τkpolaron, λ = 1.2, t*= 0.148tfree electron, λ = 0, t*= tFigure 3.7: (color online) Same data as shown in Fig. 3.6(a) but withrescaled axes; 1/(t∗τk) vs ∆/t∗ for momentum k2 = (2π/9, 0, 0) (emptysquares) is compared with the free electron lifetime (triangles) at the samemomentum, for λ = 1.2 where t∗ = 0.148t. Other parameters are Ω = t andη/t = 10−3.on-site energies. This allowed us to find an analytic expression for the lowest-order contribution from disorder to the polaron self-energy.The disorder-averaged spectral function was used to extract the lifetimeand energy shift of various polaron states. For weak electron-phonon cou-pling, we found that the MA results are in reasonable quantitative agreementwith those predicted by Fermi’s golden rule for a free particle with a properlyrenormalized mass.At intermediate and larger electron-phonon coupling where a small po-laron forms, however, the MA results quantitatively disagree with Fermi’sgolden rule estimate everywhere the lifetime is finite and for all levels ofdisorder. The reason for this is the fact that the scattering of the electronin the presence of its (robust) phonon cloud is quite different from the scat-tering of a bare particle with renormalized mass. This is the same physicsthat leads to a significant renormalization of the disorder potential seen bya polaron as compared to the bare disorder, as discussed in the previous533.4. Summary and conclusionschapter [1, 32]. This demonstrates again that, in the small polaron limit, itis wrong to assume that the only effect of the polaron cloud is to renormalizethe polaron’s mass.It is important to note that this calculation is only valid for weak disor-der. It is based on perturbation theory, and in principle it can be improvedby going to higher orders along the same lines we used to calculate thelowest-order contribution. However, one should remember that as disorderbecomes stronger, Anderson localization will eventually occur, and that thiscannot be captured within perturbation theory. Also, once disorder is largeenough to lead to localization, the disorder-averaged GF loses its meaningand usefulness. Instead, here the signature of localization becomes manifestin the distribution of various quantities such as the local density of states,not in their average value.A surprise, at least at first sight, is the fact that this calculation predictsan infinite lifetime for a range of energies at the bottom of the polaronband. This interval can include a significant fraction of the polaron states,especially at stronger electron-phonon coupling and larger disorder. As wealready mentioned, this is in fact similar to what happens for a free particle,which also is predicted, within this level of perturbation theory, to have aninfinite lifetime for all momenta for which |Ek| > 6t. The difference is onlyquantitative: the energy shift for a free particle is tiny compared with itsbandwidth of 12t, whereas for a small polaron this shift can be comparablewith its significantly narrower bandwidth even for rather weak disorder.A likely reason for this can be inferred from the fact that, for a freeparticle, states at the band edge become localized immediately upon in-troduction of disorder. In other words, we already know that there is afinite range of energies (which, for weak disorder, falls outside the free par-ticle bandwidth) where treating disorder perturbationally and calculatingthe disorder-averaged GF is meaningless. It is then reasonable to concludethat the states for which this perturbational scheme predicts infinite life-times are, in fact, already localized. If this is correct and generalizes tothe polaron case, it suggests that, unlike for a free particle, for a polaronlocalization sets in differently at the lower vs the upper polaron bandedge.We have already speculated that this difference may be due to the influenceof the higher-energy states that exist in the polaron spectrum. Confirma-tion of these conclusions will require a study going beyond a perturbationaltreatment of disorder.54Chapter 4Binding carriers to anon-magnetic impurity in atwo-dimensional square Isinganti-ferromagnet4.1 IntroductionIn chapters 2 and 3, I studied some of the effects caused by disorder ona Holstein polaron, which is a quasiparticle resulting from the dressing ofa charge carrier by lattice vibrations. I change gears in this chapter tofocus on a different type of charge carrier-boson coupling problem, namely ahole in a two-dimensional anti-ferromagnet (AFM) [64–68]. This is in partmotivated by the physics of high-Tc superconductivity in cuprates. In theirparent compound, the strong hybridization between copper 3d and oxygen 2porbitals drives their CuO2 layer into a correlated insulating state in which theholes on neighboring copper atoms align their spins anti-ferromagnetically inorder to gain the superexchange energy, J . Upon doping these AFM layerswith charge carriers, superconductivity emerges [69, 70].A major setback in the search for an analytic description of the behaviorof these charge carriers is the lack of a simple wave function for the groundstate of the undoped AFM. The semi-classical Ne´el state breaks spin rota-tion symmetry and is therefore smeared out by quantum spin fluctuationsto a significant degree that is hard to capture with simple wave functions[71]. This leaves numerical calculations as the only way to make quanti-tative predictions [72]. While implementing such numerical calculations isalready a complicated task even for a clean system, a further complicationcomes from the presence of disorder and imperfections in the real materials,introduced during the sample growth and preparation. Given the low di-mensionality, even weak disorder may have dramatic effects on the motionof charge carriers in the CuO2 layers.554.1. IntroductionImpurities have been shown to be responsible for a range of phenomenain low-dimensional correlated electron systems, and they can be also uti-lized for probing correlations which are otherwise difficult to observe in theground state [73]. For the undoped parent compound, mean-field analysisof the disordered Hubbard model predicts the emergence of an inhomoge-neous metallic phase in which the Mott gap is locally closed wherever thedisorder is strong enough to do so [74]. However, it is not always the casethat impurities destroy the order in the underlying system. For instance,impurities induce local magnetic order in one-dimensional (1D) quantummagnets [75], and long-range antiferromagnetism is predicted upon dopingsome quantum spin liquids with nonmagnetic impurities [76]. In any event,a complete understanding of the interplay between disorder and AFM corre-lations and especially of their role in controlling the carrier dynamics awayfrom half-filling is still lacking.Here, I consider a much simpler variant of this problem where, at zerotemperature, a hole is created in a 2D Ising AFM on a square lattice, andis also subject to the on-site attractive potential of an impurity that canbe visited by the hole. Thus, this model is very different from previousmodels of an impurity in a 2D Heisenberg AFM, which assumed that thehole cannot visit the impurity site, and is coupled to it at most throughexchange [77, 78]. As we discuss in the following, our results have somesimilarities but also considerable differences from those obtained numericallyin these other models.I investigate the local density of states (LDOS) near this impurity tostudy the appearance of bound states, focusing specifically on the relevanceof the magnetic sublattice on which the impurity is located. The advantageof our approach is that the wave function of the undoped 2D AFM is the sim-ple Ne´el state, and this allows us to study the problem (quasi)analytically.Of course, spin fluctuations are completely absent, but as we argue in ourdiscussion, our results allow us to speculate about (at least some of) theirlikely effects.A single hole in an Ising AFM was initially believed to be localized evenin the absence of impurities, because when the hole hops it reshuffles thespins along its path, thereby it creates a string of wrongly oriented spins.Finite mobility was believed to arise only due to spin fluctuations which canremove pairs of such defects [79, 80], but they are absent from the IsingHamiltonian.However, it is discussed in detail in the following sections that the holeis actually delocalized even in the Ising AFM by going twice around closedloops. The string of misaligned spins that are created in the first round is564.2. The Modelremoved when spins are reshuffled again during the second round. When thelast one is removed, the hole ends up at a different site from where it started,and by repeating this process it can move anywhere on its original sublattice(spin conservation ensures that the hole propagates on one sublattice). Thisraises the question of how the hole’s motion will be affected by an attractiveimpurity, especially by one located on the other sublattice than the oneon which the hole resides. While one expects the hole to become boundto the impurity if they are on the same sublattice, if they are on differentsublattices, one may expect the hole not to be sensitive to the presence of theimpurity and therefore remain unbound. We investigate this problem usinga variational method introduced in Ref. [2] to study the clean case, whichwe generalize here to systems that are not invariant under translations.I introduce the model and its variational solution in Sec. 4.2. The gener-alization of this variational method to inhomogeneous systems is discussedin Sec. 4.3, followed in Sec. 4.4 by results for a single impurity located (i) onthe same, and (ii) on the other sublattice than the quasiparticle. I concludeby giving a summary and discussing possible further developments of thiswork in Sec. 4.5.4.2 The ModelWe consider the motion of a single hole doped into a spin-12 Ising antifer-romagnet on a 2D square lattice. The Hamiltonian of the undoped systemisHAFM = J∑<i,j>[Szi Szj +14] = J¯∑<i,j>[σzi σzj + 1], (4.1)where σz is the Pauli matrix and J¯ = J/4 (> 0). Since the square latticeis bipartite, the vacuum state of Hamiltonian (4.1), |0〉, is the Ne´el-orderedstate where all spins on one sublattice point up and those on the other sub-lattice point down. Excitations are gapped spin-flips that are like localizedmagnons. The creation operator for a magnon is written in terms of thespin raising and lowering operators, σ± = σx ± iσy:d†i ={σ−i if i ∈↑ sublattice,σ+i if i ∈↓ sublattice.(4.2)Consider now the doped case, Fig. 4.1. Creating a hole on a particularlattice site corresponds to removing the spinfull particle located at that site.574.2. The ModelTherefore the hole creation operators areh†i ={ci↑ if i ∈↑ sublattice,ci↓ if i ∈↓ sublattice,(4.3)where ciσ is the fermion annihilation operator. Once the hole is created (h†i ),it can move in the lattice via nearest-neighbor (nn) hopping. As each site’snn belong to the other sublattice, when the hole switches sites with a nearbyspin, the latter will find itself surrounded by parallel spins. Therefore, amagnon is created on this site (d†i ). The same hop may annihilate a magnonfrom the hole’s arrival site (dj), if there was one already there. Therefore,the Hamiltonian can be written as [2]H = HAFM + P{−t∑<i,j>[h†jhi(d†i + dj) + H.c.]}P − Uh†0h0, (4.4)where P is the projection operator enforcing no double occupancy: at anysite, there is a hole or there is a spin which is either properly oriented or isflipped, h†ihi+d†idi+did†i = 1. Thus, the second term describes the hoppingof the hole which is accompanied by either magnon creation or annihilation.In addition, there is an attractive potential of strength U centered at theorigin r = 0, which changes the on-site energy of the visiting hole (variationsof the local hoppings and exchanges can be trivially included in the modeland our solution, but should not lead to any qualitative changes if theyare small or moderate in size). Physically, such a potential can be due toan attractive non-magnetic impurity located above the origin, in a differentlayer, which modulates the on-site energy at the origin. Another possibilitycomes from replacing the atom at the origin by an impurity atom with thesame valence, but whose orbitals lie at lower energies than those of thebackground atoms. This is very different from the impurity models studiedin previous work where the impurity is an inert site that can not be visited bycarriers [77], and there is at most exchange between the spin of the impurityand that of carriers located on neighboring sites [78].Note that bosons in this model actually help the hole to move in thelattice by modulating its kinetic energy. This contrasts with the Holsteinmodel, where the bosons (phonons) affect the potential energy of a particlethat was already mobile, thereby increase its effective mass. Hamiltonian(4.4) is similar to the limiting form (tf → 0) of the Edwards fermion-bosonmodel which describes the modulation of a particle’s hopping via coupling584.3. Propagation of the hole in the clean systemto dispersionless bosons (b) [81–83]:HEdw. = −tf∑<i,j>(c†i cj + H.c.) + Ω∑ib†i bi − tb∑<i,j>c†i cj(b†j + bi). (4.5)Unlike in our model, however, there is no constraint in this model regardingthe number of bosons at any site which can be simultaneously occupied bythe particle as well. Furthermore, in a Ne´el AFM the energy of neighboringmagnons is not additive as in the Edwards model.Figure 4.1: Square lattice of Nee´l-ordered Ising spins. Removing one of thespin-up particles creates a hole and results in a configuration with sz = −12 .4.3 Propagation of the hole in the clean systemIn this section, we construct the equations of motion for the zero-temperatureretarded Green’s function (GF) of the hole moving through the lattice in theabsence of impurity, U = 0. Due to the translational invariance in the cleansystem, one can do this in the momentum-space and calculate G(k, ω), whichis done in Ref [2]. Here, we present a real-space derivation which is neededonce the translational invariance is broken when the impurity is added. Thereal-space single hole retarded GF is defined asG0,R(ω) = 〈0|h0Gˆ(ω)h†R|0〉, (4.6)594.3. Propagation of the hole in the clean systemwhere Gˆ(ω) = limη→0+ 1/(ω −H + iη) is the resolvent of Hamiltonian (4.4)when U = 0. By separating the Hamiltonian as H = HAFM + Ht where Htis the second term responsible for hopping, equations of motion for G0,R(ω)can be generated by repeated use of the Dyson identityGˆ(ω) = GˆAFM(ω) + Gˆ(ω)HtGˆAFM(ω),in which GˆAFM(ω) = limη→0+ 1/(ω − HAFM + iη). Using this, Eq. (4.6)becomesG0,R(ω) = g0(ω)[δ0,R − t∑uF1(R,u, ω)], (4.7)where g0(ω) = 1/(ω − 4J¯ + iη) and 4J¯ is the cost of breaking four AFMbonds when the hole is created in the lattice. Here, the lattice constant isset to unity, a = 1, u = ±x,±y are the four nn vectors. F1 is a new GFwhich describes a process where the hole starts together with a magnon nextto it,F1(R,u, ω) = 〈0|h0Gˆ(ω)d†Rh†R+u|0〉. (4.8)To simplify the notation, from now on we do not explicitly write the depen-dence on ω for all these GFs.The equation of motion for F1 can be similarly generated. Upon applyingthe Dyson identity, the hole can hop back to R and remove the magnon,or it can hop further away and create a second magnon, with an associatedGF F2, and so on. However, states with many magnons are less likely tooccur due to the energy cost of creating the magnons. In order to avoid therise in the number of magnons, the hole can trace back its path to removethem, but this effectively confines the hole to the vicinity of its creation site,and results in an immobile quasiparticle which is surrounded by a cloudof magnons. It turns out that the hole is freed to move on the lattice byexecuting the so-called Trugman loop processes in which it goes twice arounda closed path. By doing so, magnons that are created at the first pass areannihilated when the hole arrives at those sites for the second time. It isnot hard to see that when the very last magnon is annihilated, the hole endsup two hops away from its starting point. This is equivalent to a second orthird nn hop on the main lattice which translates into a first and second nnhop on the hole’s sublattice, respectively.Longer loops involve more costly intermediate states with many magnons.This suggests that we can proceed within a variational approach in which alimit is set for the maximum number of magnons that can be generated asthe hole propagates. We choose to work with up to three magnons, which isthe minimum number necessary for the hole to complete a loop. Moreover,604.3. Propagation of the hole in the clean systemFigure 4.2: Effective first- and second-nearest-neighbor hoppings of the hole(the blue circle) achieved with loops involving only up to three magnons.The latter is realized when the hole starts a second loop before removingthe last magnon it created during the first loop. The magnons are shownby red circles. The properly oriented spins are not shown.we only keep configurations consistent with these short closed loops (i.e., weexclude, for example, configurations where all three magnons are collinear).Figure 4.2 shows how both types of effective hoppings can be generated withthe three magnons types of configurations that we keep in our variationalcalculation. One can include more configurations in numerical simulations,but this was shown to result in only quantitative differences as long as t/Jis not too large [2, 84].Coming back to the equation of motion for F1, it relates F1 to G0,R andalso to three GFs with two magnons, F2. One of these, with the two magnonscollinear with the hole, cannot lead to a closed loop without generating morethan three magnons. Therefore we exclude it from the variational space, asdiscussed. Hence, we are left with only three termsF1(R,u) = −tg1[G0,R +∑v⊥uF2(R,u,v)], (4.9)where F2(R,u,v) = 〈0|h0Gˆ(ω)d†Rd†R+uh†R+u+v|0〉, v = ±x if u = ±y andvice versa, g1 = 1/(ω − 10J¯ + iη) and 10J¯ is the energy of the hole and amagnon next to each other. Within our variational space, the equation ofmotion for F2 isF2(R,u,v) = −tg2[F1(R,u) + F3(R,u,v,−u)], (4.10)whereF3(R,a,b, c) = 〈0|h0Gˆ(ω)d†Rd†R+ad†R+a+bh†R+a+b+c|0〉614.3. Propagation of the hole in the clean systemand g2 = 1/(ω − 14J¯ + iη). The other pair of three-magnon configurationsthat can be reached starting from d†Rd†R+uh†R+u+v|0〉 do not belong to ourvariational space and are hence discarded. Finally, in this variational spaceF3 relates to F2 onlyF3(R,u,v,−u) = −tg3[F2(R,u,v) + F2(R + u + v,−v,−u)], (4.11)with g3 = 1/(ω − 16J¯ + iη), 16J¯ being the energy of the allowed three-spin-defect configurations.These equations can be used to eliminate all F3, F2, F1 unknowns andbe left with equations involving only G0,R. The details are presented inAppendix E. The final results isG0,R(ω) = g¯0(ω)[δR,0 − t1(ω)∑δG0,R+δ(ω)−t2(ω)∑ξG0,R+ξ(ω)], (4.12)in which g¯0(ω) = 1/(ω − 4J¯ + 4tζ1(ω) + iη), t1(ω) = 2tζ3(ω), t2(ω) = tζ2(ω)and δ = ±u± v and ξ = ±2u are all the second and third nn vectors of themain lattice, respectively. The explicit expressions of ζ1,2,3(ω) functions aregiven in Appendix E.Equation (4.12) shows that the hole moves as a quasiparticle with theeffective second and third nn hoppings t1(ω) and t2(ω), respectively, and aneffective on-site energyε(ω) = 4J¯ − 4tζ1(ω). (4.13)The hole, in this quasiparticle, is accompanied by a cloud of magnons whichare constantly created and annihilated, helping it to freely propagate in thelattice. Note that sites R, R + δ, and R + ξ belong to the same sublattice.Therefore, the quasiparticle effectively propagates in the sublattice the holeis originally created on, and whose first and second nn vectors are δ andξ, respectively. The constraint that forces the quasiparticle to move onone sublattice is very general, being due to the spin-conserving nature ofthe Hamiltonian (4.4). It prevents the hole from ending up on the othersublattice without leaving an odd number of magnons behind: if the holestarts on one sublattice and ends up on the other one, the z component ofthe total spin angular momentum of the system changes from Szi = ±1/2 toSzf = ∓1/2, therefore there needs to be an odd number of magnons aroundto compensate for the change of spin Szf − Szi = ∓1.624.3. Propagation of the hole in the clean systemBefore presenting the real-space solution of Eq. (4.12), note that we havenow enough information to identify the momentum-space Green’s functionG(k, ω) = 〈0|hkGˆ(ω)h†k|0〉 =1ω − ǫ(ω,k) + iη , (4.14)where hk = 1√N¯∑r exp(−ik · r)hr and the sum runs over N¯ → ∞ sites ofthe hole’s sublattice. ǫ(ω,k) is the self-energy encoding the effect of beingdressed by magnons, that is responsible for the dynamical generation of thehole’s energy dispersion:ǫ(ω,k) = ε(ω) − 2t1(ω)[cos(kx + ky) + cos(kx − ky)]−2t2(ω)[cos(2kx) + cos(2ky)]. (4.15)As required, this is identical to the solution derived using a momentum spaceformalism in Ref. [2] and [85]. The spectral weight A(k, ω) = − 1pi ImG(k, ω)is then used to identify the quasiparticle excitations and their various prop-erties such as energy dispersion, effective mass, etc [2].Figure 4.3: The choice of coordinate systems for the lattice with impurity.The impurity, shown as purple square, is at the origin of the xy axes thatspan the original lattice with unit vectors x,y. The XY axes are rotated by45 ◦ and span the sublattice (black dots) on which the quasiparticle propa-gates via the elementary vectors y ± x.Equation (4.12) can be solved directly in real space by the method ofcontinued fractions detailed in Ref. [86]. For completeness, we briefly outlineit here. Let n and m be the x and y components of R 6= 0 on the coordinatesaxes XY which is rotated by 45 ◦ with respect to the lattice. It spans the634.3. Propagation of the hole in the clean systemsublattice on which the quasiparticle moves, marked by dots in Fig. 4.3(a),whose elementary vectors are y ± x. In this coordinate system, Eq. (4.12)can be written asGn,mg¯0= −t1(Gn+1,m +Gn−1,m +Gn,m+1 +Gn,m−1)−t2(Gn+1,m+1 +Gn+1,m−1 +Gn−1,m+1 +Gn−1,m−1), (4.16)where Gn,m ≡ G0,R, R = n(y + x) + m(y − x), is a shorthand notation.Eq. (4.16) can be expressed as a single-index recursive relation by groupingdistinct GFs with n ≥ m ≥ 0 into column vectors VM according to theirManhattan distance, defined as M = n+mVM=2r =G2r,0G2r−1,1...Gr,r, VM=2r−1 =G2r−1,0G2r−2,1...Gr,r−1.These are the only distinct GFs at Manhattan distance M and all others canbe related to these using symmetries: Gn,m = Gm,n = Gn,−m = G−n,m. Interms of these vectors, Eqs. (4.16) can be grouped into the following matrixformλrVr = α˜rVr−2 + αrVr−1 + βrVr+1 + β˜rVr+2 (4.17)for r ≥ 2 andV0 = g¯0(ω) + β0V1 + β0V2V1 = α1V0 + β1V2 + β˜1V3(4.18)for the GFs with M = 0, 1. Here, λ, α˜, α, β and β˜ are extremely sparsematrices whose elements can be read from Eq. (4.16). Combining two copiesof Eq. (4.17) for r = 2s− 1 and r = 2s givesΓsWs = AsWs−1 +BsWs+1, (4.19)where Ws =(V2s−1V2s)andΓs =(λ2s−1 −β2s−1−α2s λ2s), As =(α˜2s−1 α2s−10 α˜2s), Bs =(β˜2s−1 0β2s β˜2s).Since Eq. (4.19) links three consecutive terms, its solution can be expressedas a continued fraction (of matrices). Assuming a solution as Ws = ΩsWs−1and using it in Eq. (4.19) givesΩs = (Γs −BsΩs+1)−1As, (4.20)644.3. Propagation of the hole in the clean systemwhich can be evaluated for Ωs starting from a cutoff c such that Ωc+1 = 0.This results in the continued fraction solution for Ωs. In particular, thisgives Ω2 which relates W2 (set of V3 and V4) to W1 (set of V1 and V2).Finally, the diagonal element of Green’s function G00(ω) is found by usingthese in Eqs. (4.18) and (4.17) with r = 2 to solve for V0 = G0,0(ω). It givesthe hole’s local density of states (LDOS):ρ(r, ω) = − 1π Im〈0|hrGˆ(ω)h†r|0〉 (4.21)= − 1π ImG0,0(ω),which is same as the total density of states in the clean system. Othercomponents of the GF, G0,R 6=0(ω), can then be calculated from G0,0(ω)using these recursive relations: using Ω2, V3 is written in terms of V1 andV2. This is used in Eq. (4.18) which is then solved for V1 and V2, and so on.In practice, the calculation is done for a finite lattice which is chosensufficiently large that the GFs become negligible beyond its boundaries (thebroadening η introduces an effective lifetime 1/η that prevents the quasi-particle from going arbitrarily far away from its original location). Notethat the equations are modified for the lattice sites close to the boundary:if the hole can not hop outside the boundary, some of the generalized GFsF1, F2, F3 must be set to zero for sites close to the boundary, ensuring thatthe hole does not go beyond the border. This modifies effective hoppingsand on-site energy near the boundary. If the cutoff is large enough, however,the solution becomes insensitive to these details.The top panels in Fig. 4.4 show the hole’s total density of states (DOS)for two moderate values of t/J , for which this variational approximationwas shown to be in good agreement with the numerical results [2]. Thequasiparticle bandwidth for t = 6 is considerably larger than that for t = 3,showing the rapid decrease of the quasiparticle’s effective mass with increas-ing hopping. In the lower panels, we plot its effective hoppings t1(ω), t2(ω)and on-site energy ε(ω) at the same energy range. It shows that their energydependence is relatively weak and that t2(ω), which would make the DOSasymmetric, is vanishingly small. This explains why the quasiparticle, inspite of being dressed with magnons, has a DOS similar to that of a feature-less bare particle with only a renormalized first nn hopping. However, itsweak effective hopping (compared to t) and low quasiparticle weight, Fig.(4.5), suggests that the polaron is actually quite different from the bare holeas it is dressed by magnons.654.4. The effect of the impurity0.01.53.04.50.000.020.040.06t1t2 , x100.0000.005teff(ω)t1t2, x300102030ρ(ω)-0.12 -0.10 -0.08ω/J-0.12-0.10ε(ω)-7.0 -6.9 -6.8ω/J-7.2-7.0-6.8t = 6J = 6t = 3J = 3Figure 4.4: The total density of states (top panels), effective hoppingst1(ω), t2(ω) (middle panels) and on-site energy ε(ω) (bottom panels) in theclean system for two different values of t/J¯ . The effective parameters arerelatively constant within the energy band, explaining why the DOS has thegeneric form expected for a bare particle with nearest-neighbor hopping ona square lattice. Here J¯ = 1, η = 10−3 and t = 3 (left panels) and t = 6(right panels), respectively.4.4 The effect of the impurityIn the previous section it was confirmed that the hole’s motion in the cleansystem is described by an effective tight-binding Hamiltonian with secondand third nearest-neighbor hoppings which keep the quasiparticle on thesame sublattice at all times. In this section I investigate the effect of anattractive impurity on the spectrum of the quasiparticle. The impurity canbe on the sublattice in which the quasiparticle moves, or it can be on theother sublattice. In the former case, one expects the quasiparticle to bind tothe impurity. As mentioned in the introduction, when they are on different664.4. The effect of the impurity-0.12-0.1-0.08E qp (k)k0.65450.6580.6615Z(k)(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)kxky (pi,pi)Figure 4.5: (Top): the energy dispersion and the quasiparticle weight forthe lowest polaron band. The ground state is at k = (0, 0), while k = (pi2 , pi2 )is a saddle point. The top-right quadrant of the full Brillouin zone (BZ) andthat of the corresponding magnetic BZ (dashed green) is shown at the inset.The band folding due to AFM order leads to symmetric dispersion along(0, 0)− (π, π). The quasiparticle weight has only minor variations which canbe related to the fact that the polaron bandwidth (≈ 0.3) is considerablysmaller than the energy of a typical spin defect (few J¯) in the magnon cloud.t = 3J¯ = 3 and η = 10−3.sublattices one might naively expect the quasiparticle to remain mobile andinsensitive to the presence of impurity. However, we will see that this is notthe case.4.4.1 Quasiparticle and impurity on the same sublatticeThe translational invariance of the clean system requires the equal spreadingof the hole’s wave function over the entire lattice. This is expected to changewhen introducing an attractive impurity and, in particular, there may existlow-energy bound states where it is energetically more favourable for thehole to stay close to the impurity. This tendency can be studied using674.4. The effect of the impurity020406080100ρ(r = 0, ω)U = 0, x4U = -0.1JU = -0.25JU = -0.5JU = -J-7.4 -7.2 -7.0 -6.8ω/J-1.0-0.50.0Ueff(r = 0, ω)/J-7.050 -7.049ω/J050100150ρ(r = 0, ω)η = 10−4η = 10−3ρ0, x 30U = 0U = -0.5JU = -JU = -0.1JFigure 4.6: (Top) LDOS at the impurity site for various values of U . Thedashed line is the DOS in the clean system, times 4. At finite U , a singlebound state splits from the continuum and its binding energy increases withU . Curves are shifted vertically to help visibility. (bottom) The effectiveon-site energy at the impurity site is essentially equal to U . Parameters aret = 6, J¯ = 1 and η = 10−3.the hole’s Green’s function, G0,R(ω), where R and the impurity site, r = 0belong to the same sublattice. This can be calculated as we explained inthe previous section, while keeping track of the position of hole with respectto the impurity in order to include the energy gain U whenever they meet.This modifies some of the equations of motion. For example, Eq. (4.7) nowreads asG0,R(ω) = g0(ω,R)[δ0,R − t∑uF1(R,u, ω)], (4.22)684.4. The effect of the impuritywhere g0(ω,R) = 1/(ω+ iη+UδR,0− 4J¯). The coefficients in the equationsof motion for F2 also become position-dependent, reflecting the possibilityof the hole being at the impurity site. The equations for F1 and F3, forwhich the hole is on the sublattice without the impurity, remain the same astheir counterparts in the clean system. Tracking these changed coefficientsand their effects on the effective hoppings and on-site energies, we now findG0,R(ω) = g˜0(ω,R)[δR,0 −∑δt˜1(R, δ, ω)G0,R+δ(ω)−∑ξt˜2(R, ξ, ω)G0,R+ξ(ω)], (4.23)which is similar to Eq. (4.12), but now t˜1 and t˜2 depend both on the locationand on the direction of hopping, if R has the impurity within the range ofits second or third nn. If R is further away, the effective parameters takethe same values as in the clean system.Equation (4.23) can be solved similar to Eq. (4.12), that is, by groupingGFs according to their Manhattan distance. Because the problem is stillsymmetric about the origin where the impurity is located, G0,R(ω) continuesto have the same symmetries as in the clean system. Therefore, only theGFs corresponding to n ≥ m ≥ 0 need to be calculated.Given the almost constant values of ε, t1 and t2 in this range of energiesand the fact that the problem is two dimensional, bound states are expectedto appear for any finite U . The top panel in Fig. 4.6 shows the LDOSat the impurity site r = 0 for various values of the on-site attraction U .The peaks that appear below the DOS of the clean system (shown by thedashed line) are proportional to Dirac delta functions which are broadenedinto Lorentzians by the finite η. They signal the appearance of quasiparticlebound states, characterized by exponential decay of the quasiparticle’s wavefunction ψb(r) away from the impurity. The inset verifies that this is trueeven for the smallest U : the height of the ”shoulder”-like feature appearingat the bottom of the band in the main figure scales like 1/η and evolves intoa separate Lorentzian for small enough η, showing the presence of a boundstate below the continuum.The bottom panel shows the effective attraction at the impurity site de-fined as the difference between the effective on-site potential, Eq. (4.13), atr = 0 and that of sites far away from the impurity (or in the clean system),Ueff(r = 0, ω) = Re [ε(r = 0, ω) − ε(ω)]. Not surprisingly, Ueff(r = 0, ω) ≈U , although a small energy dependence, associated with the hole’s fluctua-tion to nearby sites, is observed if the scale is significantly expanded.694.4. The effect of the impurity0.0 0.5 1.0 1.5 2.0r/√20.00.20.40.60.81.0|ψ b(r)/ψb(0)|U = -0.25JU = -0.50JU = -1.00Jr = x + y r = 2y r = 2x+2yFigure 4.7: Relative amplitude of the wave functions corresponding to threeof the bound states shown in Fig. 4.6, at various distances from the impuritysite. Lines are exponential fits. States with bigger binding energies haveshorter decay lengths.The exponential decay of the quasiparticle’s wave function away fromthe impurity can be verified explicitly by calculating the amplitude of thesebound states at various distances from the impurity, r. This can be doneby approximating G0,r(ω), at the bound state energy Eb, to the term in theLehmann representation which is dominated by the bound state,G0,r(ω = Eb) = limη→0∑ψψ(0)∗ψ(r)ω − Eψ + iη|ω=Eb ≈1iηψb(0)ψb(r)∗. (4.24)Figure 4.7 shows the ratio |ψb(r)/ψb(0)| = |G0,r(Eb)/G0,0(Eb)|. The dotsare the numerical values, while the lines are exponential fits. Those cor-responding to larger binding energies (more negative Eb) are more tightlybound to the impurity and therefore decay faster, as expected. This agreeswith the larger quasiparticle weight of these states at r = 0, (see Fig. 4.6).704.4. The effect of the impurity020406080ρ(r = x, ω)U = 0, x4U = -0.5JU = -JU = -2JU = -4J-7.4 -7.2 -7.0 -6.8ω/J-1.0-0.50.0Ueff(r = x, ω)/JU = 0U = -JU = -2JU = -3JU = -4JFigure 4.8: (Top) LDOS ρ(r = x;ω) with curves shifted vertically to helpvisibility; (Bottom) Ueff(r = x;ω) at the quasiparticle’s sublattice site lo-cated nearest to the impurity. Up to three bound states split from thecontinuum upon increasing U . The presence of the impurity at r = 0 in-duces a finite effective on-site attraction at r = x, whose value is significantlysmaller than U (Bottom). Parameters are t = 6, J¯ = 1, and η = 10−3.4.4.2 Quasiparticle and impurity on different sublatticesWe now investigate the more interesting case with the impurity and thequasiparticle located on different sublattices. To this end, we constructthe Green’s function Gx,R(ω) in which x and R are on the quasiparticle’ssublattice (the rotated frame XY is centered to the right of the impurity, seeFig. 4.3(b)). In particular, we are interested in the LDOS on this sublatticeclosest to the impurity ρ(r = x, ω) = − 1pi ImGx,x(ω), where the amplitude ofpossible bound states is strongest.The equations of motion for Gx,R(ω) are derived as before, however here714.4. The effect of the impuritythe equations for F1 and F3 are modified by the presence of the impurityif R is close enough to it. This leads to equations of motion for Gx,R(ω)that are similar to those in Eq. (4.23), but with different values for theeffective hoppings and on-site energies close to the impurity. We solve theseequations using the same method, but note that the number of distinctGFs is considerably higher compared to the previous case due to the lowersymmetry here.In Fig. 4.8, we plot ρ(r = x, ω) for various values of U . The appearanceof Dirac delta peaks shows that bound states exist in this case as well. Acomparison with ρ(r = 0, ω) in Fig. 4.6 for the same value of U showsthat these peaks have different energies. Therefore, they are distinct boundstates and we are not just looking at the tails of the same states. This isfurther confirmed by the fact that up to three bound states appear here forsufficiently large U , as opposed to only one when the quasiparticle and theimpurity were on the same sublattice.1.0 1.5 2.0 2.5 3.0 3.5 4.0r0.000.250.500.751.00|ψ b(r)/ψb(r = x)|upper peak (d-wave)lower peak (s-wave) r = 2x + y r = 3x r = 3x + 2yU = -2JFigure 4.9: Relative amplitude of the upper and lower bound states forU = −2J¯ , at various distances from the impurity site. Lines are exponentialfits.These bound states exist in spite of the fact that the impurity is not lo-724.4. The effect of the impuritycated on the sublattice to which the quasiparticle belongs. As noted above,within a naive picture one does not expect this: the quasiparticle shouldnot be trapped by an on-site impurity located on the other sublattice. Thisshows that the quasiparticle does not interact with the bare disorder, butwith a renormalized one which extends to the quasiparticle’s sublattice, inspite of the fact that the bare impurity is local. This is reminiscent of therenormalization of the disorder potential by electron-phonon coupling whichwe encountered in the last two chapters. It comes about because the quasi-particle’s effective motion on one sublattice is made possible via hoppingof the hole through the other sublattice, where the hole and impurity caninteract. We define the effective on-site attraction asUeff(x, ω) = Re[ε(x, ω) − ε(ω)]which again compares the effective on-site energy near the impurity to thatof sites far away from the impurity (or the clean system). This quantity isplotted in the lower panel of Fig. 4.8 for various values of U . It is finiteeven though there is no bare disorder at r = x. However, Ueff(x, ω) ismuch weaker than U since it is an indirect effect and this explains why thebinding energies for these peaks are much smaller than in the previous case.Retardation (dependence on ω) is now clearly visible, especially for the largervalues of U . It is caused by the magnons accompanying the hole: in orderto interact with the impurity, the hole must hop onto the other sublattice,however its ability to do so depends on the structure of the surroundingcloud of magnons. At low energies, the probability for the hole to visit theimpurity is further suppressed due to the energy cost associated with themagnons generated during hopping. This explains why Ueff becomes weakerat these energies. A similar effect has been predicted for hole-doped CuOladders with non-magnetic impurities that affect the propagating holes evenif they do not lie in their path [87].Perturbation theory to zero order in t suggests that there should be afinite threshold for U in order for bound states to appear. It can be estimatedby comparing the hole’s energy at any other site in the lattice, 4J¯ + O(t2),to its minimum energy at the impurity site, 10J¯ −U +O(t2) (the increasedenergy is due to the presence of at least one spin defect). If U < 6J¯ , thisimplies that it should not be energetically favorable for the hole to be at theimpurity site. Including corrections proportional to t2 does not change thisand a finite threshold value is still predicted. However, we do not find anysuch threshold in the full calculation. This emphasizes again the importanceof the (higher-order) loop processes in describing the actual behavior.734.4. The effect of the impurity0 2 4 6 8t/ J0.000.030.06E b/Jdifferent sublattice0.000.250.500.751.00E b/Jsame sublatticeFigure 4.10: Binding energy of the s-wave bound state at U = −J¯ vs. t/J¯ ,when the quasiparticle and the impurity are on the same sublattice (toppanel) and different sublattices (bottom panel). The smaller binding energyat strong hopping is due to the reduction in the quasiparticle’s effective mass,which makes it harder to trap. When the quasiparticle and the impurity areon different sublattices, the enhancement of Ueff at small t dominates overthe effective mass decrease, explaining the growth of the binding energyhere.Fig. 4.8 shows that a total of three bound states emerge upon increasingthe impurity attraction U . Further increase of U increases their bindingenergy, but it does not give rise to more bound states. One can identifythe nature of each bound state by comparing its amplitudes on the fourneighboring sites of the impurity, 〈r = u|ψb〉. These are extracted fromGx,u(ω = Eb), just as we did in Eq. (4.24). For the lower peak, we findthe same value of 〈u|ψ1b 〉 for all u, implying s-wave symmetry. A statewith s-wave symmetry is expected to have the strongest binding to the744.4. The effect of the impurityimpurity since, to the leading order in hopping, all of its four segments meetconstructively on the impurity. For the upper peak, 〈x|ψ3b 〉 = −〈y|ψ3b 〉 =〈−x|ψ3b 〉 = −〈−y|ψ3b 〉, i.e., this state has d-wave symmetry. The middle statehas px symmetry: 〈x|ψ2b 〉 = −〈−x|ψ2b 〉 and 〈y|ψ2b 〉 = 〈−y|ψ2b 〉 = 0. It has adegenerate twin bound state with py symmetry, which has zero amplitudeat r = x and therefore it does not appear in ρ(x, ω) ∝ ImGx,x(ω). Sincethe full lattice has rotational symmetry about the impurity, the resultingbound states are expected to mirror this symmetry as well. The spatialprofile of s- and d-wave states is presented in Fig. 4.9. It shows that theyhave very similar decay lengths, consistent with their fairly similar bindingenergies and with the fact that their corresponding peaks in Fig. 4.8 havesimilar quasiparticle weights. The px state, however, is expected to haveabout twice larger weight as it is distributed between only x and −x lobes,whereas the s and d states have equal weights in all four directions. Again,this is consistent with its spectral weight shown in Fig. 4.8.1 2 3 4 5t/J0.000.010.020.030.04| Es,d - E px|/JU = -JFigure 4.11: The gap between the px and either of the s or d states, for afixed U . Its enhancement as a function of t/J¯ reflects the rotational kineticenergy gain of the quasiparticle as it becomes lighter with increasing t.Fig. 4.10 shows the binding energy Eb of the s-wave state as a function ofhopping t, when U = −J¯ . It exhibits quite different trends in the two cases.754.5. Summary and conclusionsAs the hopping becomes stronger, the kinetic energy of the quasiparticle isincreased (its effective mass decreases). A lighter quasiparticle is harder totrap and this explains why its binding energy at fixed U gets weaker whenit is on the same sublattice with the impurity (top panel).When the quasiparticle is on the other sublattice (bottom panel), itinteracts with the impurity by virtue of Ueff which is dynamically generatedand therefore strongly depends on t. At t = 0, the hole is locked at a latticesite and is unaware of the presence of impurity, therefore Eb = 0. As t isincreased, Ueff is enhanced as the hole is able to visit the impurity, whereasthe quasiparticle’s effective mass is reduced as it gains more kinetic energy.The former tends to increase the binding energy while the latter reduces it,and it is their competition that sets the dependence of binding energy on thehopping, t. The initial growth of Eb implies that the enhancement of Ueffdominates over the reduction of effective mass at small t. However, sinceUeff is weaker than U for all t (Fig. 4.8), further increase of the hoppingmakes the hole too light to be easily trapped by Ueff and the binding energyeventually starts to decrease. While only the binding energy of the s-wavestate is shown in the lower panel of Fig. 4.10, all three peaks exist for smallt, although they are energetically very close to each other. With increasingt they move closer to and eventually merge into the continuum such that,at the highest value of t in Fig. 4.10, the s-wave state is the only existingbound state.The energy gaps between the three bound states (when all are present)are nearly identical. Figure 4.11 shows its evolution with t when U = −J¯ .Since this must be due to differences in the rotational kinetic energy, it isexpected to increase with t, as the quasiparticle’s effective mass decreases.This is indeed the observed behavior.4.5 Summary and conclusionsWe investigated the effect of a non-magnetic impurity on the motion of a holein a 2D square Ising AFM. The resulting quasiparticle, which propagates onone sublattice, is confirmed to form bound states around the impurity. Thisis true both when the hole and impurity are on the same sublattice and whenthey are on different ones. The latter occurs because of the renormalizationof the effective on-site energy which results in finite effective attraction atthe sites next to the impurity that can be visited by the quasiparticle. Thisalso explains why a total of (up to) three bound states with s, p, and dsymmetries were found in this case, as opposed to only one s-wave state in764.5. Summary and conclusionsthe case when the quasiparticle is on the same sublattice with the impurity.In this latter case, the impurity is located in the node of p and d symmetrystates, therefore such states do not see it and can not be bound to it (Inreality, Ueff at sites other than those occupied by the impurity is non-zero,but given their longer distance to the impurity site, it is not large enoughto bind new states).Bound states with s, p, and d symmetries have also been observed nearan inert vacancy in a Heisenberg AFM. However, in that case, it is thedistortion of the magnetic environment around the vacancy that binds thehole [77, 78]. Such a distortion is only possible in a Heisenberg model andcomes from a local modification of the spin fluctuations. In an Ising AFM,an inert site would have no effect on the AFM order of the other sites.Moreover, if the hole is not allowed to visit this inert impurity site, thereare no Trugman loops including it so the hole loses kinetic energy whenlocated in that neighborhood. As a result, we expect that in an Ising AFM,an inert impurity like that of Ref. [77] would repulse the hole. Bound statescould only appear if a sufficiently strong exchange was turned on betweenthe hole and the inert spin, so that the exchange energy gained through itcompensates for the loss of kinetic energy. Such a model was analyzed inRef. [78], although for the Heisenberg model it was found that bound statespersist only if this exchange with the inert site is rather weak. All thesedifferences show that the underlying reasons for the appearance of boundstates are very different in these other models. This is further substantiatedby the fact that while a sublattice dependence is observed in Refs. [77, 78],it consists of a variation of the spectral weight but this is associated withthe same eigenstates. By contrast, in our model, the two sublattices showdifferent spectra of bound states.This result is important because it suggests that two very different pat-terns of bound states should be observed with scanning tunneling microscopy(STM) in such systems, even if only one type of impurity is present. Notethat we assumed that the impurity is located directly at (or above) a latticesite. If, on the other hand, the impurity was located either half-way betweentwo sites or in the center of the plaquettes (such as in the case of strontiumdoppings in La2−xSrxCuO4 superconductor), then it would not break thesymmetry between the two sublattices and only one pattern of bound statesshould appear. These cases can be studied by similar means as presentedhere.A major simplifying factor of this problem was the assumption of anIsing AFM. If spin fluctuations are turned on, in a Heisenberg AFM, amajor difference is that the hole no longer needs to go twice around closed774.5. Summary and conclusionsloops in order to become delocalized: spin fluctuations can remove pairs ofneighboring magnons, thus cutting the string short and releasing the hole.As a result, one expects a significant decrease in the effective mass of thequasiparticle, which is indeed observed [88]. However, it is interesting to notethat if there is true long-range AFM order in the plane (as is the case incuprates, due to coupling between planes), the resulting quasiparticle shouldcontinue to primarily reside on one sublattice, because spin fluctuations canonly remove/create pairs of magnons and spin conservation would continueto make the two sublattices inequivalent. This suggests that the results wepresent here, which are directly traceable to the fact that the quasiparticlelives on one sublattice, could be relevant for the Heisenberg AFM as well,although it is impossible to say a priori if the effective attraction generatedwhen the quasiparticle and impurity are on different lattices would sufficeto bind states (we would still expect s-symmetry bound states to appear ifthe quasiparticle and impurity are on the same sublattice).In the broader context, these results confirm the view that coupling tobosonic degrees of freedom renormalizes not just a quasiparticle’s disper-sion, but also the effective disorder it sees. If the latter were not the case,no bound states could arise when the quasiparticle lives on a difference sub-lattice than the impurity. Similar situation for the Holstein polarons hasbeen thoroughly discussed in the previous two chapters, where coupling tolattice phonons results in a non-trivial renormalization of the disorder po-tential seen by the quasiparticle.78Chapter 5Accurate variational solutionfor a hole doped in a CuO2layer5.1 IntroductionAfter more than two decades from their discovery, cuprate high-Tc supercon-ductors have eluded a comprehensive explanation. All these layered materi-als contain CuO2 layers which, in the absence of doping, are antiferromag-netic Mott insulators. It is now widely believed that the superconductivityin cuprates arises when the additional charge carriers are doped into theselayers [69, 70]. The hole-doped superconductivity is the more interestingone: it leads to a robust superconductivity that extends over a wide rangeof dopings, it has the highest Tc at optimal doping, while on the underdopedside a pseudo-gap phase clearly shows up at higher temperatures, with in-triguing properties. The first step towards deciphering the mechanism ofsuperconductivity in these compounds is a proper description of the motionof a hole in CuO2 layer [64–68]. This elucidates the nature of quasiparticlesthat, at finite doping concentration, will bind together and condense intothe superconducting state.Modern theoretical studies of this problem have been generally basedon a three-band p − d model [89] whose unit cell includes two O 2px/y anda Cu 3dx2−y2 orbital, Fig. 5.5(a). Because of reasons that are outlined inthe next section, the doping hole tends to mainly reside on O’s and onlyvisits the Cu sites by doing virtual hops. At low energies, however, it wasproposed that this hole is locked into a singlet state with Cu holes. It canreduce its energy by making virtual hops onto the Cu’s provided that it isin a certain linear combination of surrounding O 2p orbitals that respectsthe 3dx2−y2 symmetry of the central Cu orbital. This effectively reducesthe multi-band problem to a one-band t-J model for the hole-Cu singlet,also known as the Zhang-Rice singlet [36, 90]. We studied the Ising limit795.2. Three-band model for a cuprate layer and its t-J analogueof this one-band model in chapter 4. The full problem, including the spinfluctuations, can only be solved numerically for finite size systems and thiscounteracts the apparent simplification associated with the smaller Hilbertspace of the effective one-band treatment [91–96].In this chapter I return to the multi-band model, and devise a variationalscheme for calculating the spectral properties of a single hole doped into aCuO2 layer. This will be based on an effective model introduced in Ref. [3],which was solved via exact diagonalization for finite size clusters. Instead,the variational principle, with some approximations, enables studying aninfinite system and is therefore free of finite-size effects.Sec. 5.2 summarizes the electronic configuration and discusses the multi-band model of the hole-doped CuO2 layer. The variational solution of thatmodel is then constructed in Sec 5.3. The quasiparticle dispersion and ef-fect of various terms in the Hamiltonian on it, along with the quasiparticleweights are studied in Sec. 5.4 and 5.5. A comparison of our polaron solutionand the Zhang-Rice singlet is the subject of Sec. 5.6. Finally, a five-bandgeneralization of this model which includes both in-plane 2p orbitals, andits effect on the quasiparticle’s spectral properties, are discussed in Sec. 5.7.5.2 Three-band model for a cuprate layer and itst-J analogueThe first step in describing the motion of a hole doped into a CuO2 layeris to understand the electronic configuration of the undoped system. Inthis compound, the copper cations and oxygen anions are arranged into twointertwined square lattices, Fig. 5.5(a), and have the following electronicconfigurationCu2+ : 1s2 2s22p6 3s23p63d9,O2− : 1s2 2s22p6.O’s have therefore a closed shell, while there is a single unpaired 3d-electron(or a hole) on each Cu ion. Within the three-dimensional crystal structure,coppers are located at the center of an octahedral arrangement of six O ions,with the two apical ones belonging to nearby layers, Fig 5.1(a). The energylevels of Cu in this geometry are shown in Fig. 5.1(b). The crystal field ofthe surrounding O ligands splits five Cu 3d orbitals into a lower t2g groupincluding 3dxy, 3dyz and 3dxz orbitals, and a higher eg group of 3dx2−y2and 3d3z2−r2 orbitals that directly point towards the ligands. With nine 3d-electrons, this results in an energy degenerate configuration for Cu where805.2. Three-band model for a cuprate layer and its t-J analoguethree electrons are left to occupy four eg states. In fact, the octahedronundergoes a tetragonal distortion in order to reduce its energy by removingthis degeneracy, with the in-plane ligands moving towards the Cu. Thisopens up a gap between the two eg (and within the three t2g) orbitals [97]and makes 3dx2−y2 the only partially occupied 3d orbital, Fig. 5.1(b).Figure 5.1: (a): Cu, shown in red, is surrounded by six O ligands in octa-hedral geometry. The static electric field of these ligands partially removesthe degeneracy of Cu 3d orbitals into eg and t2g manifolds. This degener-ate configuration distorts to remove the remaining degeneracies and endsup having 3dx2−y2 as its highest partially-occupied orbital, (b). The valuesof splittings for La2CuO4 are given. In the case of Sr2CuO2Cl2, splittingwithin eg is such strong that 3d3z2−r2 orbital is pushed below the entire t2gmanifold.Due to strong hybridization between 2px,y and 3dx2−y2 orbitals, thereis a considerable hopping amplitude, tpd, between neighbouring O’s andCu’s. The electronic structure of O and Cu ions in the CuO2 layer is such815.2. Three-band model for a cuprate layer and its t-J analoguethat the energy of electrons in doubly-occupied 2p orbitals of O’s is higherby Epd = Ep − Ed > 0 compared to that in the singly-occupied 3dx2−y2orbitals4, Fig. 5.2(a). Therefore, the lowest-energy charge excitation in thelayer happens between O’s and Cu’s and costs what is known as the charge-transfer energy ∆pd = Udd − Epd, where Udd is the Hubbard interactionin 3dx2−y2 orbital. Since ∆pd > 0, charge fluctuations are suppressed atlow energies and the system becomes an electrical insulator, also known ascharge-transfer insulator [98] to distinguish them from Mott insulators wherethe lowest charge excitation happens between Cu’s themselves (Epd < 0),Fig. 5.2(b). Nevertheless, virtual hopping of the Cu holes onto O’s gives riseto antiferromagnetic superexchange interaction between nearest-neighbour(nn) Cu spinsUddEpd(a) (b)UddCu O Cu O Cu O Cu Cu O CuFigure 5.2: (a) The ground state of charge-transfer insulator (left) and itslowest excitation (right) ∆pd higher in energy. ↑ and ↓ represent electronsand © stands for a hole (electron representation). (b) In a Mott insulator,however, charge excitation with the lowest energy corresponds to an electronmoving from one Cu to another at an energy cost of Udd.HJdd = Jdd∑<i,j>Si.Sj , (5.1)in which Si is the spin-12 operator for the hole at Cu site i and Jdd is thestrength of the superexchange coupling. It is a fourth-order process in tpdas illustrated in Fig. 5.3,Jdd ∼t4pd∆2pd(2∆pd + Upp),4The difference in the Madelung potential of Cu and O sites is another factor con-tributing to this.825.2. Three-band model for a cuprate layer and its t-J analogueFigure 5.3: One set of processes giving rise to the superexchange coupling,Jdd, between two neighbouring Cu holes. Arrows represent the holes andempty O levels corresponds to them being doubly-occupied with electrons(hole notation). The energies of intermediate configurations are given.where Upp is the Hubbard repulsion between two Cu holes on the same O2p orbital5.Note that the ground state of this Hamiltonian is not the simple clas-sical Ne´el-ordered state. The quantum fluctuations, Jdd/2(S−i S+j + S+i S−j ),smear out the Ne´el state by flipping pairs of spins, giving rise to spin fluctu-ations that are sometimes modelled using linear approximations [99]. Exacttreatment is only possible numerically for finite size systems [71]. Here, wemake a key simplification of replacing the full Heisenberg Hamiltonian withits Ising counterpart whose ground state, |N〉 = | ↑, ↓, ↑, ...〉, is Ne´el-orderedand has no spin fluctuationsHJdd = Jdd∑<i,j>Szi Szj . (5.2)This will enable us to set up a semi-analytic calculation of a doped hole’s5Note that there is a weak correction to this Jdd (∝ 1/Udd) due to processes involvingdouble-occupancy of Cu, which are neglected in this work.835.2. Three-band model for a cuprate layer and its t-J analoguedynamics in the layer for an essentially infinite system. We will justify thevalidity of this approximation based on the results it leads to.Figure 5.4: Energy diagram of the undoped CuO2 layer, illustrating acharge-transfer insulator. The filled bands are shaded, above them lies theFermi energy in the gap. Since ∆pd < Udd, removing electrons created holeson O’s rather than on Cu’s and this moves the Fermi energy into the O2p-band.5.2.1 Doping the CuO2 layer with a holeWhereas the CuO2 layer is an insulator at half-filling, it becomes conductingupon doping which moves the Fermi energy into either the upper Hubbardband (electron doping - inverse photoemission) or the O 2p-band (hole dop-ing - photoemission), Fig. 5.4. This picture is supported by various spec-troscopic studies [100] suggesting that holes primarily reside on O’s ratherthan on Cu’s. By doing so, the system can avoid the large Hubbard repul-sion Udd (> Upp > ∆pd) between two holes on the same Cu 3dx2−y2 orbital.Assuming Udd →∞ and projecting out configurations involving charge fluc-tuation between Cu’s and O’s (of order ∆pd), one arrives at the following845.2. Three-band model for a cuprate layer and its t-J analogueeffective Hamiltonian for the CuO2 layer with a single hole [3]H = Jdd∑<i,j>′Szi Szj + Tˆpp + HJpd + Tˆswap≡ HJdd + H1, (5.3)where the prime indicates the lack of exchange coupling between the pair ofcoppers that have the doping hole on the bridging O and the spin fluctuationsof Cu holes are neglected. Note that in reality the first term has to be theHeisenberg interaction, but we have applied the simplification suggestedbefore in Eq. (5.2). Tˆpp is the tight-binding Hamiltonian with the nearest-neighbour (nn), tpp, and Cu-bridged next nn hopping, t′pp, and describes thekinetic energy of the hole as it propagates among O 2px/y orbitals, whoserelative phases determine the sign of these hoppingsTˆpp = tpp∑R,δ,σrδp†R,σpR+δ,σ − t′pp∑R,σp†R,σ(pR−ε,σ + pR+ε,σ). (5.4)Here, p†R,σ creates a hole with spin σ on the O ligand 2p orbital locatedat R. δ and ε are the nn and next nn vectors. rδ = ±1 sets the signsof nn hoppings in accordance with the phases of 2p orbitals involved, seeFig 5.5(a). The next term, HJpd , describes the antiferromagnetic exchangecoupling between the doping hole, sR, and its two neighbouring coppers atsite i:HJpd = Jpd∑<i,R>sR.Si . (5.5)Unlike HJdd , we treat this term as a full Heisenberg Hamiltonian. Finally,Tˆswap is an additional effective kinetic energy term which describes the mo-tion of the hole between O sites that is mediated by the transfer of a Cuhole, without involving double-occupancy of Cu. Fig. 5.5(c) shows how thiscan give rise to nn and next nn hopping between O, with the same effectivehopping ts ∼ t2pd/∆pd: first, the Cu hole hops onto one of its empty neigh-bouring O orbitals, then the doping hole fills the empty Cu orbital. Notethat this results in swapping the spin of the two holes. The explicit form ofthis Hamiltonian is then as followsTˆswap = −ts∑R,ηsηp†R,σpR+η,σ′ |σ′〉〈σ|, (5.6)where |σ〉 refers to the state of Cu spins and prime reflects the possibilityof spin flip during the hop; sη = ±1 determines the sign of each nn or next855.2. Three-band model for a cuprate layer and its t-J analogueOCuOCu O1212xytt tkktt(a)(b)(c1)(c2)(d)(t <0)(t >0)a = 1in.fin.pd pdpdpdsssxyyyxxFigure 5.5: (a): The CuO2 layer with the spin-up hole occupying one of theO 2p orbitals. nn (δ) and next nn hops (ǫ) between O’s that are consideredin our model are shown. For the hole, nn hopping is negative along y = xdirection where 2p orbitals have their lobes with similar phases pointingtowards each other and positive along y = −x. (b): Cu and O ions thatmake the unit cell are shown in green. (c): Processes that give rise tonn (c1) and next nn (c2) spin-swap. (d): The full Brillouin zone of theCuO2 lattice (the outer square) which encloses the magnetic Brillouin zone(shaded). We study the quasiparticles with momenta lying along the dashedred line, starting from (0, π) and ending at (π, 0).nn hop according to the sign of 3d and 2p orbitals involved, Fig. 5.5(c). If865.3. Green’s functionsthey had opposite spin directions, swapping flips the Cu spin and creates alocalized magnon (a spin defect). Note that Jpd can also generate a magnon,but it differs from Tˆswap in that it does not involve the motion of the O hole.Hamiltonian (5.3) is reminiscent of the strong coupling limit of Hubbardmodel, namely the t-J Hamiltonian. The values of hopping parametersare determined via tight-biding fittings to LDA band structure calculations,whereas the Hubbard repulsions are estimated from cluster calculations com-pared with various spectroscopic measurements [101]. As a result, thereare no free fitting parameters. Given in units of Jdd (≈ 0.10 − 0.15 eV),the dimensionless parameters are tpp = 4.13, ts = 2.98, Jpd = 2.83 andt′pp = 0.58tpp. For complete technical details regarding the derivation ofHamiltonian (5.3) and further discussions, the reader is referred to BayoLau’s thesis [102].5.3 Green’s functionsIn this section, we calculate various Green’s functions (GF) describing thepropagation of the hole between O sites at zero temperature. The unit cellis shown in Fig. 5.5(b). Because of the underlying Ne´el order, it containstwo distinct Cu atoms (spin-up, spin-down) and therefore four different O.For each O, it is common to take into account only the ligand 2p orbitalthat hybridizes with the 3dx2−y2 orbitals of Cu and this is what we also doin this section. Although this proves to be sufficient as far as the dispersionrelations are concerned, we later extend our analysis to include the otherin-plane 2p orbital when we study the hole’s quasiparticle weight.The single-hole Green’s functions are defined asGαβ(k, ω) = 〈N|hk,α,↑Gˆ(ω)h†k,β,↑|N〉,where Gˆ(ω) = limη→0+ 1/(ω −H + iη) is the resolvent of the Hamiltonian.h†k,α,↑ creates a spin-up hole in Bloch states made out of O 2pα orbitalsh†k,α,↑ =1√N∑Rαeik.Rαp†Rα,↑, (5.7)in which N is the number of unit cells, p†Rα,↑ creates a spin-up hole at O2pα orbital located at Rα and α = px1 , px2 , py1 , py2 are the distinct O 2porbitals in the unit cell, see Fig. 5.5(b). The quasi-momentum k belongs tothe magnetic Brillouin zone (MBZ), the shaded area in Fig. 5.5(d), which875.3. Green’s functionsincludes all unique momenta in the reciprocal lattice. These sublattice GFsare useful since their poles correspond to the hole’s energy dispersion, En(k),and the associated residues give the overlap between the true eigenstates andthe single-particle Bloch states, Z(k) = |〈N|hk,α,↑|ψn,k〉|2.5.3.1 One-magnon approximationThe equations of motion for Gαβ(k, ω) can be written using the Dyson iden-tity based on the separation suggested in Eq. (5.3),Gˆ(ω) = Gˆdd(ω) + Gˆ(ω)H1Gˆdd(ω). (5.8)With four different p orbitals, there are a total of 16 different Green’s func-tion. For β = px1 , we haveGα,px1 (k, ω) = g0(ω)[δα,px1 + 〈k, α, ↑|Gˆ(ω)H1|k, px1 , ↑〉], (5.9)where g0(ω) = 1/(ω−Jdd/4+iη), Jdd/4 is energy cost of breaking a superex-change bond by the doping hole and |k, px1 , ↑〉 = h†k,px1 ,↑|N〉. H1 takes thehole from the 2px1 orbital to other 2p orbitals at its vicinity, during whichits spin can also be flipped and a magnon is createdH1p†Rpx1 ,↑|N〉 = ts[−p†Rpx1+x/2+y/2,↓|S′〉 + p†Rpx1+x/2−y/2,↓|S′〉−p†Rpx1+x,↓|S′〉] + Jpd2 p†Rpx1 ,↓|S′〉+ts[p†Rpx1−x/2+y/2,↑|N〉 − p†Rpx1−x/2−y/2,↑|N〉−p†Rpx1−x,↑|N〉] + tpp[p†Rpx1+x/2+y/2,↑|N〉+p†Rpx1−x/2−y/2,↑|N〉 − p†Rpx1−x/2+y/2,↑|N〉−p†Rpx1+x/2−y/2,↑|N〉] − t′pp[p†Rpx1+x,↑|N〉+p†Rpx1−x,↑|N〉], (5.10)where the distance between neighbouring copper atoms is taken to be unity(x and y are unit vectors) and s.S = szSz + (s+S− + s−S+)/2 is used. |S′〉is the state of Cu spins where the spin-down Cu next to the hole is flipped:885.3. Green’s functions|S′〉 = S+Rpx1+x/2|N〉. Using this in the last matrix element in Eq. (5.9) gives〈k, α, ↑ |Gˆ(ω)H1|k, px1 , ↑〉 = ts[−V 31 + V 11 − V 21 ] +Jpd2 V01 + ts[ei(kx−ky)/2×Gα,py1 (k, ω) − ei(kx+ky)/2Gα,py2 (k, ω)−eikxGα,px2 (k, ω)] + 2tpp[cos((kx + ky)/2)×Gα,py2 (k, ω) − cos((kx − ky)/2)Gα,py1 (k, ω)]−2t′pp cos(kx)Gα,px2 (k, ω). (5.11)V lM are called one-magnon Green’s functions and they describe processesFigure 5.6: Numbering various one-magnon Green’s function VMi at M = 1,3: the one with the hole lying to the far left is i = 0 and the index i increasesin the counterclockwise direction. The current configuration corresponds toV 21 .where the spin-down hole starts at any of 4M O sites (l = 0, 1, ..., 4M − 1)that are located at Manhattan distance of M from the magnon (M = 1here). The numbering scheme is shown in Fig. 5.6. The explicit expressioncorresponding to V 31 , for example, is the followingV 31 =1√N∑Rpx1eik.Rpx1 〈N|hk,α↑Gˆ(ω)p†Rpx1+x/2+y/2,↓S+Rpx1+x/2|N〉. (5.12)895.3. Green’s functionsTo generate the equation of motion for these propagators, we again use theDyson’s identity to findV 31 =g1(ω)√N∑Rpx1eik.Rpx1 〈N|hk,α↑Gˆ(ω)H1p†Rpx1+x/2+y/2,↓S+Rpx1+x/2|N〉,(5.13)where g1(ω) = 1/(ω − 7Jdd/4 + iη) and 7Jdd/4 is the superexchange energyof state with a magnon next to the doping hole.The swap and Cu-hole exchange interaction parts of H1 can annihilatethe magnon and flip the hole’s spin back to its initial direction. This relatesV1 to various Gα,γ(k, ω). They can also link V1 to new GFs W with twomagnons, which are then related to GFs with three magnons and so on.However, creating each additional magnon increases the energy of that stateby about 2Jdd as up to four Cu-Cu bonds are disturbed. Such many-magnonstates are therefore energetically expensive and cannot have a significant con-tribution to the hole’s low-energy eigenstates. Based on this, we constructa variational approximation in which we only keep states with up to a cer-tain maximum number of magnons, nm, and discard those including moremagnons. Doing so at the one-magnon level, nm = 1, gives the equation ofmotion for V1 as followsV 01 /g1(ω) = ts[−e−i(kx+ky)/2Gα,py2(k, ω) + ei(ky−kx)/2Gα,py1(k, ω)−e−ikxGα,px2(k, ω)] +Jpd2 [Gα,px1(k, ω) − V01 ]+tpp[V 13 + V 31 − V 11 − V 113 ] − t′pp[V 03 + V 21 ], (5.14)where several components of V3 are introduced. Similar equations can bederived for the other three components of V1V 11 /g1(ω) = ts[Gα,px1 (k, ω) − e−ikxGα,px2 (k, ω) − e−i(kx+ky)/2Gα,py2 (k, ω)]+Jpd2 [ei(ky−kx)/2Gα,py1 (k, ω) − V11 ] + tpp[V 23 + V 21 − V 01 − V 43 ]−t′pp[V 33 + V 31 ], (5.15)V 21 /g1(ω) = ts[−Gα,px1(k, ω) + e−i(kx+ky)/2Gα,py2(k, ω)−ei(ky−kx)/2Gα,py1(k, ω)] +Jpd2 [e−ikxGα,px2(k, ω) − V 21 ]+tpp[V 11 + V 73 − V 31 − V 53 ] − t′pp[V 01 + V 63 ], (5.16)905.3. Green’s functionsV 31 /g1(ω) = ts[e−ikxGα,px2 (k, ω) −Gα,px1 (k, ω)−ei(ky−kx)/2Gα,py1 (k, ω)] +Jpd2 [e−i(kx+ky)/2×Gα,py2 (k, ω) − V31 ] + tpp[V 01 + V 83 − V 21−V 103 ] − t′pp[V 11 + V 93 ]. (5.17)One can then generate the equation of motion for V il , l ≥ 3, in a similarway. Without allowing more magnons to be created, V il ’s are linked to otherone-magnon GFs at Manhattan distances of M = l − 2, l, l + 2. Note thatGFs with no magnons, Gα,β(k, ω), do not arise here as the hole is not animmediate neighbour of the magnon and cannot annihilate it by the actionof H1. The exact form of these equations depends on the 2p orbital fromwhich the hole starts and it is the state of nearby Cu spins that determineswhich Tswap hoppings are allowed. For example, for V 43 we haveV 43g2(ω)= tpp[V 53 + V 33 − V 11 − V 75 ] − t′pp[V 23 + V 85 ]+ts[V 75 − V 53 − V 85 ], (5.18)where g2(ω) = 1/(ω − 9Jdd/4 + iη) and 9Jdd/4 is the superexchange energyof Cu holes with one magnon which is not an immediate neighbour of thedoping hole. Eq. (5.18) is similar to the two-component recursive relationswe encountered in the last chapter and can be solved similarly using themethod of continued fractions. By collecting all one-magnon GFs at Man-hattan distance l into vectors Vl, the set of corresponding 4l equations canbe expressed as a single-indexed recursive relation as followsλlVl = ξlVl−2 + ζlVl+2, l ≥ 3 (5.19)where λl, ξl and ζl are sparse matrices. Because this relates three consecutiveterms, its solution can be expressed as a continued fraction. Assuming asolution of the form Vl = AlVl−2 and using it in Eq. (5.19) results inAl = (λl − ζlAl+2)−1ξl. (5.20)This can be solved by setting Al to zero beyond a large enough cutoff c,Ac+2 = 0. In particular, it gives A3 which relates V3 to V1. This is thenused in equations for V1, Eq. (5.14), to reduce them into equations relatingV i1 to Gα,β(k, ω) onlyPV1 = ΛGα, (5.21)915.3. Green’s functionswhere V1 =V 01V 11V 21V 31, Gα =Gα,px1Gα,px2Gα,py1Gα,py2and P and Λ are some 4×4 matrices.Using this in Eq. (5.11) we can eliminate its one-magnon GFs in favour ofvarious Gα,β(k, ω). This gives the matrix element 〈k, α, ↑ |Gˆ(ω)H1|k, px1 , ↑〉as a linear combination of Gα,β(k, ω):〈k, α, ↑ |Gˆ(ω)H1|k, px1 , ↑〉 =∑βUpx1 ,βGα,β(k, ω),where Upx1 ,β is the first row of a matrix of coefficients, U . Therefore, ouroriginal equation of motion for Gα,px1 (k, ω), Eq. (5.9), takes the followingform∑βUpx1 ,βGα,β(k, ω) −1g0(ω)Gα,px1 (k, ω) = −δα,px1 , (5.22)where α = px1 , px2 , py1 , py2 . This is a set of four linear equations betweensixteen Green’s function, Gα,β(ω). The other twelve equations are givensimilarly by the hole starting from 2px2 , 2py1 and 2py2 orbitals. These arethen trivially solved to give all Gα,β(ω).Figure 5.7: Numbering various two-magnon GFs at M = 1, 3. The configu-ration shown corresponds to W 23 .925.3. Green’s functions5.3.2 Two-magnon approximationThe one-magnon approximation can be improved by allowing for the creationof a second magnon, i.e. working in the variational space with nm = 2.In doing so, we let the second magnon be created only next to the firstone and not further away. This is a good approximation for low-energyquasiparticles whose magnons are bound to the cloud. Also, this disturbsfewer antiferromagnetic bonds and therefore costs less energy compared toconfigurations with magnons apart. The equation of motion for V 01 , forexample, is then modified as followsV 01 /g1(ω) = ts[−e−i(kx+ky)/2Gα,py2(k, ω) + ei(ky−kx)/2Gα,py1(k, ω)−e−ikxGα,px2(k, ω)] +Jpd2 [Gα,px1(k, ω) − V01 ]+tpp[V 13 + V 31 − V 11 − V 113 ] − t′pp[V 03 + V 21 ]+ts(W 113 −W 13 −W 03 ) +Jpd2 W01 , (5.23)where W jl are two-magnon GFs in which the hole starts at Manhattan dis-tance l from the first magnon (spin-up defect). Fig. 5.7 shows their num-bering scheme. The other components of V1 are also similarly modifiedV 11 /g1(ω) = {...} + ts(W 13 −W 113 −W 03 ) +Jpd2 W01 ,V 21 /g1(ω) = {...} + ts(W 113 −W 13 −W 03 ) +Jpd2 W01 ,V 31 /g1(ω) = {...} + ts(W 13 −W 113 −W 03 ) +Jpd2 W01 , (5.24)where {...} refers to those terms existing at one-magnon approximation, Eq.(5.14).For Manhattan distances l ≥ 5, similar continued-fraction solutions existfor both V and W : V5 = A5V3 and W5 = A′5W3; the prime reflects thedifferent superexchange energy of the Cu spins with one (9Jdd/4) vs twomagnons (13Jdd/4). For l = 1 and l = 3, however, various V and W aremixed together due to proximity of the hole and magnons. These linearrelations,V3 = V3(V1,W1,W3)W1 = W1(V1,V3,W3)W3 = W3(V1,V3,W1), (5.25)935.4. Resultsare solved for W1, W3 and V3 in terms of V1. This is then used to reduceEq. (5.23) to one relating V1 and Gα,β only. The rest of the solutionproceeds similar to the one-magnon case.We note that the restriction of only keeping two-magnon configurationswith neighbouring magnons can be easily removed, and indeed this was donein the work described in Ref. [103]. For low energy quasiparticle states,however, this leads to very minor quantitative differences.5.4 ResultsIn this section, we use these Green’s functions in order to extract the disper-sion relation and quasiparticle weight of the hole propagating in the CuO2layer as a magnetic polaron. We will then gauge the importance of var-ious terms in the Hamiltonian in controlling the generic features of thesequasiparticle properties.5.4.1 Dispersion relationsHaving calculated the Green’s functions, we can extract the spectral prop-erties of the hole quasiparticle by studying the associated sublattice spectralfunctions,Aα(k, ω) = −1π ImGα,α(k, ω).We let the quasi-momentum k run over the full Brillouin zone along variouscuts shown in Fig 5.5(d). The spectral weight is plotted in Fig. 5.8 incolor scale as a function of k and energy, ω, showing coherent bands ina wide range of energies. It is obtained within nm = 2 approximation,where only the two-magnon configurations with both magnons next to eachother are retained. One expects a continuum of states starting from anenergy E1,GS + 2Jdd, i.e., the ground state of nm = 1 quasiparticle plusthe energy of an unbound magnon. This is missing from Fig. 5.8, but itappears if one allows for configurations where the second magnon is awayfrom the one-magnon polaron. This is also why the energy and location ofhigher energy features is tentative since multiple-magnon states contributeheavily to these higher energy eigenstates, and only after they are includedis their convergence expected. As a result, from now on we focus on thelowest energy eigenstates (the quasiparticle band) which have the fastestconvergence and correspond to the highest electron-removal states probed945.4. Resultsby angle-resolved photoemission spectroscopy (ARPES)6.In Fig 5.9 we compare the energy dispersion of the lowest energy bandfor nm = 1 vs that of nm = 2 approximations. These are extracted from themap plots similar to Fig. 5.8. They show essentially the same features, mostnotably a deep minimum at k = (pi2 , pi2 ) which disperses in both kx = ±kydirections, as observed in experiments [28]. For nm = 2 the quasiparticlebandwidth is considerably narrower compared to the nm = 1 case. This isindeed expected since the hole becomes an effectively heavier quasiparticlewhen it is bound to a bigger cloud of magnons. This also shows that con-figurations with two magnons have significant contribution to the polaronstates as including them changes the energy dispersion quantitatively. Theoverall shift of nm = 2 energy dispersion to lower energies is a consequenceof the bigger variational space for this level of approximation.In order to gauge the accuracy of our variational approximation, wecompare it to available exact-diagonalization results for a cluster of 32 Cuof the same model [3], symbols in Fig. 5.9. It shows that the energy disper-sion given by nm = 2 approximation is almost converged to the exact result.This implies that configurations with more that two magnons have negligiblecontribution to the low-energy polaron states. Furthermore, spin fluctua-tions seem to have almost no effect on the hole’s dispersion since they arefrozen out in our approximation, but fully included in the ED calculation. Infact, spin fluctuations can bind to the doping hole to form a spin-32 polaron.In some parts of the BZ, ED results suggest that this spin-32 polaron liesenergetically below our spin-12 polaron state [3]. However, such high-spinpolaron does not appear in our model which ignores spin fluctuations.The lack of importance of spin fluctuations in the three-band model canbe understood as follows. Since the doping hole resides primarily on O lig-ands, it can move freely on the O sublattice by Tpp without being bound tomagnons on the Cu’s. The polaron can therefore move easily in the layerwhile the hole absorbs the magnons and re-creates them at other locations.Spin fluctuations act on longer time scales (Jdd < tpp) and cannot thereforehave a considerable effect on the quasiparticle that already has a fast dy-namics controlled by Tpp. This has to be contrasted with one-band modelswhere the hole is bound to magnon cloud that it creates by moving throughthe magnetic lattice [33]. In the absence of spin fluctuations, the only wayfor the hole to propagate is by executing the so-called Trugman loops [88]where the hole spends most of its time going around closed paths in order tofree itself from the magnons. This results in a very heavy quasiparticle with6Note that ARPES energies correspond to ω → −ω, because we use a hole formulation.955.4. ResultsFigure 5.8: The spectral weight at nm = 2 approximation, Apx1 (k, ω) in colorscale as a function of momentum and energy (in units of Jdd ≈ 0.1 − 0.15eV), showing several coherent bands. For the rest of this chapter, only thelowest energy band will be studied.the ground state at k = 0 [2]. Even when the spin fluctuations are included,the dispersion given by the one-band model with only nn hopping disagreeswith experiments by being almost flat along (π, 0)-(0, π). Long range hop-ping which do not disturb the antiferromagnetic order must be introducedin order to obtain an isotropic quasiparticle dispersion [104, 105]. The ex-istence of such longer range hopping terms in one-band models is justifiedusing cluster calculations [106].We also tried to mimic spin fluctuations by introducing them locally,which links two-magnon GFs to those with no magnons, Gα,β(k, ω), andvice versa. The effect on the energy dispersions, Fig. 5.9(b), was barelyvisible in that energy scale.965.4. Resultsk-21-18-15E qp(k)nm=1nm=2ED - ST=1/2-10-8E qp(k)(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)(b)(a) nm= 0Figure 5.9: (a): Energy dispersion (in units of Jdd ≈ 0.1 − 0.15 eV) of thehole without allowing for the creation of magnons, nm = 0. The groundstate is still at k = (pi2 , pi2 ). (b): Comparison between one- and two-magnonapproximation results for the hole’s energy dispersion along various cuts inthe full BZ. The dots are relevant exact-diagonalization calculations for afinite size sample of the same model, shifted in energy in order to ease thecomparison.5.4.2 Effect of various terms in the HamiltonianNow that we have established the convergence of our approximation, a nat-ural question is whether the essential features of the quasiparticle dispersionEqp(k) are robust or they depend on the specific setting of parameters. Toanswer this question, we explore the role of various terms in the Hamiltonian(5.3) by removing them one at a time. Fig. 5.10 shows the quasiparticledispersion with various terms on and off. Setting Jpd = 0 has essentiallyno major effect and results in a similar dispersion which is only shifted inenergy as some exchange energy is lost. On the other hand, Tˆswap is thecrucial term and without it (ts = 0) the dispersion changes completely. In975.4. Resultsparticular, k = (pi2 , pi2 ) ceases to be the ground state in this case. Raising tsfrom zero to half its actual value, Fig. 5.10(d), restores the proper shape ofthe dispersion and this suggest that tpp and ts have to be comparable forthe correct shape to appear. Fig. 5.9(a) shows that, including all terms,one actually gets a deep minimum at k = (pi2 , pi2 ) even without allowingfor the creation of any magnons, nm = 0. This is very different from one-band models where the dispersion along (0, 0)-(π, π) is mainly controlled byspin fluctuations, in the absence of which the bare dispersion is completelyanisotropic [107, 108]. Note that in other studies of three-band model Tˆppis only treated as a small perturbation [109–111], whereas here it has thelargest energy scale in the Hamiltonian.k-12-10E qp(k)-12-11-10E qp(k)-18-16E qp(k)-15-14E qp(k)(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)tpp = 0Jpd = 0ts = 0ts = 0.5tsact.(a)(b)(c)(d)Figure 5.10: The lowest quasiparticle dispersion (nm = 2 and in units ofJdd ≈ 0.1 − 0.15 eV) when various terms in the Hamiltonian (5.3) are onand off. Jpd has almost no effect on the dispersion (a), whereas without ts(c) the model fails to properly predict the expected ground state k = (pi2 , pi2 ).For comparable values of ts and tpp, the correct shape of energy dispersionre-appears (d).985.5. Quasiparticle weight5.5 Quasiparticle weightIn this section, we study the polaron’s sublattice quasiparticle weights Zα(k)along various cuts in the full BZ. It corresponds to residues at the poles of thesublattice spectral weight, i.e. Aα(k, ω) = Zα(k)δ(ω −Eqp(k)). In practice,delta peaks are broadened into Lorentzians due to the finite value used forη. For ω ≈ Eqp(k), we can write:Aα(k, ω) = −1π ImGα,α(k, ω)|η =ηπZα(k)(ω − Eqp(k))2 + η2→ Zα(k) = πηAα(k, Eqp(k)), (5.26)where α is any of the four distinct O 2p orbitals. Here, the quasi-momentumk is unique only inside the magnetic BZ. The upper panel of Fig. 5.11 showsthe quasiparticle weight Zpx1 (k) at the lowest band for nm = 1 and nm = 2approximations. As it was the case for the dispersions, the weights aresimilar in both levels of approximation. They are also symmetric about theground state along (0, 0) to (π, π) cut. This is a result of invariance undertranslation by reciprocal lattice vector, k → k − (π, π), and time-reversalsymmetry, k → −k. For nm = 2, part of the quasiparticle wavefunctiondescribes its two-magnon component and this explains the overall reductionof the quasiparticle weight in this case. Note that the weights follow atrend that is opposite to that of the quasiparticle dispersion, repeated atthe lower panel of Fig. 5.11. In particular, both curves show that thequasiparticle weight is enhanced near local energy minima, although thenm = 2 approximation does this slightly better by predicting the largestweight at the ground state k = (pi2 , pi2 ).Unlike the dispersion relations, quasiparticle weights depend on thechoice of orbitals on which the Green’s functions are projected. Fig. 5.12compares the quasiparticle weights Zpx1 (k) and Zpy2 (k) which correspond to|〈k, px1 , ↑ |ψk〉|2 and |〈k, py2 , ↑ |ψk〉|2, respectively. Differences in the weightsalong the first (third) cut suggests that the associated wavefunctions |ψk〉have more py (px) character for these momenta.5.5.1 Comparison with ARPESIt was mentioned at the beginning of this chapter that the problem of asingle hole in CuO2 layer is closely related to angle-resolved photo-emissionspectroscopy (ARPES) of cuprates, where the outgoing photo-electron isanalyzed in energy and momentum in order to gain information about thedispersion of the hole left behind in the sample during photo-emission. The995.5. Quasiparticle weightk00.050.10.150.2Z px 1 (k)one-magnontwo-magnons-21-19.5-18E qp( k)(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)0Figure 5.11: top: quasiparticle weights extracted from the sublattice GFApx1 (k, ω) in one- and two-magnon approximation. The correspondingenergy dispersion (for the two-magnon approximation, in units of Jdd ≈0.1 − 0.15 eV) is presented in the lower panel.latter is what we extract from our GFs. Therefore, it would be of greatinterest to compare our calculations to results available from ARPES.When a photon with an arbitrary polarization arrives in the CuO2 layer,it does not make a distinction between different O sites and is equally likelyto remove an electron from any one of the four 2p orbitals in the unit cell.It is shown in Appendix F that, after averaging over the polarization of theincident beam, the intensity of the outcoming photo-electrons in ARPES isgiven by the following expressionI¯(K, ω) ∝ A¯arpes(k, ω) ∼∑k,GδK||+k,G∑α,βeiG.RαβηαβAαβ(k, ω), (5.27)where K = (K||,Kz) is the momentum of outgoing photo-electron and ω isthe energy transferred to it during the photo-emission. Rαβ = Rα−Rβ are1005.5. Quasiparticle weightk0.030.060.09Z p(k)px1py2(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)Figure 5.12: The quasiparticle weights in nm = 2 approximation extractedfrom spectral weights that correspond to projections onto orbitals with dif-ferent polarizations, 2px and 2py.the relative positions of various 2p orbitals in the unit cell, k belong to theMBZ, G are the corresponding reciprocal lattice vectors and Aαβ(k, ω) =− 1pi ImGαβ(k, ω). ηαβ = 1 if α and β are both either 2px or 2py orbitals andzero otherwise and it reflects a selection rule resulting from averaging overthe polarization of the incoming beam. When the coupling to spins is turnedoff, the dispersion relation of the remaining Tˆpp can be analytically calcu-lated. We checked to confirm that A¯arpes(k, ω) gives the properly unfoldeddispersion relation which agrees with the analytic result.The ARPES quasiparticle weight is given in Fig. 5.13. Note that, unlikethe sublattice spectral weights, Zarpes(k) is not symmetric about the groundstate k = (pi2 , pi2 ) along (0, 0) − (π, π). Whereas all Aαβ(k, ω) show theperiodicity of the MBZ as expected, Aarpes(k, ω) does not do so becauseof the phases eiG.Rαβ : if k is inside the magnetic BZ, G = 0 and thereforecontributions from, for example, 2px1 and 2px2 orbitals add up in Eq. (5.27),eiG.Rx2x1 = 1. If k enters the second magnetic BZ, then G = (±π,±π) andthose contributions subtract since Rx2x1 = y, therefore eiG.Rx2x1 = −1.Although the sublattice quasiparticle weights are finite, the interferencebetween various spectral weights contributing to Aarpes(k, ω) actually resultsin vanishing ARPES quasiparticle weight at the gamma point. However, itsevolution along (0, 0)− (π, π) disagrees with the ED result [3] and also with1015.6. Spin polaron vs the Zhang-Rice singletk00.050.10.150.2Zarpes(k)one-magnontwo-magnon(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)Figure 5.13: Comparison between one- and two-magnon approximation re-sults for ARPES quasiparticle weight. This shows disagreement with exper-iments for momenta along (0, 0) − (π, π). See text for further discussions.experiments, which measure a large weight near k =(pi2 , pi2 ) which decreasesfast in both directions [28]. It also shows that the quasiparticle weight issensitive to the structure of magnon cloud, despite the fact that the shapeof energy dispersion is almost independent of it. We will return to addressthis issue later in section 5.7.5.6 Spin polaron vs the Zhang-Rice singletThe three-band model can be compared with the effective one-band modelin terms of their quasiparticle states, that is the spin polaron vs the Zhang-Rice singlet (ZRS). We can do so by calculating their overlap in variousparts of the BZ.The ZRS is a coherent superposition of the singlet states of Cu andthe doping hole on four surrounding O’s whose choice of phases reflects the3dx2−y2 symmetry of central Cu,|k,ZRS〉 = 1√8N∑Reik.R∑uru[p†R+u/2,↑ − p†R+u/2,↓S+R]|N〉, (5.28)where R is the position of spin-down Cu and u/2 = ±x/2,±y/2 connectany spin-down Cu to its neighbouring O, those within the dashed circle in1025.6. Spin polaron vs the Zhang-Rice singletCuOsinglet−−++upyxyxppp1221/2Figure 5.14: Inside the dashed circle are the copper and oxygens involvedin Zhang-Rice singlet. Labelling of one-magnon GFs are such that once thecopper and hole’s spin are flipped, the current configuration becomes V 21 .The sign of each singlet contribution to ZRS is given.Fig. 5.14; ru = ±1 imposes the 3dx2−y2 symmetry. If |φ〉 is a polaron statewith energy Eφ, belonging to the complete set of eigenstates |ψ〉, for ω ≈ Eφwe haveGpα,ZRS(k, Eφ) = 〈k, pα, ↑ |Gˆ(ω ≈ Eφ)|k,ZRS〉 = limη→0∑ψ〈k, pα, ↑ |ψ〉〈ψ|k,ZRS〉Eφ − Eψ + iη≈ 〈k, pα, ↑ |φ〉〈φ|k,ZRS〉iη , (5.29)where ∑ψ |ψ〉〈ψ| = 1 is used. Similarly,Gpα,pα(k, Eφ) = 〈k, pα, ↑ |Gˆ(ω ≈ Eφ)|k, pα, ↑〉≈ 〈k, pα, ↑ |φ〉〈φ|k, pα, ↑〉/iη. (5.30)Eq. (5.29) and (5.30) can be solved together to find the magnitude of overlap1035.6. Spin polaron vs the Zhang-Rice singletk-20-16-12E qp(k)0.850.90.95|<φ | k,ZRS>|2(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)Figure 5.15: top: overlap between the states in the lowest polaron band |φ〉whose energy dispersion (in units of Jdd ≈ 0.1 − 0.15 eV) is shown at thelower panel, and the corresponding Zhang-Rice singlet states.between the polaron and ZRS states, |〈φ|k,ZRS〉|2:|〈φ|k,ZRS〉|2 = |ηGpα,ZRS(k, ω)2Gpα,pα(k, ω)|ω≈Eφ . (5.31)One has to note that O’s in the ZRS are referenced with respect to thecentral Cu at R and this differs with the way we defined our GFs in whichthe Bloch states are defined as a sum over O 2p orbitals, Eq. (5.7). Forthe one-magnon GFs, Vl, the sum runs over the position of 2p orbital onwhich the hole is initially created (for example, Rpx1 in Eq. (5.12)). Theseresult in additional phases factors ϕpi = exp(−ik.ui/2). Gpα,ZRS(k, ω) is acombination of Gαβ(k, ω) and some of one-magnon GFs, V1:√8Gpα,ZRS(k, ω) = ϕpx1Gpα,px1 (k, ω) + ϕpy1Gpα,py1 (k, ω)−ϕpx2Gpα,px2 (k, ω) − ϕpy2Gpα,py2 (k, ω)−ϕpx1 (V01 + V 11 − V 21 − V 31 ), (5.32)1045.7. From three-band to five-band model: effect of other in-plane 2p orbitalswhere V i1 are those one-magnon GFs arising when the hole starts from px1 ,see Fig. 5.14.Fig. 5.15 shows the overlap along with the dispersion relation for theone-magnon approximation, nm = 1. The overlap is generally large, sug-gesting that the spin polaron has considerable ZRS character. However, itis stronger close to the ground state and this implies that the ZRS is a moreaccurate description at low energies. As one moves away from the GS, theoverlap decreases and this signals a deviation from effective one-band modelat higher energies.5.7 From three-band to five-band model: effectof other in-plane 2p orbitalsBy investigating the effect of various terms in the Hamiltonian (5.3), weshowed that the hybridization between O 2px/y and Cu 3dx2−y2 orbitalsis essential for recovering the correct quasiparticle dispersion in the CuO2layer. Out of two in-plane 2p orbitals of each O in the layer, it thereforeseems sufficient to only consider the one along Cu-Cu bonds as it directlyhybridizes with Cu orbitals. Indeed, this has been the common approachand is also what we have so far considered. However, one can show thatthere is a considerable hopping amplitude between the other 2p orbitals andthose directly hybridizing with Cu and occupied by the hole. Therefore thequasiparticle wavefunction actually spreads to all 2p orbitals and this affectsits quasiparticle properties. In this section, we extend our analysis to takethese other orbitals into account. This results in a five-band model.We start by evaluating the hopping amplitudes involving these new or-bitals, shown in blue in Fig. 5.16. The red dots represent Cu ions. The nnhopping between new 2p orbitals is same as that between the original onestpp, but the next nn hopping t′pp is different: whereas t′pp is boosted by the4s orbital of the bridging Cu, there is no such state between new orbitalsfor the hole to tunnel through and the next nn hopping between new or-bitals t˜′pp is therefore expected to be weaker than t′pp. Furthermore, the nnhopping tpp is given as sum of the hopping for π and σ bonds between 2porbitals: tpp = tpp,σ + tpp,pi = 5tpp,σ/4 where in the last equality the relationtpp,pi = tpp,σ/4 is used [112]. Since tpp,σ scales as 1/d2 with distance anddnnn =√2dnn, it follows that t˜′pp = tpp,σ/2 = 2tpp/5. This is indeed weakerthan next nn hopping between original orbitals t′pp = 0.58tpp. In any case,we observe little dependence on the precise value of t˜′pp.Having identified the new hopping amplitudes, equations of motion can1055.7. From three-band to five-band model: effect of other in-plane 2p orbitalsFigure 5.16: The CuO2 layer with two in-plane 2p orbitals per O. The nnhopping between new (blue) and original (green) orbitals, t˜pp, and the nextnn hopping between the new orbitals, t˜′pp, are shown.be modified to include the contribution of new orbitals. Inside the unit cell,O’s are labelled by a, b, c and d as indicated in Fig. 5.16 such that pax/yrefers to 2px/y orbital associated with Oa, and so on. The calculation issimilar to the three-band model but with 8 × 8 sublattice GFs Gα,β(k, ω)instead. The details are given in Appendix G.In Fig. 5.17, I plot the lowest energy band and the associated ARPESquasiparticle weight of the five-band model for nm = 1 and nm = 2 approx-imations. The energy dispersion is quite similar to that of the three-bandmodel and including the second set of 2p orbitals does not reveal new fea-tures. Indeed, the expectation that this other set of orbitals has little effecton the quasiparticle dynamics is correct. The bandwidths are slightly largerand this reflects the increase of the bare kinetic energy of the hole due toadditional degrees of freedom provided by the new orbitals.On the other hand, introducing these new orbitals has significant effecton the ARPES quasiparticle weight specially along (0, 0)-(π, π). Zarpes(K||)maintains the asymmetry about (pi2 , pi2 ), but it shows a more consistent be-1065.7. From three-band to five-band model: effect of other in-plane 2p orbitals(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)k0.00.10.2Zarpes(K||)-22-20-18-16E qp(k)one-magnontwo-magnon(a)(b)Figure 5.17: Energy dispersion (in units of Jdd ≈ 0.1−0.15 eV) and ARPESquasiparticle weight of the five-band model that includes two 2p orbitals perO. Whereas the energy dispersions are qualitatively similar to those givenby the three-band model, the quasiparticle weight is now less sensitive tothe structure of magnon cloud.haviour for both one- and two-magnon approximations: quasiparticle weightis large near (pi2 , pi2 ) and drops fast in both directions in agreement with ex-periments data [28, 105]. This has to be contrasted with the three-bandmodel, Fig. 5.13, where the quasiparticle weight changes behaviour depend-ing on how many magnons are kept in the cloud, in spite of the fact thatthe quasiparticle dispersions were robust. The ARPES quasiparticle weightis indeed expected to be different in the five band vs the three-band model.Upon introducing the new orbitals, the hole’s wavefunction re-distribute it-self by partly occupying these orbitals. This adds new terms to, and hencemodifies, the ARPES spectral function Eq. (5.27) which is sensitive to theinterference between similar 2p orbitals.The other interesting feature is the asymmetry in the quasiparticle weight1075.7. From three-band to five-band model: effect of other in-plane 2p orbitalsbetween the inside and outside of MBZ. In Fig. 5.18, we plot the spectralfunction in the first quadrant of the FBZ for a constant energy close toground state (low doping), Aarpes(k, ω ≈ EGS). It shows that the weight islarger inside MBZ along all cuts, a feature that persists for generic values ofthe parameters. A similar drop of the quasiparticle weight is measured byARPES to occur outside the MBZ [28].ππ/20ππ/20kykx0.000.500.99Figure 5.18: Constant-energy map of ARPES spectral weight for the five-band model. It clearly shows the asymmetry of the quasiparticle weightbetween the inside and outside of the magnetic BZ, which is a robust featurefor the generic values of parameters.The results of this section suggest that a complete understanding ofthe evolution of the spectral weight at low dopings, currently still missing,may require extending the theoretical models to include these additionalorbitals. Since exact numerical approaches become even more challenging1085.8. Conclusionto implement in larger Hilbert spaces, this will only increase the need foraccurate approximations like the one we proposed here.5.8 ConclusionWe developed a simple variational approximation for studying the spectralproperties of a hole doped in three- and five-band models of CuO2 layer. Weconfirmed the accuracy of our approximation by increasing the size of thevariational space and comparing the quasiparticle’s dispersion with availableexact diagonalization results for a small cluster of the same model. Since thevariational solution ignores the spin fluctuations of Cu, while they are fullyincluded in the ED study, the agreement implies that these fluctuations donot play the significant role generally attributed to them in the quasiparticledynamics. This is a significant finding because a proper description of spinfluctuations is quite a difficult task and it has hindered studying the two-hole case which provides information about the interaction between (andpossible bound states of) quasiparticles.In the three-band model, the hole can freely move on the O sublatticevia Tˆpp and therefore it is not localized near spin flips in the absence of spinfluctuations. In the one-band model, however, the ZRS moves in the mag-netic sublattice and this creates strings of spin flips whose energy increasesalmost linearly with the length of the strings [88]. The only way for the holeto move then is by executing Trugman loops and this results in a very heavyquasiparticle. Spin fluctuations are hence needed to annihilate pairs of spinflips and release the hole, but the resulting energy dispersion is wrong asbeing almost flat along (0, π) − (π, 0). This problem is overcome by intro-ducing second and third nn hopping in the one-band model which allow theZRS to move away from the spin flips created by the nn hops [106]. In thethree-band model, this is accomplished by Tˆpp and this justifies the lack ofimportance of the spin fluctuations.We also extended our study into a five band model by including thosein-plane 2p orbitals that do not directly hybridize with Cu. This resultedin a better agreement experiments for the quasiparticle weights. One hasto note that doing so with exact numerical methods is computational verycostly, whereas the efficiency of our method enabled us to readily implementthis with minimum computational effort.109Chapter 6Summary and development6.1 Summary of this workThis thesis was an effort toward understanding the spectral properties ofpolaronic quasiparticles and their interaction with external potentials fromimpurities or more extended disorder. In particular, the Holstein polaronand spin polaron in antiferromagnetic insulators were considered.In chapter 2, we extended the Momentum-Average approximation forcalculating the Green’s function of the Holstein polaron in systems withvarious types of disorder. MA is unique in that it sums all of the termsin the full diagrammatic expansion, but approximates each by discardingexponentially small contributions. It is exact in both zero coupling and zerobandwidth limit, and is also accurate everywhere is the parameter space. Ithas a variational interpretation where the number of lattice sites over whichthe polaron’s phonon cloud is extended can be varied in order to improvethe accuracy of the approximation.Dyson’s identity was used to calculate the polaron Green’s function inreal space. The main finding was that coupling to phonons, in additionto modifying the effective mass of the electron, renormalizes the disorderpotential seen by the polaron in a non-trivial and retarded fashion. Weshowed that this renormalization is the result of the binding of the electronto its phonon cloud. It affects the electron’s ability to interact with theoutside potential in a way that depends on the structure of the cloud.This formalism was used to investigate the binding of a Holstein polaronto a single impurity in three dimensions. Various kinds of impurities thatmodify the on-site energy and the strength of electron-phonon coupling, andisotope substitutions that change the coupling and phonon frequency wereconsidered. The phase diagram of polaron trapping was constructed in eachcase and differences and similarities were discussed.In chapter 3, we used the same formalism (suitably generalized) to studythe lifetime of Holstein polaron due to weak and uncorrelated disorder thatis extended over the entire system, also known as the (weak) Anderson dis-order. When the strength of the electron-phonon coupling is weak, we found1106.1. Summary of this workthat the polaron’s inverse lifetime in units of its effective mass agrees withFermi’s golden rule prediction. This suggested that the effect of weak cou-pling to phonons is to enhance the effective mass of the electron. For strongcoupling, significant deviation from the prediction of Fermi’s golden rule wasobserved. We argued that this was due to the fact that the renormalizationof the disorder potential becomes increasingly important at strong coupling.Therefore, the polaron scatters from a modified disorder potential and thisdifference is not captured by Fermi’s golden rule.In chapters 4, we switched to magnetically ordered systems where wefocused on the spin polaron description of hole-like excitations in two dimen-sional antiferromagnetic (AFM) insulators. The nearest-neighbour hoppingof a hole in an AFM on a square lattice creates a string of spin flips whoseenergy increases almost linearly with the length of the string. In the ab-sence of the spin fluctuations, the hole may therefore become localized inthe vicinity of spin flips. In fact, the hole can still propagate in the lattice byexecuting the so-called Trugman loop processes during which it goes almosttwice around closed loops. The spin flips created during the first round arecleaned during the second round such that, when the last one of them isremoved, the hole arrives at a different site of the same sublattice it startedon. An analytic solution was obtained for the hole’s Green’s function bylimiting the maximum number of such spin flips. It confirmed that the holewith the accompanying cloud of spin flips moves on one sublattice and theeffective on-site energy and hopping amplitudes characterizing its motionwere identified.We then introduced an attractive impurity that changed the hole’s on-site energy. Both cases of the impurity and the polaron being on the samevs on different sublattices were considered. Bound states were observed inboth cases, and this was in spite of the fact that the polaron in its ownsublattice was not subject to a bare impurity potential in the latter case. Itwas argued that this is a manifestation of the renormalization of impuritypotential due to coupling to bosons, the concept that was introduced in thesecond chapter in the context of Holstein polarons. The hole can be found atvarying distances from its boson cloud and this includes the other sublatticewhere it can interact with the impurity. Whereas a single bound state wasfound when the polaron and impurity were on the same sublattice, a total ofthree bound states were identified for the other case and their symmetry andthe dependence of binding energies on the bare hole hopping were analyzed.In the last chapter a more realistic model for the motion of a hole in aCuO2 layer, in the absence of disorder, was constructed which included theO 2p orbitals explicitly. The model involved two types of spin fluctuations:1116.2. Outlook and further developmentsthat of the AFM ordered Cu spins, and those caused by the doping hole. Weneglected the first type, and kept the latter. A variational scheme was thenconstructed where we allowed for the creation of up to a certain maximumnumber of spin flips as the hole moved in the lattice. The convergence ofour approximation was established by comparing the energy dispersion ofthe lowest polaron band with that given by the exact diagonalization studyof the same model, as well as with the experimental results. In this modelthe hole is able to move freely within the O sublattice and, as a result, spinfluctuations are not essential for releasing the hole as a light polaron (asthey are in one-band models with only nn hopping).Our approach is computationally quite efficient and this allowed us togeneralize it to investigate the effect of other in-plane O 2p orbitals on thehole’s dynamics. Whereas including these other orbitals resulted essentiallyin similar quasiparticle dispersions, the quasiparticle weights were found tobe considerably modified upon including new orbitals and better agreementwith the experiments was obtained.6.2 Outlook and further developmentsOne of the common themes throughout this work was the use of variationaltechniques as a way of studying many-body problems. They basically in-volved choosing a subspace, from the entire Hilbert space of the problem,that includes the energetically most relevant states. This enabled us to per-form semi-analytic and non-perturbative calculations. Although the resultsin principle depend on the choice of variational space, a desired level ofaccuracy can be reached by varying the size of the space.This variational scheme has been originally suggested and successfullyapplied for studying Holstein polarons [50]. We succeed in verifying thatthese methods are not limited to lattice polarons and can be used for study-ing polarons in magnetic systems as well. Similar ideas have been appliedby others in the context of orbital polarons.In chapters 2 and 3, we limited our study to either a single impurity orextended but weak disorder. A natural extension is to consider extendeddisorder of arbitrary strength. A new phenomenon that may occur in thisregime is the Anderson localization of Holstein polarons. It is well knownthat the disorder averaged descriptions, as the one used in chapter 3 for weakdisorder, fail to capture localization. Instead, the signature of localizationcan be traced to the distribution of quantities such as the local densityof states, P[{ρi(ω)}], which is suggested to change from normal to log-1126.2. Outlook and further developmentsnormal when the energy state at ω becomes localized [42]. The formalismdeveloped in chapter 2 allows us to calculate the local densities of states forany disorder configuration, but a challenge here is to be able to considerlarge enough system sizes in order to have enough number of LDOS valuesfor a meaningful statistical analysis.Given the success of our variational scheme in describing the spectralproperties of a single hole in the CuO2 layer, extending it to study a pairof holes can be rewarding. In particular, one can address the question ofwhether there are any coherent bound states of two holes, and if so, whattheir pairing symmetry is. The advantage of our Green’s function formalismfor this problem is that since the Green’s function of the pair gives the rela-tive wave function of the two holes, it will enable one to explicitly map outthe pairing symmetry in the real space (provided that bound states exist).This will provide a microscopic picture of the pre-formed pairs in cuprate su-perconductors and it contrasts with common macroscopic descriptions wherethe pairing symmetry is studied not in the real space, but it refers to thestructure of superconducting energy gap in the momentum space.The analysis in chapter 4 can be generalized for studying other variantsof the t-Jz model. Within the three-magnon variational space, t-J modelis mimicked by locally including spin fluctuations which relates zero- andone-magnon Green’s functions to two- and three-magnon ones, respectively,and vice versa. Further range hoppings can also be treated similarly withinthe same variational space. This enables studying t-t′-t′′-J model as well asthe effect of three-site hopping which is among the terms appearing in thestrong coupling limit of the Hubbard model [113]. We found that, unlikethe three-band model, spin fluctuations are actually essential for recoveringthe correct quasiparticle dispersion in one-band models. We argued thatthe different role played by spin fluctuations suggests that quasiparticlesdescribed by one- and three-band models are qualitatively different. Thiscasts doubt on whether one-band t-J like models are suitable for studyingcuprate physics. This work [107], which is under review for publication, doesnot appear in this thesis.113Bibliography[1] M. Berciu, A. S. Mishchenko and N. Nagaosa, EPL 89, 37007 (2010).[2] M. Berciu and and H. Fehske, Phys. Rev. B 84, 165104 (2011).[3] B. Lau, M. Berciu, and G. A. Sawatzky, Phys. Rev. Lett. 106, 036401(2011).[4] A. S. Mishchenko, N. Nagaosa, A. Alvermann, H. Fehske, G. De Fil-ippis, V. Cataudella, and O. P. Sushkov, Phys. Rev. B 79, 180301(R)(2009).[5] C. Kittle, Introduction to Solid State Physics, 8th edition (John WileySons, Inc., 2005).[6] H. J. Schulz, in ”Proceedings of Les Houches Summer School LXI”,[http://arxiv.org/abs/cond-mat/9503150] (1995).[7] L. D. Landau, Phys. Zh. Sowjet 3, 664 (1933).[8] H. Ba¨ssler and A. Ko¨hler, Top. Curr. Chem. 312, 1 (2012).[9] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175(1957).[10] C. P. J. Adolphs and M. Berciu, EPL 102, 47003 (2013).[11] C. P. J. Adolphs and M. Berciu, Phys. Rev. B 89, 035122 (2014).[12] L. Covaci and M. Berciu, Europhys. Lett. 80, 67001 (2007).[13] T. Holstein, Ann. Phys. (N.Y.) 8, 325 (1959); 8, 343 (1959).[14] M. Berciu, Phys. Rev. Lett. 97, 036402 (2006); G. L. Goodvin, M.Berciu, and G. A. Sawatzky, Phys. Rev. B 74, 245104 (2006).[15] I. G. Lang and Y. A. Firsov, Sov. Phys. JETP 16, 1301 (1963).[16] F. Marsiglio, Phys. Lett. A 180, 280 (1993).114Bibliography[17] A. A. V. Kabanov, and D. Ray, Phys. Rev. B 49, 9915 (1994).[18] E. V. L. de Mello, and J. Ranninger, Phys. Rev. B 55, 14872 (1997).[19] P. Kornilovitch, E. Pike, Phys. Rev. B (Rapid) 55, 8634 (1997).[20] P. Kornilovitch, Phys. Rev. Lett. 81, 5382 (1998).[21] P. Kornilovitch, Phys. Rev. B 60, 3237 (1999).[22] H. Fro¨hlich, H. Pelzer, and S. Zienau, Philos. Mag. 41, 221 (1950).[23] C. Slezak, A. Macridin, G. A. Sawatzky, M. Jarrell, and T. A. Maier,Phys. Rev. B 73, 205122 (2006).[24] B. Lau, M. Berciu, and G. A. Sawatzky, Phys. Rev. B 76, 174305(2007).[25] A. S. Alexanrov, and J. T. Devreese, Advances in Polaron Physics,Springer Series in Solid-State Science, Vol. 159 (2010).[26] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698(1979).[27] A. J. Heeger, S. Kivelson, J. R. Schrieffer, and W. -P. Su, Rev. Mod.Phys. 60, 781 (1988).[28] A. Damascelli, Z. Hussain, and Z. X. Shen, Rev. Mod. Phys. 75, 473(2003).[29] P. W. Anderson, Phys. Rev. 109, 1493 (1958).[30] B. Kramer, and A. MacKinnon, Rep. Prog. Phys. 56, 1469 (1993).[31] P. A. Lee and D. S. Fisher, Phys. Rev. Lett. 47, 882 (1981).[32] H. Ebrahimnejad and M. Berciu, Phys. Rev. B 85, 165117 (2012).[33] H. Ebrahimnejad and M. Berciu, Phys. Rev. B 88, 104410 (2013).[34] D. Marchand, PhD Thesis, University of British Columbia (2011).[35] G. L. Goodvin, PhD Thesis, University of British Columbia (2009).[36] F. C. Zhang, and T. M. Rice, Phys. Rev. B 37, 3759(R) (1988).115Bibliography[37] A. Lanzara, P. V. Bogdanov, X. J. Zhou, S. A. Kellar, D. L. Feng, E.D. Lu, T. Yoshida, H. Eisaki, A. Fujimori, K. Kishio, J. -I. Shimoyama,T. Noda et al., Nature 412, 510 (2001).[38] D. N. Basov, B. Dabrowski, and T. Timusk, Phys. Rev. Lett. 81, 2132(1998).[39] J. Lee, K. Fujita, K. McElroy, J. A. Slezak, M. Wang, Y. Aiura, H.Bando, M. Ishikado, T. Masui, J.-X. Zhu, A. V. Balatsky, H. Eisaki,S. Uchida et al., Nature 442, 546 (2006).[40] S. M. Girvin and M. Jonson, Phys. Rev. B 22, 3583 (1980).[41] O. Gunnarsson, M. Calandra, and J. E. Han, Rev. Mod. Phys. 75,1085 (2003).[42] F. X. Bronold, and H. Fehske, Phys. Rev. B 66, 073102 (2002); F. X.Bronold, A. Alvermann, and H. Fehske, Philos. Mag. 84, 673 (2004).For LDA, see A. Alvermann and H. Fehske: Local Distribution Ap-proach, Lecture Notes in Phys. 739, 505-526 (Springer, 2008).[43] S. Ciuchi, F. de Pasquale, S. Fratini, and D. Feinberg, Phys. Rev. B56, 4494 (1997).[44] H. Sumi, J. Phys. Soc. Jpn. 36, 770 (1974).[45] F. X. Bronold, A. Saxena and A. R. Bishop, Phys. Rev. B 63, 235109(2001).[46] J. P. Hague, P. E. Kornilovitch, and A. S. Alexandrov, Phys. Rev. B78, 092302 (2008).[47] L. Covaci and M. Berciu, Phys. Rev. Lett. 102, 186403 (2009).[48] M. Berciu and G. L. Goodvin, Phys. Rev. B 76, 165109 (2007).[49] A. N. Das and S. Sil, Phys. Lett. A 348,266 (2006); in a somewhatdifferent context, see also A. N. Das and S. Sil. J. Phys.: Condens.Matter 5, 8265 (1993) or A. Macridin, G. A. Sawatzky and MarkJarrell, Phys. Rev. B 69, 245111 (2004).[50] M. Berciu and G. L. Goodvin, Phys. Rev. B 76, 165109 (2007); M.Berciu, Phys. Rev. Lett. 98, 209702 (2007).[51] O. S. Barisic, Phys. Rev. Lett. 98, 209701 (2007).116Bibliography[52] A. M. Clogston, Phys. Rev. 125, 439 (1962).[53] P. A. Wolff, Phys. Rev. 124, 1030 (1961).[54] M. Berciu and A. M. Cook, Europhys. Lett. 92, 40003 (2010).[55] M. Berciu and H. Fehske, Phys. Rev. B 82, 085116 (2010); D. J. J.Marchand, G. De Filippis, V. Cataudella, M. Berciu, N. Nagaosa, N.V. Prokof’ev, A. S. Mishchenko, and P. C. E. Stamp, Phys. Rev. Lett.105, 266605 (2010).[56] J. Bonca, S.A. Trugman, and I. Batistic, Phys. Rev. B 60, 1633 (1999).[57] A. S. Mishchenko and N. Nagaosa, Phys. Rev. B 73, 092502 (2006).[58] C. Bernhard, T. Holden, A. V. Boris, N. N. Kovaleva, A. V. Pimenov,J. Humlicek, C. Ulrich, C. T. Lin, and J. L. Tallon, Phys. Rev. B 69,052502 (2004).[59] P. W. Anderson, Phys. Rev. 109, 1492 (1958).[60] A. Miller and E. Abrahams, Phys. Rev. 120, 745 (1960).[61] A. A. Abrikosov, L. P. Gorkov and I. Ye. Dzyaloshinski, Methods ofquantum field theory in statistical physics (Dover, New York, 1975).[62] Strictly speaking, the lifetime is defined with an additional factor of2.[63] M. Berciu, Phys. Rev. Lett. 98, 209702 (2007).[64] C. L. Kane, P. A. Lee, and N. Read, Phys. Rev. B 39, 6880 (1989).[65] F. Marsiglio, A. E. Ruckenstein, S. Schmitt-Rink, and C. M. Varma,Phys. Rev. B 43, 10882 (1991).[66] S. Sachdev, Phys. Rev. B 39, 12232 (1989).[67] A. V. Chubukov and D. K. Morr, Phys. Rev. B 57, 5298 (1998).[68] B. I. Shraiman and E. D. Siggia, Phys. Rev. Lett. 60, 740 (1988).[69] P. A. Lee, N. Nagaosa, and X. G. Wen, Rev. Mod. Phys. 78, 17 (2006).[70] P. Phillips, Rev. Mod. Phys. 82, 1719 (2010).[71] E. Manousakis, Rev. Mod. Phys. 63, 1 (1990).117Bibliography[72] E. Dagotto, R. Joynt, A. Moreo, S. Bacci and E. Gagliano, Phys. Rev.B 41, 9049 (1990); H. Fehske, V. Waas, H. Ro¨der and H. Bu¨ttner,ibid. 44, 8473 (1991); E. Dagotto, Rev. Mod. Phys. 66, 763 (1994).[73] H. Alloul, J. Bobroff, M. Gabay, and P. J. Hirschfeld, Rev. Mod. Phys.81, 45 (2009).[74] D. Heidarian and N. Trivedi , Phys. Rev. Lett. 93, 126401 (2004); YunSong, R. Wortis and W. A. Atkinson, R. Wortis and W. A. Atkinson,Phys. Rev. B 77, 054202 (2008); Phys. Rev. B 82, 073107 (2010).[75] A. Suchaneck, V. Hinkov, D. Haug, L. Schulz, C. Bernhard, A. Ivanov,K. Hradil, C. T. Lin, P. Bourges, B. Keimer, and Y. Sidis, Phys. Rev.Lett. 105, 037207 (2010).[76] S. Wessel, B. Normand, M. Sigrist, and S. Haas, Phys. Rev. Lett. 86,1086 (2001).[77] D. Poilblanc, D. J. Scalapino, and W. Hanke, Phys. Rev. Lett. 72, 884(1994).[78] D. Poilblanc, D. J. Scalapino, and W. Hanke, Phys. Rev. B 50, 13020(1994).[79] W. F. Brinkman, and T. M. Rice, Phys. Rev. B 2, 1324 (1970).[80] D. Poilblanc, H. J. Schulz, and H. T. Ziman, Phys. Rev. B 47, 3268(1993).[81] D. M. Edwards, Physica B 133, 378 (2004).[82] A. Alvermann, D. M. Edwards, and H. Fehske, Phys. Rev. Lett. 98,056602 (2007).[83] D. M. Edwards, Phys. Rev. B 82, 085116 (2010).[84] J. Riera and E. Dagotto, Phys. Rev. B 47, 15346 (1993); E. La-houd, O. Nganba Meetei, K. B. Chaska, A. Kanigel, and N. Trivedi,arXiv:1303.0649v1.[85] Note that Eq. (18) in Ref. [2] has a typo in the denominator, wherethe prefactor 4 behind t3(ω) has to be replaced with 2.[86] M. Mo¨ller, A. Mukherjee, C. P. J. Adolphs, D. J. J. Marchand, andM. Berciu, J. Phys. A: Math. Theor. 45, 115206 (2012).118Bibliography[87] P. Chudzinski, M. Gabay, and T. Giamarchi, New J. Phys. 11, 055059(2009).[88] S. A. Trugman, Phys. Rev. B 37, 1597 (1988).[89] V. J. Emery, Phys. Rev. Lett. 58, 2794 (1987).[90] H. Eskes, and G. A. Sawatzky, Phys. Rev. Lett. 61, 1415 (1988).[91] T. Giamarchi, and C. Lhuillier, Phys. Rev. B 47, 2775 (1993).[92] E. Dagotto, A. Moreo, R. Joynt, S. Bacci, and E. Gagliano, Phys. Rev.B 41, 2585(R) (1990).[93] E. Dagotto, R. Joynt, A. Moreo, S. Bacci, and E. Gagliano, Phys. Rev.B 41, 9049 (1990).[94] V. Elser, D. A. Huse, B. I. Shraiman, and E. D. Siggia, Phys. Rev. B41, 6715 (1990).[95] M. Boninsegni, and E. Manousakis, Phys. Rev. B 43, 10353 (1991).[96] D. Poilblanc, H. J. Schulz, and T. Ziman, Phys. Rev. B 46, 6435(1992).[97] M. M. Sala et al., New J. Phys. 13, 043026 (2011).[98] J. Zaanen, G. A. Sawatzky, and J. W. Allen, Phys. Rev. Lett. 55, 418(1985).[99] G. Martinez, and P. Horsch, Phys. Rev. B 44, 317 (1991).[100] A. Fujimori, E. Takayama-Muromachi, Y. Uchida, and B. Okai, Phys.Rev. B 35, 8814 (1987); J. M. Tranquada, S. M. Heald, A. R. Mood-enbaugh, and M. Suenaga, ibid. 35, 7187 (1987).[101] M. Ogata, and H. Fukuyama, Rep. Prog. Phys. 71, 036501 (2008); P.A. Lee, ibid. 71, 012501 (2008).[102] B. Lau, PhD Thesis, University of British Columbia (2011).[103] H. Ebrahimnejad, G. A. Sawatzky, and M. Berciu, Nature Physics 10,951-955 (2014).[104] P. W. Leung, B. O. Wells, and R. J. Gooding, Phys. Rev. B 56, 6320(1997).119[105] B. O. Wells, et al., Phys. Rev. Lett. 74, 964 (1995).[106] J. H. Jefferson, H. Eskes, and L. F. Feiner, Phys. Rev. B 45, 7959(1992).[107] “The t-J and the Emery model have qualitatively different quasipar-ticles”, H. Ebrahimnejad and M. Berciu, submitted (2014).[108] T. Xiang and J. M. Wheatley, Phys. Rev. B 54, R12653 (1996).[109] V. J. Emery, and G. Reiter, Phys. Rev. B 38, 4547 (1988).[110] D. M. Frenkel, R. J. Gooding, B. I. Shraiman, and E. D. Siggia, Phys.Rev. B 41, 350 (1990).[111] Y. Petrov, T. Egami, Phys. Rev. B 58, 9485 (1998).[112] W. A. Harrison, Elementary Electronic Structure (World Scientific,Singapore, (1999).[113] J. Ba la, Eur. Phys. J. B 16, 495500 (2000).[114] R. J. Elliot, J. A. Krumhansl and P. L. Leath, Rev. Mod. Phys. 46,465 (1974).[115] E. N. Economou, Green’s Functions in Quantum Physics, 3rd edition(Springer-Verlag, New York, 2006), or A. Gonis, “Green Functions forOrdered and Disordered Systems” (Elsevier, Amsterdam, 1992).120Appendix AMA0 solution for the cleansystemUsing the notation of chapter 2, the Holstein Hamiltonian for a clean systemisHH = H0 + Vˆel−ph, (A.1)where we use the real-space notation for electron-phonon coupling term:Vˆel−ph = g∑ic†ici(bi + b†i ).H0 is the sum of kinetic energy of particle and lattice vibrational energy.Our goal is to calculate the retarded single polaron Green’s functionGH(k, ω) = 〈0|ckGˆH(ω)c†k|0〉,using the Dyson’s identity GˆH(ω) = Gˆ0(ω)+ GˆH(ω)Vˆel−phGˆ0(ω). Projectingthis onto single electron Bloch states givesGH(k, ω) = G0(k, ω)[1 + g∑ieik·Ri√NF (1)ki (ω)], (A.2)where F (n)ki (ω) = 〈0|ckGˆH(ω)ci†b†ni |0〉. Note that F(0)ki (ω) = GH(k, ω) exp(−ik·Ri)/√N . Using the Dyson identity again, we find that for n ≥ 1F (n)ki (ω) = g∑j 6=iG0ji(ω − nΩ)〈0|ckGˆH(ω)c†jb†jb†ni |0〉+gG0ii(ω − nΩ)[nF (n−1)ki (ω) + F(n+1)ki (ω)], (A.3)whereG0ji(ω) =∑keik·(Rj−Ri)N G0(k, ω)121Appendix A. MA0 solution for the clean systemis the free particle propagator in real space and the sum is over the momentainside the first Brillouin zone. This propagator decays exponentially with thedistance |Rj −Ri| for energies outside the free particle continuum, |ω| > 6t.Since we are interested in energies ω−nΩ ≈ EP,GS−nΩ, where the polaronGS energy Ep,GS < −6t, all these propagators become exponentially smallfor j 6= i. It is therefore a reasonable first approximation to ignore j 6= iterms in Eq. (A.3). This is what the MA0 approximation does (higherflavors include j 6= i terms in a certain progression [48]).Within MA0, then, we have for any n ≥ 1F (n)ki (ω) = gg0(ω − nΩ)[nF (n−1)ki (ω) + F(n+1)ki (ω)], (A.4)where g0(ω) = G0ii(ω); see Eq. (3.15). This recurrence relation is solved interms of continued fractions [14]:F (n)ki (ω) = An(ω)F(n−1)ki (ω) (A.5)for any n ≥ 1, whereAn(ω) =ngg0(ω − nΩ)1 − gg0(ω − nΩ)An+1(ω). (A.6)Finally, using F (1)ki (ω) = A1(ω)F(0)ki (ω) in Eq. (A.2) leads to the MA0 solu-tion:GH(k, ω) =1ω − ε(k) − gA1(ω) + iη.ΣMA(ω) = gA1(ω) is therefore identified as the self energy associated withHolstein type electron-phonon coupling.122Appendix BIMA1At the IMA1 level, one also allows processes in which one phonon is awayfrom the phonon cloud. These are described by the propagators Sn(i, l, j;ω) =〈0|ciG(ω)c†l b†n−1l b†j |0〉 with j 6= l. In terms of these, Eq. (2.3) can be writtenasGij(ω) = Gdij(ω) +∑lglS1(i, l, l;ω)Gdlj(ω). (B.1)Once again we apply the Dyson identity to S1S1(i, l, j;ω) = gjGdjl(ω − Ωj)Gij(ω) (B.2)+∑mgmGdml(ω − Ωj)S2(i,m, j;ω).This exact equation relates S1 to the propagators S2. We can similarly findthe equation of motion of all higher, n ≥ 2, Sn(i, l, j;ω) with l 6= j and l = j,separately. For l 6= j we haveSn(i, l, j;ω) = glGdll(ω − (n− 1)Ωl − Ωj) (B.3)×[(n− 1)Sn−1(i, l, j;ω) + Sn+1(i, l, j;ω)],where we now ignore contributions from terms with a second phonon awayfrom the polaron cloud, as they are exponentially smaller than those we kept.This admits the solution Sn(i, l, j;ω) = Bn(l, j;ω)Sn−1(i, l, j;ω), whereBn(l, j;ω) =(n− 1)glGdll(ω − (n− 1)Ωl − Ωj)1 − glGdll(ω − (n− 1)Ωl − Ωj)Bn+1(l, j;ω)= An−1(l, ω − Ωj). (B.4)For l = j and n ≥ 2, Sn(i, l, l;ω) = F (n)il (ω) and we get the same solution asin MA0, i.e. Sn(i, l, l;ω) = An(l, ω)Sn−1(i, l, l;ω). The relations between S2and S1 are used in Eq. (B.2) to turn it into an equation between S1(ω) andGij(ω) onlyS1(i, l, j;ω) = gjGdjl(ω − Ωj)Gij(ω) + gjGdjl(ω − Ωj)A2(j, ω)S1(i, j, j;ω)+∑m6=jgmGdml(ω − Ωj)A1(m,ω − Ωj)S1(i,m, j;ω).123Appendix B. IMA1Together with Eq. (B.1), this can be solved to find Gij(ω). However, it isagain convenient to explicitly extract the “average” contributions, to makethese equations more efficient.We therefore remove the homogeneous part from Eq. (B.2) and includeit into a renormalized energyS1(i, l, j;ω) = gjGdjl(ω¯j)[Gij(ω) + (A2(j, ω) −A1(j, ω − Ωj))S1(i, j, j;ω)]+∑m6=jgmGdml(ω¯j)[A1(m,ω − Ωj) −A1(ω − Ω)]S1(i,m, j;ω). (B.5)where ω¯j = ω − Ωj − gA1(ω − Ω), and A1(ω − Ω) is given by Eq. (2.11)for the “average” clean system. The sum on the rhs of Eq. (B.5) againconvergences for a very small cutoff, only sites m very close to j need to beincluded. Its general solution is of the form:S1(i, l, j;ω) = xjl(ω)[Gij(ω)+(A2(j, ω)−A1(j, ω−Ωj))S1(i, j, j;ω)], (B.6)wherexjl(ω) = gjGdjl(ω¯j) +∑m6=jgmGdml(ω¯j)[A1(m,ω − Ωj) −A1(ω − Ω)]xjm(ω).In fact, using xjl(ω) = gjGdjl(ω¯j) is already a very good approximation, sincethe terms in the sum are exponentially small – but one can go beyond this.Once xjj(ω) is known, from Eq. (B.6) we find S1(i, j, j, ω) = Λj(ω)Gij(ω),whereΛj(ω) =xjj(ω)1 − xjj(ω) (A2(j, ω) −A1(j, ω − Ωj)).This can now be used in Eq. (B.1) to turn it into an equation for Gij(ω)only:Gij(ω) = Gdij(ω) +∑lglΛl(ω)Gil(ω)Gdlj(ω). (B.7)As we did in Eq. (2.13) for MA0, this can be made efficient to solve bysubtracting the MA1 self-energy and including it into the energy argumentGij(ω) = Gdij(ω˜) +∑lGil(ω)v1(l, ω)Gdlj(ω˜), (B.8)in which ω˜ = ω−ΣMA1(ω) and v1(l, ω) = glΛl(ω)−ΣMA1(ω). Here, ΣMA1(ω)is the value of glΛl(ω) in the clean, “average” system:ΣMA1(ω) =g2g0(ω¯)1 − gg0(ω¯) (A2(ω) −A1(ω − Ω))124Appendix B. IMA1where now ω¯ = ω − Ω − gA1(ω − Ω) [48]. This completes the calculation ofGF within inhomogeneous MA1 approximation.125Appendix CMA0 solution for F (n)k,si(ω) andWnm(ω)First, let’s we evaluateF (n)k,si(ω) = 〈0|ckGˆH(ω)c†i b†ns |0〉within MA0. In this level, this propagator vanishes for all i 6= s, because it isproportional to G0si(ω−nΩ). Therefore, all we need to evaluate is F(n)ki (ω) =〈0|ckGˆH(ω)c†i b†ni |0〉 which have been already calculated in Appendix A,F (n)ki (ω) = An(ω)F(n−1)ki (ω)=n∏l=1Al(ω)F (0)ki (ω)≡ Γn(ω)F (0)ki (ω), (C.1)where Γn(ω) ≡∏nl=1Al(ω) and F(0)ki (ω) = exp(−ik.Ri)GH(k, ω)/√N . Al(ω)are the continued fractions introduced in the text.Let’s now calculate theW nm(ω) = 〈0|bni ciGˆH(ω)c†i b†mi |0〉. Because W nm(ω) =Wmn(ω), we only need to find W nm(ω) for m ≤ n. Note that we alreadyknow W 00(ω) = 〈0|ciGˆH(ω)c†i |0〉 = 1N∑kGH(k, ω) = g0(ω−ΣMA(ω)). Fur-thermore,W n0(ω) = 〈0|bni ciGˆH(ω)c†i |0〉=∑ke−ik·Ri√N[F (n)ki (ω)|η→−η]∗= Γn(ω)g0(ω − ΣMA(ω)). (C.2)For any m ≥ 1, the equations of motion for W nm(ω) within the MA approx-imation (i.e., not allowing the electron to change its site when phonons are126Appendix C. MA0 solution for F (n)k,si(ω) and W nm(ω)present) leads toW nm(ω) = m!g0(ω −mΩ)δnm + gg0(ω −mΩ)×[mW n,m−1(ω) +W n,m+1(ω)]. (C.3)For m > n, the delta function vanishes and Eq. (C.3) is identical to Eq.(A.4), henceW nm(ω) = Am(ω)W n,m−1(ω). (C.4)In particular, this gives W n,n+1(ω) = An+1(ω)W nn(ω). Using this in Eq.(C.3) with n = m relates W nn(ω) to W n,n−1(ω),W nn(ω) = An(ω)W n,n−1(ω) +(n − 1)!g An(ω). (C.5)This is taken together with the EOM for 1 ≤ m ≤ n− 1W nm(ω) = gg0(ω −mΩ) × [mW n,m−1(ω) +W n,m+1(ω)]to give a system of n equations with n unknowns W nm(ω), m = 1, ..., n[W n0(ω) is known; see above]. This can be solved in many ways, includingdirect numerical solution. A nicer approach is to use the linearity of this sys-tem of equations to split it into two different systems, with (n−1)!g An(ω) andgg0(ω−Ω)W n0(ω) as inhomogeneous parts. These can be solved analyticallyto giveW nm(ω) = Γm(ω)W n0(ω) + Γ˜m(ω)(n− 1)!An(ω)g[1 −An(ω)Bn(ω)],where Γ˜m(ω) = Bm+1(ω)Bm+2(ω) · · ·Bn(ω) for m < n where Γ˜n(ω) = 1,andBm+1(ω) =gg0(ω −mΩ)1 − (m− 1)gg0(ω −mΩ)Bm(ω)are continued fractions ending at B2(ω) = gg0(ω − Ω).127Appendix DDisorder self-energy inAverage T-matrixapproximationAverage T-matrix approximation (ATA) [114, 115] relates the disorder partof the self-energy of a single particle to the disorder average of its transfermatrix through a single impurity whose on-site energy is ǫ:ΣATA(ω) =t¯1 + t¯g0(ω), (D.1)where t = ǫ/(1 − ǫg0(ω)) is the sum over all single impurity scattering con-tributions. g0(ω) is the momentum average of the free particle propagator,see Eq. (3.15). For Anderson-type disorder of width 2∆ we findt¯ = 12∆∫ ∆−∆ǫdǫ1 − ǫg0(ω)= − 1g0(ω)+ 12∆g20(ω)ln 1 + ∆g0(ω)1 − ∆g0(ω).Expanding to lowest order in ∆ regains the perturbational result, Eq. (3.14):ΣATA(ω) ≈∆23 g0(ω) = σ2g0(ω).Therefore, differences between ATA and FGR show that disorder is so largethat multiple scattering processes off the same impurity cannot be ignoredanymore.To extend ATA to the Holstein model, we note that difference betweenthe MA clean polaron’s GF and that of the free electron is the appearanceof ΣMA(ω). This simply modifies g0(ω) → g0(ω − ΣMA(ω)) in all the aboveequations. As discussed above, this approximation implies that there is nocrossing between phonon lines and scattering lines, in diagramatic terms.Our results show that this is a bad approximation for medium to strongelectron-phonon coupling.128Appendix EThe equations of motion forG0,R(ω)Here we present the details of the calculations that lead to Eq. (4.12), whichrelates the various G0,R(ω) GFs. Eq. (4.11) enables us to eliminate F3 fromEq. (4.10) to obtain:F2(R,u,v) − t2g¯2g3F2(R + u + v,−v,−u) = −tg¯2F1(R,u), (E.1)andF2(R + u + v,−v,−u) − t2g¯2g3F2(R,u,v) = −tg¯2F1(R + u + v,−v),(E.2)where g¯2 = 1/(ω − 14J¯ − t2g3 + iη) and Eq. (E.2) results from Eq. (E.1)after changing the coordinates R → R + u + v, u → −v, v → −u. Solvingthe coupled equations (E.1) and (E.2), we findF2(R,u,v) = γ1F1(R,u) + γ2F1(R + u + v,−v), (E.3)in which γ1 = −tg¯2/[1 − (t2g¯2g3)2] and γ2 = t2g¯2g3γ1. Using this in Eq.(4.9) givesF1(R,u) = − tg¯1[G0,R + γ2∑v⊥uF1(R + u + v,−v)], (E.4)where g¯1 = 1/(ω − 10J¯ + 2tγ1 + iη) and the sum includes the two nearest-neighbor vectors, ±v, along the direction perpendicular to u. With a properchange of coordinates, each F1 on the right-hand side of Eq. (E.4) can beexpressed in term of a component of G and new F1’s. For example,F1(R + u + v,−v) + tg¯1G0,R+u+v = −tg¯1γ2[F1(R + 2u,−u)+F1(R,u)], (E.5)andF1(R + u− v,v) + tg¯1G0,R+u−v = −tg¯1γ2[F1(R + 2u,−u)+F1(R,u)], (E.6)129Appendix E. The equations of motion for G0,R(ω)which results after applying either of R→ R + u± v, u → ∓v, v → u toEq. (E.4), respectively. The additionally introduced F1 can be written interms of the existing ones by doing R → R + 2u, u → −u on Eq. (E.4)F1(R + 2u,−u) + tg¯1G0,R+2u = −tg¯1γ2[F1(R + u + v,−v)+F1(R + u− v,v)]. (E.7)The four equations (E.4) to (E.7) can be simultaneously solved for the fourF1’s in terms of the existing components of G. In particular, we find:F1(R,u) = ζ1G0,R + ζ2G0,R+2u+ζ3[GR+u+v +GR+u−v], (E.8)where ζ1 = −tg¯1[1 − 2(tg¯1γ2)2]/[1 − 4(tg¯1γ2)2], ζ2 = −2tg¯1(tg¯1γ2)2/[1 −4(tg¯1γ2)2] and ζ3 = −tg¯1γ2(ζ1 + ζ2). Finally, using this in Eq. (4.7) resultsin the equation of motion for the GFG0,R(ω) = g¯0(ω)[δR,0 − t1(ω)∑δG0,R+δ(ω)−t2(ω)∑ξG0,R+ξ(ω)], (E.9)and its various coefficients are given in the text following Eq. (4.12).These effective hoppings and on-site energies are identical to those de-rived for the clean system in Ref. [2][85]. In the presence of disorder, thesolution proceeds similarly but now various g functions acquire dependenceon the location since their argument is shifted by U if R = 0. This leadsto dependence on location (and even direction of hopping) for the effectivehopping and on-site energies, at sites close enough to the impurity.130Appendix FDerivation of ARPESintensityIn ARPES, an incoming photon kicks out an electron from the sample andthis leaves a hole behind. For an energy transfer of ω, the intensity ofoutgoing electrons with momentum K can be written asI(K, ω) ∼∑n,k|〈K|A(r).p|n,k〉|2δ(ω − En(k)), (F.1)where A(r) = ǫA0eiq.r is the electromagnetic potential of the incomingphoton that results in scattering of electrons from polaron states |n,k〉 insidethe CuO2 layer to free electron state |K〉 outside. As the momentum of thephoton q is much smaller than that of the electrons inside the layer (λph ≫lattice spacing), the dipole approximation q ≈ 0 can be usedI(K, ω) ∼∑n,k|〈−Khole|ǫ.p|n,k〉|2δ(ω − En(k)), (F.2)and Khole = qph −K ≈ −K. This matrix element M = 〈−Khole|ǫ.p|n,k〉corresponds to the following integralM = −i~∫dr eiK.r ǫ.∇∑pcp(n,k)1√N∑Rpeik.Rpφp(r−Rp)= −i~∑pcp(n,k)1√N∑Rpeik.Rp∫dr eiK.r ǫ.∇φp(r−Rp) (F.3)in which φp(r−Rp) is the wavefunction of a 2p orbital located at Rp, andcp(n,k) relates |n,k〉 to the bare-hole eigenstates |p,k〉|n,k〉 =∑pcp(n,k)|p,k〉 + ... , (F.4)where ... includes contributions with at least one magnon. In Eq. (F.3), thelowest level of approximation is used where the terms describing the possi-bility for creation of magnons when the electron is removed are neglected.131Appendix F. Derivation of ARPES intensityThe integral is written aseiK.Rp∫dρ eiK.ρ ǫ.∇φp(ρ) (F.5)where ρ ≡ r −Rp . Using this in Eq. (F.3) and doing the summation overRp we getM = −i~√N∑pcp(n,k)∑Gδk+K||,G eiG.δp∫dρ eiK.ρ ǫ.∇φp(ρ). (F.6)Here, G’s are the reciprocal lattice vectors of the magnetic sublattice R,eiG.R = 1, and δp’s define the location of the four 2p orbitals inside the unitcell, Rp = R + δp. The integral is evaluated using φp(ρ) ∼ xpe−λρ wherexp = x(y) if p is an x(y) orbital, for examplemx ∼∫dρ eiK.ρ ǫ.∇(xe−λρ)= ǫx∫dρ eiK.ρ−λρ − λ∫dρ eiK.ρ−λρxǫ.ρ/ρ. (F.7)The relevant limit here is when the wavelength of the outgoing electrons 1/Kis much greater than the size of 2p orbitals, i.e. λ ≫ K. mx is dominatedby the first integral in this limit and it is given asmx ≈ ǫx∫dρ eiK.ρ−λρ = 8πλǫx(λ2 + K2)2 ≃8πλ3 ǫx. (F.8)ǫy and ǫz components are also given similarly. Hence,M = −i~√N∑pcp(n,k)∑Gδk+K||,G eiG.δp8πλ3 ǫp (F.9)and|M |2 = 64π2~2Nλ6∑Gδk+K||,G∑p,p′cp(n,k)c∗p′(n,k) eiG.(δp−δp′)ǫpǫp′ . (F.10)The intensity I(K, ω) is then written asI(K, ω) ∼∑G,kδk+K||,G∑p,p′eiG.(δp−δp′)ǫpǫp′−1π ImGp,p′(k, ω) (F.11)132Appendix F. Derivation of ARPES intensitywhere we have used the followingδ(ω − En(k)) = −1π Im〈n,k|Gˆ(ω)|n,k〉 ≈ −1π Im∑p,p′Gp,p′(k, ω)in which only the first term of Eq. (F.4) is retained as before. Eq. (5.27)is finally recovered by averaging Eq. (F.11) with respect to the polarizationof the incoming photon 〈ǫpǫp′〉 = 〈ǫ2p〉ηp,p′,I¯(K, ω) ∼∑G,kδk+K||,G∑p,p′eiG.(δp−δp′)ηp,p′Ap,p′(k, ω), (F.12)where ηp,p′ = 1 if p and p′ are of the same type (both px or both py) and zerootherwise. The approximation K ≪ λ is valid for small K, but for larger-K||00.080.160.24Zarpes (K||)λ = ∞λ = 2pi(0,pi) (0,0) (pi,pi) (0,pi) (pi,0)Figure F.1: The effect of finite K/λ (lowest order correction) on the ARPESquasiparticle weight of the five-band model. The difference between thequasiparticle dispersions is invisible on the same energy scale.momenta one has to keep the second integral in Eq. (F.7). Including theleading order correction to mα givesmα ≈8πλ3 ǫα +32πKαλ5 (Kβǫβ +Kγǫγ). (F.13)133Appendix F. Derivation of ARPES intensitywhere β and γ are the other two indices. This modifies the averages as〈mpmp′〉 ∼ [ηpp′ +8KpKp′λ2 (1 − ηpp′)]〈ǫ2p〉.This results in corrections that involves interference between 2p orbitals ofdifferent types and modifies the ARPES quasiparticle weight accordingly,Fig. F.134Appendix GSolving for Gα,β(k, ω) in thefive-band modelFor Gα,pax (k, ω), the second term inGα,pax (k, ω) = g0(ω)[δα,pax + 〈k, α, ↑|Gˆ(ω)H1|k, pax , ↑〉],becomes as follows〈k, α, ↑ |Gˆ(ω)H1|k, pax , ↑〉 = ts[−V 31 + V 11 − V 21 ] +Jpd2 V01 + ts[ei(kx−ky)/2×Gα,pcy (k, ω) − ei(kx+ky)/2Gα,pdy (k, ω)−eikxGα,pbx (k, ω)] + 2tpp[cos((kx + ky)/2)×Gα,pdy (k, ω) − cos((kx − ky)/2)Gα,pcy (k, ω)]−2t′pp cos(kx)Gα,pbx (k, ω) − 2t˜pp[cos((kx + ky)/2)×Gα,pdx (k, ω) + cos((kx − ky)/2)Gα,pcx (k, ω)].(G.1)Equation of motion for GFs where the hole starts from the rest of original2p orbitals is〈k, α, ↑ |Gˆ(ω)H1|k, pbx , ↑〉 = {...} − 2t˜pp[cos((kx + ky)/2)Gα,pcx (k, ω)+ cos((kx − ky)/2)Gα,pdx (k, ω)], (G.2)〈k, α, ↑ |Gˆ(ω)H1|k, pcy , ↑〉 = {...} − 2t˜pp[cos((kx + ky)/2)Gα,pby (k, ω)− cos((kx − ky)/2)Gα,pay (k, ω)], (G.3)〈k, α, ↑ |Gˆ(ω)H1|k, pdy , ↑〉 = {...} − 2t˜pp[cos((kx + ky)/2)Gα,pay (k, ω)+ cos((kx − ky)/2)Gα,pby (k, ω)], (G.4)where {...} refers to those terms that existed before introducing the neworbitals, as in Eq. (5.10). V il (implicit within {...}) are now one-magnon135Appendix G. Solving for Gα,β(k, ω) in the five-band modelGFs where the hole occupies one of the original 2p orbitals. In one-magnonapproximation, their equation of motion links to various Gαβ(k, ω), otherone-magnon GFs of the same type V and also to a new class of one-magnonGFs, U , where the hole occupies one of the new orbitals in the presence ofa magnon. For example,V 01g1(ω)= ts[−e−i(kx+ky)/2Gα,pdy (k, ω) + ei(ky−kx)/2Gα,pcy (k, ω)−e−ikxGα,pbx (k, ω)] +Jpd2 [Gα,pax (k, ω) − V01 ]+tpp[V 13 + V 31 − V 11 − V 113 ] − t′pp[V 03 + V 21 ]−t˜pp[U13 + U31 + U11 + U113 ]. (G.5)When the hole start from the new 2p orbitals, due to lack of hybridizationwith Cu, there is no Tˆswap and Hpd and the equations of motion are thereforesimpler:Gα,pay (k, ω)g˜0(ω)− δα,pay = 2tpp[cos((kx + ky)/2)Gα,pdx (k, ω)− cos((kx − ky)/2)Gα,pcx (k, ω)]−2t˜pp[cos((kx + ky)/2)Gα,pdy (k, ω)+ cos((kx − ky)/2)Gα,pcy (k, ω)]−2t˜′pp cos(ky)Gα,pby (k, ω), (G.6)and so on. Note that the doping hole does not inhibit the superexchangecoupling between its neighbouring Cu in this case as it does not occupy the2p orbital hybridizing with them, hence g˜0(ω) = 1/(ω + iη). Similarly, wehaveGα,pby (k, ω)g˜0(ω)− δα,pby = 2tpp[cos((kx + ky)/2)Gα,pcx (k, ω)− cos((kx − ky)/2)Gα,pdx (k, ω)]−2t˜pp[cos((kx + ky)/2)Gα,pcy (k, ω)+ cos((kx − ky)/2)Gα,pdy (k, ω)]−2t˜′pp cos(ky)Gα,pay (k, ω), (G.7)136Appendix G. Solving for Gα,β(k, ω) in the five-band modelGα,pcx (k, ω)g˜0(ω)− δα,pcx = 2tpp[cos((kx + ky)/2)Gα,pby (k, ω)− cos((kx − ky)/2)Gα,pay (k, ω)]−2t˜pp[cos((kx + ky)/2)Gα,pbx (k, ω)+ cos((kx − ky)/2)Gα,pax (k, ω)]−2t˜′pp cos(kx)Gα,pdx (k, ω), (G.8)Gα,pdx (k, ω)g˜0(ω)− δα,pdx = 2tpp[cos((kx + ky)/2)Gα,pay (k, ω)− cos((kx − ky)/2)Gα,pby (k, ω)]−2t˜pp[cos((kx + ky)/2)Gα,pax (k, ω)+ cos((kx − ky)/2)Gα,pbx (k, ω)]−2t˜′pp cos(kx)Gα,pcx (k, ω). (G.9)Equations of motion for U il link them to various other one-magnon GFs only,U01g˜1(ω)= tpp(U31 + U13 − U11 − U113 ) − t˜pp(V 31 + V 13+V 11 + V 113 ) − t˜′pp(U103 + U23 ), (G.10)U11g˜1(ω)= tpp(U21 + U23 − U01 − U43 ) − t˜pp(V 21 + V 23+V 01 + V 43 ) − t˜′pp(U13 + U53 ), (G.11)U21g˜1(ω)= tpp(U11 + U73 − U31 − U53 ) − t˜pp(V 11 + V 73+V 31 + V 53 ) − t˜′pp(U43 + U83 ). (G.12)U31g˜1(ω)= tpp(U01 + U83 − U21 − U103 ) − t˜pp(V 01 + V 83 (G.13)+V 21 + V 103 ) − t˜′pp(U73 + U113 ),where g˜1(ω) = 1/(ω − 2Jdd + iη) since all four AFM bonds are turnedferromagnetic. The equations of motion for one-magnon GFs at M ≥ 3 are137Appendix G. Solving for Gα,β(k, ω) in the five-band modelsimilar to those in the three-band model, Eq. (5.19), with equations for Vand U coupled together:λlVl = ξlVl−2 + ζlVl+2 + ξ′lUl−2 + ζ ′lUl+2 − λ′lUl,λ˜lUl = ξ˜lUl−2 + ζ˜lUl+2 + ξ′′l Vl−2 + ζ ′′l Vl+2 − λ′′lVl. (G.14)This can expressed as a single-index recursive equation for Zl =[VlUl]:ΛlZl = ΩlZl−2 + ΓlZl+2,where Λl =(λl λ′lλ˜l λ′′l)and so on. This results in a continued-fraction solu-tion relating V3 and U3 to both V1 and U1, similar to Eq. (5.20). Therest of the calculation proceeds as in Sec. 5.3.1, solving the remaining linearequations to yield all 64 sublattice Green’s function Gαβ(k, ω).138
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Variational studies of dressed quasiparticles' properties...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Variational studies of dressed quasiparticles' properties and their interactions with external potentials Ebrahimnejad Rahbari, S. Hadi 2014
pdf
Page Metadata
Item Metadata
Title | Variational studies of dressed quasiparticles' properties and their interactions with external potentials |
Creator |
Ebrahimnejad Rahbari, S. Hadi |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | In this thesis, we investigated the spectral properties of polaronic quasiparticles resulting from the coupling of a charge carrier to the bosonic excitations of ordered environments. Holstein polarons, describing an electron locally coupled to dispersionless phonons, and spin polarons in hole-doped antiferromagnets are considered. The Green's function of the Holstein polaron in a lattice with various kinds of impurities is calculated using the momentum average approximation. The main finding is that the scenario where the mere effect of the coupling to phonons is to enhance the electron's effective mass is incomplete as the phonons are found to also renormalize the impurity potential the polaron interacts with. This formalism is applied first to the case of a single impurity and the range of parameters where the polaron's ground state is localized is identified. The lifetime of the polaron due to scattering from weak but extended disorder is then studied and it is shown that the renormalization of the disorder potential leads to deviation of the strong coupling results from Fermi's golden rule's predictions. Motivated by the hole-doped cuprate superconductors, the motion of a hole in a two-dimensional Ising antiferromagnet and its binding to an attractive impurity is studied next, based on a variational scheme that allows for configurations with a certain maximum number of spin flips. The role of the magnetic sublattices in determining the symmetry of the resulting bound states is discussed. Next, a more realistic model describing a hole in a CuO₂ layer which retains the O explicitly, is considered. By neglecting the fluctuations of the Cu spins and using a variational principle similar to that of the previous chapter, a semi-analytic solution for the Green's function of the hole in an infinite 2D lattice is constructed. The resulting quasiparticle dispersion shows the proper ground state and other features observed in experiments. The lack of importance of the background spin fluctuations is justified based on the hole's ability to move on the O sublattice without disturbing the Cu spins. Finally, the model is generalized to gauge the importance of other relevant O orbitals. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2014-12-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0166096 |
URI | http://hdl.handle.net/2429/51589 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_february_ebrahimnejadrahbari_seyedhadi.pdf [ 2.67MB ]
- Metadata
- JSON: 24-1.0166096.json
- JSON-LD: 24-1.0166096-ld.json
- RDF/XML (Pretty): 24-1.0166096-rdf.xml
- RDF/JSON: 24-1.0166096-rdf.json
- Turtle: 24-1.0166096-turtle.txt
- N-Triples: 24-1.0166096-rdf-ntriples.txt
- Original Record: 24-1.0166096-source.json
- Full Text
- 24-1.0166096-fulltext.txt
- Citation
- 24-1.0166096.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0166096/manifest