Modeling Polarons in Transition-Metal Oxides by Bayo Lau B.A.Sc, University of British Columbia, 2004 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Physics) The University Of British Columbia (Vancouver) April 2011 c Bayo Lau, 2011 Abstract This work is a series of reports on progress in the description of electron-electron and electron-lattice interactions in transition-metal oxides, with an emphasis to the class of high-TC superconducting cuprates. Combinations of numerical and ana- lytical approaches were devised and developed to study large-scale models, which distinguish cation and anion sites of the realistic lattice structure. The many results range from incremental deviation to significant difference compared to those of the widely-accepted simple models without such distinction. A previously proposed numerical scheme and a perturbation approach were adapted to study the one-dimensional breathing-mode model, which describes a charge carrier interacting with vibrating anions. The effort yielded the first accurate benchmark result for all parameters values of the model. Based on a physical insight about the spin-1/2 Heisenberg antiferromagnet on a two-dimensional square lattice, an octapartite approach was devised to model the low-energy states. The efficiency of the resulting numerical approach was show- cased with the explicit solution to a record-breaking 64-spin torus. A spin-polaron model was derived to model holes injected into cuprate’s quasi two-dimensional copper-oxygen layer. Total-spin-resolved exact diagonalization was performed for a single fermionic hole in a record-breaking cluster with 32 copper and 64 oxygen sites. The solutions revealed important physics missed by previous studies. The octapartite approach was then developed to solve the spin-polaron model with two extra holes in the same cluster. The accuracy and efficiency of the method were established. Enhanced singlet correlation between two holes was observed. The preliminary results justify the need for an in-depth study. ii Preface A version of Chapter 2 has been published as B. Lau, M. Berciu, and G. A. Sawatzky. “Single-polaron properties of the one-dimensional breathing-mode Hamil- tonian.” Physical Review B 76 174305 (2007). The effort was initiated by Prof. M. Berciu and Prof. G. A. Sawatzky. I was responsible for all computational and theoretical work. Further physical interpretations and the writing of the manuscript were carried out with suggestions from the co-authors. A version of Chapter 3 has been published as B. Lau, M. Berciu, and G. A. Sawatzky. “Computational scheme for the spin- 12 Heisenberg antiferromag- net based on an octapartite description of the square lattice.” Physical Review B 81 172401 (2010). I was responsible for the initial effort as well as all computa- tional and theoretical work. Further physical interpretations and the writing of the manuscript were carried out with suggestions from the co-authors. A version of Chapter 4 has been published as B, Lau, M. Berciu, and G. A. Sawatzky. “High-Spin Polaron in Lightly Doped CuO2 Planes.” Physical Review Letter 106 036401 (2011). I was responsible for the initial effort as well as all com- putational and theoretical work. Further physical interpretations and the writing of the manuscript were carried out with suggestions from the co-authors. A manuscript for a refereed journal is being prepared from the results in Chap- ter 5. I was responsible for the initial effort as well as all computational and theoret- ical work. M. Berciu and G. A. Sawatzky have been involved in the interpretation of the results. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Glossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Hole-Doped Cuprates . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 The t-J Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Towards Realistic Oxide Structures . . . . . . . . . . . . . . . . . 11 1.4 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2 Lattice Polaron on the Oxide Chain . . . . . . . . . . . . . . . . . . 15 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 The Holstein and Breathing-Mode Models . . . . . . . . . . . . 16 2.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.1 Strong-Coupling Perturbation Results . . . . . . . . . . . 18 2.3.2 Matrix Computation . . . . . . . . . . . . . . . . . . . . 21 2.3.3 Green’s Function Computation . . . . . . . . . . . . . . . 23 iv 2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.1 Low-Energy States . . . . . . . . . . . . . . . . . . . . . 24 2.4.2 Spectral Function . . . . . . . . . . . . . . . . . . . . . . 33 2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3 Octapartite Description of the Antiferromagnetic Square Lattice . . 38 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 The spin- 12 Heisenberg Antiferromagnet on the Square Lattice . . 39 3.2.1 The Neel Picture . . . . . . . . . . . . . . . . . . . . . . 40 3.2.2 The Resonating Valence Bond Picture . . . . . . . . . . . 42 3.3 Octapartite Computational Scheme . . . . . . . . . . . . . . . . . 43 3.3.1 Physical Insights . . . . . . . . . . . . . . . . . . . . . . 43 3.3.2 The Octapartite Approach . . . . . . . . . . . . . . . . . 44 3.3.3 Numerical Implementation . . . . . . . . . . . . . . . . . 46 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4 Spin Polaron on the Copper-Oxygen Plane . . . . . . . . . . . . . . 52 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2 The Microscopic Model . . . . . . . . . . . . . . . . . . . . . . . 54 4.2.1 Derivation of the Spin Polaron Model . . . . . . . . . . . 54 4.2.2 Polaron Coherence and Other Considerations . . . . . . . 56 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.3.1 Energetics . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3.2 Wavefunction Analysis . . . . . . . . . . . . . . . . . . . 59 4.3.3 Further Experimental Implications . . . . . . . . . . . . . 62 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5 Two Holes in the Spin Polaron Model . . . . . . . . . . . . . . . . . 64 5.1 An Intuition about the Doped Antiferromagnet . . . . . . . . . . . 64 5.2 The Spin Polaron Model for Two Holes . . . . . . . . . . . . . . 67 5.3 Hilbert Space Formulation . . . . . . . . . . . . . . . . . . . . . 71 5.3.1 Translational and Total-Spin Symmetries . . . . . . . . . 72 5.3.2 Enumeration . . . . . . . . . . . . . . . . . . . . . . . . 74 v 5.4 Low-Energy Two-Hole Solutions . . . . . . . . . . . . . . . . . . 75 5.4.1 Energetics . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.4.2 Wavefunction Analysis . . . . . . . . . . . . . . . . . . . 78 5.5 Speculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2 Opportunities for Further Developments . . . . . . . . . . . . . . 93 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 A Software Integrity Checks . . . . . . . . . . . . . . . . . . . . . . . . 103 B Computation Development . . . . . . . . . . . . . . . . . . . . . . . 104 B.1 Software Optimizations . . . . . . . . . . . . . . . . . . . . . . . 104 B.2 Details of the Octapartite Approach . . . . . . . . . . . . . . . . 105 B.2.1 General Observations . . . . . . . . . . . . . . . . . . . . 106 B.2.2 Implementations . . . . . . . . . . . . . . . . . . . . . . 108 vi List of Tables Table 1.1 Types of electronic models for the CuO2 plane, their assump- tions, and the Fork space of n extra holes injected into a half- filled system with N unit cells. . . . . . . . . . . . . . . . . . 7 Table 2.1 Most-relaxed cut-off condition for the 1D breathing-mode po- laron Green’s function computation. . . . . . . . . . . . . . . 23 Table 2.2 an ;n+ vs. n+;n for the ground state . . . . . . . . . . . . . . 29 Table 2.3 an ;n+ vs. n+;n for the first-excited state . . . . . . . . . . . 29 Table 3.1 Size of the Cs = 14 subspace as compared to that of the com- monly used basis. . . . . . . . . . . . . . . . . . . . . . . . . 47 Table 5.1 (N=32) The S= 12 k = ( p 2 ; p 2 ) single-hole ground state’s weight in subspaces of particular SA N SB. Numbers are percentage adding up to 100%. Two lines are drawn to highlight theCS = 14 truncation which discards states with SA;SB < 6 and yield an energy within 3.6% of the exact value. . . . . . . . . . . . . . 66 Table 5.2 Single-hole eigenstates of HJpd . p † s creates an oxygen hole and the arrows in the ket indicate the spins of the two copper sites neighboring the oxygen hole. . . . . . . . . . . . . . . . . . . 81 vii List of Figures Figure 1.1 Phase diagram of hole-doped cuprates (Temperature vs electron- removal per unit-cell). . . . . . . . . . . . . . . . . . . . . 3 Figure 1.2 (top) A CuO2 layer embedded in certain chemical host.(bottom) Lattice structure of single-layer cuprate HgBa2CuO4. . . . 5 Figure 1.3 Two adjacent unit cells of the CuO2 plane. The orbitals kept in the 3-band model of Eq. (1.1) are shown, with white/shaded for positive/negative signs. The two e vectors (solid arrow) and the four d vectors (dashed arrow) are also shown. . . . . . 6 Figure 1.4 Mean-field phase diagram for the t-J model. . . . . . . . . 10 Figure 2.1 (a) GS energy and (b) GS qp weight as a function of the corre- sponding dimensionless coupling parameter. The dashed line corresponds to the Holstein model, while the breathing-mode results are shown by circles (line is a guide to the eye). Param- eters are t = 2;W= 1. . . . . . . . . . . . . . . . . . . . . . 25 Figure 2.2 (a) Energy of the first excited K = 0 state, measured from EGS, and (b) its qp weight. The dashed line corresponds to the Hol- stein model, while the breathing-mode results are shown by circles (line is a guide to the eye). Parameters are t = 2;W= 1. 26 Figure 2.3 P(0) for t = 2;W = 1 for the breathing-mode (symbols) and Holstein (line) models, respectively. The solid and dashed lines correspond to ground state and first excited state, respec- tively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 viii Figure 2.4 The k-dependent (a) energy and (b) qp weight of the GS (cir- cles) and first bound state (square) for the 1D breathing mode model, for t = 2;W= 1;lB = 1:07. . . . . . . . . . . . . . . 30 Figure 2.5 Ratio of effective polaron mass to that of the free electron. Cir- cles and squares show breathing-mode results for GS and first bound state, respectively. The other lines correspond to the Holstein model GS (full) and second bound state (dashed). Pa- rameters are t = 2;W= 1. . . . . . . . . . . . . . . . . . . . 31 Figure 2.6 GS qp weights for the breathing-mode (circles) and Holstein (squares) models. Full symbols correspond to a = 8. For com- parison purposes, the empty symbols show the a = 4 results of Fig. 2.1(b). (W= 1) . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 2.7 The spectral function of the 1D breathing-mode model, on lin- ear (left) and logarithmic (right) scales. The curves are shifted for better viewing, and correspond (top to bottom) to lB = 1:5;1:25;1:0 and 0:75. Vertical lines indicates EGS and EGS+ W. K = 0; t = 2;W= 1;h = 0:004W. . . . . . . . . . . . . . 33 Figure 2.8 A(k;w) vs w , for K=p=0; 0.25; 0.5; 0.75 and 1 (top to bottom) for intermediate coupling t = 2;W= 1;lB = 1:25. h=4=DE = 0:001W. The height of the spectral weight is plotted in linear scale on the left and logarithmic scale on the right. The two vertical lines indicate EGS and EGS+W. . . . . . . . . . . . . 34 Figure 3.1 Square lattice divided in octads (shaded areas). Sites posi- tioned similarly inside octads have the same label. Each is surrounded by neighbors with different labels. . . . . . . . . . 45 Figure 3.2 (a) Overlap jhGS;CsjGS;Cs+ 2N ij2 between ground-state wave functions corresponding to adjacent values ofCs forN= 16;32;64. Note that for N = 64, the overlap is already 95% forCs = 14 ; (b) Fractional change 1EGS dEGS dCs in ground-state energy correspond- ing to differentCs cutoffs. The lines are linear fits to the data. 48 ix Figure 3.3 GS energy computed at specific Cs values (symbols) and the extrapolation to Cs = 1 (solid lines) for N = 16;32;64. The horizontal dashed lines are the lowest values from Ref. In- set: expectation values of the non-commuting operators p hm2i and hSr Sr+( L2 ; L2 )i. . . . . . . . . . . . . . . . . . . . . . . . 49 Figure 4.1 (a) Two adjacent unit cells of the CuO2 plane. The orbitals kept in the 3-band model of Eq. (4.1) are shown, with white/shaded for positive/negative signs. The two e vectors (solid arrow) and the four d vectors (dashed arrow) are also shown. (b) Sketch of a virtual process of Tswap. . . . . . . . . . . . . . . . . . . 54 Figure 4.2 a) Energy and b) quasiparticle weight (bottom) for the lowest eigenstates with ST = 12 and 3 2 vs. momentum. Different sets are shifted so as to have the same GS energy. . . . . . . . . . 58 Figure 4.3 hCx(d ;D)i for the lowest energy state at (a) (p2 ; p2 )with ST = 12 , and (b) at (p ,p) with ST = 32 . The darkly-shaded bullet denotes the oxygen position at l+ ex. Each bullet shows the correla- tion value between the two sandwiching Cu sites. The cen- tral 12 Cu sites are shown; the correlations between the other 20 Cu spins converge fast towards the AFM value of -0.33. hCy(d ;a)i is the P̂x$y reflection. . . . . . . . . . . . . . . . . 60 Figure 4.4 hCx=y(d ;a)i for lowest state at k=(0,p) with ST=1/2. . . . . . . 62 Figure 5.1 Convergence of the one-hole ground state. (left) GS energy calculated for increasing CS. The value approaches rapidly to the exact value marked by the dotted line. (right) fractional change in GS energy for the next increment of Cs. Solid lines are extrapolation and linear fits. . . . . . . . . . . . . . . . . 66 Figure 5.2 Lowest two-hole states in the CS = 14 subspace. . . . . . . . . 76 Figure 5.3 Convergence of the ST = 0 lowest energy states. The two- hole energy level corresponding to DEb = 0 is indicated by the horizontal dash line. . . . . . . . . . . . . . . . . . . . . . . 77 x Figure 5.4 (top) Probability of charge separation, hCcharge(R)i for the low- est states in theCS = 12 subspace. (bottom) The difference from a random distribution. . . . . . . . . . . . . . . . . . . . . . 80 Figure 5.5 (top) Probability of singlet three-spin polaron (3SP) pair, hC3SP;singlet(R)i for the lowest states in theCS = 12 subspace. (bottom) The dif- ference from a randomly distributed paramagnetic configuration. 83 Figure 5.6 c3SP;singlet(d ) 104 for one of the two K = (p;p) GS in the CS = 12 subspace. Cu and O are marked by circles and rectan- gles, respectively. Values for the other state are obtained by 90o rotation. (top, black!red) 28 x-y, (bottom, black!black) 16 x-x and (bottom,red!red) 16 y-y 3SP pair configurations are shown; 6 O-O configurations, which involve the two O holes sharing the same copper site, have low probability and are omitted. For example, the bottom figure shows that c3SP;singlet(2a;0)= 0:0109 for two 3SPs both centered on x-rung O sites. . . . . . 84 Figure 5.7 c3SP;singlet(d )104 for the lowest K= (p;0) state in theCS = 12 subspace. Cu and O are marked by circles and rectangles, re- spectively. Values for the K = (0;p) state are obtained by 90o rotation. (top, black!red) 28 x-y, (bottom, black!black) 16 x-x and (bottom,red!red) 16 y-y 3SP pair configurations are shown; 6 O-O configurations, which involve the two O holes sharing the same copper site, have low probability and are omitted. For example, the bottom figure shows that c3SP;singlet(2a;0)= 0:0265 for two 3SPs both centered on x-rung O sites. . . . . . 86 Figure 5.8 ccharge(d ) 104 for the lowest K = (p;0) state in the CS = 12 subspace. Cu and O are marked by circles and rectangles, re- spectively. Values for the K = (0;p) state are obtained by 90o rotation. (top, black!red) 32 x-y, (bottom, black!black) 17 x-x and (bottom,red!red) 17 y-y 3SP pair configurations are shown. For example, the bottom figure shows that ccharge(2a;0)= 0:0289 for two 3SPs both centered on x-rung O sites. . . . . . 87 xi Figure 5.9 c3SP;singlet(d ) 104 for the lowest K = (0;0) ST = 0 state in the CS = 12 subspace. Cu and O are marked by circles and rectangles, respectively. (top, black!red) 28 x-y, (bottom, black!black) 16 x-x and (bottom,red!red) 16 y-y 3SP pair configurations are shown; 6 O-O configurations, which involve the two O holes sharing the same copper site, have low prob- ability and are omitted. For example, the bottom figure shows that c3SP;singlet(0; 2a) = 0:0367 for two 3SPs both centered on x-rung O sites. . . . . . . . . . . . . . . . . . . . . . . . . . 88 xii Glossary ZRS Zhang-Rice singlet AFM antiferromagnet ARPES Angular-resolved photoemission spectroscopy RVB resonating valence bond BCS Bardeen, Cooper, and Schrieffer 3SP three-spin polaron xiii Acknowledgments I would like to thank my supervisor, George. A. Sawatzky, who guided my growth as a scientist and a problem solver, particulary with many messages which seem trivial at first but proved to be very insightful and stimulating. I am grateful for him always making himself available, looking past my blunders, and most importantly, showing tremendous sympathy and support during the hardest time. I would also like to thank my PhD committee for making themselves available. In particular, M. Berciu has played a far greater role than required, as is evident in her co-authorship in my publications. I would like to thank her for the many insights and for going out of her way to support me during hardships. Before my graduate studies, A. Damascelli mentioned “One thing [I] can’t do is to say ‘It’s too hard!’ and give up.”, which has since become a motto of my life. While our research never overlapped, I. Elfimov has taught me many pecu- liar aspects about scientific computation and has also spent many extra hours in supporting the computing platform I used. M. Choptuik, M. Feeley, C. Greif, O. Toader, and A. Wagner have also provided computational advices for my work. The foundation for my research was laid by the excellent Canadian education system. I would like to acknowledge Kitsilano Secondary School’s teachers, Mr. Sarna and Ms. Kolson, without whom I would have never received a higher edu- cation. I would also like to praise the UBC Engineering Physics program, which provided me the multidisciplinary training that made this work possible. Last but not least, I would like to thank my family who has been supportive. My parents and my brother have always been there, praying for me, during the ups and the downs. My other half Kimberley has been undyingly supportive and has also tolerated many missed moments such as birthdays and valentines. xiv Chapter 1 Introduction When the temperature drops below the melting point of a material, atoms break the continuous spatial symmetries of the gaseous and liquid phase and condense into an energetically favorable solid state. Electrons are more than a thousand times lighter than the nuclei which form the lattice, so the electronic structure is a prominent characteristic. The most primitive theoretical description starts with the scenario of an electron placed in a periodic lattice potential with the potential of the inner core electrons screening part of the bare nuclear potential. For many materi- als such as normal metals and band insulators, a good description can be obtained by correcting, mostly by perturbation theory, the independent-particle picture with the effects of electron-electron and electron-lattice interactions[1, 2]; however, this way of treating many-body correlation cannot explain the properties of all mate- rials, and the community has put forth continuous effort in the understanding of the anomalous cases such as high-temperature superconductors[3–7] and colossal magneto resistors. The single-particle band theory predicts that the E k2 energy dispersion of an electron in free space is broken into energy bands with different E-vs-k dis- persion periodic in the Brillouin zone. The full many-body quantum state is then constructed by assigning the many electrons into these bands of single-particle lev- els according to Fermi statistics which forbids multiple electrons to have the same set of quantum numbers. If all bands are either empty or completely occupied, the material is predicted to be an insulator because electrons must excite to higher 1 energy band to acquire net momentum for charge flow and electrical conduction. If the valence band is only partially occupied, its electrons can acquire momentum by exciting to an unoccupied level in the same band with only an infinitesimal cost of energy. These low-lying unoccupied levels provide a channel for scattering due to electron-electron interaction, resulting in states different from the single-particle states in the band theory. The Fermi liquid theory considers such scattering and predicts a one-to-one correspondence between the original single electron state and the “quasiparticles” formed as a result of electron-electron interactions; that is, the single-particle description remains intact, but the states near the Fermi level are quasiparticles, or dressed electrons since they are accompanied by a polarization cloud, instead of free electrons. The Fermi liquid theory predicts such systems with strong electron-electron interaction to be metallic with a resistivity that scales as T 2 at low temperature. The electrons can also interact with phonons, quantized modes of lattice vibrations, which induces an effective attraction among the quasi- particles near the Fermi level. This introduces an instability against the Fermi sea, leading to a superconducting groundstate, from which excitation of certain fermion pairs of momenta k is characterized by a Bardeen, Cooper, and Schrieffer (BCS) energy gap D(k) = D0 which is constant. The preceding picture is taught in standard condensed matter courses, but such description of electronic structure is not universal. One of the most striking coun- terexamples is the class of cuprate compounds, whose properties contradict with all of the above statements. The undoped “parent” cuprate compound is an insula- tor with an odd number of electrons per unit cell. This is not the insulating scenario predicted by band theory, which requires an even number of electrons to fill the up- and down-spin states of a full band for insulation. Removing about one electron per five unit cells creates the “optimally doped” compounds, which exhibit a super- conducting transition temperature TC. Above TC, the materials were found to have a resistivity that scales as T instead of T 2 as predicted by the Fermi-liquid theory. Below TC, the system superconducts with a k-dependent BCS gap D(k) of dk2x k2y symmetry, which is different from the s-symmetric constant predicted due to electron-phonon interaction in conventional theory. For the “underdoped” com- pounds with electron concentration between that of the “parent” and “optimally doped” compounds, an extra excitation energy scale DPG > DSC has been measured 2 Figure 1.1: Phase diagram of hole-doped cuprates (Temperature vs electron- removal per unit-cell). by various experiments. This is another deviation from the phonon-driven super- conductor theory, which lacks an extra higher energy scale separating the normal and superconducting state. 1.1 Hole-Doped Cuprates Since their discovery in 1986 , cuprates have been classified as a high-TC super- conductor, but their importance to condensed matter theory lies in the many anoma- lous properties, including the aforementioned contradictions to the conventional theory, when the material is tuned away from the superconducting phase. The class of hole-doped cuprates, which allows electron removal from the parent com- pound, is of particular interest because of the existence of pseudogap phenomena, in addition to antiferromagnetism, superconductivity, Fermi liquid and non-Fermi liquid physics in different regions of the phase diagram as shown in Fig. 1.1. The undoped antiferromagnet (AFM) insulating phase is well understood in terms of the superexchange physics discussed below. The relatively low 300K magnetic or- dering temperature compared to a 3D antiferromagnet is due to the weak coupling 3 between the strongly correlated 2D AFM planes. The superconducting phase cen- tered at 0.15% doping is characterized by a BCS gap D(k) of d-symmetry, which is zero for kx = ky. Between the AFM phase and superconducting phase lies the region of pseudogap phenomena, which is a cross-over region instead of a definite phase. Spin-sensitive measurements showed a gapped excitation behavior but in- plane conductivity showed no gapped behavior. Angular-resolved photoemission spectroscopy (ARPES) measurements show excitations DPG(k) of the same trend as D(k) with maximum at (0;p) and zero for kx = ky. Breaking of spatial symme- tries have also been observed. The connection between different regions has not been elucidated. The mysteries stem from cuprates’ unique structure, which differs greatly from the three-dimensional isotropic lattice assumed by conventional theories. The com- mon feature among cuprates is the copper-oxygen layer embedded in different chemical hosts as shown in Fig. 1.2. One of the simplest examples is the HgBa2CuO4 compound with tetragonal crystal structure and one copper-oxygen layer per unit cell. Although in general TC is higher for complicated materials with multi- ple copper-oxygen planes and chains per unit cell, a full understanding of a single layer has not been achieved. While some compounds, such as La1 xSrxCaO4, exhibit lattice buckling at low temperature, HgBa2CuO4 does not exhibit lattice distortion. Lattice dis- tortion is not a universal property, and the layer structure as shown in Fig. 1.2 is believed to be an adequate representation of the system. For the undoped case, each Cu2+ ([Ar]d9) is surrounded by an octahedron of six O2 ([Ne]2s22p6). Four “in-plane” oxygens are linked to other Cu2+’s. The other two apical out-of-plane oxygens are further away from the copper than the in-plane ones are. Crystal field theory predicts the splitting of the ten 3d states into four energy levels with dif- ferent spatial symmetries: E(dxy) < E(dxz) = E(dyz) < E(dz2 r2) < E(dx2 y2). The nine d-electrons thus fully fill all but the dx2 y2 orbital. With one electron (or hole) per dx2 y2 orbital, this configuration is the so-called half-filled limit because the other four d-orbitals are commonly treated as inert. Doping is achieved by changing the chemical host outside of the CuO2 plane shown in Fig. 1.2. An electron removal does not simply introduce a hole in the dx2 y2 band because the oxygen 2p level is closer to the chemical potential than the 4 Figure 1.2: (top) A CuO2 layer embedded in certain chemi- cal host.(bottom) Lattice structure of single-layer cuprate HgBa2CuO4. 5 Figure 1.3: Two adjacent unit cells of the CuO2 plane. The orbitals kept in the 3-band model of Eq. (1.1) are shown, with white/shaded for pos- itive/negative signs. The two e vectors (solid arrow) and the four d vectors (dashed arrow) are also shown. copper dx2 y2 orbital. The material is in the charge-transfer regime in the Zaanen- Sawatzky-Allen classification scheme, and a hole would have strong character of the oxygen ligands. Although there are three orthogonal orientations per 2p oxygen level, the dx2 y2 wavefunction overlaps predominantly with the s orienta- tion of the four surrounding in-plane oxygens. The most widely studied model with one copper and two oxygen sites per unit cell is the p d model whose unit cell consists of one dx2 y2 and two 2ps orbitals[14, 15]. For copper sites located at l, some integral multiples of the lat- tice parameters ax = (a;0) and ay = (0;a), the oxygen sites are located at l+ ex=y, where ex=y = 12ax=y. Since there are more electrons than vacancies in the problem, the model is specified with fermion operators p†l+e;s and d † l;s creating holes of spin s on a vacuum with one hole per unit cell. For the orbital setting illustrated in Fig.1.3, the model reads H3B = Tpd +Tpp+Dpdånl+e;s +Uppånl+e;"nl+e;#+Uddånl;"nl;#; (1.1) 6 Model Approximations Number of Fermions 3-band dx2 y2 on Cu and two ps bands N+n 2-band dx2 y2 and one ps band N+n Ch. 3-5 hnd"+nd#i= 1 on each Cu and two ps bands N+n 1-band dx2 y2 band N+n t-J each ps hole is locked to a Cu spin as a ZRS N n Table 1.1: Types of electronic models for the CuO2 plane, their assumptions, and the Fork space of n extra holes injected into a half-filled system with N unit cells. with Tpd = tpdå( p†l+e;s + p†l e;s )dl;s +h:c: Tpp = tppåsd p†l+e+d ;s pl+e;s t 0ppå p†l e;s + p†l+3e;s pl+e;s nl;s = d † l;sdl;s nl+e;s = p † l+e;s pl+e;s : Tpd is the direct oxygen-copper hopping. Tpp contains the direct tpp oxygen-oxygen nearest-neighbor hopping as well as the t 0pp oxygen-oxygen next-nearest-neighbor hopping mediated by the copper 4s orbitals. The sign of the matrix elements are determined by the phases of the overlapping orbital wavefunctions in Fig.1.3a. In particular, for upper-right/lower-left (upper-left/lower-right) O-O hopping, sd = 1 (sd = 1). Dpd is the charge transfer energy, which is the energy cost of the d9 ! d10L ¯ Ligand hole excitation. Udd=pp are repulsive on-site Hubbard interac- tions, neglecting interaction of longer range. The magnitude of parameters have been established to be tpp < tpd < Dpd < Upp < Udd. In practice, neighbor Coulomb interactions such as Updd†dp†p should be considered; however, such considerations do not add any new matrix elements into the effective Hamiltonian of interest (Sec. 4.2.1). While ab initio density functional theory made important contributions in de- termining physical parameters, its applicability is based on groundstate energy and electron density, and the description of strongly correlated many-body physics would require models such as the aforementioned three-band model. In this 7 context, microscopic modeling has been carried out using models with one [5– 7, 16–19], two [20–28], three [14, 15], or more [29–32] bands. Table 1.1 gives a comparison of the models. While exact analytical solutions seem to be out of reach, numerical studies are also carried out with compromises such as the use of small clusters or variational approaches, classical Néel state plus spin-waves, etc [33–35]. Given these difficulties and the drive to find the simplest model, the t-J model is unsurprisingly the most studied [5–7]. 1.2 The t-J Model Recently, the two-dimensional t-J model has emerged as the most well understood microscopic model[5–7]. This route is based on the super-exchange physics and the formation of Zhang-Rice singlet (ZRS) upon electron removal. In the half-filled limit with no hopping, the groundstate of Eq. (1.1) has one hole per copper orbital. When hopping increases from zero, a pair of neighboring copper holes with opposite spins can hop onto the oxygen orbital shared by both sites, but such process is blocked by the Pauli principle if the two holes have the same spin. The net effect leads to superexchange, described by the Heisenberg quantum AFM coupling spins of neighboring copper sites of a 2D square lattice: HJ = Jå hi; ji SziS z j+ 1 2 S+i S j +S i S + j = Jå hi; ji Si S j (1.2) Although the 2D Hamiltonian has not been solved exactly, intensive studies over the years have shown that the undoped AFM insulating phase in Fig. 1.1 is well- described by the super-exchange mechanism. The ZRS mechanism was proposed to extend the undoped AFM description away from half filling. Although the lowest energy electron removal happens at the oxygen orbitals, Zhang and Rice pointed out that a specific linear combination of single oxygen hole states can form a singlet, the ZRS, with a copper hole to sig- nificantly lower local energy. The injection of an extra hole can be effectively described by removing a spin from HJ . It was shown that the resulting vacancy can propagate, which is equivalent to having the electron at the destination site hopping to the originally vacant site. The effective physics leads to the t-J model, which can 8 also be derived from theU ! ¥ limit of the one-band Hubbard model. HtJ = P̂ " å i; j ti j c†j;sci;s +h:c: # P̂+ Jå hi; ji Si S j (1.3) where the c†i;s creates a spin-s electron at site i and P̂ projects out states with double electron occupancy. In the literature, t-J model refers to the original case of ti; j = 0 for hopping range longer than the nearest neighbor. The community soon realized that further hopping is needed for better agreement with experiments, especially the dispersion observed by ARPES for the undoped parent insulating compound[38, 39]. The case with next-nearest-neighbor hopping is known as t-t’- J model, and further corrections in hopping are denoted as t-t’-t”-J and so forth. The t-J model is more manageable compared to other electronic models be- cause each additional hole removes a spin from the problem instead of adding one and because the effective Hamiltonian acts on a square lattice (Tab. 1.1). The (N+n)-fermion problem is reduced to a (N-n)-spin problem with a no-double- occupancy constraint. Equation 1.3 is deceptively simple, as it took many years of effort from the community to gain a satisfactory understanding of the model. Although exact analytical and numerical solutions are still lacking, the com- bination of mean-field treatment and variational Monte-Carlo calculation makes the t-J model the most understood microscopic model[5–7]. The most prominent mean-field description evolved from the connection between the resonating va- lence bond (RVB) ansatz to an undoped AFM and the BCS wavefunction. In the slave-boson approach to the t-J model, an electron creation is equivalent to a simultaneous creation of a fermionic spin and destruction of a bosonic vacancy c†i;s = f † i;sbi (1.4) if the restriction of f †i;" fi;"+ f † i;# fi;#+b † i bi = 1 (1.5) is enforced at each site. The boson b†i takes on the meaning of a ZRS, which is conveniently a S=0 singlet. The t-J model can thus be viewed as spinons f † and holons b† coupled via Eq. 1.3 plus the constraint in Eq. 1.5. Treating the constraint 9 Figure 1.4: Mean-field phase diagram for the t-J model. approximately, the mean-field theory yields the phase diagram shown in Fig. 1.4. Holons Bose condense below TB. Spinons acquire coherent motion below TD and pairs into singlets to form a RVB state with d-symmetry below TRVB. In the left region for TB < T < TRVB, the RVB state is formed without Bose condensation, and the pseudogap region in Fig. 1.1 is explained in terms of the formation of spin singlets. The superconducting region for T < TRVB;TB is explained as a RVB spin singlet state with Bose condensation of the charge carrying holons. To the right for TRVB < T < TB;TD, the spinons no longer form into singlets, and the remaining coherent motion leads to Fermi-liquid like behavior. In the region above the super- conducting area for TRVB;TB < T < TD, the non-Fermi-liquid region is attributed to spinons coherence without singlet formation. The highest TC predicted by mean- field theory is about a factor of two to four higher than the highest known TC for single-layer compounds[5, 6]. Remarkably, the non-Fermi-liquid T resistivity was found in accordance with experiment. The most significant difference between the mean-field phase diagram and the general one in Fig. 1.1 is the superconducting region which extends all the way to half-filling as well as the lack of separate AFM phase in that region. In this 10 low-doping regime, variational Monte-Carlo calculation showed that the RVB-type ansatz does not yield the lowest energy. A combination of RVB and AFM character in the variational ansatz can yield a lower energy state, signaling the coexistence of superconductivity and antiferromagnetism at low-doping and low temperature[41, 42]. Another inconsistency between the t-J model and experiment at low doping is that the one-hole solution predicts a band of coherent quasiparticles with sharp peak in the spectral weight; however, ARPES measurement at 300K reported a wide peak of 300meV in width, which is as broad as the band dispersion, decreasing lin- early with temperature. It has been suggested that coupling to phonons might be the cause of this broadening, and that phonons might be relevant to cuprate physics, but in a way different from conventional superconductor. This sparked the study of carrier-phonon interactions for models relevant to the cuprate lattice structure. For an electron-phonon coupling parameter much larger than that predicted by density functional theory, a wide peak broadening that decreases lin- early with temperature was found. However, the issue has not been elucidated because phonon theories predict that the linear decrease ceases at around 200K, but no experimental data is available below that temperature due to technical lim- itations. It was also pointed out that phonons can at least enhance spin-mediated superconductivity. 1.3 Towards Realistic Oxide Structures Holistically speaking, a full understanding of the physics of a CuO2 layer doped with a few holes has still not been achieved despite continuous effort[11, 46, 47]. In particular, the role of phonons[44, 48], time-reversal symmetry breaking, and spatial symmetry breaking have not been elucidated. Recent experiments have raised further questions regarding the understand- ing of cuprates, in addition to the ARPES line shape broadening [43, 51–54] dis- cussed in the last section. Neutron experiments have been performed in the pseu- dogap phase and reported magnetic response throughout the Brillouin zone, not restricted to the region of the much discussed magnetic resonance[55, 56]. The neutron results have been used to rule out one-band models such as the t-J model. 11 These and other issues including the broken local fourfold symmetry, which is taken for granted by the t-J model in Eq. 1.3, seen in scanning tunneling probe (STM) and x-ray scattering  remain either open questions or are contro- versial. While oxygen-bond-specific characteristics observed by XAS, EELS, and STM [10, 59, 60] cannot be described using one band, the significance of omitting other bands cannot be quantified without a comparison to unbiased solutions of more detailed models. Aside from adding more parameters such as long range hopping and carrier- phonon coupling for better agreement with experiments, a logical step is to study models more faithful to the realistic structure of the CuO2 plane to verify the main concepts of the t-J paradigm. For example, even though the RVB nature is com- patible with the BCS wavefunction, the robustness of spin-charge separation c†i;s = f † i;sbi, where the holon b † takes on the ZRS’s identity, should be verified in more detailed scenarios. In particular, cuprates exhibit charge-transfer band-gap behavior with mobile holes located mainly on anion ligands and unpaired electrons on cation d orbitals . One trade-off for the t-J model’s elegance is the use of momentum-independent effective parameters, even though it is well known that the ZRS state has a strong k-dependent renormalization . The impact of such approximations must be verified with models that distinguish anion and cation sites. Carrier-phonon interactions have been studied almost exclusively in the frame- work of the Holstein model. The Holstein coupling is compatible with the t-J model because it does not distinguish anion and cation sites. The effect of such distinction should be studied as well 1.4 Scope Given how the conventional understanding of electron-electron and electron-lattice interactions fails to transcribe into cuprate, and how the most developed micro- scopic description relies on significant assumptions, it is important to carefully examine these interactions in detail. The description of fermions interacting with other degrees of freedom is challenging, especially whenmodels retain the specifics of realistic material structures. One powerful concept is the formation of polaron— 12 a quasiparticle composed of the propagating carrier accompanied by some changes in the otherwise undisturbed background. While the quasiparticle description is certainly robust in the dilute limit, intuition about the many-body background is crucial in the actual solution process and is also the central theme of this disserta- tion. Physical insights about the many-body background were exploited to advance the treatment of models relevant to the transition metal oxide structures. Chapter 2 discusses the modeling of lattice polarons. A previously proposed numerical scheme is adapted to study the breathing-mode Hamiltonian which de- tails the interaction between a carrier and a background of vibrating oxygen ions. Results are compared with those of the Holstein model which is the most primitive starting point of carrier-phonon interaction. The “parent” cuprate compound is an AFM insulator (Eq. 1.2). The many elec- trons order antiferromagnetically with an energy scale of J 150meV . The AFM behavior is expected to be robust at low doping at low temperature, so an injected charge carrier can be thought of as a spin polaron of the carrier with a local dis- turbance to the AFM background. This spin polaron problem does not enjoy an exactly known background like the bosonic b j0i = 0 in a carrier-phonon inter- action scenario. Chapter 3 discusses progress in dealing with this issue. A novel octapartite numerical approach is devised to capture the two-dimensional Heisen- berg antiferromagnetic spin background with breakthrough efficiency. A surprising feature is revealed during the application process of the aforemen- tioned scheme to the doped copper-oxygen plane. Contrary to the common belief of ZRS formation and to the intuition that spin- 12 quasiparticles would result from a spin- 12 hole interacting with a spin-0 anti-ferromagnetic background, Chapter 4 shows that, when details are treated properly, spin- 12 and spin- 3 2 polarons become energetically favorable in different regions of the Brillouin zone. The chapter cov- ers the derivation of a spin polaron model relevant to cuprates, the solution process, and the nature of the low-energy wavefunctions. Chapter 5 studies two fermionic holes in the aforementioned spin-polaron frame- work. The computational approach from Ch. 3 is utilized to bypass various tech- nical barriers in modeling large systems. The result serves as a proof of concept, demonstrating the scheme’s prowess in capturing the AFM background disturbed by multiple holes. The results show enhanced singlet correlation. 13 Lastly, the conclusion is presented in Ch. 6. It summarizes the progress achieved in this series of works and presents various opportunities for further developments. Appendix A outlines the extensive efforts spent in ensuring the integrity of these novel numerical approaches. Computational optimizations that enabled these calculations are listed in Appendix B. 14 Chapter 2 Lattice Polaron on the Oxide Chain 2.1 Introduction In a solid-state system, the interaction between a charge carrier and phonons (quan- tized lattice vibrations) leads to the formation of polarons. This mechanism is a key ingredient in the physics of the manganites, Bechgaard salts, and, pos- sibly, of the cuprates. There are various model Hamiltonians describing the coupling of the particle and bosonic degrees of freedom. The asymptotic lim- its of weak or strong coupling can be investigated analytically using perturbation theory; however, the intermediate-coupling regime generally requires numerical simulations. Recently, investigations of basic model Hamiltonians have progressed rapidly thanks to the development of efficient analytical and computational tools, and we are now able to begin to study more and more realistic models. We studied numerically various low-energy properties and the spectral func- tion of the single polaron in the 1D breathing-mode Hamiltonian which captures some details of the transition metal oxide structure. The results are compared with the relevant results for the single Holstein 1D polaron, allowing us to contrast the behavior of the polarons in the two models. The chapter is organized as follows: Sec. 2.2 introduces the two models of interest. Sec. 2.3 discusses relevant asymp- totic results and describes the numerical methods we use to calculate low-energy 15 properties and the spectral functions for both models. Sec. 2.4 presents the results. Conclusions are made in Sec. 2.5. 2.2 The Holstein and Breathing-Mode Models The simplest electron-phonon coupling is described by the Holstein Hamiltonian. It is essentially the tight-binding model with an on-site energy proportional to the lattice displacement Xi = 1p2MW(b † i +bi): HH = t å <i j> c†i c j+ c † jci +Wå i b†i bi+gå i niXi (2.1) Here ci is the annihilation operator for an electron at site i (since we only consider the single electron case, the spin is irrelevant and we drop its index. We also set h̄ = 1). t is the hopping matrix, ni = c † i ci. For the Einstein phonons, bi is the annihilation operator at site i, W is the frequency, M is the atomic mass, and g is the electron-phonon coupling strength. The model has been widely studied numer- ically by Monte-Carlo calculations[65–76], variational methods[78–89], and exact diagonalization[90–95]. Analytic approximations, such as the momentum- average approximation, have also progressed over the years[96–99]. For some materials, a more appropriate model is provided by the breathing- mode Hamiltonian. For example, consider a half-filled 2D copper-oxygen plane of a parent cuprate compound. Injection of an additional hole should fill an oxygen 2p orbital. Due to hybridization between oxygen 2p and copper 3dx2 y2 orbitals, the hole resides in fact in a ZRS with a binding energy proportional to 8t2dp, where tdp is the hopping between neighboring O and Cu orbitals. The dynamics of the ZRS can be described by an effective one-band model with orbitals centered on the copper sub-lattice. If lattice vibrations are considered, the motion of the lighter oxygen ions, which live on the bonds connecting Cu sites, are the most relevant. The hopping integral tdp and charge-transfer gap between Cu and O orbitals are now modulated as the oxygen moves closer or further from its neighboring Cu atom. Both the on-site energy and hopping integral are modulated in the effective one-band model, but the former has been shown to be dominant. The breathing-mode Hamiltonian describes the physics of the linear modulation of on- 16 site energy. While this breathing-mode model is motivated as a 2D model, here we investi- gate numerically only its 1D version, relevant, e.g. for CuO chains. In 1D we can investigate accurately and efficiently not only ground-state (GS) properties, but also some excited state properties. For the Holstein model, it was found that po- laron properties are qualitatively similar in different dimensions,[83, 98] but with a sharper large-to-small polaron crossover in higher dimensions. We will show that a sharp crossover is already present in the 1D model, and we expect less of a dimensionality effects in the breathing-mode Hamiltonian. Quasi-1D systems, such as the Bechgaart salts, also involve electron-phonon coupling which is more complicated than the Holstein model; however, these sys- tems have a rather complicated structure and involve both strong electron-electron and electron-phonon interaction. The electron-phonon coupling also modulates the intermolecular hopping integrals in addition to the on-site energy. This adds a con- siderable degree of complication to the calculations and will be a subject of future studies. The 1D breathing-mode Hamiltonian that we investigate here is described by: HB = tå i c†i ci+1+ c † i+1ci +Wå i b† i+ 12 bi+ 12 + gp 2MWåi ni b† i+ 12 +bi+ 12 b † i 12 bi 12 (2.2) The notation is the same as before, except that now the phonons live on an inter- laced lattice. The difference between the two models is more apparent in momen- tum space. The Holstein model has constant coupling to all phonon modes VH = gp N p 2MWåkq c†k qck b†q+b q ; whereas the breathing-mode model has a coupling strength that increases mono- tonically with increasing phonon momentum VB = 2igp N p 2MWåkq sin q 2 c†k qck b†q+b q : 17 Here N is the number of lattice sites, and becomes infinite in the thermodynamic limit. The momenta k;q are restricted to the first Brillouin zone ( p;p] (we take the lattice constant a= 1). While numerical and analytical studies of the Holstein polaron abound, there is much less known about models with g(q) coupling. In particular, there is no detailed numerical study of the 1D breathing-mode model, apart from an exact diagonalization of a simplified t-J-breathing-mode model in restricted basis , and an investigation based on the self-consistent Born approximation , which is known to become inaccurate for intermediate and strong couplings. 2.3 Methodology 2.3.1 Strong-Coupling Perturbation Results Perturbational results for the strong-coupling limit g t provide a good intuitive picture of the problem even for the intermediate coupling regime. In the absence of hopping, t = 0, both Hamiltonians can be diagonalized by the Lang-Firsov transformation Õ= eSOe S: (2.3) Using SH = g W p 2MWåi ni(b † i bi); and SB = g W p 2MWåi ni( b†i 12 +bi 12 +b † i+ 12 bi+ 12 ) respectively, the diagonal forms of the Hamiltonians are, in terms of the original (undressed) operators: H̃H = T̃H +Wå i b†i bi g2 2MW2åi n2i (2.4) H̃B = T̃B+Wå i b† i+ 12 bi+ 12 2g2 2MW2åi ni(ni ni+1) (2.5) 18 where the kinetic energies are: T̃H = te g2 2MW3 å i c†i+1cie g( b†i+1+b † i ) W p 2MW e g( bi+1+bi) W p 2MW +h:c: T̃B = te 3 g2 2MW3 å i c†i+1cie g W p 2MW (b† i+ 32 2b† i+ 12 +b† i 12 ) e g W p 2MW (b i+ 12 2b i 12 +b i 32 ) +h:c: For a d-dimensional lattice, T̃B is modified by i) extra creation and annihilation operators of phonons in the direction transverse to hopping, and ii) change of the -3 factor in the exponent to (z+ 1);z = 2d. The third term in Eq. (2.4) and in Eq. (2.5) signifies that the mere presence of an electron would induce a lattice deformation, leading to the formation of a polaron to lower the system’s energy. For a single polaron, the lattice deformation energy is proportional to the number of nearest phonon sites (one for the Holstein model and z for the breathing-mode model). For t = 0, the ground-state energy is degenerate over momentum-space: E(0)H (k) = g2 2MW2 E(0)B (k) = z g2 2MW2 Each model has three energy scales, therefore the parameter space can be charac- terized by two dimensionless ratios. It is natural to define the dimensionless (ef- fective) coupling as the ratio of the lattice deformation energy to the free-electron ground-state energy zt: lH = g2 2MW2 1 zt (2.6) lB = g2 2MW2 1 t ; (2.7) where z is also the number of nearest neighbors in the electron sublattice. It should be noted that since W 1=pM, the l ’s do not depend on the ion mass, M. l has 19 been shown to be a good parameter to describe the large-to-small polaron crossover in the Holstein model. It will be shown in later sections that the definition also works well for the breathing-mode model. The other parameter is the adiabatic ratio which appears naturally from the perturbation in t: a = zt W : (2.8) Using standard perturbation theory, the first-order corrections to the energy of the lowest state of momentum k are: E(1)H (k) = 2te alH cos(k); (2.9) E(1)B (k) = 2te 3 2alB cos(k); (2.10) showing that the polaron bandwidth is exponentially suppressed in the strong cou- pling limit. As is well known, this is due to the many-phonon clouds created on the electron site (Holstein) or on the two phonon sites bracketing the electron site (breathing-mode model). As the polaron moves from one site to the next, the over- laps between the corresponding clouds become vanishingly small and therefore teff ! 0. To first order in t, the suppression is stronger for the breathing-mode model simply because the overlap integral involves phonon clouds on 3 sites in- stead of just 2, as in the case for Holstein. The second-order corrections are: E(2)H (k) = 2 t2 W e 2alH fH(k;a;lH); E(2)B (k) = 2 t2 W e 3alB fB(k;a;lB): The functions f can be written in the form fH;B = AH;B+BH;B cos(2k) with AH = Ei(2alH) ln(2alH) g BH = Ei(alH) ln(alH) g 20 and AB = Ei(3alB) ln(3alB) g BB = Ei(2alB) 2Ei(alB)+ ln(alB2 )+ g: Here, g is the Euler-Mascheroni constant and Ei(x) is the exponential integral with the series expansion Ei(x) = g+ ln(x)+ ¥ å n=1 xn n!n : The result can be further simplified in the limit alH ;alB 1 using limx >¥å¥n=1 x n n!n ex x . This leads to the simplified expressions: E(2)H (k) 2 t2 alH 1 2 + e alH cos(2k) (2.11) E(2)B (k) 2 t2 alB " 1 3 + e alB 2 cos(2k) # (2.12) Thus, the breathing-mode model’s ground state energy is slightly higher for any finite t. For the Holstein model, the dispersion is monotonic, since the second order cos(2k) contribution is suppressed more strongly than the first order cos(ka) contribution. However, a quick comparison between Eq. (2.10) and (2.12) shows that in the breathing-mode model, the second order cos(2k) contribution becomes dominant at large enough coupling. As a result, at strong couplings we expect the breathing-mode polaron energy to exhibit a maximum at a finite k< p , and then to fold back down. 2.3.2 Matrix Computation The computation method we use is a direct generalization of the method introduced by Bonca et al. for the Holstein model, in Ref. . This approach requires sparse matrix computations to solve the problem. It gives us a systematic way to compute excited state properties, which would be more difficult to achieve using other numerical methods. The idea is to use a suitable basis in which to represent the Hamiltonian as 21 a sparse matrix. For the Holstein model, this basis contains states of the general form jS;Ki= 1p Nåj eiK jc†j Õ mefSg b†nmj+mp nm! j0i; (2.13) where K is the total momentum and S denotes a particular phonon configuration, with sets of nm phonons located at a distance m away from the electron. For the Holstein model, m are integers, since the phonons are located on the same lattice as the electrons. The generalization for the breathing-mode is simple: here m are half-integers, since here the phonons live on the interlaced sublattice. All states in either basis can be obtained by repeatedly applying the Hamiltonian to the free electron state which has all nm = 0. The possible matrix elements are teiK , Wn, and gp 2MW p n, where n are integers related to the numbers of phonons. Since the Hilbert space of the problem is infinite, this basis must be truncated for computation. The original cut-off scheme in Ref.  was optimized for com- putation of ground-state properties of the Holstein model, by restricting the number of matrix elements between any included state and the free-electron state. Getting the higher-energy states is more involved, as it is evident from Eq. (2.3) that each state jfiLF in the Lang-Firsov basis correspond to a state jfiR = e SjfiLF (2.14) in real space. With our choice for the S operators, this reverse transformation in- duces a phonon coherent state structure at the electron site (Holstein), respectively the two bracketing phonon sites (breathing-mode). The phonon statistics of the co- herent state obeys the Poisson distribution. In the anti-adiabatic regime (zt > W), the splitting due to the hopping (off-diagonal hopping matrix elements) is signifi- cant compared to the diagonal matrix elements proportional to W. The underlying Lang-Firsov structure needs to be modeled by the hopping of the electron away from the coherent state structure created by the e S operator. To capture these characteristics, the basis is divided into subspaces with fixed numbers of phonons. Each subspace is enlarged by the addition of states with phonons further and fur- ther away from the electron site (increase of maximum value of m in Eq. (2.13)), until convergence is reached. This procedure allows for efficient generation of all 22 Subspace’s # phonons jmjmax # States 1-11 22.5 - (# phonons) 17053356 12-13 10.5 16871582 14-15 9.5 28274774 16-17 8.5 33423071 18-20 7.5 41757650 21-30 5.5 42628080 31-40 4.5 38004428 41-50 3.5 12857573 Table 2.1: Most-relaxed cut-off condition for the 1D breathing-mode polaron Green’s function computation. basis states required to model the higher-energy states. Matrices of dimension up to 106 107 were needed to compute the two lowest- energy states accurately. These two states were calculated numerically using the Lanczos method with QR shift, [105, 106] which works efficiently for the low energy bound states. For the larger n values used in the Green’s function calculation (see below), SMP machines were used for efficient computation. Table 2.1 lists the most-relaxed cut-off condition we used to calculate the Green’s function (see next section) for the 1D breathing-mode model. We relaxed the cut-off condition rather roughly until convergence was observed. As a result, the size of these matrices is certainly much larger than it has to be. 2.3.3 Green’s Function Computation Computation of higher-energy properties requires much larger matrices. The mem- ory and fflops needed for such computations are formidable, especially to exten- sively investigate the multi-dimensional parameter space (l ;a;K). Furthermore, one characteristic of the single polaron problem is a continuum of states starting at one phonon quantum above the K = 0 groundstate. Lanczos-type methods typ- ically are problematic in dealing with bands of eigenvalues with small separation. Therefore, in order to study higher energy states, we calculate directly the Green’s function: G(w;k) = hkj 1 w H+ ih jki; (2.15) 23 where jki = c†k j0i. This can be written as the solution of a linear system of equa- tions: G(w;k) = hkjyi (w H+ ih)jyi = jki One can iteratively tri-diagonalize H by the vanilla Lanczos process: H = QTQ† (w+ ih T )Q†jyi= Q†jki If the RHS is of the form [100:::]T , Crammer’s rule can be used to express G = hkjQQ†jyi as a continuous fraction in terms of the matrix elements of the tri- diagonal matrix (w + ih T ). In particular, this condition is achieved by pick- ing the initial Lanczos vector to be jki. This method is efficient because it does not require the complete solution of the linear system nor of the eigenvalue prob- lem. It is well known that this type of iterative process suffers from numerical instability, which leads to the loss of orthogonality in Q and incorrect eigenvalue multiplicity in T . We perform the vanilla Lanczos tri-diagonalization and re-orthogonalize each state in Q against the starting vector jki to validate the con- tinuous fraction expansion. Then, numerical errors may come from the fact that T may have the wrong eigenvalue multiplicity. However, this will not affect the location of poles in the spectral weight; that is, the eigenenergies are accurate. 2.4 Results 2.4.1 Low-Energy States The ground state energy EGS and quasiparticle (qp) weight Z0 = jhFGSjc†k=0j0ij2, where jFGSi is the ground-state eigenfunction, are shown in Fig. 2.1 for the 1D breathing-mode and Holstein models. For a fixed value of a , we see the expected crossover from a large polaron (at weak coupling lB;H) to a small polaron (at strong coupling lB;H), signaled by the collapse of the qp weight. 24 0 1 2 λΗ,λΒ -9 -8 -7 -6 -5 -4 EGS 0 1 2 λΗ,λΒ 0 0.2 0.4 0.6 0.8 1 Z0 (a) (b) Figure 2.1: (a) GS energy and (b) GS qp weight as a function of the cor- responding dimensionless coupling parameter. The dashed line cor- responds to the Holstein model, while the breathing-mode results are shown by circles (line is a guide to the eye). Parameters are t = 2;W= 1. The ground state energy of both models decreases monotonically with increas- ing coupling, but that of the Holstein polaron is lower. This is in agreement with the second order strong-coupling perturbation results in Eq. (2.11) and (2.12). Unlike the rather gradual decrease in the quasiparticle weight of the 1D Holstein polaron, the 1D breathing-mode polaron shows a large Z0 at weak couplings, followed by a much sharper collapse at the crossover near lB 1. The reason for the enhanced Z0 at weak couplings is straightforward to understand. Here, the wave-function is well described by a superposition of the free electron and electron-plus-one-phonon states. Given the conservation of the total polaron momentum K = 0 = k+ q and the large electron bandwidth t, states with high electron and phonon momenta have high energies and thus contribute little to the weak-coupling polaron ground-state. On the other hand, the coupling g(q) sin(q=2) to the low-energy states with low electron and phonon momenta is very small for the breathing-mode model. This explains the slower transfer of spectral weight at weak coupling for breathing-mode 25 0 1 2 λΗ,λΒ 0.5 0.6 0.7 0.8 0.9 1.0 1.1 E1-EGS 0 1 2 λΗ,λΒ 0 0.1 0.2 0.3 0.4 0.5 Z0 (a) (b) Figure 2.2: (a) Energy of the first excited K = 0 state, measured from EGS, and (b) its qp weight. The dashed line corresponds to the Holstein model, while the breathing-mode results are shown by circles (line is a guide to the eye). Parameters are t = 2;W= 1. versus the Holstein polaron. The energy (measured from EGS) and qp weight of the first excited K = 0 state are shown in Fig. 2.2. For both models, at weak-coupling this state is precisely at W above the ground-state energy, at the lower edge of the polaron-plus-one-phonon continuum. As the coupling increases above a critical value, a second bound state gets pushed below the continuum. This second bound state is absent in SCBA calculations. The separation between the two lowest energy states now first decreases and then increases back towards W as lH;B ! ¥. This behavior is well- known for the Holstein polaron. The breathing-mode polaron shows the same qualitative behavior. Note that below the critical coupling, the computed energy of the first excited state is slightly larger than W. The reason is a systematic error that can be reduced by increasing the number of one-phonon basis states, in order to better simulate the delocalized phonon that appears in this state. The qp weight of the first excited state is zero below the critical coupling, due to the crossing 26 0 0.5 1 1.5 2 λΗ , λΒ 0 0.2 0.4 0.6 0.8 1 < P( 0)> Figure 2.3: P(0) for t = 2;W= 1 for the breathing-mode (symbols) and Hol- stein (line) models, respectively. The solid and dashed lines correspond to ground state and first excited state, respectively. between on-site and off-site phonon states. The nature of these states is revealed by checking the locality of the phonon cloud. We define the projection operator P(K) = å Slocal jS;KihS;Kj where the summation is over all states with mH = 0 and mB =0:5 in Eq. (2.13). Comparison with Eq. (2.14) shows that this operator selects only basis states with phonons only on the electron site (Holstein) and only on the two phonon sites bracketing the electron site (breathing-mode). Note that the coherent state struc- ture is expressed as phonon creation operators in an exponential function, which includes states with zero number of phonons. Fig. 2.3 shows the expectation value of this operator for the two lowest eigenstates of both models. For both ground- states hP(0)i 1, indicating that here most phonons are nearest to the electron. However, at weak coupling the first excited state (which is here the band-edge of 27 the polaron + free phonon continuum) has hP(0)i ! 0, precisely because the free phonon can be anywhere in the system. When the second bound states form, hP(0)i becomes large, showing that phonons in these states are primarily localized near the electron. While there appears to be a crossing between the ground-state and first-excited state values, we emphasize that P(K) measures the locality of the phonon cloud, not its structure. For the breathing-mode model, these results suggest the possibility to describe them using the on-site coherent-state structure. That is, a Lang-Firsov state with n and n+ number of phonons excited to the left and right of the electron site, mapped to real space by Eq. (2.14). We note that we are no longer in the strong coupling regime and the transformation cannot be determined by g and W alone; therefore, we seek an effective transformation with S̃(D) = SBj gW=D The computed eigenstates jfi are projected into such structure an ;n+ by 1p Nål eiKlc†l ¥ å n ;n+=1 an ;n+ (b† l 12 )n p n ! (b† l+ 12 )n+ p n+! j0i = eS̃(D)P(K)jfi: Tables 2.2 and 2.3 show the results of such projections for the ground state and for the first excited bound state. It is clear that they can be rather well described as jGSi e S̃(t;g;W;K) 1p N ål e iKlc†l j0i, respectively j1ibound e S̃(t;g;W;K) 1pN ål eiKlc † l (e iqb† l 12 e iqb† l+ 12 )j0i for some phase q(t;g;W;K) needed to satisfy time-reversal symme- try. These states no longer have definite parity symmetry like those of the Holstein model. The symmetry is broken by the anti-symmetric coupling term in the model. If the free-electron component is non-zero for an eigenstate, its components with odd/even number of phonons should have odd/even parity. For increasing momen- tum, this description is valid as long as the excited state remains bound, with energy less than EGS+W. Fig. 2.4 shows momentum dependent energy and qp weight for the two low- est eigenstates of the breathing-mode model, for an intermediate coupling strength 28 t = 2 W= 1 g= 0 D= 0 n+n 0 1 2 3 0 1.0000 0.0000 0.0000 0.0000 1 0.0000 0.0000 0.0000 0.0000 2 0.0000 0.0000 0.0000 0.0000 3 0.0000 0.0000 0.0000 0.0000 t = 2 W= 1 g= 1:5 D= 1:05 n+n 0 1 2 3 0 0.9150 0.0584 0.1372 -0.0477 1 -0.0584 -0.1681 0.0520 -0.0536 2 0.1372 -0.0520 0.0549 -0.3300 3 0.0477 -0.0536 0.0330 -0.0255 t = 2 W= 1 g= 1:964 D= 1:964 n+n 0 1 2 3 0 0.9876 -0.0509 0.0082 -0.0025 1 0.0509 -0.0100 0.0034 -0.0020 2 0.0082 -0.0034 0.0022 -0.0019 3 0.0025 -0.0020 0.0019 -0.0019 Table 2.2: an ;n+ vs. n+;n for the ground state t = 2 W= 1 g= 1:5 D= 1:05 n+n 0 1 2 3 0 0.0228 -0.6515 0.0152 0.1151 1 0.6515 0.0062 0.1700 -0.0473 2 0.0152 -0.1700 0.0505 -0.0456 3 0.1151 -0.0473 0.0456 -0.0220 t = 2 W= 1 g= 1:964 D= 1:964 n+n 0 1 2 3 0 -0.0859 -0.6675 0.0630 -0.0164 1 0.6675 -0.0864 0.0246 -0.0120 2 0.0630 -0.0246 0.0134 -0.0099 3 0.0164 -0.0120 0.0099 -0.0089 Table 2.3: an ;n+ vs. n+;n for the first-excited state 29 0 0.2 0.4 0.6 0.8 1 k/pi -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 Ek/Ω 0 0.2 0.4 0.6 0.8 1 k/pi 0 0.1 0.2 0.3 0.4 Zk (a) (b) Figure 2.4: The k-dependent (a) energy and (b) qp weight of the GS (circles) and first bound state (square) for the 1D breathing mode model, for t = 2;W= 1;lB = 1:07. lB = 1:07. The polaron band has a maximum at k 0:4p , in qualitative agreement with the strong-coupling perturbation theory results. This behavior is not captured by SCBA which is only accurate for low coupling strength. For the Holstein polaron, the polaron dispersion is a monotonic function of momentum. Even though the qpweights remain moderately high at zero momentum, the weights col- lapse towards zero with increasing momentum, similar to the well-known Holstein case. This is due to the fact that at large total momentum, the significant contri- bution to the eigenstate comes from states with at least one or more phonons. The free electron state has a large energy for large momentum, and contributes very little to the lowest energy eigenstates, so indeed Z! 0. The effective masses for the two lowest eigenstates of both models are shown in Fig. 2.5, as a function of lH;B. These were calculated from the second deriva- tive of the energy at momentum K = 0. For the GS of both models, the effective mass increases monotonically with lH;B. At weak couplings, the breathing-mode polaron is lighter than the Holstein polaron. As already discussed, this is due to 30 0 0.5 1 1.5 2 λΗ,λΒ 1 10 100 1000 m/m0 Figure 2.5: Ratio of effective polaron mass to that of the free electron. Cir- cles and squares show breathing-mode results for GS and first bound state, respectively. The other lines correspond to the Holstein model GS (full) and second bound state (dashed). Parameters are t = 2;W= 1. the vanishingly weak coupling to low-momentum phonons. At strong coupling, however, the effective mass is larger for the breathing-mode polaron. This is in agreement with predictions of the strong-coupling perturbation theory, and results from the fact that the hopping of a breathing-mode polaron involves phonon clouds on 2z 1 phonon sites, whereas hopping of a Holstein polaron involves phonons at only two sites. The effective mass of the first excited state can only be defined once this state has split-off from the continuum. It has non-monotonic behavior, first decreasing and then increasing with increasing lH;B. This can be understood through the link of the effective mass and the qp weight. In terms of derivatives of the self-energy S(k;w), the effective mass m is given by: m m = 1 ¶S ¶w 1+ m h̄2 ¶ 2S ¶k2 1 ; 31 0 0.5 1 1.5 2 λΗ , λΒ 0 0.2 0.4 0.6 0.8 1 Zk=0 Figure 2.6: GS qp weights for the breathing-mode (circles) and Holstein (squares) models. Full symbols correspond to a = 8. For compari- son purposes, the empty symbols show the a = 4 results of Fig. 2.1(b). (W= 1) where derivatives are evaluated at K = 0 and at the corresponding eigenenergy. The first term is linked to the qp weight, Z = (1 ¶S¶w ) 1, so that m 1=Z. As shown in Fig. 2.2(b), the qp weight of the first excited state has non-monotonic behavior, leading to the non-monotonic behavior of the effective mass. All the results shown so far were for a = 4. For higher a (lowerW and/or larger t), the difference between the two models can be grasped from Fig. 2.6. Similar to the Holstein model, the large-to-small polaron cross-over occurs at lower l for increasing a[83, 98]. At weak and moderate coupling, the qpweight and the effec- tive mass (not shown) in the breathing-mode model are much less sensitive to an increase in a than is the case for the Holstein polaron. This suggests that breathing- mode polarons should be better charge carriers than the Holstein polarons in this regime. 32 -6 -4 -2 ω/Ω A(k=0,ω) -7 -6 -5 -4 -3 ω/Ω ln A(k=0,ω) Figure 2.7: The spectral function of the 1D breathing-mode model, on lin- ear (left) and logarithmic (right) scales. The curves are shifted for bet- ter viewing, and correspond (top to bottom) to lB = 1:5;1:25;1:0 and 0:75. Vertical lines indicates EGS and EGS+W. K = 0; t = 2;W= 1;h = 0:004W. 2.4.2 Spectral Function The spectral function is proportional to the imaginary part of the Green’s function: A(k;w) = 1 p ImG(k;w) (2.16) In terms of single electron eigenstates and eigenfunctions Hjni= Enjni, we obtain the Lehmann representation: A(k;w) =å n jhnjc†k j0ij2d (w En) Of course, since we use a finite small h in numerical calculations, the d -functions are replaced by Lorentzians of width h [see Eq. (2.15)]. We calculate the Green’s functions as discussed in the previous section. Fig. 2.7 shows the spectral function for zero momentum as a function of en- 33 -6 -4 -2 0 2 4 ω/Ω A(k=0,ω) -6 -5 ω/Ω ln A(k=0,ω) Figure 2.8: A(k;w) vs w , for K=p=0; 0.25; 0.5; 0.75 and 1 (top to bottom) for intermediate coupling t = 2;W= 1;lB = 1:25. h=4=DE = 0:001W. The height of the spectral weight is plotted in linear scale on the left and logarithmic scale on the right. The two vertical lines indicate EGS and EGS+W. ergy. Results corresponding to four different coupling strengths lB from near the crossover region are shown for the breathing-mode polaron. We note that there is always a continuum starting at one phonon quantum above the ground state energy. This is more clearly visible in the right panel, where the spectral weight is shown on a logarithmic scale, and vertical lines mark the position of the ground-state energy EGS, respectively of EGS+W. As the coupling lB increases, we see the appearance of the second bound state below the continuum. We only find at most one extra bound state in this energy range for all coupling strengths. As lB increases, the spectral weight of the first continuum decreases dramatically. Other bound states form above it, followed by higher energy continua whose weight is also systemati- cally suppressed. This is qualitatively similar to the behavior exhibited by Holstein polaron. Fig. 2.8 illustrates the momentum dependence of the breathing-mode polaron’s spectral weight. The results correspond to a coupling above the critical value, 34 where there is a second bound state. The majority of the spectral weight is trans- ferred to much higher energies as the momentum increases, and a broad feature develops at roughly the position of the free-electron energy for that momentum. This spectral weight transfer is also qualitatively similar to what is observed for Holstein polarons. Our results have a high-enough resolution to clearly show the continuum at EGS +W for all values of K. This is part of the kink-like structure reported in Ref. . The logarithmic plot clearly reveals a non-monotonic dis- persion of the ground-state like in Fig. 2.4, characteristic for the breathing-mode polaron. Fig. 2.8 shows only one peak located between the ground-state and the polaron- plus-one-phonon continuum, even though in fact we believe that there are more eigenstates within this region. We found, from eigenvalue computation, additional energy states below the continuum; however, computation of exact energy values requires prohibitively long computation time due to the clustering of eigenvalues. By observing the convergence behavior due to increasing basis size, we can con- clude that additional bound states do exist below the continuum. The lack of their contribution to the spectral function can be understood by the fact that the single particle Green’s function only contains information about eigenstates with finite qp weight, jhf jc†K j0ij2 > 0, see Eq. (2.15). These eigenstates must have com- ponents corresponding to some Lang-Firsov eigenstate with no off-site phonons (Eq. (2.14)). Also, the wave-function of these states must have a peculiar space inversion symmetry: S-symmetric for all even-phonon-number components and P-symmetric for all odd-phonon-number components. The ground-state always satisfies this requirement, but above a critical coupling, only one other state below the phonon threshold satisfies this requirement. 2.5 Conclusions In summary, we have reported here the first accurate numerical study of the 1D breathing-mode polaron. A previous study based on the self-consistent Born approximation proves to be inadequate to describe correctly the behavior for medium and strong couplings, as expected on general grounds. Comparison with the Holstein model results, which correspond to a coupling 35 g(q) = const:, reveals some of the similarities and differences of the two models. The breathing-mode polaron is much more robust (has a much larger qp weight, and less variation with parameters) at weak couplings. This is a direct consequence of the fact that coupling to low-momentum phonons, which is relevant here, be- comes vanishingly small g(q) sin(q=2)! 0. Similar behavior is expected for any other g(q) model if limq!0 g(q)! 0. On the other hand, at strong couplings the breathing-mode polaron is much heavier and has a lower qp weight than the Holstein polaron. This also results from strong-coupling perturbation results, and is due to the fact that in order to move from site i to site i+ 1, a small breathing- mode polaron must (i) create a new polaron cloud at site i+ 32 ; (ii) rearrange the polaron cloud at site i+ 12 , so that its displacement is now pointing towards site i+1, not towards site i; and (iii) relax (remove) the phonon cloud at site i 12 . This results in a suppressed polaron kinetic energy, and an enhanced effective mass. Because of the larger Z at weak coupling, and the lower Z at strong couplings, the crossover from large to small polaron is much sharper for the breathing-mode polaron. Another interesting observation is that the polaron dispersion becomes non-monotonic with momentum k for medium and large couplings. This can be understood in terms of strong-coupling perturbation theory, which shows that the second-order cos(2k) correction is larger than the first-order cos(k) correction for large-enough couplings. Similarities with the Holstein behavior regard the appearance of the polaron- plus-free-phonon continuum at EGS +W, and the appearance of a second bound state with finite qp weight for large-enough couplings. The convergence of numer- ics points to the existence of additional bound states, whose absence from the spec- tral function can be explained by symmetry or missing free-electron components in the wavefunction; however, this issue is not fully settled. Also, the importance of such states for the physical properties is not known. The general aspect of the higher energy spectral weight at strong couplings, as a succession of bound states with large spectral weight and continua with less and less spectral weight, is also reminiscent of the Holstein polaron results. This chapter considers the propagation of a single fermionic carrier interacting with a many-body phonon background. The effect of the phonon background on the realistic strongly-correlated many-fermion scenario (Ch. 1) can be different. 36 As will be shown in the rest of the thesis, treating the details of the fermionic part already presents significant challenges and leads to remarkable properties. The following chapters are devoted into the modeling the strongly-correlated fermionic part, leaving the many-fermion-many-phonon scenario for future exploration. 37 Chapter 3 Octapartite Description of the Antiferromagnetic Square Lattice 3.1 Introduction The d-wave superconducting gap in cuprates is inconsistent with the s-wave gap due solely to electron-phonon interactions in normal metals. Because of the belief that the dominant physics lies in the two-dimensional copper-oxygen plane and be- cause such planes are antiferromagnetic (AFM) in the undoped limit, much effort has been spent on understanding of charged fermions moving in a two-dimensional background of antiferromagnetically coupled spins. In spite of recent efforts that consider carriers interacting with both phonon- and spin-degrees of freedom, an exact analytical solution for the spin part of the problem — doped or not — re- mains elusive to the best of my knowledge. This adds another layer to the puzzle as compared to the carrier-phonon problem, which takes for granted an exactly known bosonic background, b j0i = 0 and b+jni =pn+1jn+ 1i. Numerical modeling of large-scale, strongly correlated spin systems is thus of extreme importance and has been pursued in many different ways, each with their own advantages and dis- advantages. For example, the powerful Monte Carlo (MC) methods give extremely accu- rate information for the undoped AFM, providing the benchmark for the ground state (GS) energies and correlation functions of systems with up to thousands 38 of spins[110–114]. However, they are limited by the sign problem if additional fermions are present. This led to work on alternative optimized ways of treating the undoped AFM, with the latest development in MC sampling being the use of vari- ational RVB-type ansätze[5, 40–42, 115–117]. In the same context, density matrix renormalization group (DMRG) has also progressed over the years[118, 119]. Since the interactions with the fermions are often comparable or bigger than the AFM’s energy scale, most such AFM-specialized methods require very non-trivial modifications to adapt to doped systems, due to technicalities of the RVB ansatz or the DMRG special boundary conditions. In fact, recent studies still perform exact diagonalization (ED) of the full Hilbert space to complement these some- what biased schemes[17, 120, 121]. ED studies have the huge advantage that they provide the GS wavefunction, from which any GS properties, including all corre- lation functions, can be calculated. This would make ED the numerical method of choice, were it not for the extreme restriction on the sizes of doped systems that can be currently treated[17, 121]. The motivation behind this chapter is the belief that some additional physical insights about the AFM background can help in the modeling of the aforemen- tioned unresolved issues. The accuracy and efficiency of the phonon basis, Eq. 2.13 hinges on the insights about the phonon coherent state background Eq. 2.14. Al- though Bonca et. al. did successfully adapt the idea for a carrier in the spin back- ground , the method has left something to be desired because its Néel starting point cannot capture a real AFM background. The strong quantum fluctuations leading to strong deviations from a Néel-like background are believed to be an essential part of the low-energy physics of the doped systems. It is therefore nec- essary to find accurate yet efficient ways to describe the AFM background, which can also be straightforwardly extended to the doped problem. 3.2 The spin-12 Heisenberg Antiferromagnet on the Square Lattice The dominant physics in an undoped copper-oxide plane is due to the superex- change mechanism, which leads to AFM Heisenberg exchange among nearest neighbor S= 12 copper spins arranged in the square lattice. 39 HAFM = å hi; ji SziS z j+ 1 2 S+i S j +S i S + j = å hi; ji Si S j: (3.1) Here, a 2N4 constant is omitted and the characteristic energy scale Jdd 150meV has been set to unity. The “natural” basis has a dimension of 2N and a wavefunction can be written as a superposition of different spin configurations with jslil specifying the z-projection of spin at the copper site l. å s csÕ l jslil (3.2) The total spin ST 2 [0; N2 ] is conserved with degenerate z-projections so the full spectrum can thus be obtained from the SzT = ål sl = 0 basis, which contains N!=(N2 !) 2 states without exploiting spatial symmetries. Efficient matrix-vector multiplication can be performed using Algorithm 151  as discussed in Ap- pendix B. Although the exact ground state wavefunction remains unknown, Marshall showed very early that the groundstate wavefunction must have ST = 0 and that its cs in Eq. 3.2 obeys cs = ( 1)ns jcs j (3.3) where ns is the number of up spins in one sub-lattice. This is now known as the “Marshall sign rule” in the community. Among the many routes in treating an undoped AFM, two intuitions have been prominently linked to the doped case. 3.2.1 The Neel Picture The model is bi-partite because the square lattice can be divided into two interlaced sub-lattices such that spins within each sub-lattice are not coupled to each other by HAFM . The bi-partite nature and the dot-product interpretation in Eq. 3.1 motivated the popular “up-down-up-down” classical Neel state starting point. In this picture, the spin-wave theory justified in the S 1 limit worked surprisingly well for S = 1 2 . Starting from a Neel background, the z-projection of a spin at site la in the spin-up or at lb in the spin-down sublattice is mapped to bosonic number operators. 40 nla = a + laala = S Szla nlb = b + lbblb = S+S z lb (3.4) Eq. 3.1 can then be expanded in powers of nl2S due to the p n matrix elements of bosonic operators. The collection of O(1) zeroth-order terms is known as the linear spin-wave theory, which can be diagonalized as HLSW = 0:658N+2å q r 1 1 4 (cos(qx)+ cos(qy)) 2 a+q aq+b+q bq (3.5) where a+ and b+ are S=1 “magnon” excitations related to a+ and b+ by a unitary transformation. The prediction of gapless excitation is correct. The groundstate energy-per-site of -0.658 is a bit higher than the best known value of -0.6692. The magnon bandwidth is grossly underestimated to 2 compared to the actual bandwidth of 2:5. The next order 1/S correction in spin-wave theory lowers the groundstate energy to -0.6705, increases the bandwidth to 2.36 but also in- troduces three magnon-magnon interaction terms: a+p apa+q aq, b+p bpb+q bq, and a+p apb+q bq. Even though the results are excellent, the major challenge against the spin-wave theory is that the S 1 bosonic representation in Eq. 3.4 simply fails for a S = 12 system because the Hilbert space allows only zero or one occupation. For the undoped case, this occupation constraint can be treated by less-transparent bi-partite methods which cannot be easily generalized to the doped case. The classical Neel state has a finite overlap with some low-energy states of the model, and correcting a Neel state starting point, for example with off-diagonal terms in Eq 3.1, can eventually lead to a decent approximation of the low-energy sector; however, strictly speaking, the zero-temperature ground state do not have a finite Neel characteristic. After all, each S S term acting on a Neel state yields 14 as the diagonal matrix element but 12 as the off-diagonal matrix element. An upper-bound of the overlap between the classical Neel state and the true ground state wavefunction is derived here. In the bipartite picture, the Neel state has two sub-lattices of opposite spins 41 SA = N4 ;SzA =N4 SB = N4 ;SzB =N4 : (3.6) Knowing that the ground state has ST = 0, one can immediately conclude that the Neel state overlaps with only one ST = 0 basis state from the Clebsch-Gordon coefficients: hSA;mAjhSB;mBjS= 0;Sz = 0i= ( 1) SA mA p 2SA+1 dSA;SBdmA; mB (3.7) The ground state’s weight in a Neel state is thus at most 1N=2+1 which vanishes in the N ! ¥ limit. Even though the solution works out in the undoped case, when deviations from the Neel background are equal everywhere, one should be careful when injecting extra fermions into the otherwise homogeneous system. To be specific, the idea of modeling “local flipped spins” without correcting all spins from the Neel background  can yield an incomplete solution because a few holes can hardly change the low-range behavior in the N! ¥ limit. For example, the high-spin polarons in Chapter 4 have not been reported until we performed calculation with specific ST , which cannot be properly differentiated if some spins remain uncorrected from the Neel background. 3.2.2 The Resonating Valence Bond Picture Anderson has pointed out a vastly different starting point in dealing with Eq. 3.1 by the class of resonating valence bond (RVB) singlet liquid wave functions. The original observation was that a superposition of states with N2 pairs of singlets can be energetically competitive in low-dimensional and/or frustrated systems. The crucial idea is that a RVB state can be constructed by first removing all N spins from the lattice, and then projecting a mean-field, BCS wavefunction back to the lattice with the N2 pairs and no-double-occupancy constrains. jRVBi=Õ l (1 nl"nl#)PN=2jBCSi (3.8) Over the years, the RVB framework has been implemented for an AFM torus by two approaches. The bottom-up route consists of the explicit optimization of 42 the bond amplitudes. The top-down approach is by energy minimiza- tion by tweaking the unprojected wavefunction jBCSi in Eq. 3.8 to other related ansatz. Because the optimal variational ansatz is determined as the mean-field solution to some Hamiltonian with optimal values for parameters such as superconducting, antiferromagnetic, and density-wave orders, the real strength of this route is to test a priori postulations about the nature of the solution. The RVB starting point can indeed lead to the correct ground state wavefunction, but the projective, variational nature of the solution procedure hampers its neutrality and transparency. 3.3 Octapartite Computational Scheme One message from section 3.2 is that, in the lightly-doped scenario, the Neel start- ing point eases the description of local interactions at the cost of accuracy in dealing with the spin background, whereas RVB excels in implementing holistical insights with little details about the local picture. Here, I discuss how a physical insight can be exploited to bridge the gap between these two starting points. 3.3.1 Physical Insights Knowing that the total spin of the ground state must be zero, one can start with the bipartite division and write down all relevant basis states as Clebsch-Gordon series jS= 0;Sz = 0i=å ( 1) SA mA p 2SA+1 dSA;SBdmA; mB jSA;mAijSB;mBi (3.9) for SA=B 2 [0; N4 ];mA=B 2 [ SA=B;SA=B]. Note that the Marshall sign rule Eq. 3.3 is automatically satisfied so all of Marshall’s observations have been exploited. The Hilbert space contains a huge number of such singlets. They can be labeled by an index a , and the GS wave function can be expressed as a superposition jGSi=å a ca j0;0ia (3.10) I then noticed one characteristic shared by all a spin states with high values of 43 jca j. This can be illustrated by the staggered magnetization, which is a well-studied observable of the ground-state wave function. For a bi-partite lattice in a staggered magnetic field, this quantity is formally defined as the difference between hmAi and hmBi in the zero-field limit. In the absence of a staggered field, this quantity has alternative definitions, such as: m̂2 = 1 N2 å r ( 1)jrjŜr 2 = ŜA ŜB 2 N2 ; (3.11) where ŜA=B are the total sublattice spin operators. In the AFM GS, m = hm̂2i 12 has been extrapolated to 0:3 as N ! ¥, and increases like 1=pN to 0:45 for N = 32. In terms of the total spin Ŝ= ŜA+ ŜB, we have: m̂2 = 1 N2 (2Ŝ2A+2Ŝ 2 B Ŝ2) (3.12) The ground state is a singlet, S= 0. It follows that the wave function must allocate significant weight to sectors with high values of SA = SB in order to yield such large values of m. In fact, for the RHS of Eq. 3.12 to be larger than the expected value, the sublattice spins SA=B have to be within N16 of their maximum values. The physical meaning is that N8 spins add to a total spin of zero while the rest adds to the maximum 3 N16 . 3.3.2 The Octapartite Approach There are many ways to add spins quantum mechanically, but not all are viable because the truncation error needs to be bounded systematically and the compu- tational effort must be small. There are two extremes: if the original single-site basis of Eq. 3.1 is used, no basis transformation is needed but enforcing truncation based on Eq. 3.12 is costly. However, if a random basis tabulated according to values of SA is used, transforming Eq. 3.1 into the new basis would be costly due to the many Clebsch-Gordon series needed. We propose a parametrization which is a good compromise between these two extremes. Starting from the case where all N2 sublattice spins add to the maximum of N 4 , if we want to include states with spin down to 3N 16 = N 4 N16 , the basis must 44 Figure 3.1: Square lattice divided in octads (shaded areas). Sites positioned similarly inside octads have the same label. Each is surrounded by neighbors with different labels. allow many new configurations, including ones where N8 spins have a total spin of 0 while the other 3N8 spins have a total spin of 3N 16 . This suggests that groups of N 8 or fewer spins must be allowed to take all possible spin quantum numbers. Since a larger group would overshoot the sublattice spin below the 3N16 “threshold”, while a smaller group would introduce extra “enforcement” costs as discussed above, we divide the full lattice into groups of 8 sites, see Fig. 3.1. The resulting octads repeat periodically with translational vectors 2a(1;1). Each spin is identified by the octad it belongs to, and by its position inside the octad. A sublattice is composed from 4 groups of spins indexed with the same label, e.g. ŜA = Ŝ0+ Ŝ1+ Ŝ2+ Ŝ3. With this arrangement, each spin interacts with spins from all the 4 groups of the other sublattice, e.g. a group 0 spin is always to the west/east/south/north of a spin in group 4/5/6/7. This partition is the minimum division that permits such a description and is a fundamental building block for the wave function in the model. Even though we arrived at it by looking to optimize the computation, our partition has been used in other contexts. 45 We can now generalize Eq. 3.9 and write singlets of the total lattice as: j0;0; f (fSig)i=åcfSi;mig 7 Õ i=0 jSi;mii (3.13) where Si;mi are the quantum numbers for Ŝi, i.e. the total spin for group i = 0;7. Here, cfSi;mig is the product of the appropriate Clebsch-Gordon coefficients for the 8 pairs of quantum numbers (Si;mi). The f (fSig) index on the LHS keeps track of the many ways in which these spins can be added to form a singlet. If all these singlets are kept into the computational basis, the GS calculation is exact. The total spin of each group takes values Si 2 [0; N16 ], and there is no reason to restrict them. However, our previous discussion suggests that while all possible Si values must be allowed, the singlets with highest GS weight are those whose total sublattice spin is within N16 of the maximum value. We therefore parameterize the singlet Hilbert space in terms of a completeness parameter,Cs 2 [0;1], and include in the singlet basis only states for which SA=B N 4 (1 Cs). For Cs = 1, the calculation is thus exact. For Cs = 14 , this means that the maximum number of anti-aligned spins in the sublattice is N8 . Based on the discussion for the staggered magnetization, we expect this to already be a good variational basis for the GS. For any Cs 6= 1, this formulation allows all possible values of Si for each i = 0;7, but restricts the ways in which they can combine to give the total sublattice spin. The number of discarded states is combinatorially large due to the large degeneracy of states with low sublattice spins. Equation 3.13 is invariant under spin rotations; thus, the anti-alignment is en- forced here in terms of angular momentum addition, very different from the clas- sical Néel picture. This approach also maintains the full translational symmetry of the Hamiltonian. If the RVB bond-amplitude optimization is a bottom-up way of adding AFM order to a singlet[115, 116], our method is a compatible top-down approach without the need of projecting an ansatz. 3.3.3 Numerical Implementation The cost of this truncation scheme is the one-time computation of matrix elements between these basis states. This is acceptable because conventional diagonalization 46 N Cs = 14 S= 0 Sz = 0 Full 16 50 1430 12870 216 32 11042 35357670 601080390 232 64 9:8108 5:551016 1:831018 264 Table 3.1: Size of theCs = 14 subspace as compared to that of the commonly used basis. schemes are limited by storage, not processor speed, and most importantly, the time for matrix element computations was less than the amount of time required to develop the actual numerical solver. The octapartite division of Fig. 3.1 is valuable for greatly minimizing the num- ber of non-zero matrix elements of the Hamiltonian in the basis states of Eq. 3.13. In particular, the unique neighboring property dictates that each Si S j term in Eq. 3.1 changes only 2 out of 8 kets in Eq. 3.13. This relatively small change in the wavefunction can by exploited by building the N-spin S=0 basis by mixing two identical N2 -spin subspaces, which are in turn built from two identical N 4 -spin sub- spaces, which are formed by two identical N8 -spin subspaces. The matrix element, hajSi S jja 0i, between two N-spin configurations a and a 0 can be calculated from the Szl , S + l , and S l matrix elements within the N 8 subspace and the Clebsch-Gordon coefficients in the mix hierarchy. The viability of the computation hinges on the fact that Szl , S + l , and S l yield extremely sparse matrices among the N 8 subspace. This efficiency was complemented by tracing the spin-mix hierarchy according to the delta functions in the Clebsch-Gordon coefficients. While the full Sz = 0 basis can be enumerated by, for example, algorithm 151, enumerating states of the form Eq. 3.13 is tricky, more so when the Cs condition is enforced and when the above computation needs to be efficient. Increasing Cs increases the basis by lowering the minimum SA = SB in the N2 -spin enlarges so the basis is enumerated as consecutive blocks of decreasing SA. The number of states within each block is nSA 2SA+1 2 , where nSA is the number of unique N 2 -spin SA sectors (with 2SA+1 z-projections) obtained by mixing lower level sub- spaces. To speed up the matrix element computation, the N2 basis is also sorted in blocks of decreasing SA=B. Each of those blocks is further divided into sectors of 2SA+ 1 containing consecutive z-projection of each unique SA sector. Details of 47 0.0 0.5 1.0 CS 10-5 10-4 10-3 10-2 10-1 100 dE /E dC 0 0.5 1 CS 0.6 0.7 0.8 0.9 1 O ve rla p N=16 N=32 N=64 a) b) Figure 3.2: (a) Overlap jhGS;CsjGS;Cs + 2N ij2 between ground-state wave functions corresponding to adjacent values of Cs for N = 16;32;64. Note that for N = 64, the overlap is already 95% for Cs = 14 ; (b) Frac- tional change 1EGS dEGS dCs in ground-state energy corresponding to different Cs cutoffs. The lines are linear fits to the data. the implementation is provided in Appendix B. The overall matrix can be constructed column-by-column and term-by-term and is extremely sparse. The practical implication is that the grand problem is decomposed into many small independent computations which can be massively parallelized (See Appendix B). Table 3.1 compares the Cs = 14 basis size to those of other bases commonly used by non-iterative methods. At Cs = 14 , the N = 32 systemwas solved by a single CPUwithin seconds using a generic Lanczos routine. To convince the critics, the single-processor computation was modified to cal- culate an explicit wavefunction for N = 64 to break the N=40 limit which required computer clusters. The numerical vector was about 8GB in size, and the matrix was too large for the physical memory of a single workstation. The remedy was the use of one processor core in fetching data from fast hard drive storages while another core performs floating-point operations. Starting from theCS = 18 solution, the GS vector was obtained after 15 matrix-vector products, each taking less than 30 minutes without any optimizations. 48 0.00 0.25 0.50 0.75 1.00 CS -0.75 -0.70 -0.65 -0.60 -0.55 -0.50 E G S/N N=16 N=32 N=64 0 0.5 1 CS 0.12 0.16 0.2 0.24 < S r S r + (L /2, L/ 2)> 0 0.5 1 CS 0.48 0.52 0.56 m Figure 3.3: GS energy computed at specific Cs values (symbols) and the ex- trapolation to Cs = 1 (solid lines) for N = 16;32;64. The horizontal dashed lines are the lowest values from Ref. Inset: expectation val- ues of the non-commuting operators p hm2i and hSr Sr+( L2 ; L2 )i. 3.4 Results We have performed computations for the full basis for N = 16 as a benchmark for larger N values. For N = 32 (64), we computed up to Cs = 12 (Cs = 1 4). While the computations with Cs < 1 are variational, the goal here is to demonstrate that Cs 14 is a good rule-of-thumb value to capture well the AFM ground-state. The stability of this formulation is demonstrated in Fig. 3.2(a) which shows the overlap between the GS wave functions for consecutive values ofCs. The deviation decreases exponentially with increasing Cs, and for Cs = 14 the overlap is 95%. We conclude that the GS vector is already pointed in the correct direction even for small Cs, and subsequent increments of Cs merely result in minor improvements. This is reinforced by the convergence of the GS energy demonstrated in Fig. 3.2(b), which supports the scaling law dEGSdCs e aCs . The error decays exponentially, with a rate that increases with N. The N = 32;64 lines cross at Cs 14 where the frac- tional error is 10%. This exponential efficiency is amplified by the fact that a linear decrease in Cs leads to a combinatorial decrease of the basis size because of the numerous low-spin combinations removed from the basis. These GS energies are not as accurate as for established iterative methods; 49 however, a single-pass computation at Cs = 1=4 already has less than 2% error. Using a linear fit of dydx e x, an estimate of EGS at Cs = 1 is obtained by simple integration (see Fig. 3.3). From the Cs = 1 estimates and the N 3 2 scaling, we extrapolate limN!¥ EGSN = 0:6671, within 0:5% of the best published value of -0.6692, even though it is achieved at a significantly reduced computational cost. We have thus far demonstrated the exponential convergence of the GS eigen- value and eigenvector with increasingCs, which validates our proposed variational basis. While a systematic approximation of the wavefunction can provides esti- mates of all quantities – more accurately so for those commuting with the Hamil- tonian – derivable from wavefunctions, our approach is not directly competitive with established methods in evaluating particular correlation functions of interest. For example, Monte-Carlo sampling is inherently one of the best, if not the best, way of evaluating numerical integrals. The undoped system is predominantly char- acterized by correlation operators which do not commute with the Hamiltonian. The insets of Fig. 3.3 show expectation values of two such operators (the staggered magnetization of Eq. 3.12 and the spin-spin correlation at maximum distance) eval- uated from our eigenvector for a given value of Cs. Although these operators are non-commuting, their matrix representation is diagonal in our basis and therefore trivial to compute. The convergence of these quantities is essentially linear for Cs < 14 . For greater values ofCs, d dCs hSr Sr+( L2 ; L2 )i is suppressed faster than e Cs so the post-threshold convergence is even better than for the GS energy. On the other hand, because m is exploited to discard low-spin parts of the Hilbert space, its approximated value is always higher than the true value; convergence is only asymptotical for Cs > 14 . Thus, the computational discount for the GS wavefunc- tion at a small Cs is achieved at the cost of lowered accuracy for one of many definitions of staggered magnetization. As discussed in Sec. 3.1, the polaronic de- scription of doped cuprate systems is difficult due the non-bosonic nature of the undoped background. The ultimate goal of the approach is to model the wave- function of the AFM background with and without doped hole. The application to lightly doped AFM system will be demonstrated in Ch. 4 and 5. 50 3.5 Conclusions In summary, one dominant characteristic of the AFM background is that, within each bipartite sublattice, N8 spins add to a total spin of zero while the other 3 8N add to the maximum value of 316N. The octa-partition provides a natural, efficient means of exploiting this knowledge in modeling the background. In particular, Cs was introduced to systematically parameterize the resulting Hilbert space, and the N8 -spin idea was used as an apriori argument in postulating that Cs 14 is the minimal value that captures the GS accurately. The resulting basis is combinatorially smaller than in any other schemes, apriorily known, and computationally very efficient. Numerical results confirmed the postulation. The stability of the formulation was demonstrated in Fig. 3.2(a) and the sources of error were identified. Using a single thread on a commodity computer, explicit calculation of the AFM wavefunction is feasible for systems with up to N = 64 spins, significantly exceeding the current full ED record of N = 40. This breakthrough was achieved without transforming the problem away from the real-space basis of the square lattice with periodic boundary conditions, in which other interactions are defined naturally. Neither translational nor point- group symmetries are exploited. Thus, it can be used as an excellent starting point for studying doped models where these symmetries are broken, or as an ef- ficient kernel for iterative methods, such as renormalization and quantum cluster theory. Of course, this insight might also be adapted in devising more effi- cient variational Monte-Carlo formulations. From a holistic point of view, the N8 -spin observation contributes to the de- scription of the AFM background in analogy to the b j0i = 0 statement taken for granted by a phonon problem. Excited states analogous to b+j0i are indeed im- portant when carriers are added to the system. Because the total spin of the AFM plus carriers is conserved while AFM order is expected to persist for low doping at low temperature, arguments such as Eq. 3.12 can still shed light on the AFM phase space spanned by the solution. This will be explored in Ch. 4 and 5. 51 Chapter 4 Spin Polaron on the Copper-Oxygen Plane 4.1 Introduction It is widely believed that a proper description of holes in a spin- 12 2D antiferro- magnet (AFM) with full quantum fluctuations could provide the answers to many unsolved puzzles in cuprates. Consideration of more exotic scenarios [44, 45, 48, 126] are indeed exciting developments; however, a detailed modeling of the hole+AFM is a crucial first step to understanding the significance of such addi- tions. This problem is difficult because of the complicated nature of the 2D AFM background, whose quantum fluctuations in the presence of extra holes were never fully captured for a large CuO2 lattice. As discussed in Ch. 1, the simplest model required to describe a CuO2 plane is the three-band p d model whose unit cell consists of one dx2 y2 and two 2ps orbitals[14, 15]. For copper sites located at l, some integral multiples of the lattice parameters ax = (a;0) and ay = (0;a), the oxygen sites are located at l + ex=y, where ex=y = 12ax=y. Since there are more electrons than vacancies in the problem, the model is specified with fermion operators p†l+e;s and d † l;s creating holes of spin s on a vacuum with one hole per unit cell. For the orbital setting illustrated in Fig.4.1a, the model reads 52 H3B = Tpd +Tpp+Dpdånl+e;s +Uppånl+e;"nl+e;#+Uddånl;"nl;#; (4.1) with Tpd = tpdå( p†l+e;s + p†l e;s )dl;s +h:c: Tpp = tppåsd p†l+e+d ;s pl+e;s t 0ppå p†l e;s + p†l+3e;s pl+e;s nl;s = d † l;sdl;s nl+e;s = p † l+e;s pl+e;s : Tpd is the direct oxygen-copper hopping. Tpp contains the direct tpp oxygen-oxygen nearest-neighbor hopping as well as the t 0pp oxygen-oxygen next-nearest-neighbor hopping mediated by the copper 4s orbitals. The sign of the matrix elements are determined by the phases of the overlapping orbital wavefunctions in Fig.4.1a. In particular, for upper-right/lower-left (upper-left/lower-right) O-O hopping, sd = 1 (sd = 1). Dpd is the charge transfer energy, which is the energy cost of the d9 ! d10L¯ Ligand hole excitation. In reality, this energy decreases with the distance between the d10 site and the L ¯ Ligand hole, and a nearest neighbor Hubbard interaction should be added to the three-band model to capture this deviation. However, the following procedure is sensitive only to the local charge transfer energy, and such inclusion would merely lead to minor variations of effective parameters with no significant effects on the results, at least in the single hole scenario. Udd=pp are repulsive Hubbard interactions. Although one issue in cuprate physics is the lack of significantly small parameters, the magnitude of parameters have been established to be tpp < tpd < Dpd <Upp <Udd. In practice, neighbor Coulomb interactions such asUpdd†dp†p should be considered; however, such considerations do not add any new matrix elements into the effective Hamiltonian of interest (Sec. 4.2.1). This three-band model has not been solved exactly. Unbiased numerical solu- tion has also been limited due to the combinatorially large Hilbert space and the sign problem in quantum Monte-Carlo integration. To make progress, this chapter derives an effective spin polaron model which 53 δ ε Cu O O (a) swt Cu O tpd tpd O (b) Cu O Cu Figure 4.1: (a) Two adjacent unit cells of the CuO2 plane. The orbitals kept in the 3-band model of Eq. (4.1) are shown, with white/shaded for pos- itive/negative signs. The two e vectors (solid arrow) and the four d vectors (dashed arrow) are also shown. (b) Sketch of a virtual process of Tswap. captures the anion-cation nature of the CuO2 plane in the context of superexchange. The model is then solved exactly in a 32-unit-cell CuO2 cluster with the aid of the optimization in Ch. 3. 4.2 The Microscopic Model 4.2.1 Derivation of the Spin Polaron Model Noting that all Tpd virtual processes in Eq. 4.1 increase energy by eitherU and/or D, and that high-order terms in tpp always come with an even power of tpd , we derive an effective Hamiltonian for the states p†l+e;sÕ jsli using degenerate Rayleigh- Schrodinger perturbation theory Heff = Pgs[H+T 1 Pgs E0 H0T + ]Pgs Pgs = Õ l å s d†l;sdl;s using Pgs as the projector into this subspace of degenerate states with all N Cu spins intact. The resulting Hamiltonian has the form Heff = Tpp+Tswap+HJpd +HJdd (4.2) 54 where Tpp is the oxygen-oxygen hopping in Eq. 4.3. The most straight-forward derivation can be given by considering the Udd ! ¥ limit. The second-order tpd virtual processes shown in Fig 4.1b lead to an additional spin swap O-O hopping term: Tswap = tswåsh p†l+e+h ;s pl+e;s 0 js 0le;h ihsle;h j where tsw = t2pd=Dpd and le;h = l+ e +(h e=jh ej)e . This comes about when a hole from a neighboring Cu site located at le ;h first hops to one of its other 3 hole- free neighbor O sites, followed by the original O hole falling into the now-empty Cu site. The end result is either NN O-O hopping if h = d (Fig 4.1b), or NNN O-O hopping if h = 2e . In the latter case sh = 1. For h = d , the sign sh is opposite to that of the direct NN O-O hopping, due to the phases of the dx2 y2 orbital. As shown below, such processes do not simply renormalize tpp or t 0pp. If the original Cu and O holes have antiparallel spins, it is also possible for the Cu hole to hop onto the occupied O orbital, followed by one of the O holes falling back on the empty Cu site. This leads to Heisenberg exchange: HJpd = JpdåSl Sle where Jpd = 2t2pd=[Upp+Dpd ] and Sle ;Sl are respectively the spin operators at O and Cu sites. Finally, there is superexchange between NN Cu holes, except for the bonds “blocked” by p-holes: HJdd = JddåSl2e SlPs (1 nle;s ) where Jdd = 8t4pd=[D 2 pd(Upp+2Dpd)]. There are also processes involving a Cu hole hopping to an adjacent empty O and then back. These are analogous to lollipop diagrams and renormalize all the above terms by roughly equal amounts when all orders are included, hence we ignore them. We also discard a third-order tppt2pd “Kondo-hopping” term for reasons discussed below. While we find that the three-spin polaron (3SP) plays an important role, our approach is different from previous work [20–24] by recognizing (i) Tpp’s signif- 55 icant contribution to polaron coherence (as discussed in the next session), (ii) its complementing process Tswap (Fig. 4.1b), and (iii) suppression of superexchange along the bond inhabited by the hole. 4.2.2 Polaron Coherence and Other Considerations To illustrate the importance of the terms we kept, consider the GS of HJpd , which has energy Jpd . It consists of two out-of-phase terms, each with the O hole plus a triplet of two d9 Cu spin neighbors. This is the 3SP j3SPi= r 1 6 (j "#i+ j #"i) jsih r 2 3 jssij sih; (4.3) first studied by Emery and Reiter . The other eigenvalues are 0 (Cu singlet plus O hole) and Jpd=2 (in-phase terms of Cu triplets plus O hole). The Cu-Cu exchange in HJdd is different from that of other models as it ac- counts for the lack of exchange between Cu spins bridged by an O hole. This lowers the 3SP by a further Jdd . Basically, the O hole activates an effective Cu-Cu FM coupling to frustrate the bi-partite description (which, however, is not a requirement for long-range AFM order [115, 128]). This is qualitatively different from the t J model, where a ZRS breaks 4 AFM bonds without frustrating the bi-partite lattice. Due to the sign difference between Tpp and Tswap, the effective hopping of the O hole is reduced to tpp tsw for a Cu-O triplet (with the “central” Cu at le;h ), but is enhanced to tpp+ tsw for a singlet. The interference is completely kinetic in nature, and thus conceptually different from the ZRS mechanism which is due to intra- plaquette phase coherence . The 3SP can also be written as a superposition of two Cu-O singlets: j3SPi (j ";sij #ih j #;sij "ih)+(js ;"ij #ih js ;#ij "ih). Thus, the 3SP already stabilized by Jpd and Jdd is also expected to have an increased bandwidth because of Tpp+Tswap. We have considered corrections to Heff numerically and have found that they do not affect our results. The reasons for this can be understood as follows. A third-order term in the previous section and finite Udd yield an additional 56 ”Kondo-hopping” term TKondo = tkå( 1)nh p†l+e+h ;a pl+e;a(Sl+e Sle;h 1 4 ) tk = tppt2pd D(D+Upp) + t2pd (Udd Dpd) ! considered previously [19–24]. TKondo is non-zero only for a Cu-O singlet so it is merely a O(t2pd=Udd) reinforcement to the above terms. A finite Udd also modifies matrix elements by some t2pd=(Udd Dpd) correc- tions. The effects are marginal because we scale the model in units of Jdd which itself contains such a correction. There are also fourth-order corrections to Heff. These are either reinforcement to HJpd or S +S terms. These terms are small, es- pecially against the coherence discussed above, due to the lack of the factor of 4 as in superexchange, and because their denominator involvesUdd and D+Upp. 4.3 Results Using tpd = 1:3eV , tpp = 0:65eV , t 0pp = 0:58tpp, Dpd = 3:6eV , andUpp = 4eV , we scale the parameters in units of Jdd to find their dimensionless values to be tpp = 4:13, tsw = 2:98, and Jpd = 2:83. We push the computational limit to perform total-spin-resolved exact diagonal- ization (ED) of a topologically superior  cluster of N = 32 CuO2 unit cells, treating the AFM background exactly. ED provides the transparency, flexibility, and neutrality to support new results. The price for a systematic mapping of the excited states is the limited cluster size, denying us access to higher k-space reso- lution. N = 16 results are provided to check for size dependence. All low-energy eigenstates have a total spin of either ST = 12 or 3 2 . The z pro- jections for each ST are degenerate. The ST = 12 subspace is due to the s= 1 2 hole mixing with various S = 0 background states, including the AFM GS, or mixing with the S= 1 background states, including the “single-magnon” states. The s= 12 carrier can also mix with S = 1 or 2 background states to yield the ST = 32 sub- space. The partition of the SzT subspace into separate ST sectors was managed by the optimizations from Ch.3. Unlike there, no basis truncation was employed here 57 -44 -43 -42 E(k) a) ST=1/2 32A ST=3/2 32A ST=1/2 16B ST=3/2 16B t-t’-t’’-J 32A (0,1) (0,0) (1/2,1/2) (1,1) (0,1) (1,0)(1/2,1/2) (k x ,ky)/pi 0 0.1 0.2 0.3 Z(k) b) Figure 4.2: a) Energy and b) quasiparticle weight (bottom) for the lowest eigenstates with ST = 12 and 3 2 vs. momentum. Different sets are shifted so as to have the same GS energy. for rigorous results. The (k;ST = 32 ;S z T = 1 2 ) sector contains 0.44109 states. 4.3.1 Energetics Fig. 4.2(a) shows the lowest eigenenergies. The GS has k = (p2 ; p 2 ) and ST = 1 2 , and is consistent with the 3-spin polaron (Eq. 4.3). Remarkably, we find similar dispersion along (0,0)!(p ,p) and (0,p)!(p ,0) without having to add longer range hopping or fine-tune parameters as is needed in one band models. The biggest surprise, though, are the low-lying ST = 32 states which go below the 1 2 states near (0,0) and (p ,p). The solution’s robustness is supported by the difference E 3 2 -E 1 2 which decreases with increasing N at each k. We note that the magnon spectrum is gapped as 1=N for finite lattices . Going from N = ¥! 32, excitation energies increase at all q, e.g. from 0!0.28J at q = 0 and 2.4J!2.6J at the BZ 58 edge . Given this finite-size scaling of the “free magnon” energy, we find that ST = 32 states at (0,0), (p4 ; p 4 ), ( 3p 4 ; 3p 4 ) and (p;p) have much lower energy than minq[E 12 (k q)+Wq] for a ST = 12 state plus a free magnon. At other k-points, the ST = 3 2 states are within 0:2J of this value, so we cannot draw definite conclusions without ac- cess to larger N data. It follows that the ST = 32 states are stable polarons at least in the regions marked by thick solid lines in Fig. 4.2(a). Thus, a ST = 12 quasiparticle cannot describe the low-energy states throughout the BZ. 4.3.2 Wavefunction Analysis To compare the lowest energy states on both sides of the crossing, we note that the lowest kx = ky eigenstates have odd parity upon a P̂x$y reflection (Fig.4.1a) so they can be expressed as 2 1=2(1+ P̂x$y)åeikl p†l+ex;s js ; lix. The band-crossing results in noticeable change in the expectation values of the correlation function: Ĉx(d ;a) = 2å l;s Sl+d Sl+d+anl+ex;s (4.4) which measures the correlation between two Cu sites separated by a lattice constant at a certain distance away from the hole. hĈyi is a reflection with P̂x$y for kx = ky. hCi ranges from 34 for singlet, to 0:33 for 2D AFM GS, to 14 for triplet. Fig. 4.3a shows the hĈxi correlations between neighbor Cu spins, when the hole is located at the darkly shaded bullet, in the GS: k= (p2 ; p 2 ), ST = 1 2 . The hole affects the AFM order in its vicinity. Because of the hole-spin exchange HJpd and the blocked superexchange between the two Cu spins neighboring the hole, these “central” spins have triplet correlations, of0.13. Also, hHJpd i -0.9Jpd , showing that locally this is consistent with the 3SP solution (Eq. 4.3). More interesting are the correlations with the other 3 neighbors of each of these central Cu spins: with two of them, there are robust AFM correlations of 0:22, while with the third the correlation nearly vanishes (lightly shaded bullet). This is counterintuitive if one views the system as a fluctuating Neél background, where a spin-flip would change the spin-spin correlation to all four neighbors. Although the two central Cu spins have 23 weight in triplet configuration which is hardly bi-partite, long-range 59 (b)(a) −0.28 −0.28 −0.26 −0.29 −0.30 − 0. 26 − 0. 05 − 0. 22 − 0. 05 − 0. 31 −0.33 −0.26 −0.30 − 0. 26 +0.13 −0.33 −0.29 − 0. 31 − 0. 22 +0.17 − 0. 07 −0.21 −0.14 −0.21 −0.21 −0.14 −0.21 − 0. 23 − 0. 23 − 0. 23 − 0. 23 − 0. 07 − 0. 07 − 0. 07 Figure 4.3: hCx(d ;D)i for the lowest energy state at (a) (p2 ; p2 ) with ST = 12 , and (b) at (p ,p) with ST = 32 . The darkly-shaded bullet denotes the oxy- gen position at l+ ex. Each bullet shows the correlation value between the two sandwiching Cu sites. The central 12 Cu sites are shown; the correlations between the other 20 Cu spins converge fast towards the AFM value of -0.33. hCy(d ;a)i is the P̂x$y reflection. AFM order cannot be automatically discounted [115, 128]. Indeed, the correlations we find are consistent with such order, except for the zigzag of 3 bonds shown by shaded bullets. This strange shape is dictated by the hopping mechanism. For a Bloch wave, oxygen-oxygen hole hopping in the upper-left/lower-right direction yields a phase shift of ei0=ei(kx ky) and hence constructive interference if kx = ky. In contrast, hopping in the upper-right/lower-left direction yields a phase shift of eikx=e iky , so the interference is scaled down by cos(kx=y). For kx = ky = p2 , having a mixture of singlets and triplets upper-left/lower-right to the O hole lowers energy with the least disturbance to AFM order. Thus, the two outside zigzag bonds are triplet “disturbance tails” pointing orthogonal to the momentum direction. This is different from the ZRS, which freezes a copper spin by the intraplaquette coherence of four oxygen sites. Fig. 4.3b shows the correlation values when the ST = 32 polaron becomes the lowest energy state at (p ,p). The results look similar at (0,0). hHJpd i remains 0:9Jpd , but there are now four more heavily disturbed bonds. This further supports this being a stable polaron with an extra magnon locally bound close to the O hole. We stress here that this 32 polaron is formed by a spin disturbance around the 3SP. This is very different from the S = 32 excitation local to HJpd with energy + Jpd2 . 60 Before we proceed further, let us note that the T=0 single-hole Green’s function hAFMjpeiHt p†jAFMi contains information only of the states in the Krylov sub- space which is the subspace spanned by the class of wavefunction Hnp†jAFMi, for all positive integers n. It follows that the spectral weight of first electron re- moval Z(k), being the imaginary part of G(k;w), is exactly zero if the lowest hole state is not in the Krylov space of the electron removal state. Fig. 4.2b shows the quasiparticle weight Z(k) for the first electron removal state. The major difference from other models is that Z(k) = 0 in three regions: a) Z(0;0) = Z(p ,p)=0 because here the lowest eigenstate has ST = 32 which due to spin-conservation is not in the Krylov space of any ST = 12 state, and b) Z(0;p)=0 because the lowest eigenstate is not in the Krylov space Hnp†jAFMi. The t-t’- t”-J model restricts a spin degree of freedom, leading to Z(0;0) 0:1 and a fi- nite Z(p ,p)  instead of zero as dictated by a 32 polaron. Unlike the t-t’-t”-J model , Z(k) is concave down along (p/2,p/2)!(0,p) in closer agreement with ARPES which measured a maximum between (p2 ; p 2 ) and ( p 4 ; 3p 4 ) [52, 53]; however, this observation might be sensitive to finite-size effects. Our Z(k) is smaller than that of the t-J model, suggesting less “free particle” nature of the polaron. The lowest energy state at k = (0;p) has ST = 12 , but the spectral weight of first electron removal state is exactly zero. This is confirmed by solving for the lowest energy state in the full basis, the Krylov subspace of åeikl p†l+ex jAFMi, and Krylov subspace of åeikl p†l+ey jAFMi. The k = (0;p) lowest state of the two latter cases are at least 0.006J higher than the first. The lowest energy state is not in the Krylov space of electron removal states so Z is exactly zero. This seems to be due to symmetry, although we do not yet fully understand this. Their close existence above the Z=0 lowest state may be related to pseudogap phenomena in this region; however, this needs to be investigated in more detail. Fig. 4.4 shows the correlation for the lowest eigenstate. Compared to the GS (see Fig. 4.3a), there are two more disturbed bonds as required by the reflection parity about the momentum direction. This larger disturbance range is accompanied by more negative (AFM) correlation values. 61 −0.29 −0.24 −0.19 −0.24 −0.24 − 0 .2 3 − 0 .2 3 −0.29 − 0 .1 0 − 0 .1 0 − 0 .1 0 −0.24 − 0 .2 3 −0.19 − 0 .2 3 +0.16 − 0 .1 0 −0.09 − 0. 28 −0.24 − 0. 28 − 0. 19 − 0. 28 −0.24−0.24 − 0. 26 − 0. 28 − 0. 19 −0.09 −0.09 −0.09 − 0. 26 −0.24 + 0. 15 Σ(−1)ly [ p † l+x,σ |σ, l〉x + p † l+y,σ |σ, l〉y ] Figure 4.4: hCx=y(d ;a)i for lowest state at k=(0,p) with ST=1/2. 4.3.3 Further Experimental Implications Although we are restricted to rather low momentum resolution, more can be said about the E 3 2 E 1 2 = 0 band crossings in the nodal direction by looking at the k points between which the difference switch signs. The observation in Fig. 4.2a is that, going away from the (p2 ; p 2 ) GS, E 32 E 12 is larger towards (0;0) than towards (p;p). The 12 ! 32 band crossing would induce an abrupt change in Z(k) from non- zero to exactly zero, irrespective of the Z(k) value on the finite side. The larger E 3 2 E 1 2 towards k = (0;0) would suggest the non-zero region is larger towards k= (0;0) than towards k= (p;p). Fig. 4.2a also shows that the 32 states get pushed further down as N ! ¥ so the crossing is expected to be closer to the (p2 ; p2 ) GS. This is consistent with ARPES which indeed observed an abrupt peak suppression in the nodal direction as well as the peaks surviving longer towards k= (0;0) [52, 53]. Even when the ST = 32 states are not lowest in energy, they hug the ST = 1 2 band. This provides a . Jdd /2 energy scale for spin excitations. At finite T , as magnons become thermally activated, these 32 states become “visible” to ARPES. This sug- gests a T -dependent broadening mechanism of . Jdd /2 scale. Coincidentally, this is the same as the energy scale linked to phonons [43, 45, 54]. Recent neutron experiments on samples at higher doping reveal 50 meV magnetic response centered at q = 0, away the AFM resonance momentum [55, 56]. The bottom of the single-particle band structure in Fig 2(a) indeed has a 62 q = 0 12 -to- 3 2 excitation of this energy scale (the computed energy gap is larger than the 1N S=1 gap for N=32 lattice). It has also been pointed out that the q = 0 magnetic excitation can be explained by ordering multiple spins on oxygen sites in the higher-doping scenario . In addition to the 32 band, there are internal excitations of the local 3SP since HJPD also has a S = 1 2 doublet and a S = 3 2 quartet separated in energy by Jpd and 3Jpd /2 from the 3SP. Magnetic excitations at these energy scales have been observed via inelastic resonant x-ray scattering for highly doped samples . 4.4 Conclusions We derived and solved a model which includes the O sites and takes full account of the AFM quantum fluctuations, for large N = 32 clusters. The phases of the p and d orbitals lead to phase coherence via Tpp+Tswap. This is reenforced by HJpd and the blocking of the AFM superexchange, making corrections such as TKondo negligible. The dispersion is similar to that measured by ARPES; however, the lifting of the Cu-O singlet restriction present in ZRS-based models leads to wave functions of a different nature, namely the 3SP where the O hole correlates with both its neighbor Cu sites. This model also provides low-energy channels for S= 1 excitations. Z(k) was found to be identically zero in certain regions of the BZ for two reasons: (1) the spin- 32 of the lowest states close to (0;0) and (p;p); and (2) around the antinodal region because of the lowest state there being exactly orthogonal to the single electron removal state. 63 Chapter 5 Two Holes in the Spin Polaron Model 5.1 An Intuition about the Doped Antiferromagnet The work presented in Ch. 3 and Ch. 4 was actually motivated by the drive to study multiple holes injected into the CuO2 layer. The Hilbert space of the doped cuprate is combinatorially large. Brute-force Monte-Carlo simulations are limited by the Fermion sign problem, and zero-temperature computations have been carried out with variational Monte-Carlo and NRG, which employ various approximation or special boundary conditions. The largest two-hole system studied by an unbiased method is the N=32 t-t’-J model. Because each additional hole removes a spin instead of adding one in this one-band model, Leung et. at. managed to use exact diagonalization to show that a realistic value of t’ destroys short-range hole-hole attraction. The result contradicts other studies using variational or mean- field approaches which predict binding. The discrepancy between this exact N=32 solution and approximated solutions on larger lattices has been attributed to a finite size issue. As discussed in Ch. 4, the solution of the spin polaron model (Eq. 4.2) solved in a cluster with 32 Cu and 64 O orbitals reveals properties completely missed by previous studies which employed various approximation or small clusters. In fact, 12 unit cells were required to illustrate the difference between the ST = 12 and 64 ST = 32 polarons in Fig. 4.3. The two-hole solution for a larger cluster could be different from previous studies limited to N=16. I realized early that, even after exploiting translational and spin-projection sym- metries, the two holes injected into the N=32 spin polaron model would have a Hilbert space with 0:1541012 states, and larger N would involve literally an as- tronomical number of states. A system of this size challenges the capability of unbiased methods, and such an endeavor is extremely high-risk, especially consid- ering the many unknowns that come with the novelty of the spin polaron model. With the drive to solve for the zero-temperature ground state with 2-dimensional periodic boundary condition, I came up with a rather radical plan. In the N ! ¥ limit, two additional holes should not be able to drastically change the undoped AFM groundstate; therefore, the compact yet systematic way of modeling AFM in Ch. 3 should be a stable, systematic way of modeling the two-hole scenario. The octapartite scheme with CS truncation discussed in Ch. 3 was formulated under this belief. Its first application has not been publicized, but the CS = 14 truncation was crucial in the development of the work in Ch. 4. The formulation enabled an approximate solution to Eq. 4.2 for N=32 within seconds, allowing me to scan the parameter space extensively and note various properties missed by previous studies. Nevertheless, full exact diagonalization was performed to prove the new findings are not mere numerical artifacts. The convergence of the one-hole ground state is illustrated in Fig. 5.1. The energy computed for CS = 14 and 1 2 is within 3.6% and 0.5%, respectively, of the exact value. Table 5.1 shows the exact wavefunction’s probability in specific sub- spaces SA N SB of sub-lattice total-spin SA=B. While theCS = 14 truncation captures 95% of the undoped wavefunction as discussed in Ch. 3, Table 5.1 shows that CS = 14 captures 73% of the one-hole groundstate. The next two increments of CS contains 17% and 6%, respectively, of remaining weight. The addition of a hole couples the spin background to adjacent values of SA N SB, and states added by increasing CS have decreasing importance in the wavefunction; therefore, it is reasonable to expect an increment of dCS 1N from CS = 14 to suffice. At least in the very dilute limit, my postulate holds, and the formulation in Ch. 3 does provide a good starting point to systematically capture the low energy state. These observations provide reasonable merits for the application of the scheme to 65 0 0.5 1 CS -50 -45 -40 -35 -30 -25 E G S/J dd 0 0.5 1 CS 10-4 10-3 10-2 10-1 100 ∆E /E Figure 5.1: Convergence of the one-hole ground state. (left) GS energy cal- culated for increasing CS. The value approaches rapidly to the exact value marked by the dotted line. (right) fractional change in GS energy for the next increment ofCs. Solid lines are extrapolation and linear fits. SAnSB 0 1 2 3 4 5 6 7 8 0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1 0.00 0.02 0.03 0.00 0.00 0.00 0.00 0.00 0.00 2 0.00 0.03 0.14 0.16 0.00 0.00 0.00 0.00 0.00 3 0.00 0.00 0.16 0.69 0.68 0.00 0.00 0.00 0.00 4 0.00 0.00 0.00 0.68 2.52 2.21 0.00 0.00 0.00 5 0.00 0.00 0.00 0.00 2.21 7.06 5.28 0.00 0.00 6 0.00 0.00 0.00 0.00 0.00 5.28 14.24 8.44 0.00 7 0.00 0.00 0.00 0.00 0.00 0.00 8.44 18.06 6.85 8 0.00 0.00 0.00 0.00 0.00 0.00 0.00 6.85 9.97 Table 5.1: (N=32) The S= 12 k = ( p 2 ; p 2 ) single-hole ground state’s weight in subspaces of particular SA N SB. Numbers are percentage adding up to 100%. Two lines are drawn to highlight the CS = 14 truncation which discards states with SA;SB < 6 and yield an energy within 3.6% of the exact value. 66 the two-hole scenario, whose groundstate should be captured in the subspace of CS 14 +O(dCS) for N ! ¥. To provide a conservative error analysis though, I would aim for a capability of up to CS 12 , which in turn limits the cluster size to N=32. This chapter reports the application of the octapartite scheme discussed in Ch. 3 in solving the two-hole scenario in the spin polaron framework discussed in Ch. 4. In particular, two-body interactions are derived for the two-hole scenario. The Hilbert space formulation and truncation is then documented in detail. Numerical results are presented. Due to the novelty of the model and the numerical scheme as well as the lack of any exact solution, many details remain open questions and this effort is on-going at the time of writing. Nevertheless, the numerical approach proves to be successful in capturing doped antiferromagnet and the solution shows enhanced singlet correlations. 5.2 The Spin Polaron Model for Two Holes An effective model for two holes interacting with an antiferromagnetic background can be derived using the same procedure discussed in Sec. 4.2.1. When hopping vanishes in the three-band model (Eq. 4.1), the two-hole GS is p†l+e;s p † l0+e 0;s 0Õ jsl00i, with a degeneracy of 2N 2 2N+2. The Rayleigh-Schrodinger expansion in Sec. 4.2.1 operates in the N+2 hole sector of the Fock space, and the outcome is slightly dif- ferent from the single-hole scenario. In particular, the projector used in the expan- sion becomes P2h =Õ(1 nl+e;"nl+e;#)Õ(nl;"+nl;# 2nl;"nl;#); (5.1) which allows only states with a full lattice of copper spins and no double-occupancies (due toUpp=dd). The Hamiltonian can be written as two parts He f f = P2hH1P2h+H2; (5.2) where H1 is Eq. 4.2 derived for the single-hole scenario and H2 is the correction when the two holes are close to one another. The idea here is to start with the established single-hole scenario and correct 67 P2hH1P2h from second- to forth-order with H2 =H (2) 2 +H (3) 2 +H (4) 2 in the presence of a second hole. In the Rayleigh-Schrodinger expansion with projector P2h , the appearance of (1 P2h) dictates that all terms in H2 must have an even power of tpd because the final states must have one copper spin per unit cell. Because H2 accounts for only 2-hole corrections, one can deduce that all terms in the second- order H(2)2 share a factor of t 2 pp. All terms in the third-order H (3) 2 share a t 2 pdtpp factor because the second step of any three-step t3pp process would yield states projected out by (1 P2h) in the perturbation. The forth-order H(4)2 can have terms with a prefactor of t4pp, t 2 ppt 2 pd and t 4 pd . The second order corrections are the bare hole-hole correlations. H(2)2 = 2t2pp Upp å( 1)d (d1d2) S S 1 4 a 0b 0 ab p†l+e+d1;a 0 pl+e+d1;a p † l+e+d1+d2;b 0 pl+e ;b (5.3) Here, d1;2 sum over (ex;ey). The matrix element is non-zero only if the two holes are jd1j = ap2 apart. The sign of the matrix element is positive when d1 d2 = 0. There are 8 non-zero matrix elements in this situation. Two of these correspond to d1+ d2 = 0 so the static AFM exchange is 2 2t 2 pp Upp . There are also 6 other ways for one hole to “skip” over the other with such a Heisenberg factor. The factor of 14 plays an important role here in setting the matrix element to zero whenever the two oxygen holes are in a triplet configuration. There is an abundance of terms in third- and forth-order corrections which are a sub-set of terms studied in the literature of 2-band models[20–24]. To provide a simple physical picture, we consider terms that are greater or equal to Jdd , which is roughly 816t 4 pp 3U3pp 1681 tpp for tpd 2tpp and Dpd Upp tpp 6 . The basis of this approximation is the observation that the dominant short-range physics should be already captured by the numerous low-order terms with an energy scale of tpp and Tswap;Jpd ;Jpp 0:66tpp, which are already greater than the long-range physics at the scale of Jdd 0:2tpp. Adding the relevant short-range corrections with mag- nitude greater than 0:2tpp should be more than adequate to capture the physics. We also discards terms that can be factored into tpp;Tswap;Jpd ;Jpp and the identity processes, which merely contribute some combinatorial constants scaling in the 68 infinite-order perturbation. As discussed above, all third-order corrections have a prefactor of tppt2pd . The perturbation goes through two virtual states so the largest possible magnitude is tppt2pd DpdUpp tppt2pd DpdDpd 19 tpp. The splitting due to such a pair of Hermitian matrix ele- ments is 29 tpp which is the same as the Jdd splitting. Evidently, all these processes can be factored into consecutive low-order processes and are discarded. Because Upp Dpd , other processes would have denominators that are at least a factor of two greater; that is, constructive quantum interference is required for any terms to be non-negligible compared to the superexchange. Constructive interference among tppt2pd processes requires the same orbital occupancy in the initial and final state, and this can happen only if the two oxygen holes are d = (ex;ey) apart. The transition of interest is then p†lex;a p † ley;bd † l;g ! p†lex;a 0 p † ley;b 0d † l;g 0 . The ef- fective matrix elements have the form ha 0;b 0;g 0jH(3)2 ja;b ;gi, a local three-spin ring involving the copper spin sandwiched by the oxygen holes. The correction term can be expressed as a summation over each copper spin with two additional vector Dx=y summed over ex=y. Noting that all 3-step hopping would give an overall phase of 1 (Fig. 4.1) and virtual states have double-occupancy, it’s not surprising that all operators contain a shifted Heisenberg factor. 69 H(3)2 = 2t2dptpp Upp(Dpd +Upp)å S S 1 4 b 0a 0 ba p†l+Dx;b 0 pl+Dx;a p † l+Dy;g pl+Dy;b jl;a 0ihl;gj + S S 1 4 a 0b 0 ab p†l+Dx;g pl+Dx;a p † l+Dy;a 0 pl+Dy;b jl;b 0ihl;gj + S S 1 4 g 0a 0 ga p†l+Dx;g 0 pl+Dx;a p † l+Dy;a 0 pl+Dy;b jl;b ihl;gj + S S 1 4 b 0g 0 bg p†l+Dx;b 0 pl+Dx;a p † l+Dy;g 0 pl+Dy;b jl;aihl;gj + 2t2dptpp (Dpd +Upp)2å dab S S 1 4 b 0g 0 bg p†l+Dx;a pl+Dx;a p † l+Dy;b 0 pl+Dy;b jl;g 0ihl;gj + dgb S S 1 4 a 0b 0 ab p†l+Dx;g pl+Dx;a p † l+Dy;a 0 pl+Dy;b jl;b 0ihl;gj + dga S S 1 4 b 0a 0 ba p†l+Dx;b 0 pl+Dx;a p † l+Dy;g pl+Dy;b jl;a 0ihl;gj + dba S S 1 4 a 0g 0 ag p†l+Dx;a 0 pl+Dx;a p † l+Dy;b pl+Dy;b jl;g 0ihl;gj(5.4) These terms are finite if there is a 2-spin singlet formation among the 3 spins. The first four terms give an eigenvalue of +4 t2dptpp Upp(Dpd+Upp) 0:2tpp for oxygen- oxygen singlet pair and 0 for triplet pair. The last four terms give an eigenvalue of t2dptpp (Dpd+Upp)2 0:02tpp for singlet pairs and 0; 3 t 2 dptpp (Dpd+Upp)2 0:06tpp for triplets. Therefore we can approximate this term by raising the energy of oxygen-oxygen singlet pair appropriately. In the forth-order correction, all t4pd processes can be factored into two Tswap or Jpd processes. The t2pdt 2 pp processes are smaller than Jdd by a factor of 44 and no constructive interference is possible. 2t4pp processes are even smaller. Therefore, we set H(4)2 0. In summary, the 2-hole correction is in essence 70 H2 + J(2)pp å( 1)d (d1d2) S S 1 4 a 0b 0 ab p†l+e+d1;a 0 pl+e+d1;a p † l+e+d1+d2;b 0 pl+e;b J(3)pp å S S 1 4 a 0b 0 ab p†l+e+d ;a 0 pl+e+d ;a p † l+e;b 0 pl+e;b (5.5) with J(2)pp = 2t2pp Upp and J(3)pp = 4t2dptpp Upp(Dpd+Upp) for hole-hole interaction due to second and third order corrections. 5.3 Hilbert Space Formulation The Hilbert space of the two-hole spin polaron model (Eq. 5.2 with Eq. 4.2 and Eq. 5.5) contains 2N 2 2N+2 states of the form p†l+e;s p † l0+e 0;s 0Õ jsl00i, with l + e 6= l0 + e 0. To perform computations for large N, this basis must be reformu- lated to take advantage of translational symmetry, total-spin symmetry, total-spin- projection symmetry, and, most importantly, the truncation scheme discussed in Chapter 3. Point-group symmetries were not exploited to reduce the basis. The Hilbert space formulation discussed below is non-trivial; the matrix transformation re- quired to implement the basis and the truncation is a convoluted process. The code written for the project totaled 0.5MB in size, containing 65536 charac- ters. Software development of this scale is prone to minor but fatal mistakes such as a typo between “-” and “+” of adjacent keys on the keyboard. The point-group operators are used as a consistency check at various stages of the computation as discussed in Appendix A. Nevertheless, point-group symmetries are not compat- ible with all directions of total momentum K and the basis reduction is negligible to the combinatorial reduction by the octapartite scheme as shown in Table 3.1. Leaving this redundancy is a small price to pay for a consistency check. 71 5.3.1 Translational and Total-Spin Symmetries First, I’ll define singlet and triplet creation operators for the two oxygen holes. s†l;l0 = 1p 2 (p†l"p † l0# p†l#p†l0") t† 1;l;l0 = p † l#p † l0# t†0;l;l0 = 1p 2 (p†l"p † l0#+ p † l#p † l0") t†1;l;l0 = p † l"p † l0": (5.6) According to quantum mechanical angular momentum addition, N+2 spins can add up to a total spin ST by mixing a two-hole singlet to a background with total spin ST;N = ST and also a two-hole triplet mixing with a background with total spin ST;N = ST ;ST 1. Taking SzT = 0, any spin background can be specified orthonor- mally with respect to a position l. †l;l0 jail ( s†l;l0 ja;Sz = 0il å1z= 1 c(z;ST;N ;ST )t † z;l;l0 ja;Sz = zil (5.7) a denotes a particular group of 2ST;N + 1 spin configurations related by the total spin raising and lowering operators S+= T = ål S += l , summed over all Cu sites. The total spin of this a group can be ST;N = ST for two-hole singlet and ST;N = ST ;ST 1 for two-hole triplets. Due to the choice of the overall projection SzT = 0, the two-hole singlet would mix only with backgrounds with SzT;N = 0. The two- hole triplets would mix with three different projections ja;Sz = 1il , ja;Sz = 0il , and ja;Sz = 1il from the group a . The weight c(z;ST;N ;ST ) is the Clebsch-Gordon coefficients for mixing these three states with the three two-hole triplets to achieve a state with total spin ST and SzT = 0. Exploitation of the translational symmetry is performed by the use of a Fourier series of the form åeiKl†l+e;l0+e 0 jail; however, care must be taken to ensure orthonormality. Due to the commutation relation of the triplet and singlet (Eq. 5.6), the Fourier series is not straightforward 72 for e = e 0 = ex=y when both oxygen holes occupy x- or y-rung oxygen orbitals (see Fig. 4.1a). The specification of these two-hole configurations requires the two orthogonal periodic lattice vector L0 and L1 of length p N for the 2D, N-unit-cell lattice. Defining d l = (l00 l0; l01 l1); (5.8) most hole-hole configurations can be classified in the region 0 d l0 < L02 ; 0< d l1 < L1 2 0< d l0 < L0 2 ; d l1 = L1 2 0< d l0 < L0 2 ; L1 2 < d l1 0 d l0 = L0 2 ; L1 2 < d l1 < 0: (5.9) These states can be expressed using N-term Fourier series jxx=yy;d l;a;Ki= 1p Nåe iKl†l+ex=y;l+d l+ex=y jail (5.10) Because the two oxygen holes are indistinguishable fermions (Eq. 5.6), there are three remaining d l values which require special attention due to the periodic bound- ary condition. d l0 = L0 2 ; d l1 = 0 d l0 = 0 ; d l1 = L1 2 d l0 = L0 2 ; d l1 = L1 2 (5.11) The number of terms in the Fourier series depends on the spin background trans- lated by d l: Td ljail . If such a translation yields an orthogonal state, lhajTd ljail = 73 0, the Fourier series still has N terms. For the example of d l = (L02 ;0), jxx=yy;d l;a;Ki= 1p Nåe iKl†l+ex=y;l+d l+ex=y jail = 1p N L1 1 å l1=0 eiK1l1 L0 2 1 å l0=0 eiK0l0†l+ex=y;l+d l+ex=y 1+ seiK0 L0 2 TL0 2 jail (5.12) where s is the sign change due to hole-swapping in the singlet or triplet (†a;b = s†b;a). For the case of hajTd ljail = 1, the above expansion makes clear that there can be only N2 terms in the Fourier series due to the term 1+ seiK0 L0 2 TL0 2 . For this case the series has the form jxx=yy;d l;a;Ki= r 2 N L1 1 å l1=0 eiK1l1 L0 2 1 å l0=0 eiK0l0†l+ex=y;l+d l+ex=y jail (5.13) The formulation for the scenario in which one hole occupies a pl+ex and the other occupies pl+ey is straight forward. All values of d l are unique and the Fourier series has the form jxy;d l;a;Ki= 1p Nåe iKl†l+ex;l+d l+ey jail: (5.14) In summary, a full orthonormal Hilbert space can be specified by the states jxy;d l;a;Ki and jxx=yy;d l;a;Ki. Translational symmetry is specified by the quantum number K. Total spin and its projection are specified by the singlet/triplet nature of the two oxygen holes () (singlet-triplet) in conjunction to the spin back- ground a . 5.3.2 Enumeration The reason behind the formulation in the previous section is to take advantage of the combinatorial truncation of the octapartite scheme discussed in Chapter 3. For a given completeness value CS (see section 3.3), the number of different spin configurations (;a) (Eq. 5.7) is reduced to an integer nSpinRow. I would also define a number n2h for the total number of oxygen-oxygen hole configura- 74 tions. For N=32, there are N=32 distinct d l’s for †l+ex;l+d l+ey , 17 distinct d l’s for †l+ex;l+d l+ex , and 17 distinct d l’s for † l+ey;l+d l+ey . So n2h=66. The total number of states jxy;d l;a;Ki and jxx=yy;d l;a;Ki is thus n2hnSpinRow. Because the matrix elements of the Cu-Cu superexchange å<l;l0> Sl Sl0 is the same irrespective of the location of oxygen holes, the numerical basis was chosen such that an arbitrary wavefunction can be expressed in terms of a vector c of complex numbers. jYi= (n2h 1) å i=0 (nSpinRow 1) å j=0 cnSpinRowi+ jji; ji (5.15) where j specifies a particular combination of singlet/triplet and background (Eq. 5.7) and i specifies the oxygen-oxygen configuration, i.e. x-y, x-x, y-y occupancies and separation d l. Each ji; ji corresponds to a Fourier series of a definite set of (;d l;K;ST ) as discussed in the previous section. This is an orthonormal numeri- cal basis hi; jji0; j0i = di;i0d j; j0 (5.16) nSpinRown2h 1 å m=0 jcmj2 = 1 (5.17) Of course, the values of vector elements are the probability amplitudes of a quan- tum mechanical description. For example, the probability of the two oxygen holes in one of the n2h configurations is P(i) = nSpinRow 1 å j jcnSpinRowi+ jj2 (5.18) 5.4 Low-Energy Two-Hole Solutions Exact diagonalization in the subspace of CS = 0 to 12 was performed for the two- hole spin polaron model (Eq. 5.2 with Eq. 4.2 and Eq. 5.5), with parameters tpp = 4:1293, t 0pp = 0, tsw = 2:9822, Jpd = 2:8253, J (2) pp = 1:3420, and J (3) pp = 0:9182. t 0pp was set to zero to simplify the model. Another reason is that the Jpp corrections 75 (0,1) (0,0) (1/2,1/2) (1,1) (0,1) (1/2,1/2) (1,0) (K x ,Ky)/pi -67 -66 -65 -64 E( K) G S/J dd ST=0, CS=1/4 ST=1, CS=1/4 Figure 5.2: Lowest two-hole states in the CS = 14 subspace. due to t 0pp would hop a hole half-way across the N = 32 cluster, worsening finite size effects. The one-hole GS was recalculated with t 0pp = 0 for comparison. 5.4.1 Energetics Fig. 5.2 shows the dispersion of the lowest two-hole states for different total- momentum K and total-spin ST . The K = (p;p) ST = 0 groundstate is doubly degenerate. The K = (p;0) and K = (0;p) states are 0:17Jdd higher in energy. The S = 1 band crosses below the S = 0 band in the region around K = 0. These states are higher than the ground state by at least 0:3Jdd , which is roughly the smallest difference between the 12 and 3 2 single-hole states as shown in Fig. 4.3. The convergence of the lowest states at high-symmetry points is shown in Fig. 5.3. The trend of exponential convergence is similar to the undoped (Fig. 3.3) and sin- gle hole case (Fig. 5.1), suggesting that the dominant part of the Hilbert space has been captured. The binding energy of two holes is the energy difference between the two- hole energy minus twice the one-hole energy, shifted by the energy of the undoped 76 0 0.5 1 CS -70 -65 -60 -55 E/ J d d K=(pi,pi) ST=0 K=(pi,0) ST=0 K=(0,0) ST=0 ∆Eb=0 0 0.5 1 CS 10-4 10-3 10-2 10-1 100 ∆E /E Figure 5.3: Convergence of the ST = 0 lowest energy states. The two-hole energy level corresponding to DEb = 0 is indicated by the horizontal dash line. system. DEb = (E2h E0h) 2(E1h E0h) (5.19) The energy of the K = (p;p) groundstate was extrapolated to be 69:167Jdd , giving DEb = +0:18Jdd which is non-binding. Of course, this value is sensitive to the extrapolation. If the result of the smallest CS = 0 subspace is excluded from the linear fit of DEE in Fig. 5.3, the extrapolated groundstate energy lowers to 70:1228Jdd . This lowest estimate gives a binding of DEb = 0:0252Jdd . This is roughly 44K, which is three times smaller than the extrapolated pseudo gap temperature of T 145K at 2=32 = 6:25% doping. The finding here is that the groundstate of two holes does not have significant, if any, binding energy in the spin polaron model for a N=32 cluster. However, one should note that for a N=32 cluster, the undoped AFM has a 1N 0:3Jdd q=0 magnon excitation, so a rigorous conclusion cannot be drawn because jDEbj is smaller than this finite-size effect. 77 5.4.2 Wavefunction Analysis The K = (p;p) GS was found to be degenerate. It is well known that Lanczos- type diagonalization suffers greatly from degeneracy due to numerical noise. In this case, neither the generic ARPACK solver or my custom solver managed to produce a pair of orthogonal eigenvectors. As a quick fix, I modified my plain Lanczos solver to first calculate one GS vector jGS0i, then extract the GS of (I jGS0ihGS0j)H(I jGS0ihGS0j). The crucial point here is to make sure the starting vectors of the two runs are not the same. If the two starting vectors are the same, for example because the pseudo-random number generator uses the same seed, the projected Hamiltonian would not have the other GS in the Krylov space and the computation would yield the lowest excited state instead. To get an orthogonal jGS1i, I start with the same pseudo-random vector in both runs, but before starting the second Lanczos procedure, I maximize orthogonality of the starting vectors by rotating pairs of vector components by the block-diagonal matrix of 22 rotation0BBBBB@ 0 1 1 0 ! 0 0 . . . 0 ... 0 . . . 1CCCCCA (5.20) The resulting pair of groundstate wavefunctions satisfies hGS0jGS0i= hGS1jGS1i = 1 hGS0jGS1i = 0 hGS0jPxyjGS0i= hGS1jPxyjGS1i = cos2q sin2q jGS+i jGS i ! = cosq sinq sinq cosq ! jGS0i jGS1i ! (5.21) where Pxy and Pxy are reflections about the diagonal directions (a;a) or ( a;a) in terms of lattice parameters (Fig. 4.1a). The total momentum K = (p;p) is a high 78 symmetry point invariant to Pxy and Pxy. The eigenvalues are hGSjPxyjGSi= hGSjPxyjGSi=1 (5.22) which indicate p-symmetry. The K = (p;0) and K = (0;p) states are merely 0:17Jdd higher in energy. These two momenta are invariant to Px and Py, the reflection about (a;0) and (0;a), respectively. The lowest K = 0 state is a doublet of total spin of S = 1 (6 states in total counting all z-projections) and should be an S=1 excitation from the K = (p;p)GS doublet. Disregarding such S= 1 excitations, one can consider the K = 0 S = 0 state. This energy level is not degenerate, and the wavefunction can be classified using Pxy and Pxy. The expectation value of both operators are -1 so this state has a d symmetry. This is different from the t-t’-J model diagonalized in the N=32 cluster, wherein s-state is lower than the d-state. The charge correlation as function of hole-hole separation can be illustrated by the expectation value of ĉcharge(d ) =å p†l+e;s pl+e;s p†l+e+d ;s 0 pl+e+d ;s 0 (5.23) Ĉcharge(R) = å jd j=R ĉcharge(d ) (5.24) 1 = å R Ĉcharge(R) (5.25) In a finite cluster, the number of hole-hole configuration separated by a distance R, is limited by the periodic boundary condition to be smaller than that of L2 =N!¥ limit if R > L2 . Therefore, such correlations should be compared to åRP(R) = 1, the probability of two randomly distributed holes to be separated by R in the same finite cluster. For this purpose, the correlations can be gauged by DĈcharge(R) = Ĉcharge(R) P(R) (5.26) where the correlation is the same as that of a random distribution if DC(R) = 0. The expectation values as a function of hole-hole separation are shown in 79 0 1 2 3 4 R/a 0 0.1 0.2 C c ha rg e(R ) K=(pi,pi) ST=0 K=(pi,0) ST=0 K=(0,0) ST=0 Random 0 1 2 3 4 R/a -0.05 0 0.05 ∆C ch ar ge (R ) Figure 5.4: (top) Probability of charge separation, hCcharge(R)i for the lowest states in the CS = 12 subspace. (bottom) The difference from a random distribution. Fig. 5.4. DCcharge(R) indicates charge repulsion at K = (p;p), a small peak of correlation at R=2 for (p;0), and attraction for K = (0;0). However, one should note that two holes in a N=32 lattice corresponds to 6:25% doping, which could be in the pseudogap regime without the simultaneous charge and phase coherence needed for superconductivity. It is thus important to also study the spin correlation between the two charges. For the single-hole scenario in Chapter 4, the mobile carrier was found to be the S = 12 3SP, which is a composite of the oxygen spin and two neighboring copper spins. The two-hole solutions yielded hHJpd i 2 0:9Jpd , meaning HJpd can still provide a good characterization. Determining the spin-spin correlation between the two polarons requires phase information among the four copper spins. This is not a straight forward exercise because annihilating copper spins would take the wave- function outside of the håd†l+sdl+s i=N numerical basis formulated in Sec. 5.3. It is instructive to first consider the eight single-hole states of HJpd listed in Tab. 5.2. In the two-hole scenario, the oxygen holes can neighbor the same copper spin to form a 5-spin object, but Ccharge(R) shown in Fig. 5.4 shows that the probability is low. Ignoring these common-copper configurations, the wavefunction contains two 80 Wavefunction Total Spin hHJpd i Jpd j *i= q 1 3 j"#i+j#"ip 2 p†"
UBC Theses and Dissertations
Modeling polarons in transition-metal oxides Lau, Bayo 2011
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
- 24-ubc_2011_spring_lau_bayo.pdf [ 1.32MB ]
- JSON: 24-1.0071682.json
- JSON-LD: 24-1.0071682-ld.json
- RDF/XML (Pretty): 24-1.0071682-rdf.xml
- RDF/JSON: 24-1.0071682-rdf.json
- Turtle: 24-1.0071682-turtle.txt
- N-Triples: 24-1.0071682-rdf-ntriples.txt
- Original Record: 24-1.0071682-source.json
- Full Text