UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Modeling polarons in transition-metal oxides Lau, Bayo 2011

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2011_spring_lau_bayo.pdf [ 1.32MB ]
Metadata
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
24-1.0071682-fulltext.txt
Citation
24-1.0071682.ris

Full Text

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 analytical 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 showcased 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 Hamiltonian.” 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 antiferromagnet 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 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 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 computational 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 Chapter 5. I was responsible for the initial effort as well as all computational and theoretical 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  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  35  3 Octapartite Description of the Antiferromagnetic Square Lattice . .  38  2.5  3.1 3.2  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  38  Heisenberg Antiferromagnet on the Square Lattice . .  39  3.2.1  The Neel Picture . . . . . . . . . . . . . . . . . . . . . .  40  3.2.2  The Resonating Valence Bond Picture . . . . . . . . . . .  42  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  3.3  The  spin- 12  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  Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  57  4.3.1  Energetics . . . . . . . . . . . . . . . . . . . . . . . . . .  58  4.3.2  Wavefunction Analysis . . . . . . . . . . . . . . . . . . .  59  4.3.3  Further Experimental Implications . . . . . . . . . . . . .  62  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  63  5 Two Holes in the Spin Polaron Model . . . . . . . . . . . . . . . . .  64  4.3  4.4  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 assumptions, and the Fork space of n extra holes injected into a halffilled system with N unit cells. . . . . . . . . . . . . . . . . .  Table 2.1  7  Most-relaxed cut-off condition for the 1D breathing-mode polaron Green’s function computation. . . . . . . . . . . . . . .  23 29  Table 2.3  αn− ,n+ vs. n+ , n− for the ground state . . . . . . . . . . . . . . αn− ,n+ vs. n+ , n− for the first-excited state . . . . . . . . . . .  Table 3.1  Size of the Cs =  Table 2.2  1 4  subspace as compared to that of the com-  monly used basis. . . . . . . . . . . . . . . . . . . . . . . . . Table 5.1  29  47  (N=32) The S= 12 k = ( π2 , π2 ) single-hole ground state’s weight in subspaces of particular SA  ⊗  SB . Numbers are percentage  adding up to 100%. Two lines are drawn to highlight the CS =  1 4  truncation which discards states with SA , SB < 6 and yield an energy within 3.6% of the exact value. . . . . . . . . . . . . . Table 5.2  Single-hole eigenstates of HJpd .  pσ†  66  creates an oxygen hole and  the arrows in the ket indicate the spins of the two copper sites neighboring the oxygen hole. . . . . . . . . . . . . . . . . . .  vii  81  List of Figures Figure 1.1  Phase diagram of hole-doped cuprates (Temperature vs electronremoval per unit-cell)[3]. . . . . . . . . . . . . . . . . . . . .  Figure 1.2  (top) A CuO2 layer embedded in certain chemical host[4].(bottom) Lattice structure of single-layer cuprate HgBa2 CuO4 [12]. . . .  Figure 1.3  3 5  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 ε vectors (solid arrow) and the four δ vectors (dashed arrow) are also shown. . . . . .  6  Figure 1.4  Mean-field phase diagram for the t-J model[5]. . . . . . . . .  10  Figure 2.1  (a) GS energy and (b) GS qp weight as a function of the corresponding 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). Parameters are t = 2, Ω = 1.  Figure 2.2  . . . . . . . . . . . . . . . . . . . . .  25  (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, Ω = 1.  Figure 2.3  26  P(0) for t = 2, Ω = 1 for the breathing-mode (symbols) and Holstein (line) models, respectively. The solid and dashed lines correspond to ground state and first excited state, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  viii  27  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, Ω = 1, λB = 1.07.  Figure 2.5  . . . . . . . . . . . . . .  30  Ratio of effective polaron mass to that of the free electron. Circles 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, Ω = 1. . . . . . . . . . . . . . . . . . . .  Figure 2.6  31  GS qp weights for the breathing-mode (circles) and Holstein (squares) models. Full symbols correspond to α = 8. For comparison purposes, the empty symbols show the α = 4 results of Fig. 2.1(b). (Ω = 1) . . . . . . . . . . . . . . . . . . . . . . .  Figure 2.7  32  The spectral function of the 1D breathing-mode model, on linear (left) and logarithmic (right) scales. The curves are shifted for better viewing, and correspond (top to bottom) to λB = 1.5; 1.25; 1.0 and 0.75. Vertical lines indicates EGS and EGS + Ω. K = 0,t = 2, Ω = 1, η = 0.004Ω.  Figure 2.8  . . . . . . . . . . . . .  33  A(k, ω ) vs ω , for K/π =0; 0.25; 0.5; 0.75 and 1 (top to bottom) for intermediate coupling t = 2, Ω = 1, λB = 1.25. η /4 = ∆E = 0.001Ω. 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 + Ω. . . . . . . . . . . . .  Figure 3.1  34  Square lattice divided in octads (shaded areas). Sites positioned similarly inside octads have the same label. Each is surrounded by neighbors with different labels. . . . . . . . . .  Figure 3.2  (a) Overlap |⟨GS,Cs |GS,Cs +  2 2 N ⟩|  45  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 = 41 ; (b) Fractional change  1 dEGS EGS dCs  in ground-state energy correspond-  ing to different Cs cutoffs. The lines are linear fits to the data.  ix  48  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. [37]In√ set: expectation values of the non-commuting operators ⟨m2 ⟩ and ⟨Sr · Sr+( L , L ) ⟩. . . . . . . . . . . . . . . . . . . . . . . . 2 2  Figure 4.1  49  (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 ε vectors (solid arrow) and the four δ vectors (dashed arrow) are also shown. (b) Sketch of a virtual process of Tswap . . . . . . . . . . . . . . . . . . .  Figure 4.2  a) Energy and b) quasiparticle weight (bottom) for the lowest eigenstates with ST =  1 2  and  3 2  vs. momentum. Different sets  are shifted so as to have the same GS energy. . . . . . . . . . Figure 4.3  54  58  ⟨Cx (δ , ∆)⟩ for the lowest energy state at (a) ( π2 , π2 ) with ST = 21 , and (b) at (π ,π ) with ST = 32 . The darkly-shaded bullet denotes the oxygen 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. ⟨Cy (δ , a)⟩ is the Pˆx↔y reflection. . . . . . . . . . . . . . . . .  60  Figure 4.4  ⟨Cx/y (δ , a)⟩ for lowest state at k=(0,π ) 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. . . . . . . . . . . . . . . . . 1 4  Figure 5.2  Lowest two-hole states in the CS =  Figure 5.3  Convergence of the ST = 0 lowest energy states. The two-  subspace. . . . . . . . .  66 76  hole energy level corresponding to ∆Eb = 0 is indicated by the horizontal dash line. . . . . . . . . . . . . . . . . . . . . . .  x  77  Figure 5.4  (top) Probability of charge separation, ⟨Ccharge (R)⟩ for the lowest states in the CS = 21 subspace. (bottom) The difference from a random distribution. . . . . . . . . . . . . . . . . . . . . .  Figure 5.5  80  (top) Probability of singlet three-spin polaron (3 SP) pair, ⟨C3SP,singlet (R)⟩ for the lowest states in the CS =  1 2  subspace. (bottom) The dif-  ference from a randomly distributed paramagnetic configuration. 83 Figure 5.6  c3SP,singlet (δ ) × 104 for one of the two K = (π , π ) GS in the CS =  1 2  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 3 SP 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 3 SPs both centered on x-rung O sites. . . . . . Figure 5.7  c3SP,singlet  (δ )×104  for the lowest K = (π , 0) state in the CS =  84  1 2  subspace. Cu and O are marked by circles and rectangles, respectively. Values for the K = (0, π ) 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 3 SP 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 3 SPs both centered on x-rung O sites. . . . . . Figure 5.8  ccharge  (δ ) × 104  for the lowest K = (π , 0) state in the CS =  86  1 2  subspace. Cu and O are marked by circles and rectangles, respectively. Values for the K = (0, π ) 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 3 SP pair configurations are shown. For example, the bottom figure shows that ccharge (2a, 0) = 0.0289 for two 3 SPs both centered on x-rung O sites. . . . . .  xi  87  Figure 5.9  c3SP,singlet (δ ) × 104 for the lowest K = (0, 0) ST = 0 state in the CS =  1 2  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 3 SP 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 (0, −2a) = 0.0367 for two 3 SPs both centered on x-rung O sites. . . . . . . . . . . . . . . . . . . . . . . . . .  xii  88  Glossary ZRS  Zhang-Rice singlet  AFM  antiferromagnet  ARPES  Angular-resolved photoemission spectroscopy  RVB  resonating valence bond  BCS  Bardeen, Cooper, and Schrieffer  3 SP  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 peculiar 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 education. 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 materials 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 materials, 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[8]. 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 dispersion 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 levels 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 quasiparticles 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 ∆(k) = ∆0 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 counterexamples is the class of cuprate compounds, whose properties contradict with all of the above statements[4]. The undoped “parent” cuprate compound is an insulator 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 upand 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 superconducting 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 ∆(k) of  dkx2 −ky2 symmetry, which is different from the s-symmetric constant predicted due to electron-phonon interaction in conventional theory. For the “underdoped” compounds with electron concentration between that of the “parent” and “optimally doped” compounds, an extra excitation energy scale ∆PG > ∆SC has been measured 2  Figure 1.1: Phase diagram of hole-doped cuprates (Temperature vs electronremoval per unit-cell)[3]. by various experiments. This is another deviation from the phonon-driven superconductor theory, which lacks an extra higher energy scale separating the normal and superconducting state.  1.1  Hole-Doped Cuprates  Since their discovery in 1986 [9], cuprates have been classified as a high-TC superconductor, but their importance to condensed matter theory lies in the many anomalous 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 compound, 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 ordering temperature compared to a 3D antiferromagnet is due to the weak coupling 3  between the strongly correlated 2D AFM planes. The superconducting phase centered at ∼0.15% doping is characterized by a BCS gap ∆(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 inplane conductivity showed no gapped behavior. Angular-resolved photoemission spectroscopy (ARPES) measurements show excitations ∆PG (k) of the same trend as ∆(k) with maximum at (0, π ) and zero for kx = ky [6]. Breaking of spatial symmetries have also been observed[10]. The connection between different regions has not been elucidated[11]. The mysteries stem from cuprates’ unique structure, which differs greatly from the three-dimensional isotropic lattice assumed by conventional theories. The common 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 HgBa2 CuO4 compound with tetragonal crystal structure and one copper-oxygen layer per unit cell[12]. Although in general TC is higher for complicated materials with multiple 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−x Srx CaO4 , exhibit lattice buckling at low temperature[4], HgBa2 CuO4 does not exhibit lattice distortion[12]. Lattice distortion 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]2s2 2p6 ). 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 different spatial symmetries: E(dxy ) < E(dxz ) = E(dyz ) < E(dz2 −r2 ) < E(dx2 −y2 )[4]. 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 chemical host[4].(bottom) Lattice structure of single-layer cuprate HgBa2 CuO4 [12].  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 ε vectors (solid arrow) and the four δ vectors (dashed arrow) are also shown. copper dx2 −y2 orbital. The material is in the charge-transfer regime in the ZaanenSawatzky-Allen classification scheme, and a hole would have strong character of the oxygen ligands[13]. Although there are three orthogonal orientations per 2p oxygen level, the dx2 −y2 wavefunction overlaps predominantly with the σ orientation 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 2pσ 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 + εx/y , where εx/y = 12 ax/y . Since there are more electrons than vacancies in the problem, the model is specified with fermion operators p†l+ε ,σ and dl,†σ creating holes of spin  σ on a vacuum with one hole per unit cell. For the orbital setting illustrated in Fig.1.3, the model reads H3B = Tpd + Tpp + ∆ pd ∑ nl+ε ,σ +Upp ∑ nl+ε ,↑ nl+ε ,↓ +Udd ∑ nl,↑ nl,↓ , 6  (1.1)  Model 3-band 2-band Ch. 3-5 1-band t-J  Approximations dx2 −y2 on Cu and two pσ bands dx2 −y2 and one pσ band ⟨nd↑ + nd↓ ⟩ = 1 on each Cu and two pσ bands dx2 −y2 band each pσ hole is locked to a Cu spin as a ZRS  Number of Fermions N +n N +n N +n N +n 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 = t pd ∑(−p†l+ε ,σ + p†l−ε ,σ )dl,σ + h.c.  ′ Tpp = t pp ∑ sδ p†l+ε +δ ,σ pl+ε ,σ − t pp ∑ p†l−ε ,σ + p†l+3ε ,σ pl+ε ,σ  nl,σ nl+ε ,σ  = dl,†σ dl,σ = p†l+ε ,σ pl+ε ,σ .  Tpd is the direct oxygen-copper hopping. Tpp contains the direct t pp oxygen-oxygen ′ oxygen-oxygen next-nearest-neighbor nearest-neighbor hopping as well as the t pp  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, sδ = 1 (sδ = −1). ∆ pd is the charge transfer energy, which is the energy cost of the d 9 → d 10 L 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 t pp < t pd < ∆ pd < Upp < Udd [5]. In practice, neighbor Coulomb interactions such as Upd d † d p† 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 determining 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[4]. 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´eel 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[36], described by the Heisenberg quantum AFM coupling spins of neighboring copper sites of a 2D square lattice: [ )] 1( + − z z − + HJ = J ∑ Si S j + S S + Si S j = J ∑ Si · S j 2 i j ⟨i, j⟩ ⟨i, 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[37]. 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 significantly lower local energy[16]. 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 the U → ∞ limit of the one-band Hubbard model[5]. [  ] ( ) HtJ = Pˆ ∑ ti j c†j,σ ci,σ + h.c. Pˆ + J ∑ Si · S j  (1.3)  ⟨i, j⟩  i, j  where the c†i,σ creates a spin-σ 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  for the undoped parent insulating  ARPES  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 because 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-doubleoccupancy 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 combination 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 valence bond (RVB) ansatz to an undoped  AFM  and the  BCS  wavefunction[40]. 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,σ = fi,†σ bi  (1.4)  † † fi,↑ fi,↑ + fi,↓ fi,↓ + b†i bi = 1  (1.5)  if the restriction of  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[5]. 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  region for TB < T < TRV B , the  state with d-symmetry below TRV B . In the left  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 < TRV B , TB is explained as a  RVB  spin  singlet state with Bose condensation of the charge carrying holons. To the right for TRV B < 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 superconducting area for TRV B , TB < T < TD , the non-Fermi-liquid region is attributed to spinons coherence without singlet formation. The highest TC predicted by meanfield 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 10  AFM  phase in that region. In this  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 linearly with temperature. It has been suggested that coupling to phonons might be the cause of this broadening[43], 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[44]. For an electron-phonon coupling parameter much larger than that predicted by density functional theory, a wide peak broadening that decreases linearly with temperature was found[45]. 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 limitations. 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[49], and spatial symmetry breaking[50] have not been elucidated. Recent experiments have raised further questions regarding the understanding 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 pseudogap 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)[57] and x-ray scattering [58] remain either open questions or are controversial. 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 carrierphonon 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 patible with the c†i,σ  =  fi,†σ bi ,  BCS  RVB  nature is com-  wavefunction[40], the robustness of spin-charge separation  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 [13]. 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 [16]. 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 framework of the Holstein model[61]. 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 microscopic 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 when models 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 dissertation. 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 details 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 electrons 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 disturbance to the  AFM  background. This spin polaron problem does not enjoy an  exactly known background like the bosonic b− |0⟩ = 0 in a carrier-phonon interaction scenario. Chapter 3 discusses progress in dealing with this issue. A novel octapartite numerical approach is devised to capture the two-dimensional Heisenberg antiferromagnetic spin background with breakthrough efficiency. A surprising feature is revealed during the application process of the aforementioned 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- 23 polarons become energetically favorable in different regions of the Brillouin zone. The chapter covers 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 framework. The computational approach from Ch. 3 is utilized to bypass various technical 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 (quantized lattice vibrations) leads to the formation of polarons. This mechanism is a key ingredient in the physics of the manganites[8], Bechgaard salts[62][63], and, possibly, of the cuprates[43][64]. There are various model Hamiltonians describing the coupling of the particle and bosonic degrees of freedom. The asymptotic limits 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 function 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 asymptotic 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.[61] It is essentially the tight-binding model with an on-site energy proportional to the lattice displacement Xi = HH = −t  √ 1 (b† + bi ): 2MΩ i  ∑  <i j>  (  ) c†i c j + c†j ci + Ω ∑ b†i bi + g ∑ ni Xi i  (2.1)  i  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, Ω is the frequency, M is the atomic mass, and g is the electron-phonon coupling strength. The model has been widely studied numerically by Monte-Carlo calculations[65–76], variational methods[77][78–89], and exact diagonalization[90–95]. Analytic approximations, such as the momentumaverage approximation, have also progressed over the years[96–99]. For some materials, a more appropriate model is provided by the breathingmode 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 −8td2p , where td p 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.[16] 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 td p 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[100][101]. The breathing-mode Hamiltonian describes the physics of the linear modulation of on16  site energy. While this breathing-mode model is motivated as a 2D model, here we investigate 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 polaron 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 systems 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 considerable 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 ∑ c†i ci+1 + c†i+1 ci + Ω ∑ b†i+ 1 bi+ 1 i  +√  g 2MΩ  i  2  ( ) † † ∑ ni bi+ 1 + bi+ 1 − bi− 1 − bi− 1 i  2  2  2  2  2  (2.2)  The notation is the same as before, except that now the phonons live on an interlaced lattice. The difference between the two models is more apparent in momentum space. The Holstein model has constant coupling to all phonon modes ( ) g VH = √ √ c†k−q ck b†q + b−q , ∑ N 2MΩ kq whereas the breathing-mode model has a coupling strength that increases monotonically with increasing phonon momentum ( ) 2ig q VB = √ √ sin c†k−q ck b†q + b−q . ∑ 2 N 2MΩ kq  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 (−π , π ] (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 [102], and an investigation based on the self-consistent Born approximation [103], 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[104] O˜ = eS Oe−S . Using SH = and SB =  (2.3)  g √ ni (b†i − bi ), ∑ Ω 2MΩ i  g ∑ ni (−b†i− 12 + bi− 12 + b†i+ 12 − bi+ 12 ) Ω 2MΩ i √  respectively, the diagonal forms of the Hamiltonians are, in terms of the original (undressed) operators: H˜ H = T˜H + Ω ∑ b†i bi − i  g2 n2i 2MΩ2 ∑ i  2g2 H˜ B = T˜B + Ω ∑ b†i+ 1 bi+ 1 − ni (ni − ni+1 ) 2 2MΩ2 ∑ 2 i i  18  (2.4)  (2.5)  where the kinetic energies are: g2  T˜H = −te− 2MΩ3  g(−b†i+1 +b†i ) √ Ω 2MΩ  ∑  c†i+1 ci e  g2  ∑ c†i+1 ci e Ω  −  e  g(−bi+1 +bi ) √ Ω 2MΩ  + h.c.  i  T˜B =  −te−3 2MΩ3  √g (b† −2b† 1 +b† 1 ) 2MΩ i+ 3 i+ 2 i− 2 2  ×  i  −  e  √g (b −2bi− 1 +bi− 3 ) Ω 2MΩ i+ 12 2 2  + 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: g2 2MΩ2 g2 (0) EB (k) = −z 2MΩ2 (0)  EH (k) = −  Each model has three energy scales, therefore the parameter space can be characterized by two dimensionless ratios. It is natural to define the dimensionless (effective) coupling as the ratio of the lattice deformation energy to the free-electron ground-state energy −zt: g2 1 2MΩ2 zt g2 1 λB = , 2MΩ2 t  λH =  (2.6) (2.7)  where z is also the number of nearest neighbors in the electron sublattice. It should √ be noted that since Ω ∼ 1/ M, the λ ’s do not depend on the ion mass, M. λ 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:  α=  zt . Ω  (2.8)  Using standard perturbation theory[95], the first-order corrections to the energy of the lowest state of momentum k are: EH (k) = −2te−αλH cos(k), (1)  (1) EB (k)  = −2te  − 32 αλB  cos(k),  (2.9) (2.10)  showing that the polaron bandwidth is exponentially suppressed in the strong coupling 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 overlaps 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 instead of just 2, as in the case for Holstein. The second-order corrections are: t2 (2) EH (k) = −2 e−2αλH fH (k, α , λH ), Ω t2 (2) EB (k) = −2 e−3αλB fB (k, α , λB ). Ω The functions f can be written in the form fH,B = AH,B + BH,B cos(2k) with AH  = Ei(2αλH ) − ln(2αλH ) − γ  BH  = Ei(αλH ) − ln(αλH ) − γ 20  and AB = Ei(3αλB ) − ln(3αλB ) − γ BB = Ei(2αλB ) − 2Ei(αλB ) + ln(  αλB ) + γ. 2  Here, γ is the Euler-Mascheroni constant and Ei(x) is the exponential integral with the series expansion  ∞  xn . n=1 n!n  Ei(x) = γ + ln(x) + ∑  x The result can be further simplified in the limit αλH , αλB ≫ 1 using limx−>∞ ∑∞ n=1 n!n ∼ n  ex x.  This leads to the simplified expressions: [ ] t2 1 −αλH +e cos(2k) ∼ −2 αλH 2 ] [ t 2 1 e−αλB (2) EB (k) ∼ −2 + cos(2k) αλB 3 2  (2) EH (k)  (2.11) (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 < π , 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. [77]. 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[77] m b†n 1 j+m iK j † |S, K⟩ = √ ∑ e c j ∏ √ |0⟩, N j mε {S} nm !  (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 −te±iK , Ωn, g √ and ± √2MΩ 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. [77] was optimized for computation 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 |ϕ ⟩LF in the Lang-Firsov basis correspond to a state |ϕ ⟩R = e−S |ϕ ⟩LF  (2.14)  in real space. With our choice for the S operators, this reverse transformation induces a phonon coherent state structure at the electron site (Holstein), respectively the two bracketing phonon sites (breathing-mode). The phonon statistics of the coherent state obeys the Poisson distribution. In the anti-adiabatic regime (zt > Ω), the splitting due to the hopping (off-diagonal hopping matrix elements) is significant compared to the diagonal matrix elements proportional to Ω. 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 further 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 1-11 12-13 14-15 16-17 18-20 21-30 31-40 41-50  |m|max 22.5 - (# phonons) 10.5 9.5 8.5 7.5 5.5 4.5 3.5  # States 17053356 16871582 28274774 33423071 41757650 42628080 38004428 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 lowestenergy 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 memory and fflops needed for such computations are formidable, especially to extensively investigate the multi-dimensional parameter space (λ , α , 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 typically 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:[107] G(ω , k) = ⟨k|  1 |k⟩, ω − H + iη  23  (2.15)  where |k⟩ = c†k |0⟩. This can be written as the solution of a linear system of equations: G(ω , k) = ⟨k|y⟩ (ω − H + iη )|y⟩ = |k⟩ One can iteratively tri-diagonalize H by the vanilla Lanczos process:[108] H = QT Q† (ω + iη − T )Q† |y⟩ = Q† |k⟩ If the RHS is of the form [100...]T , Crammer’s rule can be used to express G = ⟨k|QQ† |y⟩ as a continuous fraction in terms of the matrix elements of the tridiagonal matrix (ω + iη − T ). In particular, this condition is achieved by picking the initial Lanczos vector to be |k⟩. This method is efficient because it does not require the complete solution of the linear system nor of the eigenvalue problem. 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 [109]. We perform the vanilla Lanczos tri-diagonalization and re-orthogonalize each state in Q against the starting vector |k⟩ to validate the continuous 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 = |⟨ΦGS |c†k=0 |0⟩|2 , where |ΦGS ⟩ is the ground-state eigenfunction, are shown in Fig. 2.1 for the 1D breathing-mode and Holstein models. For a fixed value of α , we see the expected crossover from a large polaron (at weak coupling λB,H ) to a small polaron (at strong coupling λB,H ), signaled by the collapse of the qp weight. 24  Z0  EGS -4  1 (a)  (b)  -5  0.8  -6  0.6  -7  0.4  -8  0.2  -9 0  1 λΗ,λΒ  2  0 0  1 λΗ,λΒ  2  Figure 2.1: (a) GS energy and (b) GS qp weight as a function of the corresponding 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). Parameters are t = 2, Ω = 1. The ground state energy of both models decreases monotonically with increasing 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 λB ≈ 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  E1-EGS 1.1  Z0 0.5 (a)  (b)  1.0  0.4  0.9 0.3 0.8 0.2 0.7 0.1  0.6 0.5 0  1 λΗ,λΒ  2  0 0  1 λΗ,λΒ  2  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, Ω = 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 Ω 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[103]. The separation between the two lowest energy states now first decreases and then increases back towards Ω as λH,B → ∞. This behavior is wellknown for the Holstein polaron[77]. 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 Ω. 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  1  <P(0)>  0.8 0.6 0.4 0.2 0 0  0.5  1 λΗ , λΒ  1.5  2  Figure 2.3: P(0) for t = 2, Ω = 1 for the breathing-mode (symbols) and Holstein (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[77]. The nature of these states is revealed by checking the locality of the phonon cloud. We define the projection operator P(K) =  ∑ |S, K⟩⟨S, K|  Slocal  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 structure 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 groundstates ⟨P(0)⟩ ∼ 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 ⟨P(0)⟩ → 0, precisely because the free phonon can be anywhere in the system. When the second bound states form, ⟨P(0)⟩ becomes large, showing that phonons in these states are primarily localized near the electron[77]. 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 Ω alone; therefore, we seek an effective transformation with ˜ S(∆) = SB | Ωg =∆ The computed eigenstates |ϕ ⟩ are projected into such structure αn− ,n+ by (b†l− 1 )n− (b†l+ 1 )n+ ∞ 1 iKl † √ ∑ e cl ∑ αn− ,n+ √ 2 √2 |0⟩ n ! n+ ! N l − n− ,n+ =1 ˜  = eS(∆) P(K)|ϕ ⟩. 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 |GS⟩ ∼ e−S(t,g,Ω,K) √1N ∑l eiKl c†l |0⟩, respectively |1⟩bound ∼ e−S(t,g,Ω,K) √1N ∑l eiKl c†l (eiθ b†l− 1 − ˜  ˜  e−iθ b†l+ 1 )|0⟩ for some phase θ (t, g, Ω, K) needed to satisfy time-reversal symme2  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 momentum, this description is valid as long as the excited state remains bound, with energy less than EGS + Ω. Fig. 2.4 shows momentum dependent energy and qp weight for the two lowest eigenstates of the breathing-mode model, for an intermediate coupling strength 28  2  t =2Ω=1g=0∆=0 0 1 2 3 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 t = 2 Ω = 1 g = 1.5 ∆ = 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 Ω = 1 g = 1.964 ∆ = 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 n+ n− 0 1 2 3  Table 2.2: αn− ,n+ vs. n+ , n− for the ground state  t = 2 Ω = 1 g = 1.5 ∆ = 1.05 0 1 2 3 0.0228 -0.6515 0.0152 0.1151 0.6515 0.0062 0.1700 -0.0473 0.0152 -0.1700 0.0505 -0.0456 0.1151 -0.0473 0.0456 -0.0220 t = 2 Ω = 1 g = 1.964 ∆ = 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 n+ n− 0 1 2 3  Table 2.3: αn− ,n+ vs. n+ , n− for the first-excited state  29  Zk  Ek/Ω -4.2 -4.4  0.4 (a)  0.3  (b)  -4.6 0.2 -4.8 0.1  -5.0 -5.2 0  0.2 0.4 0.6 0.8 k/π  1  0 0  0.2 0.4 0.6 0.8 k/π  1  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, Ω = 1, λB = 1.07.  λB = 1.07. The polaron band has a maximum at k ∼ 0.4π , 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[103]. For the Holstein polaron, the polaron dispersion is a monotonic function of momentum[77]. Even though the qp weights remain moderately high at zero momentum, the weights collapse 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 contribution 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 λH,B . These were calculated from the second derivative of the energy at momentum K = 0. For the GS of both models, the effective mass increases monotonically with λH,B . At weak couplings, the breathing-mode polaron is lighter than the Holstein polaron. As already discussed, this is due to  30  m/m0 1000  100  10  1 0  0.5  1 λΗ,λΒ  1.5  2  Figure 2.5: Ratio of effective polaron mass to that of the free electron. Circles 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, Ω = 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 λH,B . This can be understood through the link of the effective mass and the qp weight. In terms of derivatives of the self-energy Σ(k, ω ), the effective mass m∗ is given by: ( ) ( )−1 m∗ ∂Σ m ∂ 2Σ = 1− · 1+ 2 2 , m ∂ω h¯ ∂ k  31  Zk=0 1 0.8 0.6 0.4 0.2 0 0  0.5  1 λΗ , λΒ  1.5  2  Figure 2.6: GS qp weights for the breathing-mode (circles) and Holstein (squares) models. Full symbols correspond to α = 8. For comparison purposes, the empty symbols show the α = 4 results of Fig. 2.1(b). (Ω = 1) where derivatives are evaluated at K = 0 and at the corresponding eigenenergy. The first term is linked to the qp weight, Z = (1 − ∂∂ωΣ )−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 α = 4. For higher α (lower Ω 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 λ for increasing α [83, 98]. At weak and moderate coupling, the qp weight and the effective mass (not shown) in the breathing-mode model are much less sensitive to an increase in α than is the case for the Holstein polaron. This suggests that breathingmode polarons should be better charge carriers than the Holstein polarons in this regime.  32  A(k=0,ω)  -6  ln A(k=0,ω)  -4 ω/Ω  -7  -2  -6  -5 ω/Ω  -4  -3  Figure 2.7: The spectral function of the 1D breathing-mode model, on linear (left) and logarithmic (right) scales. The curves are shifted for better viewing, and correspond (top to bottom) to λB = 1.5; 1.25; 1.0 and 0.75. Vertical lines indicates EGS and EGS + Ω. K = 0,t = 2, Ω = 1, η = 0.004Ω.  2.4.2 Spectral Function The spectral function is proportional to the imaginary part of the Green’s function: 1 A(k, ω ) = − ImG(k, ω ) π  (2.16)  In terms of single electron eigenstates and eigenfunctions H|n⟩ = En |n⟩, we obtain the Lehmann representation: A(k, ω ) = ∑ |⟨n|c†k |0⟩|2 δ (ω − En ) n  Of course, since we use a finite small η in numerical calculations, the δ -functions are replaced by Lorentzians of width η [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 en33  A(k=0,ω)  -6  ln A(k=0,ω)  -4  -2 0 ω/Ω  2  4  -6  ω/Ω  -5  Figure 2.8: A(k, ω ) vs ω , for K/π =0; 0.25; 0.5; 0.75 and 1 (top to bottom) for intermediate coupling t = 2, Ω = 1, λB = 1.25. η /4 = ∆E = 0.001Ω. 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 + Ω. ergy. Results corresponding to four different coupling strengths λB 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 + Ω. As the coupling λB 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 λB 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 systematically suppressed. This is qualitatively similar to the behavior exhibited by Holstein polaron[99]. 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 transferred 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 + Ω for all values of K. This is part of the kink-like structure reported in Ref. [103]. The logarithmic plot clearly reveals a non-monotonic dispersion 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 polaronplus-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 conclude 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, |⟨ϕ |c†K |0⟩|2 > 0, see Eq. (2.15). These eigenstates must have components 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[103] 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, becomes 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 breathingmode polaron must (i) create a new polaron cloud at site i + 23 ; (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 polaronplus-free-phonon continuum at EGS + Ω, and the appearance of a second bound state with finite qp weight for large-enough couplings. The convergence of numerics points to the existence of additional bound states, whose absence from the spectral 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 because 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[6][7]. 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 — remains 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− |0⟩ = 0 and b+ |n⟩ = n + 1|n + 1⟩. 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 disadvantages. For example, the powerful Monte Carlo (MC) methods give extremely accurate 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 variational RVB-type ans¨atze[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 somewhat biased schemes[17, 120, 121]. ED studies have the huge advantage that they provide the GS wavefunction, from which any GS properties, including all correlation 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 aforementioned 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. Although Bonca et. al. did successfully adapt the idea for a carrier in the spin background [34], the method has left something to be desired because its N´eel starting point cannot capture a real AFM background. The strong quantum fluctuations leading to strong deviations from a N´eel-like background are believed to be an essential part of the low-energy physics of the doped systems. It is therefore necessary 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 superexchange mechanism[36], which leads to AFM Heisenberg exchange among nearest neighbor S =  1 2  copper spins arranged in the square lattice. 39  [ )] 1( + − z z − + HAFM = ∑ Si S j + S S + Si S j = ∑ Si · S j . 2 i j ⟨i, j⟩ ⟨i, j⟩  (3.1)  Here, a −2 N4 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 |σl ⟩l specifying the z-projection of spin at the copper site l.  ∑ cσ ∏ |σl ⟩l σ  (3.2)  l  The total spin ST ∈ [0, N2 ] is conserved with degenerate z-projections so the full spectrum can thus be obtained from the STz = ∑l σl = 0 basis, which contains N!/( N2 !)2 states without exploiting spatial symmetries. Efficient matrix-vector multiplication can be performed using Algorithm 151 [122] as discussed in Appendix B. Although the exact ground state wavefunction remains unknown, Marshall showed very early that the groundstate wavefunction must have ST = 0 and that its cσ in Eq. 3.2 obeys cσ = (−1)nσ |cσ |  (3.3)  where nσ is the number of up spins in one sub-lattice[123]. This is now known as the “Marshall sign rule” in the community. Among the many routes in treating an undoped AFM[37], 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 [37].  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  z = a+ la ala = S − Sla  nlb  z = b+ lb blb = S + Slb  Eq. 3.1 can then be expanded in powers of  nl 2S  due to the  (3.4) √ 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  ) ( 1 1 − (cos(qx ) + cos(qy ))2 αq+ αq + βq+ βq 4  (3.5)  where α + and β + 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 introduces three magnon-magnon interaction terms: α p+ α p αq+ αq , β p+ β p βq+ βq , and  α p+ α p βq+ βq . 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[37]. 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  N N SA = , SAz = ± 4 4  ⟩  N N SB = , SBz = ∓ 4 4  ⟩ .  (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: (−1)SA −mA ⟨SA , mA |⟨SB , mB |S = 0, Sz = 0⟩ = √ δS ,S δm ,−m 2SA + 1 A B A B The ground state’s weight in a Neel state is thus at most  1 N/2+1  (3.7)  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 [34] 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  N 2  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  N 2  pairs and no-double-occupancy constrains[40]. |RV B⟩ = ∏(1 − nl↑ nl↓ )PN/2 |BCS⟩  (3.8)  l  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[115][116]. The top-down approach is by energy minimization by tweaking the unprojected wavefunction |BCS⟩ in Eq. 3.8 to other related ansatz[41][117][42][5]. 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 starting 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[123], one can start with the bipartite division and write down all relevant basis states as Clebsch-Gordon series (−1)SA −mA |S = 0, Sz = 0⟩ = ∑ √ δS ,S δm ,−m |SA , mA ⟩|SB , mB ⟩ 2SA + 1 A B A B  (3.9)  for SA/B ∈ [0, N4 ], mA/B ∈ [−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 α , and the GS wave function can be expressed as a superposition |GS⟩ = ∑ cα |0, 0⟩α  (3.10)  α  I then noticed one characteristic shared by all α spin states with high values of 43  |cα |. 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 ⟨mA ⟩ and ⟨mB ⟩ in the zero-field limit. In the absence of a staggered field, this quantity has alternative definitions, such as: 1 mˆ2 = 2 N  (  ∑(−1)  |r|  )2 Sˆr  r  )2 ( SˆA − SˆB = , N2  (3.11)  1 where SˆA/B are the total sublattice spin operators. In the AFM GS, m = ⟨mˆ2 ⟩ 2 √ has been extrapolated to ∼ 0.3 as N → ∞, and increases like 1/ N to ∼ 0.45 for N = 32[37]. In terms of the total spin Sˆ = SˆA + SˆB , we have:  1 mˆ2 = 2 (2SˆA2 + 2SˆB2 − Sˆ2 ) N  (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 physical meaning is that the maximum  N 8  N 16  of their maximum values. The  spins add to a total spin of zero while the rest adds to  N 3 16 .  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 computational 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 N 4,  N 2  sublattice spins add to the maximum of  if we want to include states with spin down to  44  3N 16  =  N 4  N − 16 , the basis must  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 0 while the other  3N 8  spins have a total spin of  3N 16 .  N 8  spins have a total spin of  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  3N 16  “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. SˆA = Sˆ0 + Sˆ1 + Sˆ2 + Sˆ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[124].  45  We can now generalize Eq. 3.9 and write singlets of the total lattice as: 7  |0, 0, f ({Si })⟩ = ∑ c{Si ,mi } ∏ |Si , mi ⟩  (3.13)  i=0  where Si , mi are the quantum numbers for Sˆi , i.e. the total spin for group i = 0, 7. Here, c{Si ,mi } is the product of the appropriate Clebsch-Gordon coefficients for the 8 pairs of quantum numbers (Si , mi ). The f ({Si }) 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 N exact. The total spin of each group takes values Si ∈ [0, 16 ], 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  N 16  of the maximum value.  We therefore parameterize the singlet Hilbert space in terms of a completeness parameter, Cs ∈ [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  N 8.  Based on the  discussion for the staggered magnetization, we expect this to already be a good variational basis for the GS. For any Cs ̸= 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 enforced here in terms of angular momentum addition, very different from the classical N´eel 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 16 32 64  Cs = 14 50 11042 9.8 × 108  Table 3.1: Size of the Cs = used basis.  S=0 1430 35357670 5.55 × 1016 1 4  Sz = 0 12870 601080390 1.83 × 1018  Full 216 232 264  subspace as compared to that of the commonly  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 number 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  N 2 -spin  subspaces, which are in turn built from two identical  spaces, which are formed by two identical ⟨α |Si · S j  |α ′ ⟩,  N 8 -spin  N 4 -spin  sub-  subspaces. The matrix element,  between two N-spin configurations α and α ′ can be calculated from  the Slz , Sl+ , and Sl− 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 Slz , Sl+ , and Sl− 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[122], 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  N 2 -spin  enlarges so the basis is enumerated as consecutive blocks of decreasing SA . The ( )2 nA number of states within each block is 2SAS+1 , 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  N 2  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  1 0  a)  -1  10  dE/EdC  Overlap  0.9  -2  10  0.8  N=16 N=32 N=64  0.7 0.6 0  b)  10  -3  10  -4  10  -5  0.5  CS  10 0.0  1  1.0  0.5  CS  Figure 3.2: (a) Overlap |⟨GS,Cs |GS,Cs + N2 ⟩|2 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 = 41 ; (b) FracGS tional change E1GS dE 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 =  1 4  basis size to those  of other bases commonly used by non-iterative methods. At Cs = 14 , the N = 32 system was solved by a single CPU within seconds using a generic Lanczos routine. To convince the critics, the single-processor computation was modified to calculate 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 the CS =  1 8  solution,  the GS vector was obtained after 15 matrix-vector products, each taking less than 30 minutes without any optimizations.  48  EGS/N  -0.60  0.56 0.24 0.52  0.2 0.16 0.12  m  -0.55  N=16 N=32 N=64  <Sr Sr+(L/2,L/2)>  -0.50  0.48 0  -0.65  0.5 CS  10  1  0.5 CS  -0.70 -0.75 0.00  0.25  0.50  0.75  1.00  CS 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. [37]Inset: expectation values of the non-commuting operators ⟨m2 ⟩ and ⟨Sr · Sr+( L , L ) ⟩. 2 2  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 =  1 2  (Cs = 14 ). While  the computations with Cs < 1 are variational, the goal here is to demonstrate that Cs ∼  1 4  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 of Cs . The deviation decreases exponentially with increasing Cs , and for Cs =  1 4  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  dEGS dCs  ∼ e−αCs . The error decays exponentially, with  a rate that increases with N. The N = 32, 64 lines cross at Cs ∼  1 4  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  dy dx  ∼ 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 − 2 scaling, we 3  extrapolate limN→∞ ENGS = −0.6671, within 0.5% of the best published value of -0.6692[37], even though it is achieved at a significantly reduced computational cost. We have thus far demonstrated the exponential convergence of the GS eigenvalue and eigenvector with increasing Cs , which validates our proposed variational basis. While a systematic approximation of the wavefunction can provides estimates of all quantities – more accurately so for those commuting with the Hamiltonian – 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 characterized 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) evaluated 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 of Cs ,  d dCs ⟨Sr · Sr+( L2 , L2 ) ⟩  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 wavefunction 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 description 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 wavefunction 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,  N 8  spins add to a total spin of zero while the other 38 N  add to the maximum value of  3 16 N.  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 that Cs ∼  1 4  N 8 -spin  idea was used as an apriori argument in postulating  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[121]. 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 pointgroup 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 efficient kernel for iterative methods, such as renormalization and quantum cluster theory[125]. Of course, this insight might also be adapted in devising more efficient variational Monte-Carlo formulations. From a holistic point of view, the  N 8 -spin  observation contributes to the de-  scription of the AFM background in analogy to the b− |0⟩ = 0 statement taken for granted by a phonon problem. Excited states analogous to b+ |0⟩ are indeed important 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 antiferromagnet (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 additions. 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 2pσ 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 + εx/y , where εx/y = 12 ax/y . Since there are more electrons than vacancies in the problem, the model is specified with fermion operators p†l+ε ,σ and dl,†σ creating holes of spin  σ 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 + ∆ pd ∑ nl+ε ,σ +Upp ∑ nl+ε ,↑ nl+ε ,↓ +Udd ∑ nl,↑ nl,↓ ,  (4.1)  with Tpd = t pd ∑(−p†l+ε ,σ + p†l−ε ,σ )dl,σ + h.c.  ′ Tpp = t pp ∑ sδ p†l+ε +δ ,σ pl+ε ,σ − t pp ∑ p†l−ε ,σ + p†l+3ε ,σ pl+ε ,σ  nl,σ nl+ε ,σ  = dl,†σ dl,σ = p†l+ε ,σ pl+ε ,σ .  Tpd is the direct oxygen-copper hopping. Tpp contains the direct t pp oxygen-oxygen ′ oxygen-oxygen next-nearest-neighbor nearest-neighbor hopping as well as the t pp  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, sδ = 1 (sδ = −1). ∆ pd is the charge transfer energy, which is the energy cost of the d 9 → d 10 L ¯ Ligand hole excitation. In reality, this energy decreases with the distance between the d 10 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 t pp < t pd < ∆ pd < Upp < Udd [5]. In practice, neighbor Coulomb interactions such as U pd d † d p† 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 solution 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  (a)  O  00 11 11 00 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 11 00 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11  000 111 111 000 000 111 000 111  ε  δ 000 111 111 000 000 111 000 111  Cu O  11 00 00 11 00 11 00 11 0000 1111 00 11 0000 00 1111 11 0000 00 1111 11 0000 1111 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11  11 00 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11  11 00 00 11 00 11 0000 1111 00 1111 11 0000 00 11 0000 00 1111 11 00 11 0000 1111 00 11 00 11 00 11 00 11  11 00 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11 00 11  (b) t pd  t pd  Cu O Cu  Cu O  t sw  O  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 ε vectors (solid arrow) and the four δ 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 either U and/or ∆, and that high-order terms in t pp always come with an even power of t pd , we derive an effective Hamiltonian for the states p†l+ε ,σ ∏ |σl ⟩ using degenerate RayleighSchrodinger perturbation theory[127] 1 − Pgs T + · · · ]Pgs E0 − H0 = ∏ ∑ dl,†σ dl,σ  Heff = Pgs [H + T Pgs  l  σ  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 54  (4.2)  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 t pd virtual processes shown in Fig 4.1b lead to an additional spin swap O-O hopping term: Tswap = −tsw ∑ sη p†l+ε +η ,σ pl+ε ,σ ′ |σl′ε ,η ⟩⟨σlε ,η | 2 /∆ where tsw = t pd pd and lε ,η = l + ε + (η · ε /|η · ε |)ε . This comes about when a  hole from a neighboring Cu site located at lε ,η first hops to one of its other 3 holefree 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 η = δ (Fig 4.1b), or NNN O-O hopping if η = ±2ε . In the latter case sη = 1. For η = δ , the sign sη 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 t pp or t pp  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 · Sl±ε 2 /[U + ∆ ] and S where Jpd = 2t pd pp pd l±ε , 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 ∑ Sl±2ε · Sl Πσ (1 − nl±ε ,σ ) 4 /[∆2 (U + 2∆ )]. where Jdd = 8t pd pp pd pd  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 2 “Kondo-hopping” term for we ignore them. We also discard a third-order t ppt pd  reasons discussed below. While we find that the three-spin polaron (3 SP) plays an important role, our approach is different from previous work [20–24] by recognizing (i) Tpp ’s signif55  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 d 9 Cu spin neighbors. This is the 3 SP √ |3SP⟩ =  1 (| ↑↓⟩ + | ↓↑⟩) |σ ⟩h − 6  √  2 |σ σ ⟩| − σ ⟩h , 3  (4.3)  first studied by Emery and Reiter [25]. 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 accounts for the lack of exchange between Cu spins bridged by an O hole. This lowers the 3 SP 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 t pp −tsw for a Cu-O triplet (with the “central” Cu at lε ,η ), but is enhanced to t pp +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 [16]. The 3 SP can also be written as a superposition of two Cu-O singlets: |3SP⟩ ∼ (| ↑, σ ⟩| ↓⟩h − | ↓, σ ⟩| ↑⟩h ) + (|σ , ↑⟩| ↓⟩h − |σ , ↓⟩| ↑⟩h ). Thus, the 3 SP 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 1 TKondo = −tk ∑(−1)nη p†l+ε +η ,α pl+ε ,α (Sl+ε · Slε ,η − ) 4 ) ( 2 2 t pd t ppt pd + tk = ∆(∆ +Upp ) (Udd − ∆ pd ) considered previously [19–24]. TKondo is non-zero only for a Cu-O singlet so it is 2 /U ) reinforcement to the above terms. merely a ∼ O(t pd dd 2 /(U − ∆ ) correcA finite Udd also modifies matrix elements by some t pd dd pd  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, especially against the coherence discussed above, due to the lack of the factor of 4 as in superexchange, and because their denominator involves Udd and ∆ +U pp .  4.3  Results  ′ = 0.58t , ∆ = 3.6eV , and U = 4eV [5], Using t pd = 1.3eV , t pp = 0.65eV , t pp pp pp pd  we scale the parameters in units of Jdd to find their dimensionless values to be t pp = 4.13, tsw = 2.98, and Jpd = 2.83. We push the computational limit to perform total-spin-resolved exact diagonalization (ED) of a topologically superior [129] 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 resolution. N = 16 results are provided to check for size dependence. All low-energy eigenstates have a total spin of either ST = jections for each ST are degenerate. The ST =  1 2  1 2  or 32 . The z pro-  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 = carrier can also mix with S = 1 or 2 background states to yield the ST = space. The partition of the  STz  3 2  1 2  sub-  subspace into separate ST sectors was managed by  the optimizations from Ch.3. Unlike there, no basis truncation was employed here  57  E(k)  a)  -42 -43 ST=1/2 32A ST=3/2 32A  -44  Z(k)  ST=1/2 16B  b)  0.3  ST=3/2 16B t-t’-t’’-J 32A  0.2 0.1 0 (0,1)  (0,0)  (1/2,1/2)  (1,1)  (0,1)  (1/2,1/2)  (1,0)  (kx,ky)/π Figure 4.2: a) Energy and b) quasiparticle weight (bottom) for the lowest eigenstates with ST = 12 and 32 vs. momentum. Different sets are shifted so as to have the same GS energy. for rigorous results. The (k, ST = 32 , STz = 21 ) sector contains ∼0.44×109 states.  4.3.1 Energetics Fig. 4.2(a) shows the lowest eigenenergies. The GS has k = ( π2 , π2 ) and ST = 12 , and is consistent with the 3-spin polaron (Eq. 4.3). Remarkably, we find similar dispersion along (0,0)→(π ,π ) and (0,π )→(π ,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 =  3 2  states which go below the  1 2  states near  (0,0) and (π ,π ). The solution’s robustness is supported by the difference E 3 -E 1 2  2  which decreases with increasing N at each k. We note that the magnon spectrum is gapped as ∼ 1/N for finite lattices [37]. 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 [120]. Given this finite-size scaling of the “free magnon” energy, we find that ST =  3 2  states at (0,0), ( π4 , π4 ), ( 34π , 34π ) and (π , π ) have much lower energy than minq [E 1 (k − q) + Ωq ] for a ST =  1 2  state plus a free magnon. At other k-points, the ST =  3 2  2  states  are within ∼ 0.2J of this value, so we cannot draw definite conclusions without access to larger N data. It follows that the ST =  3 2  states are stable polarons at least in  the regions marked by thick solid lines in Fig. 4.2(a). Thus, a ST =  1 2  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 ,σ |σ , l⟩x . The band-crossing results in noticeable change in the expectation values of the correlation function: Cˆx (δ , a) = 2 ∑ Sl+δ · Sl+δ +a nl+ex ,σ  (4.4)  l,σ  which measures the correlation between two Cu sites separated by a lattice constant at a certain distance away from the hole. ⟨Cˆy ⟩ is a reflection with Pˆx↔y for kx = ky . ⟨C⟩ ranges from − 43 for singlet, to ∼ −0.33 for 2D AFM GS, to 14 for triplet. Fig. 4.3a shows the ⟨Cˆx ⟩ correlations between neighbor Cu spins, when the hole is located at the darkly shaded bullet, in the GS: k = ( π2 , π2 ), ST = 12 . 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, of ∼0.13. Also, ⟨HJpd ⟩ ∼ -0.9Jpd , showing that locally this is consistent with the 3 SP 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´el background, where a spin-flip would change the spin-spin correlation to all four neighbors. Although the two central Cu spins have  2 3  weight in triplet configuration which is hardly bi-partite, long-range  59  −0.33  −0.26  −0.21  −0.30  −0.14  −0.23 −0.23  −0.28 −0.07  +0.17 −0.07  −0.26  −0.23  −0.28  −0.21 −0.07  −0.14 −0.07  −0.23  −0.22 −0.29 −0.05  +0.13 −0.22  −0.31  −0.29  (b) −0.21  −0.33 −0.31  −0.26 −0.05  −0.26  (a) −0.30  −0.21  Figure 4.3: ⟨Cx (δ , ∆)⟩ for the lowest energy state at (a) ( π2 , π2 ) with ST = 12 , and (b) at (π ,π ) with ST = 23 . The darkly-shaded bullet denotes the oxygen 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. ⟨Cy (δ , a)⟩ 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 = π2 , 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 =  3 2  polaron becomes the  lowest energy state at (π ,π ). The results look similar at (0,0). ⟨HJpd ⟩ 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  3 2  polaron is formed by a spin disturbance  around the 3 SP. This is very different from the S = energy  J + 2pd [29].  60  3 2  excitation local to HJpd with  Before we proceed further, let us note that the T=0 single-hole Green’s function ⟨AFM|peiHt p† |AFM⟩ contains information only of the states in the Krylov subspace which is the subspace spanned by the class of wavefunction H n p† |AFM⟩, for all positive integers n. It follows that the spectral weight of first electron removal Z(k), being the imaginary part of G(k, ω ), 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(π ,π )=0 because here the lowest eigenstate has ST = spin-conservation is not in the Krylov space of any ST =  1 2  3 2  which due to  state, and b) Z(0, π )=0  because the lowest eigenstate is not in the Krylov space H n p† |AFM⟩. The t-t’t”-J model restricts a spin degree of freedom, leading to Z(0, 0) ∼ 0.1 and a finite Z(π ,π ) [18] instead of zero as dictated by a  3 2  polaron. Unlike the t-t’-t”-J  model [18], Z(k) is concave down along (π /2,π /2)→(0,π ) in closer agreement with ARPES  which measured a maximum between ( π2 , π2 ) and ( π4 , 34π ) [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, π ) has ST = 21 , 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+εx |AFM⟩, and  Krylov subspace of ∑ eikl p†l+εy |AFM⟩. The k = (0, π ) 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.19 −0.09  −0.24  −0.28  −0.26  −0.28  −0.23  −0.09  −0.24  −0.28 −0.09  +0.15  −0.09  −0.24 −0.26  −0.24 y  −0.19  −0.23  −0.10 −0.19  −0.24  −0.29 −0.10  +0.16 −0.10  −0.23  −0.29  + p†l+ y ,σ |σ, l  −0.24  −0.19 −0.10  −0.23  −0.24  x  −0.28  Σ(−1)ly p†l+ x,σ |σ, l  −0.24  Figure 4.4: ⟨Cx/y (δ , a)⟩ for lowest state at k=(0,π ) 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 − E 1 = 0 band crossings in the nodal direction by looking at the k 2  2  points between which the difference switch signs. The observation in Fig. 4.2a is that, going away from the ( π2 , π2 ) GS, E 3 − E 1 is larger towards (0, 0) than towards (π , π ). The  1 2  →  3 2  2  2  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 − E 1 towards k = (0, 0) would suggest the non-zero region is larger towards 2  2  k = (0, 0) than towards k = (π , π ). Fig. 4.2a also shows that the  3 2  states get pushed  further down as N → ∞ so the crossing is expected to be closer to the ( π2 , π2 ) 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 = 12 band. This provides a  Jdd /2 energy scale for spin excitations. At finite T , as magnons  become thermally activated, these  3 2  states become “visible” to  gests a T -dependent broadening mechanism of  ARPES .  This sug-  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- 32 excitation of this energy scale (the computed energy gap is larger than the  1 N  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 [130]. In addition to the HJPD also has a S =  3 2 1 2  band, there are internal excitations of the local 3 SP since doublet and a S =  3 2  quartet separated in energy by Jpd  and 3Jpd /2 from the 3 SP. Magnetic excitations at these energy scales have been observed via inelastic resonant x-ray scattering for highly doped samples [131].  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 ;  lifting of the Cu-O singlet restriction present in  models leads to wave  ZRS -based  however, the  functions of a different nature, namely the 3 SP 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 (π , π ); 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[132]. The result contradicts other studies using variational or meanfield 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 = 64  1 2  and  ST =  3 2  polarons in Fig. 4.3. The two-hole solution for a larger cluster could be  different from previous studies limited to N=16[24]. I realized early that, even after exploiting translational and spin-projection symmetries, the two holes injected into the N=32 spin polaron model would have a Hilbert space with 0.154 × 1012 states, and larger N would involve literally an astronomical number of states. A system of this size challenges the capability of unbiased methods, and such an endeavor is extremely high-risk, especially considering 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 =  1 4  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 =  1 4  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 subspaces SA  ⊗  SB of sub-lattice total-spin SA/B . While the CS =  1 4  truncation captures  ∼ 95% of the undoped wavefunction as discussed in Ch. 3, Table 5.1 shows that CS =  1 4  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  ⊗  SB , and states added  by increasing CS have decreasing importance in the wavefunction; therefore, it is reasonable to expect an increment of ∼ dCS ∼  1 N  from CS =  1 4  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  -25  10  -30  0  -1  -35  ∆E/E  EGS/Jdd  10  10  -2  -40  10  -45  -3  -4  -50 0  10 0  1  0.5  CS  1  0.5  CS  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. SB SA \  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 = ( π2 , π2 ) single-hole ground state’s weight in ⊗ subspaces of particular SA 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 ∼  1 4  + 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+ε ,σ p†l ′ +ε ′ ,σ ′ ∏ |σl ′′ ⟩, ( ) N+2 with a degeneracy of 2N . The Rayleigh-Schrodinger expansion in Sec. 4.2.1 2 2 operates in the N+2 hole sector of the Fock space, and the outcome is slightly different from the single-hole scenario. In particular, the projector used in the expansion becomes P2h = ∏(1 − nl+ε ,↑ nl+ε ,↓ ) ∏(nl,↑ + nl,↓ − 2nl,↑ nl,↓ ),  (5.1)  which allows only states with a full lattice of copper spins and no double-occupancies (due to Upp/dd ). The Hamiltonian can be written as two parts He f f = P2h H1 P2h + 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  (2)  (3)  (4)  P2h H1 P2h from second- to forth-order with H2 = H2 + H2 + H2 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 t pd 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(2)  order H2  (3)  2 . All terms in the third-order H share a factor of t pp 2  2 t share a t pd pp  3 process would yield states factor because the second step of any three-step t pp (4)  projected out by (1 − P2h ) in the perturbation. The forth-order H2 can have terms 4 , t 2 t 2 and t 4 . with a prefactor of t pp pp pd pd  The second order corrections are the bare hole-hole correlations.  2 2t pp (2) H2 = − Upp  ∑(−1)  δ (δ1 ·δ2 )  (  1 S·S− 4  )α ′ β ′ αβ  p†l+ε +δ1 ,α ′ pl+ε +δ1 ,α p†l+ε +δ1 +δ2 ,β ′ pl+ε ,β (5.3)  Here, δ1,2 sum over (±ex , ±ey ). The matrix element is non-zero only if the two holes are |δ1 | =  √a 2  apart. The sign of the matrix element is positive when  δ1 · δ2 = 0. There are 8 non-zero matrix elements in this situation. Two of these 2 2t pp . There are also correspond to δ1 + δ2 = 0 so the static AFM exchange is 2 × Upp 6 other ways for one hole to “skip” over the other with such a Heisenberg factor. The factor of 41 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 8·16·t 4 t pp is roughly 3U 3 pp ∼ 16 81 t pp for t pd ∼ 2t pp and ∆ pd ∼ U pp ∼ 6 . The basis of this pp 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 t pp and Tswap , Jpd , Jpp ∼ 0.66t pp , which are already greater than the long-range physics at the scale of Jdd ∼ 0.2t pp . Adding the relevant short-range corrections with magnitude greater than 0.2t pp should be more than adequate to capture the physics. We also discards terms that can be factored into t pp , Tswap , Jpd , Jpp and the identity processes, which merely contribute some combinatorial constants scaling in the  68  infinite-order perturbation. 2 . The As discussed above, all third-order corrections have a prefactor of t ppt pd  perturbation goes through two virtual states so the largest possible magnitude is 2 t pp t pd ∆ pd U pp  ∼  ments is  2 t pp t pd 1 ∆ pd ∆ pd ∼ 9 t pp . The splitting due 2 9 t pp which is the same as the Jdd  to such a pair of Hermitian matrix elesplitting. Evidently, all these processes  can be factored into consecutive low-order processes and are discarded. Because Upp ∼ ∆ pd , 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 2 processes requires the same orbital occupancy in the initial and final among t ppt pd  state, and this can happen only if the two oxygen holes are δ = (±ex , ±ey ) apart. The transition of interest is then p†l±ex ,α p†l±ey ,β dl,†γ → p†l±ex ,α ′ p†l±ey ,β ′ dl,†γ ′ . The effective matrix elements have the form ⟨α ′ , β ′ , γ ′ |H2 |α , β , γ ⟩, a local three-spin (3)  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 ∆x/y summed over ±εx/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  (3)  H2 =  (  2td2pt pp  Upp (∆ pd +Upp ) ∑  S·S− ( S·S−  + (  S·S−  + ( + +  2td2pt pp  (∆ pd +Upp )2 ∑  S·S−  1 4 1 4 1 4 1 4  )β ′ α ′ βα  )α ′ β ′ αβ  )γ ′ α ′ γα  )β ′ γ ′ βγ  p†l+∆x ,β ′ pl+∆x ,α p†l+∆y ,γ pl+∆y ,β |l, α ′ ⟩⟨l, γ | p†l+∆x ,γ pl+∆x ,α p†l+∆y ,α ′ pl+∆y ,β |l, β ′ ⟩⟨l, γ | p†l+∆x ,γ ′ pl+∆x ,α p†l+∆y ,α ′ pl+∆y ,β |l, β ⟩⟨l, γ | p†l+∆x ,β ′ pl+∆x ,α p†l+∆y ,γ ′ pl+∆y ,β |l, α ⟩⟨l, γ |  ( ) ′ ′ 1 βγ † δαβ S · S − p pl+∆x ,α p†l+∆y ,β ′ pl+∆y ,β |l, γ ′ ⟩⟨l, γ | 4 β γ l+∆x ,α ) ′ ′ ( 1 αβ † + δγβ S · S − p pl+∆x ,α p†l+∆y ,α ′ pl+∆y ,β |l, β ′ ⟩⟨l, γ | 4 αβ l+∆x ,γ ) ′ ′ ( 1 βα † † ′ + δγα S · S − p ′ pl+∆x ,α pl+∆ ,γ pl+∆y ,β |l, α ⟩⟨l, γ | y 4 β α l+∆x ,β ( ) ′ ′ 1 αγ † † ′ + δβ α S · S − p (5.4) γ| ′ pl+∆x ,α pl+∆ ,β pl+∆y ,β |l, γ ⟩⟨l, y 4 αγ l+∆x ,α  These terms are finite if there is a 2-spin singlet formation among the 3 spins. t 2 t pp  p The first four terms give an eigenvalue of +4 Upp (∆d pd +U pp ) ∼ 0.2t pp for oxygen-  oxygen singlet pair and 0 for triplet pair. The last four terms give an eigenvalue of td2p t pp (∆ pd +U pp )2  ∼ 0.02t pp for singlet pairs and 0, −3 (∆  td2p t pp 2 pd +U pp )  ∼ −0.06t pp for triplets.  Therefore we can approximate this term by raising the energy of oxygen-oxygen singlet pair appropriately. 4 processes can be factored into two T In the forth-order correction, all t pd swap or 2 t 2 processes are smaller than J by a factor of ∼ 4 × 4 and Jpd processes. The t pd dd pp 4 processes are even smaller. Therefore, no constructive interference is possible. 2t pp (4)  we set H2 ∼ 0. In summary, the 2-hole correction is in essence  70  ( ) ′ ′ 1 αβ † δ (δ1 ·δ2 ) (−1) S · S − p p† ′p ′p ∑ 4 αβ l+ε +δ1 ,α l+ε +δ1 ,α l+ε +δ1 +δ2 ,β l+ε ,β ) ′ ′ ( 1 αβ † (3) − Jpp ∑ S · S − p p† (5.5) ′p ′p 4 αβ l+ε +δ ,α l+ε +δ ,α l+ε ,β l+ε ,β (2)  H2 ∼ + Jpp  (2)  with Jpp =  2 2t pp U pp  (3)  4t 2 t pp  and Jpp = Upp (∆dpdp +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 ( ) N+2 Eq. 5.5) contains 2N states of the form p†l+ε ,σ p†l ′ +ε ′ ,σ ′ ∏ |σl ′′ ⟩, with l + 2 2  ε ̸= l ′ + ε ′ . To perform computations for large N, this basis must be reformulated to take advantage of translational symmetry, total-spin symmetry, total-spinprojection 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 required to implement the basis and the truncation is a convoluted process. The code written for the project totaled ∼0.5MB in size, containing ∼ 65536 characters. 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 compatible 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,l ′ † t−1,l,l ′ † t0,l,l ′ † t1,l,l ′  1 √ (p†l↑ p†l ′ ↓ − p†l↓ p†l ′ ↑ ) 2 † † = pl↓ pl ′ ↓ =  1 √ (p†l↑ p†l ′ ↓ + p†l↓ p†l ′ ↑ ) 2 = p†l↑ p†l ′ ↑ . =  (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 STz = 0, any spin background can be specified orthonormally with respect to a position l. { † l,l ′ |α ⟩l  ≡  s†l,l ′ |α , Sz = 0⟩l  ∑1z=−1 c(z, ST,N , ST )tz,l,l ′ |α , Sz = −z⟩l †  (5.7)  α denotes a particular group of 2ST,N + 1 spin configurations related by the total +/− +/− spin raising and lowering operators ST = ∑l Sl , summed over all Cu sites. The total spin of this α 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 STz = 0, z the two-hole singlet would mix only with backgrounds with ST,N = 0. The twoz hole triplets would mix with three different projections |α , S = −1⟩l , |α , Sz = 0⟩l , and |α , Sz = 1⟩l from the group α . 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 STz = 0. Exploitation of the translational symmetry is performed by the use of a Fourier series of the form ∼ ∑ eiKl †l+ε ,l ′ +ε ′ |α ⟩l ; 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 ε = ε ′ = εx/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 N for the 2D, N-unit-cell lattice. Defining  δ l = (l0′ − l0 , l1′ − l1 ),  (5.8)  most hole-hole configurations can be classified in the region L0 2 L0 0 < δ l0 < 2 L0 0 < δ l0 < 2 L0 δ l0 = 2 0 ≤ δ l0 <  , 0 < δ l1 < , δ l1 =  L1 2  L1 2  L1 < δ l1 ≤ 0 2 L1 , − < δ l1 < 0. 2 , −  (5.9)  These states can be expressed using N-term Fourier series |  xx/yy , δ l, α , K⟩  1 = √ ∑ eiKl N  † l+εx/y ,l+δ l+εx/y |α ⟩l  (5.10)  Because the two oxygen holes are indistinguishable fermions (Eq. 5.6), there are three remaining δ l values which require special attention due to the periodic boundary condition.  δ l0 =  L0 2  , δ l1 = 0 L1 2 L1 , δ l1 = 2  δ l0 = 0 , δ l1 = δ l0 =  L0 2  (5.11)  The number of terms in the Fourier series depends on the spin background translated by δ l: Tδ l |α ⟩l . If such a translation yields an orthogonal state, l ⟨α |Tδ l |α ⟩l =  73  0, the Fourier series still has N terms. For the example of δ l = ( L20 , 0), | =  xx/yy , δ l, α , K⟩  1 L1 −1 √ ∑ eiK1 l1 N l1 =0  1 = √ ∑ eiKl N  L0 2 −1  ∑  eiK0 l0  l0 =0  † l+εx/y ,l+δ l+εx/y |α ⟩l  † l+εx/y ,l+δ l+εx/y  (  ) L0 1 + s eiK0 2 T L0 |α ⟩l (5.12) 2  where s is the sign change due to hole-swapping in the singlet or triplet (  † a,b  =  For the case of ⟨α |Tδ l |α ⟩l = ±1, the above expansion makes clear that ) ( L0 there can be only N2 terms in the Fourier series due to the term 1 + s eiK0 2 T L0 . s  † b,a ).  2  For this case the series has the form √ |  xx/yy , δ l, α , K⟩  =  2 L1 −1 iK1 l1 e N l1∑ =0  L0 2 −1  ∑  eiK0 l0  l0 =0  † l+εx/y ,l+δ l+εx/y |α ⟩l  (5.13)  The formulation for the scenario in which one hole occupies a pl+εx and the other occupies pl+εy is straight forward. All values of δ l are unique and the Fourier series has the form |  xy , δ l, α , K⟩  1 = √ ∑ eiKl N  † l+εx ,l+δ l+εy |α ⟩l .  (5.14)  In summary, a full orthonormal Hilbert space can be specified by the states |  xy , δ l, α , K⟩  and |  xx/yy , δ l, α , K⟩.  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 background α .  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 ( , α ) (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 δ l’s for † l+εx ,l+δ l+εx ,  of states |  and 17 distinct δ l’s for  xy , δ l, α , K⟩  and |  † l+εx ,l+δ l+εy ,  † l+εy ,l+δ l+εy .  xx/yy , δ l, α , K⟩  17 distinct δ l’s for  So n2h=66. The total number  is thus n2h×nSpinRow.  Because the matrix elements of the Cu-Cu superexchange ∑<l,l ′ > Sl · Sl ′ 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. (n2h−1) (nSpinRow−1)  |Ψ⟩ =  ∑  ∑  i=0  cnSpinRow×i+ j |i, j⟩  (5.15)  j=0  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 δ l. Each |i, j⟩ corresponds to a Fourier series of a definite set of ( , δ l, K, ST ) as discussed in the previous section. This is an orthonormal numerical basis ⟨i, j|i′ , j′ ⟩ = δi,i′ δ j, j′  (5.16)  nSpinRow×n2h−1  ∑  |cm |2 = 1  (5.17)  m=0  Of course, the values of vector elements are the probability amplitudes of a quantum mechanical description. For example, the probability of the two oxygen holes in one of the n2h configurations is nSpinRow−1  P(i) =  ∑  |cnSpinRow×i+ j |2  (5.18)  j  5.4  Low-Energy Two-Hole Solutions  Exact diagonalization in the subspace of CS = 0 to  1 2  was performed for the two-  hole spin polaron model (Eq. 5.2 with Eq. 4.2 and Eq. 5.5), with parameters t pp = ′ ′ = 0, t = 2.9822, J = 2.8253, J 4.1293, t pp pp = 1.3420, and J pp = 0.9182. t pp sw pd (2)  (3)  was set to zero to simplify the model. Another reason is that the Jpp corrections 75  -64 ST=0, CS=1/4 ST=1, CS=1/4  E(K)GS/Jdd  -65  -66  -67 (0,1)  (0,0)  (1/2,1/2)  (1,1)  (0,1)  (1/2,1/2)  (1,0)  (Kx,Ky)/π  Figure 5.2: Lowest two-hole states in the CS =  1 4  subspace.  ′ would hop a hole half-way across the N = 32 cluster, worsening finite due to t pp ′ = 0 for comparison. size effects. The one-hole GS was recalculated with t pp  5.4.1 Energetics Fig. 5.2 shows the dispersion of the lowest two-hole states for different totalmomentum K and total-spin ST . The K = (π , π ) ST = 0 groundstate is doubly degenerate. The K = (π , 0) and K = (0, π ) 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  1 2  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 single 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 twohole energy minus twice the one-hole energy, shifted by the energy of the undoped  76  -55  10  0  -1  10  E/Jdd  ∆E/E  -60  10  K=(π,π) ST=0 K=(π,0) ST=0 K=(0,0) ST=0 ∆Eb=0  -2  -65  10  -3  -70 -4  0  1  0.5  10 0  CS  0.5  1  CS  Figure 5.3: Convergence of the ST = 0 lowest energy states. The two-hole energy level corresponding to ∆Eb = 0 is indicated by the horizontal dash line. system. ∆Eb = (E2h − E0h ) − 2(E1h − E0h )  (5.19)  The energy of the K = (π , π ) groundstate was extrapolated to be −69.167Jdd , giving ∆Eb = +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  ∆E E  in Fig. 5.3, the extrapolated groundstate energy lowers to  −70.1228Jdd . This lowest estimate gives a binding of ∆Eb = −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  1 N  ∼ 0.3Jdd q=0 magnon excitation, so a rigorous  conclusion cannot be drawn because |∆Eb | is smaller than this finite-size effect.  77  5.4.2 Wavefunction Analysis The K = (π , π ) GS was found to be degenerate. It is well known that Lanczostype 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 |GS0⟩, then extract the GS of (I − |GS0⟩⟨GS0|)H(I − |GS0⟩⟨GS0|). 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 |GS1⟩, 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 2 × 2 rotation (  0  1  )   · · ·    .. . 0  .. . 0    −1 0    0  .. .  0  (5.20)  The resulting pair of groundstate wavefunctions satisfies ⟨GS0|GS0⟩ = ⟨GS1|GS1⟩ = 1 ⟨GS0|GS1⟩ = 0 ⟨GS0|Pxy |GS0⟩ = −⟨GS1|Pxy |GS1⟩ = cos2 θ − sin2 θ ( ) ( )( ) |GS+⟩ cos θ sin θ |GS0⟩ = |GS−⟩ − sin θ cos θ |GS1⟩ (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 = (π , π ) is a high  78  symmetry point invariant to Pxy and Pxy . The eigenvalues are ⟨GS ± |Pxy |GS±⟩ = −⟨GS ± |Pxy |GS±⟩ = ±1  (5.22)  which indicate p-symmetry. The K = (π , 0) and K = (0, π ) 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 = (π , π ) 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 cˆcharge (δ ) = ∑ p†l+ε ,σ pl+ε ,σ p†l+ε +δ ,σ ′ pl+ε +δ ,σ ′  ∑  Cˆcharge (R) =  cˆcharge (δ )  (5.23)  (5.24)  |δ |=R  1 =  ∑ Cˆcharge (R)  (5.25)  R  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 ∑R P(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 ∆Cˆcharge (R) = Cˆcharge (R) − P(R)  (5.26)  where the correlation is the same as that of a random distribution if ∆C(R) = 0. The expectation values as a function of hole-hole separation are shown in  79  Ccharge(R)  0.2  0.1  K=(π,π) ST=0 K=(π,0) ST=0 K=(0,0) ST=0 Random  0 0  2  1  3  4  3  4  R/a  ∆Ccharge(R)  0.05  0  -0.05 0  2  1  R/a  Figure 5.4: (top) Probability of charge separation, ⟨Ccharge (R)⟩ for the lowest states in the CS = 21 subspace. (bottom) The difference from a random distribution. Fig. 5.4. ∆Ccharge (R) indicates charge repulsion at K = (π , π ), a small peak of correlation at R=2 for (π , 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=  1 2  3 SP, which is a composite of the oxygen spin and two neighboring copper  spins. The two-hole solutions yielded ⟨HJpd ⟩ ∼ 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 ⟨∑ dl+ σ dl+σ ⟩ = 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 √ √ † 2 √ | ⇑⟩ = 13 |↑↓⟩+|↓↑⟩ p − | ↑↑⟩p†↓ ↑ 2 √3 √ √ p†↓ − 23 | ↓↓⟩p†↑ | ⇓⟩ = 13 |↑↓⟩+|↓↑⟩ 2 √ √ p†↑ |0+⟩ = 13 |↑↓⟩−|↓↑⟩ 2 √ √ p†↓ |0−⟩ = 13 |↑↓⟩−|↓↑⟩ 2  Total Spin  = | ↑↑⟩p†↑ √ √ † † 1 √ = 23 |↑↓⟩+|↓↑⟩ p + ↑ 3 | ↑↑⟩p↓ 2 √ √ √ | 32 , − 21 ⟩ = 23 |↑↓⟩+|↓↑⟩ p†↓ + 13 | ↓↓⟩p†↑ 2 | 32 , − 23 ⟩ = | ↓↓⟩p†↓  | 32 , | 32 ,  3 2⟩ 1 2⟩  ⟨HJ pd ⟩ J pd  1 2  −1  1 2  −1  1 2  0  1 2 3 2 3 2  0  3 2 3 2  1 2 1 2  1 2 1 2  Table 5.2: Single-hole eigenstates of HJpd . p†σ creates an oxygen hole and the arrows in the ket indicate the spins of the two copper sites neighboring the oxygen hole. polarons involving a total of six spins. Noting that ⟨HJpd ⟩ ∼ −1.8Jpd , the singlepolaron levels in Tab. 5.2 suggest that the dominant part of the wavefunction con(| ⇑⟩| ⇓⟩ ± | ⇑⟩| ⇓⟩),| ⇑⟩| ⇑⟩, { } ⟨HJ ⟩ and | ⇓⟩| ⇓⟩. Other pair configurations would have Jpdpd ∈ 0, ± 21 , ±1 . The  tains four possible 3 SP pairs with ⟨HJpd ⟩ = −2Jpd :  √1 2  numerical wavefunction can be projected into the subspace of 3 SP pairs by P3SP =  ⟨HJpd ⟩/Jpd − ∆λ −2 − ∆λ ∆λ ∈{0,± 1 ,±1}  ∏  (5.27)  2  wherein each term in the product zeros out the eigenvalue of a particular level while scaling the  ⟨HJ pd ⟩ J pd  = −2 level to unity. The probability of having two 3 SP with any  spin-spin correlation and at any separation is thus ⟨P3SP ∑δ cˆcharge (δ )P3SP ⟩. In a six-spin problem, the projected wavefunction is a superposition of the four possible 3 SP pairs |6spin⟩ = a  | ⇑⟩| ⇓⟩ − | ⇑⟩| ⇓⟩ | ⇑⟩| ⇓⟩ + | ⇑⟩| ⇓⟩ √ √ +b| ⇑⟩| ⇑⟩+c +d| ⇓⟩| ⇓⟩. (5.28) 2 2  The probability of finding a singlet among oxygen holes in the left and right ket,  81  p†L,σ and p†R,σ respectively, is ⟨  p†L↑ p†R↓ − p†L↓ p†R↑ pL↑ pR↓ − pL↓ pR↑ √ √ 2 2  ⟩ 3 2 = a2 + (b2 + c2 + d 2 ) 9 9  (5.29)  Rearranging for a2 and generalizing the expression to different hole-hole separation  δ , the nature of singlet 3 SP pair can be gauged by c3SP,singlet (δ ) = cˆsinglet (δ ) =  ( ) ⟨P3SP 9cˆsinglet (δ ) − 2cˆcharge (δ ) P3SP ⟩ ⟨P3SP ∑δ ′ cˆcharge (δ ′ )P3SP ⟩  ∑ s†l+ε ,l+ε +δ sl+ε ,l+ε +δ  (5.30) (5.31)  where s† is the singlet creation operator defined in Eq. 5.6. c3SP,singlet (δ ) is the probability, within the projected subspace, of the two oxygen holes and the  AFM  background forming a 3 SP singlet pair separated at distance δ . The value ranges from zero for no singlet nature to unity for pure singlet at hole-hole separation of δ ; however, these two extreme values are not possible due to, for example, the spatial spreading required to lower the kinetic energy of Tpp and Tswap . The measure can be summed as a function of hole-hole separation C3SP,singlet (R) =  ∑  c3SP,singlet (δ )  (5.32)  |δ |=R  The corresponding triplet measure is C3SP,triplet (R) = 1 =  ⟨P3SP cˆcharge (δ )P3SP ⟩ − c3SP,singlet (δ ) (5.33) ⟨P ˆ (δ ′ )P3SP ⟩ ′c |δ |=R 3SP ∑δ charge  ∑  ∑ C3SP,singlet (R) +C3SP,triplet (R)  (5.34)  R  The singlet correlation should be compared to ∑R P′ (R) = 1, the random distribution of two 3 SPs spreading over the 64 oxygen sites, with no common copper spin. The probability of singlet correlation is  1 4  in a paramagnetic state. The difference  of interest is thus ∆C3SP,singlet (R) = C3SP,singlet (R) − 82  P′ (R) . 4  (5.35)  C3SP,singlet(R)  0.1  K=(π,π) ST=0 K=(π,0) ST=0 K=(0,0) ST=0 Random  0.05  0 0  2  1  3  4  3  4  ∆C3SP,singlet(R)  R/a  0.05  0 0  2  1  R/a  Figure 5.5: (top) Probability of singlet 3 SP pair, ⟨C3SP,singlet (R)⟩ for the lowest states in the CS = 12 subspace. (bottom) The difference from a randomly distributed paramagnetic configuration. The low-energy states all have ⟨P3SP ∑δ cˆcharge (δ )P3SP ⟩ ∼ 0.8. The rest of the weight are distributed towards the sector with one or two |0±⟩ (Tab. 5.2). Figure 5.5 shows the resulting C3SP,singlet (R) and ∆C3SP,singlet (R). It is evident that all states have enhanced short-range singlet nature compared to the random distribution – a feature absent from previous studies using small cluster or the t-J model which does not distinguish oxygen and copper sites. This short-range nature is robust against finite-size issues because long-range nature can indeed be found as for the charge correlation Fig. 5.4. The K = (π , π ) GS has enhanced short-range singlet correlation but the overall charge distribution is repulsive as shown in Fig. 5.4. This implies that the most of the long-range charge correlations belongs to the sector of non-3 SP-pair nature. The individual c3SP,singlet (δ ) values are shown in Fig. 5.6. The local maximum occurs when the two 3 SP are R/a ∼ 1.5 − 2 apart. I stress here that, although the fluctuation of the values over space is merely about a factor of five, such fluctuation is not present in the t-J model solved for two holes propagating in a cluster with the same geometry[132]. This enhanced correlation in the singlet nature must be an anomalous signature because, for example, the singlet enhancement for R/a ≤ 2 83  18 57  63  20  108 9  63  18  92  107 9 63  128  128  92 20  24 9 107 20  63  12 57 9 108  12 63 20  24  20  8  67  97  32  13  13  27  15  68  68  47  26  13  107  20 24  16  13  33  109  8  109 23  98  46  24  110  64  23  67  Figure 5.6: c3SP,singlet (δ ) × 104 for one of the two K = (π , π ) GS in the CS = 12 subspace. Cu and O are marked by circles and rectangles, 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 3 SP 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 3 SPs both centered on x-rung O sites. 84  is observed even though the charge measure shows no such tendency (Fig. 5.4). In fact, a singlet correlation that is not confined in real space implies that the singlet correlation should be more confined in momentum space; however, k-space features have not been extracted because each 3 SP involves two background copper spin in real space, and k-space features would require Fourier transforming the local spins while leaving the rest intact. The K = (π , 0) state is slightly higher in energy. Figure 5.5 shows enhanced short-range singlet correlation. However, the charge correlation in Fig. 5.4 shows a local maximum at R/a = 2 unlike the K = (π , π ) state. The spatial distribution of the singlet and charge are shown in Fig. 5.7 and Fig. 5.8. The observation is that the charge and singlet maxima happens when both holes are aligned along the momentum direction. This alignment can potentially be pinned by external interactions as stripes, but a much larger system with more holes are required to verify this. Disregarding the S=1 excitation from the above states, the K = (0, 0), ST = 0 state is much higher in energy and has d-symmetry, unlike the s-symmetric state found in N=32 exact diagonalization of the t-t’-J model[132]. The singlet correlation is shown in Fig. 5.9. The two polarons are R/a = 1 apart without sharing a copper site.  5.5  Speculations  The results presented above are fairly unbiased interpretations of the numerical solutions. This section makes some speculations about larger systems with more holes. The goal is to list several potential features of interest for future developers. Although we have broken some technological barriers in extracting explicit multi-hole wavefunctions for large lattice at zero temperature, finite system modeling is inherently limited by cluster size. Two holes in a periodic N=32 system is often interpreted as an approximation to 6.25% doping, but the solution might or might not be representative of a real system because inhomogeneous concentration is not allowed. To be precise, if the holes have tendency to cluster, two holes in a N=32 system should be interpreted as limN→∞ N2 doping because a nominal 32unit-cell patch in a real 6.25% doping system would often not contain two holes.  85  54 54  66 43  11  43  67  45  43  51 51 35  66  51 51 35  35  45  67  67 35  45 43 66 66  45 54 54  67  8  15  9  26  48  26  14  83  103  11  104  23  104 15  17  25  12  265  49  24 183  26 17  144  8  26  25  12  47  20  24  Figure 5.7: c3SP,singlet (δ ) × 104 for the lowest K = (π , 0) state in the CS = 12 subspace. Cu and O are marked by circles and rectangles, respectively. Values for the K = (0, π ) 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 3 SP 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 3 SPs both centered on x-rung O sites. 86  222  222  222  154 230  240  184 185 144  58  72 108 141  72  72 108 141  72  108  143  154  141  141 108  154 144 185 184 143  154 230 230  222  156  165  35  156  135  153  76  225  137  139  318  58  173  99  81  172 153  289 156  146  141 87  171  187  134 265  165 98  76  224  76  138  141  Figure 5.8: ccharge (δ ) × 104 for the lowest K = (π , 0) state in the CS = 21 subspace. Cu and O are marked by circles and rectangles, respectively. Values for the K = (0, π ) 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 3 SP pair configurations are shown. For example, the bottom figure shows that ccharge (2a, 0) = 0.0289 for two 3 SPs both centered on x-rung O sites. 87  13  13  25  73 73  25  22  13  43  22  25  22  73 25  43  73  73 25  73  25  25  43 22 73 73  43 25 25  13  41  17 46  31  47  31  25  41  37  47  14  25  47  46  47 31  37  31  58 58  367  14  41  47 47  367  58  38  58  31  25  38  29 41  Figure 5.9: c3SP,singlet (δ ) × 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 3 SP 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 (0, −2a) = 0.0367 for two 3 SPs both centered on x-rung O sites. 88  Although in the last section we noted that the enhanced singlet correlations are not confined in real space and should therefore be confined in momentum space, one can speculate on the origin of the singlet tendency by considering the 6-spin wavefunctions  √1 2  (| ⇑⟩| ⇓⟩ ± | ⇓⟩| ⇑⟩). By expanding and taking the S · S  between a copper spin from the first 3 SP and one from the second 3 SP, the ex5 1 pectation value is − 18 Jdd if they are singlet and + 18 Jdd if they are triplet. There-  fore, two 3 SPs can take advantage of this nature by forming a singlet if they are separated by an empty oxygen site, which mediate an ordinary S · S Heisenberg AFM bond. The two-hole numerical solution indeed shows a local maximum when R/a ≤ 2 (Fig. 5.5). As doping increases, the average spacing between oxygen holes would decrease if they do not cluster. This can be checked by expanding the 12-spin wavefunctions obtained as the direct product of two 3 SP singlet pairs:  1 2  (| ⇑⟩| ⇓⟩ − | ⇓⟩| ⇑⟩)  ⊗  (| ⇑⟩| ⇓⟩ − | ⇓⟩| ⇑⟩). The expectation value of S · S  between a copper spin from the first pair and one from the second pair is exactly zero. Therefore, different pairs tend to be separated by more than one empty oxygen sites for some negative correlations. There is no immediate clustering tendency beyond two holes, at least according to this back-of-the-envelope calculation. This provides some merit for interpreting the finite-size results as a realistic doped scenario with little clustering. Having noted the lack of immediate clustering mechanism, I would like to point out that three features of the numerical solution are likely to prevail in larger systems with more holes: i) the kinetic and potential energy leads to the formation of 3 SPs, ii) 3 SP destroys AFM correlation in its vicinity (Fig. 4.3 and 4.4), and iii) two 3 SPs tend to have enhanced singlet correlation compared to the pure random scenario (Fig. 5.5). Noting these three points, the above back-of-the-envelope calculations sug5 gests a possible binding mechanism for a real doped system. The − 18 Jdd energy  for a 3 SP singlet pair separated by one empty oxygen site is of course higher than the ∼ −0.33Jdd copper-copper correlation of an undoped system; however, the 5 − 18 Jdd is a local AFM correlation which is independent of doping, and the cor-  relations between copper spins away from the 3 SP are raised on average due to point ii). If pair-pair clustering is indeed absent, doping would reduce the average distance between oxygen holes while increasing the average copper-copper corre89  5 lation. At certain doping, the − 18 Jdd short-range AFM correlation can become  energetically important and the holes form bipolaron pairs. This pairing tendency cannot increase indefinitely because the pairs want to be more than one empty oxygen apart. Evidently, ∼ 25% is the highest doping level which allows 3 SP singlet pairs without having pairs within one empty oxygen site. I stress again that this section is merely providing a list of speculations based on the two-hole solution. The two-hole solution has not been obtained when previous speculations of the same type were made. In particular, the numerical two-hole wavefunctions are not node-less s-symmetric as suggested by earlier approximations by Emery and Reiter[14, 25]. The temperature scale of the mechanism is likely to be much lower than Jdd because the binding tendency is determined by 5 the difference between the − 18 Jdd short-range AFM 3 SP-3 SP correlation and the  other copper-copper spin correlation, which increases from the undoped value of ∼ −0.33Jdd . The stability of the 3SP itself increases and also the pairing tendency increases as long range antiferromagnetic correlations are destroyed by further doping. The temperature cannot increase indefinitely because overpopulation would destroy the tendency of bipolaron singlet formation. If the pairing physics were indeed due to this singlet binding mechanism, the transition from the superconducting regime to the Fermi-liquid regime at ∼ 25% doping (Fig. 1.1) can be explained as the spatial saturation of 3 SP singlet pairs. The copper spin background can no longer be correlated into 3 SP with the oxygen holes, and the Fermi-liquid behavior is due to the “free” oxygen holes. This speculative picture should be explored when solutions of larger systems become available.  5.6  Conclusions  The octapartite scheme has been successfully applied to the scenario of two holes propagating in the spin-polaron model with 32 Cu and 64 O sites, breaking the technological barrier. Exponential convergence similar to the undoped and onehole scenarios has been established. While the energetics of the low energy states yield inconclusive results, one should be reminded about the danger of relying solely on ground-state energy, which is the “‘Quantum Chemist’ Fallacy No. 2” coined by by Anderson[133].  90  The perturbative derivation of the spin polaron model in Sec. 4.2.1 is bound to leave out terms that would slightly affect quantities such as binding energy. In fact, the pseudo gap temperature of T ∗ ∼ 150K is less than one-tenth of Jdd , smaller than the ∝  1 N  ∼ 0.3Jdd q = 0 magnon excitation for an N = 32 undoped AFM as  well as any energy difference in the one- or two-hole energy spectrum (Fig. 4.2 and Fig. 5.2). It would be fruitful to refine the octapartite scheme to treat larger lattices. Careful analysis of the wavefunctions at K = (0, 0), (π , 0), (π , π ) shows that all states exhibit anomalous enhanced hole-hole singlet correlation, compared to the ED results of the t-t’-J model solved on a N=32 cluster with the same geometry[132]. Some interesting speculations have been made about the behavior of larger systems, providing the merits for further developments.  91  Chapter 6  Conclusions 6.1  Summary  This work describes progress in the treatment of electron-electron and electronphonon interactions. Several numerical approaches were devised, developed, and adapted to accurately solve large models with details. These technological developments were complemented by analytical approaches to reveal new physics, stressing the importance of distinguishing anions and cations in transition-metal oxides. Chapter 2 explores the breathing-mode model which describes a charge carrier centered at cation sites interacting with quantized modes of vibration at anion sites. In this more realistic scenario, the coupling amplitude becomes momentumdependent unlike the widely-studied Holstein model. This momentum dependence directly leads to non-monotonic energy dispersion, sharper large-to-small polaron crossover as a function of coupling strength, and lighter effective mass in the lowcoupling regime compared to the Holstein model. The result serves as the first available accurate benchmark for the model for all coupling regime. Chapter 3 points out a low-energy characteristic of the much-studied spin- 12 Heisenberg antiferromagnet on a 2D square lattice, which is relevant to the undoped cuprate compounds. The octapartite approach is then developed to take advantage of this intuition. The groundstate wave function can be solved explicitly for a record breaking 64-spin torus. The error can be systematically reduced, and 92  the convergence of the approach is established to be exponential with respect to the completeness parameter CS ∈ [0, 1]. A spin polaron model is derived in Chapter 4 to capture the coherence of holes doped into the copper-oxygen layer in cuprate compounds. The solution of a single hole is solved via a novel total-spin-resolved exact diagonalization approach for a record-breaking cluster of 32 copper and 64 oxygen sites. The combination of the new model and large cluster reveals important physics, missed by previous studies, with direct experimental consequences. Chapter 5 combines Ch. 3 and 4 to bypass several computational restrictions, such as the sign problem and special periodic boundary condition, in dealing with two holes propagating in a cluster with 32 copper and 64 oxygen sites. The results prove that the octapartite scheme can indeed accurately capture the disturbed spin background. The solutions exhibit enhanced singlet correlation compared to the simplest t-J model solved in the same cluster.  6.2  Opportunities for Further Developments  The promising results in the two-hole scenario provide merits for further developments. Theoretically, the model can be tweaked by adding other terms, including coupling to phonons. Larger clusters are within reach because of various ways of improving the numerical scheme. Technologically speaking, specialized optimization can be worthwhile because the octapartite approach has proven to be reliable. In particular, the current implementation requires a pre-transformed matrix for computation. It would be fruitful to develop an algorithm fast enough to by-pass this step. Implementing this concept via variational Monte-Carlo sampling or renormalization can significantly extend its prowess. In a more general context, the octapartite approach (Ch. 3) is devised to capture an AFM background on the square lattice complicated by the presence of doped holes. Orthogonal extensions to the square lattice  AFM  problem include the frus-  trations due to non-nearest-neighbor AFM bonds and due to other lattice structures. Because the current approach is formulated to capture the essence of doped cuprates, the method is not expected to be directly applicable to these frustrated scenarios. Nevertheless, the use of nonlocal spin partitions to capture important  93  physical essences is a holistic way which could lead to powerful case-by-case solutions in dealing with these peculiar spin systems.  94  Bibliography [1] Ashcroft, N. W. and Mermin, N. D. Solid State Physics. Holt, Rinehart and Winston, (1976). [2] Taylor, P. L. and Heinonen, O. A Quantum Approach to Condensed Matter Physics. Cambridge University Press, Cambridge, (2002). [3] Maple, M. Journal of Magnetism and Magnetic Materials 177-181(Part 1), 18 – 30 (1998). International Conference on Magnetism. [4] Imada, M., Fujimori, A., and Tokura, Y. Rev. Mod. Phys. 70(4), 1039–1263 Oct (1998). [5] Ogata, M. and Fukuyama, H. Rep. Prog. Phys. 71, 03650 (2008). [6] Lee, P., Nagaosa, N., and Wen, X. Rev. Mod. Phys. 78, 17 (2006). [7] Lee, P. Rep. Prog. Phys. 71, 012501 (2008). [8] Salamon, M. B. and Jaime, M. Rev. Mod. Phys. 73, 583 (2001). [9] Bednorz, J. G. and Mller, K. A. Zeitschrift fr Physik B Condensed Matter 64, 189–193 (1986). 10.1007/BF01303701. [10] Vojta, M. Adv. Phys. 58, 699 (2009). [11] Nature Physics 2, 138 (2006). [12] Bariˇsi´c, N., Li, Y., Zhao, X., Cho, Y.-C., Chabot-Couture, G., Yu, G., and Greven, M. Phys. Rev. B 78(5), 054518 Aug (2008). [13] Zaanen, J., Sawatzky, G. A., and Allen, J. W. Phys. Rev. Lett. 55(4), 418–421 Jul (1985). [14] Emery, V. J. Phys. Rev. Lett. 58(26), 2794–2797 Jun (1987). 95  [15] Eskes, H. and Sawatzky, G. A. Phys. Rev. Lett. 61(12), 1415–1418 Sep (1988). [16] Zhang, F. C. and Rice, T. M. Phys. Rev. B 37, 3759 (1988). [17] Leung, P. Phys. Rev. B 73, 075104 (2006). [18] Leung, P. W., Wells, B. O., and Gooding, R. J. Phys. Rev. B 56(10), 6320–6326 Sep (1997). [19] Barabanov, e. a. JETP Lett. 75, 107 (2002). [20] Zaanen, J. and Ole´s, A. M. Phys. Rev. B 37(16), 9423–9438 Jun (1988). [21] Frenkel, D. M., Gooding, R. J., Shraiman, B. I., and Siggia, E. D. Phys. Rev. B 41(1), 350–370 Jan (1990). [22] Shen, J. L. and Ting, C. S. Phys. Rev. B 41(4), 1969–1972 Feb (1990). [23] Ding, H.-Q., Lang, G. H., and Goddard, W. A. Phys. Rev. B 46(21), 14317–14320 Dec (1992). [24] Petrov, Y. and Egami, T. Phys. Rev. B 58(14), 9485–9491 Oct (1998). [25] Emery, V. J. and Reiter, G. Phys. Rev. B 38(7), 4547–4556 Sep (1988). [26] Macridin, A., Jarrell, M., Maier, T., and Sawatzky, G. A. Phys. Rev. B 71(13), 134527 Apr (2005). [27] Annett, J. F. and Martin, R. M. Phys. Rev. B 42(7), 3929–3934 Sep (1990). [28] Eder, R. and Becker, K. W. Z. Phys. B 79(7), 333 (1990). [29] Klein, L. and Aharony, A. Phys. Rev. B 45(17), 9915–9925 May (1992). [30] Digor, D. F. and et. al. Theor. Math. Phys. 75(17), 024517 (2007). [31] Hozoi, L., Nishimoto, S., Kalosakas, G., Bodea, D. B., and Burdin, S. Phys. Rev. B 75(2), 024517 Jan (2007). [32] Hozoi, L., Laad, M. S., and Fulde, P. Phys. Rev. B 78(16), 165107 Oct (2008). [33] Yanagisawa, T., Koike, S., and Yamaji, K. Phys. Rev. B 64, 184509 (2001). [34] Bonca, J., Maekawa, S., and Tohyama, T. Phys. Rev. B 76, 035121 (2007). 96  [35] Tan, F. and Wang, Q. H. Phys. Rev. Lett. 100, 117004 (2008). [36] Anderson, P. W. Phys. Rev. 115(1), 2 Jul (1959). [37] Manousakis, E. Rev. Mod. Phys. 73, 1 (1991). [38] Sushkov, O. P., Sawatzky, G. A., Eder, R., and Eskes, H. Phys. Rev. B 56(18), 11769–11776 Nov (1997). [39] Jefferson, J. H., Eskes, H., and Feiner, L. F. Phys. Rev. B 45(14), 7959–7972 Apr (1992). [40] Anderson, P. Science 235, 1196 (1987). [41] Himeda, A. and Ogata, M. Phys. Rev. B 60, R9935 (1999). [42] Sorella, S. Phys. Rev. B 71, 251103 (2005). [43] Shen, K. M., Ronning, F., Lu, D. H., Lee, W. S., Ingle, N. J. C., Meevasana, W., Baumberger, F., Damascelli, A., Armitage, N. P., Miller, L. L., Kohsaka, Y., Azuma, M., Takano, M., Takagi, H., and Shen, Z.-X. Phys. Rev. Lett. 93(26), 267002 Dec (2004). [44] Mischenko, A. S. Advances in Condensed Matter Physics 2010(306106) (2010). [45] Cataudella, V., De Filippis, G., Mishchenko, A. S., and Nagaosa, N. Phys. Rev. Lett. 99(22), 226402 Nov (2007). [46] Bonn, D. . [47] Hufner, S., Hossain, M. A., Damascelli, A., and Sawatzky, G. A. . [48] Newns, D. M. and Tsui, D. . [49] Weber, C., L¨auchli, A., Mila, F., and Giamarchi, T. Phys. Rev. Lett. 102(1), 017005 Jan (2009). [50] Sangiovanni, G., Gunnarsson, O., Koch, E., Castellani, C., and Capone, M. Phys. Rev. Lett. 97(4), 046404 Jul (2006). [51] Damascelli, A., Hussain, Z., and Shen, Z. Rev. Mod. Phys. 75, 473 (2003). [52] Ronning, F., Sasagawa, T., Kohsaka, Y., Shen, K. M., Damascelli, A., Kim, C., Yoshida, T., Armitage, N. P., Lu, D. H., Feng, D. L., Miller, L. L., Takagi, H., and Shen, Z.-X. Phys. Rev. B 67(16), 165101 Apr (2003). 97  [53] Ronning, F., Shen, K. M., Armitage, N. P., Damascelli, A., Lu, D. H., Shen, Z.-X., Miller, L. L., and Kim, C. Phys. Rev. B 71(9), 094518 Mar (2005). [54] Shen, K. M., Ronning, F., Meevasana, W., Lu, D. H., Ingle, N. J. C., Baumberger, F., Lee, W. S., Miller, L. L., Kohsaka, Y., Azuma, M., Takano, M., Takagi, H., and Shen, Z.-X. Phys. Rev. B 75(7), 075115 Feb (2007). [55] Yu, G., Li, Y., Motoyama, E. M., Zhao, X., Bariˇsi´c, N., Cho, Y., Bourges, P., Hradil, K., Mole, R. A., and Greven, M. Phys. Rev. B 81(6), 064518 Feb (2010). [56] Li, Y., Baledent, V., Yu, G., Barisic, N., Hradil, K., Mole, R. A., Sidis, Y., Steffens, P., Zhao, X., Bourges, P., and Greven, M. Nature 468(283) (2010). [57] Lawler, M. J., Fujita, K., Lee, J., Schmidt, A. R., Kohsaka, Y., Kim, C. K., Eisaki, H., Uchida, S., Davis, J. C., Sethna, J. P., and Kim, E.-A. Nature 466(347) (2010). [58] Abbamonte, P., Rusydi, A., Smadici, S., Gu, G. D., Sawatzky, G. A., and Feng, D. L. Nature Physics 1(155) (2005). [59] Merz, M., N¨ucker, N., Schweiss, P., Schuppler, S., Chen, C. T., Chakarian, V., Freeland, J., Idzerda, Y. U., Kl¨aser, M., M¨uller-Vogt, G., and Wolf, T. Phys. Rev. Lett. 80(23), 5192–5195 Jun (1998). [60] Schuster, R., Pyon, S., Knupfer, M., Fink, J., Azuma, M., Takano, M., Takagi, H., and B¨uchner, B. Phys. Rev. B 79(21), 214517 Jun (2009). [61] Holstein, T. Ann. Phys. (N.Y.) 8, 325 (1959). [62] Mila, F. Phys. Rev.B 52, 4788 (1995). [63] Magna, A. L. and Pucci, R. Eur. Phys. J. B 4, 421 (1998). [64] Cuk, T., Lu, D. H., Zhou, X. J., Shen, Z.-X., Devereaux, T. P., and Nagaosa, N. phys. stat. sol. (b) 242, 1 (2005). [65] Marsiglio, F. Phys. Rev. B 42, 2416 (1990). [66] DeRaedt, H. and Lagendijk, A. Phys. Rev. Lett. 49, 1522 (1982). [67] H. DeRaedt, A. L. Phys. Rev. B 27, 6097 (1983). [68] H. DeRaedt, A. L. Phys. Rev. B 30, 1671 (1984).  98  [69] Kornilovitch, P. J. Phys.: Condens. Matter 9, 10675 (1997). [70] Kornilovitch, P. E. and Pike, E. R. Phys. Rev. B 55, R8634 (1997). [71] Kornilovitch, P. E. Phys. Rev. B 60, 3237 (1999). [72] Kornilovitch, P. E. Phys. Rev. Lett. 81, 5382 (1998). [73] M. Hohenadler, H. G. E. and von der Linden, W. Phys. Status Solidi B 242, 1406 (2005). [74] M. Hohenadler, H. G. E. and von der Linden, W. Phys. Rev. B 69, 024301 (2004). [75] Prokof’ev, N. V. and Svistunov, B. V. Phys. Rev. Lett. 81, 2514 (1998). [76] Macridin, A. PhD thesis. [77] J. Bonca, S. A. Trugman, I. B. Phys. Rev. B 60, 1633 (1999). [78] Romero, A. H., Brown, D. W., and Lindenberg, K. Phys. Rev. B 59, 13728 (1999). [79] Romero, A. H., Brown, D. W., and Lindenberg, K. Phys. Rev. B 60, 4618 (1999). [80] Romero, A. H., Brown, D. W., and Lindenberg, K. Phys. Rev. B 60, 14080 (1999). [81] Cataudella, V., Filippis, G. D., and Iadonisi, G. Phys. Rev. B 60, 15163 (1999). [82] Cataudella, V., Filippis, G. D., and Iadonisi, G. Phys. Rev. B 62, 1496 (2000). [83] Ku, L.-C., Trugman, S. A., and Bonca, J. Phys. Rev. B 65, 174306 (2002). [84] Cataudella, V., Filippis, G. D., Martone, F., and Perroni, C. A. Phys. Rev. B 70, 193105 (2004). [85] Filippis, G. D., Cataudella, V., Ramaglia, V. M., and Perroni, C. A. Phys. Rev. B 72, 014307 (2005). [86] Barisic, O. S. Phys. Rev. B 65, 144301 (2002). [87] Barisic, O. S. Phys. Rev. B 69, 064302 (2004). 99  [88] Marsiglio, F. Phys. Rev. B 42, 2416 (1990). [89] Wellein, G. and Fehske, H. Phys. Rev. B 58, 6208 (1998). [90] Alexandrov, A. S., Kabanov, V. V., and Ray, D. K. Phys. Rev. B 49, 9915 (1994). [91] G. Wellein, H. R. and Fehske, H. Phys. Rev. B 53, 9666 (1996). [92] Wellein, G. and Fehske, H. Phys. Rev. B 56, 4513 (1997). [93] de Mello, E. V. L. and Ranninger, J. Phys. Rev. B 55, 14872 (1997). [94] Capone, M., Stephan, W., and Grilli, M. Phys. Rev. B 56, 4484 (1997). [95] Marsiglio, F. Physica C 244, 21 (1995). [96] Kornilovitch, P. E. Euro. phys. Lett. 59, 735 (2002). [97] Berciu, M. Phys. Rev. Lett. 97, 036402 (2006). [98] G. L. Goodvin, M. Berciu, G. A. S. Phys. Rev. B 74, 245104 (2006). [99] Berciu, M. and Goodvin, G. L. Phys. Rev. B 76, 165109 (2007). [100] Szczepanski, K. J. and Becker, K. W. Z. Phys. B 89, 327 (1992). [101] Rosch, O. and Gunnarsson, O. Phys. Rev. Lett. 92, 146403 (2004). [102] Poilblanc, D., Sakai, T., Scalapino, D. J., and Hanke, W. Europhys. Lett. 34, 367 (1996). [103] Slezak, C., Macridin, A., Sawatzky, G. A., Jarrell, M., and Maier, T. A. Phys. Rev. B 73, 205122 (2006). [104] Mahan, G. D. Many-Particle Physics. Plenum, NewYork, (1981). [105] Lehoucq, R., Sorensen, D., and Yang, C. ARPACK Users’ Guide: Solution of Large-scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM, US, (1998). [106] I have reported a software bug in the parallel version of the ARPACK package, but the package has not been updated as of to date. [107] Dagotto, E. Rev. Mod. Phys. 66, 3 (1994).  100  [108] Demmel, J. W. Applied Numerical Linear Algebra. SIAM, Philadelphia, (1981). [109] Simon, H. Math. Comp. 42, 165 (1984). [110] Reger, J. and Young, A. Phys. Rev. B 37, 5978 (1988). [111] Evertz, H., Lana, G., and Marcu, M. Phys. Rev. Lett. 70, 875 (1993). [112] Sandvik, A. Phys. Rev. B. 59, R14157 (1999). [113] Foulkes, W., Mitas, L., Needs, R., and Rajagopal, G. Rev. Mod. Phys. 73, 33 (2001). [114] Mezzacapo, F., Schuch, N., Boninsegni, M., and Ignaciocirac, J. New J. Phys. 11, 083026 (2009). [115] Liang, S., Doucot, B., and Anderson, P. Phys. Rev. Lett. 61, 365 (1988). [116] Lou, J. and Sandvik, A. Phys. Rev. B 76, 104432 (2007). [117] Ivanov, D. Phys. Rev. B 70, 104503 (2004). [118] White, S. and Chernysheb, A. Phys. Rev. Lett. 99, 127004 (2007). [119] White, S. and Scalapino, D. Phys. Rev. B 79, R220504 (2009). [120] Luscher, A. and Lauchli, A. Phys. Rev. B 79, 195102 (2009). [121] J Richter, et al.. Quantum Magnetism, Lecture Notes in Physics. (2004). [122] Walter, H. F. Commun. ACM 6(2), 68 (1963). [123] Marshall, W. Proc. Roy. Soc. (London) A232, 48 (1955). [124] Brugger, C., Kampfer, F., Moser, M., Pepe, M., and Wiesee, U. Phys. Rev. B 74, 224432 (2006). [125] Maier, T., Jarrell, M., Pruschke, T., and Hettler, M. Rev. Mod. Phys. 88, 1027 (2005). [126] Varma, C. Nature 468, 184 (2010). [127] Lindgren, I. and Morrison, J. Atomic Many-body Theory. Springer, Berlin, (1986). [128] Eder, R. Phys. Rev. B 59(21), 13810–13823 Jun (1999). 101  [129] Betts, D. D., e. a. Can. J. Phys. 77, 353 (1999). [130] Fauqu´e, B., Sidis, Y., Hinkov, V., Pailh`es, S., Lin, C. T., Chaud, X., and Bourges, P. Phys. Rev. Lett. 96(19), 197001 May (2006). [131] Keimer, B. to be published . [132] Leung, P. W. Phys. Rev. B 65(20), 205101 Apr (2002). [133] Anderson, P. W. Basic Notions of Condensed Matter Physics. Benjamin/Cummings, California, (1984).  102  Appendix A  Software Integrity Checks • Trial computations were performed by zeroing parameters such that the nonzero parts are known analytically. • Trial computations were performed without assuming H † = H. • All computations were checked by different combinations of at least two numerical solvers and two algorithmic routines developed independently. • Solutions are verified against those computed from larger Hilbert space without reduction due to translational, total-spin, total-spin-projection, or pointgroup symmetries. • The symmetries of solutions were always explicitly characterized. • Symmetries were checked whenever possible during intermediate steps. • For all operators that preserve 2-norm, ∑ f |⟨ f |O|i⟩|2 = 1 was checked explicitly for an initial state |i⟩ and all final states | f ⟩. • For some cases, HP ± PH was verified explicitly for a transformed Hamiltonian, H, and known symmetry operators, P. • H † − H = 0 was verified explicitly. • When possible, floating-point variables were initialized as NaN, and unsigned integers were initialized as 0xFFFFFFFFFFFFFFFF. 103  Appendix B  Computation Development This appendix first lists the software optimizations employed to perform recordbreaking numerical computations on relatively modest hardwares. The technical reasons behind the novel octapartite approach (Ch. 3-5) as well as its implementation are then discussed in details.  B.1 Software Optimizations Software was developed using C++ for the x86 and, in some cases, powerPC, platorms. Some legacy F77 libraries were integrated using a C interface. Here is a list of optimizations utilized throughout the course of this work. • Multi-thread parallelization was implemented using OpenMP due to the ease of synchronization. • Multi-node parallelization was implemented using MPI. This includes the use of processor farm for pipe-line procedures. • Performance critical part of the computation was written with in-line assembly instructions. In particular, data parallelism was exploited using SSE instructions with the 128-bit registers. • Matrix-vector product can be optimized using compressed sparse row (CSR) formulation. CSR compression can be implemented efficiently using the 104  gnu cxx::hash map template with the floating-point type overloaded with a value-zeroing default constructor. • Spin configurations of the form ∏ |σz ⟩ can be represented by 0 and 1 of an unsigned integer. Translation, spin flip, and spin-spin correlations can be performed in parallel with bitwise operators such as <<, >>, &, |, and ˆ. • Before peculiar transformations and truncations, spin and phonon basis can be enumerated using Algorithm 151[122]. • Computations with precomputed matrix transformations were optimized by the parallel use of separate computation core, write-only core, and read-only core. The throughput of the read-only process was optimized by reading from a RAID00 configuration of SSD drives connected via PCIe slots. • In general, the precomputed approach is better than on-the-fly matrix-vector product only if the basis reduction allows solution to system size that is otherwise impossible. Also, the preliminery computation should be massively parallelizable. The strategy is viable only if the code does not loop over the matrix elements that are indentically zero. This can be accomplished by paying close attention to, for example, the delta functions of the Clebsch-Gordan coefficients.  B.2 Details of the Octapartite Approach As discussed in Ch. 3- 5, the “octapartite approach” is motivated by the need of efficient modeling of large AFM background with and without the presence of extra charge carriers. The idea is to transform the spin background into a suitable basis such that the essence can be captured by a manageable number of states. There might be more efficient ways of doing so, but the current implementation was the sure way of producing interesting results for a set period of time for a pioneering project. The aim was to demonstrate an overall idea which combinatorially decreases the number of states required to model the interesting part of the Hilbert space. This reduction is exactly the reason why, even without the most sophisticated implementation, technological barriers were immediately broken for 105  combinatorially large system when available computational power advances only exponentially over time according to Morse’s Law. In particular, the implementation partitions the computation into consecutive stages for two main reasons. First, the Hilbert space truncation for picking out the “ N8 spins adding to zero” charateristic (Sec. 3.3) requires a peculiar quantummechanical reformulation with no apparent way of enumeration. Developing matrixvector product within such strange basis is highly prone to errors, and stage-bystage consistnecy checks are required anyways in the development process. Second, the strategy allowed me to simultaneously run one stage of computation while coding for the next. By doing so, the “precomputation time” mentioned in the last section is effectively zero because coding and quality assurance took longer than actual run time. The rest of this appendix will point out some general observations of the problem then describe the current implementation. Future developers are encourage to consider some general observations before either duplicating the current implementation or developing their own.  B.2.1 General Observations There are several challenges in working within a basis formulated to single out the “ N8 spins adding to zero” subspace. First, the truncation of less important states should be controllable systematically and flexible enough to adapt to the many unknowns of doped systems. Second, within the truncated basis transformed away from the natural z-projected description, there should be a fast way of indexing initial and final states upon Hamiltonian operations. Third, one should have apriori knowledge about the identically-zero matrix elements in order to avoid them for sparse matrices’ performace scaling. Let’s review several general aspects of an arbitrary basis for the spin problem in Ch. 3-5. For a small number of oxygen holes and a copper spin background, an arbitrary translationally invariant basis can be written as ( 1 √ ∑ eiKl p†l,σl N  )  ∏ p†l+∆l,σ ∆l  106  ∆l  |σ ⟩l  (B.1)  with one oxygen hole and a spin background centered at position l and all other oxygen holes positioned ∆l away from the reference point l. The Hamiltonian acted upon an initial state would produce some final states scaled by matrix elements. Each term in the Hamilonian is a summation over all unit cells; therefore, the final states always come in the form of a Fourier sum, ∑ eiKl . The overlap between the final state and orthnormal basis states is thus some complex constant with a sum, ∑ δ , wherein delta functions arise from the phases, the fermion operators and the bra-ket of the spin backgrounds. On top of a particular spin background σ l , 1 N  the different configurations of oxygen holes, p†∆l,σδ l , can be enumerated as shown in Sec. 5.3. The remaining part is then the connections between different spin backgrounds. A spin polaron model of the form introduced in Ch. 4 and 5 would yield several types of connections. Direct oxygen-oxygen hopping Tpp moves either the “nonreference” oxygen holes or the oxygen hole at the reference position, l. The former is simply changing the ∆l tabulation while the later result in no p† occupation at the reference point l. This simply means that the final state can be specified in terms of a new reference point l ′ occupied by a p† hole. The matrix element is thus ′  complex with a phase eiK(l−l ) and a spin background shifted by l − l ′ : ′  eiK(l−l ) ⟨σ ′ |Tl−l ′ |σ ⟩,  (B.2)  which requires the knowledge of how one spin configuration connect to another upon translation. Noting that the oxygen holes p† are already taken care of, the Jpd exchange z,± requires knowledge of S∆L . Sz is computed using  1 † † z ⟨σ ′ |S∆L |σ ⟩ = ⟨σ ′ |d∆L,↑ d∆L,↑ − d∆L,↓ d∆L,↓ |σ ⟩. 2  (B.3)  S± requires information about kicking out an original spin, α , then injecting a new spin β = −α : † ⟨σ ′ |K∆L (α , β )|σ ⟩ = ⟨σ ′ |d∆L, β d∆L,α |σ ⟩.  (B.4)  The Tswap hopping is the combination of the direct oxygen-oxygen hopping, plus kicking process K∆L (α , β ) with a α to β spin swapping. Note that β is now 107  the original spin of the oxygen hole with no definite sign relative to α . Jdd matrix elements involve two copper spins and cannot be tabulated as easily. For nearest neighbor separation a, one would require the overlaps z z ⟨σ ′ |S∆L S∆L+a |σ ⟩  (B.5)  ⟨σ ′ |K∆L (α , −α )K∆L+a (−α , α )|σ ⟩  (B.6)  Note that for a full Hilbert space specified in terms of the z-projection of copper spins, all these connections can be computed trivially by flipping bits of an unsigned integer and using Algorithm 151[122] to enumerate the states of the basis. However, such full Hilbert space is unmanageable for the large system sizes of interest (Sec. 5.1) and an alternative way is needed. The observations made in Ch. 3 are specified in terms of sA/B , the total spin of the sublattices; therefore, a formulation in terms of total quantum mechanical spin is required to reduce the dimension of the Hilbert space. The task is to find one such formulation which also allows easy computation of the above matrix elements. The following section discusses a way of doing so.  B.2.2 Implementations The previous section pointed out that the matrix representation of an effective spin polaron model requires several “basic” operations on the spin background: z . K∆L (α , β ), Tl−l ′ , and S∆L  There are two immediate advantages of building the full spin background in terms as a Clebsch-Gordan series of total-sub-lattice-spin states (Eq. 3.9). First, it is apparently easier to truncate out the low-sub-lattice-spin part of the Hilbert space according the CS criterion (Ch. 3.3.2). Second, note that processes Tpp/swap and Jpd reqires a single operation of one of the three “basic” operators, and the Jdd z copper-copper exhange requires two consecutive operations of either K∆L or S∆L .  Jdd is bi-partite in nature so the pair of consecutive operations is decomposed into two simultaneous operations on orthogonal part of the Hilbert space, for example ⟨σ ′ |K∆L (α , −α )K∆L+a (−α , +α )|σ ⟩ → ⟨σA′ |K∆L (α , −α )|σA ⟩⟨σB′ |K∆L+a (−α , α )|σB ⟩. (B.7) 108  It was found that the choice of sub-lattice basis yields extremely sparse matrices. The consecutive matrix multiplications within the final overall basis are reduced to matrix elements look up within a sparse matrix of a much smaller basis. Having established the usefulness of bi-partite separation in dealing with truncation as well as the the evaluation of the Jdd term, the remaining question is how to formulate such spin basis. The truncation according to sA/B requires the sublattice spins to be added quantum mechanically, e.g. in terms of Clebsch-Gordan series. However, the hopping processes require knowledge of how each of the spin background relate to one another upon the 2D translation Tl−l ′ , which has no immediately apparent relationship with a spin addition exercise for a large number of spins. For example, the apparent way of determinationing the overall sign of a matrix element is to follow the signs of Clebsch-Gordan addition for both the initial and final states. The tabulation of the N2 sublattice spins into a basis of |sA/B , szA/B ⟩ can be done ) ( by adding i spins with N2 − i spins. One way to proceed is to add the group of spins neighbored by the oxygen holes with the rest, but such approach would complicate the determination of T( l, l ′ ) among spin background. For example, each oxygen configuration, ∏ p†l+δ l,σ , would have a different basis for the spin backδl ground centered at l. While this approach could potentially be useful for detailed truncation of the disturbed background, I decided to use a simple coarse-grain approach (Sec. 3.3.2) to eradicate this complication. As discussed there, the square lattice is partitioned into eight sublattices such that the sites of a particular sublattice are connected by the vectors (Fig. 3.1). For a group of N 8  N 8  √1 2  (a, ±a)  spins, I’d start with the z-projected representation  of 2 states and transform into a Clebsch-Gordan basis |s N , szN ⟩. I would repre8  8  sent each z-projected basis state using bits of an unsigned integer and store each Clebsch-Gordan series as gnu cxx::hash map<size t,double>. With a particular enumeration of |s N , szN ⟩ states, I can mix two identical enumerations to 8  8  build |s N , szN ⟩ according to Clebsch-Gordan addition. Then with two identical enu4  4  merations of |s N , szN ⟩, I’d build a single enumeration of |s N , szN ⟩. Finally, I can sim4  2  4  2  ilarly enumerate the overal background |sN , szN ⟩ from the enumeration of |s N , szN ⟩. 2  2  The enumeration can be specified by a non-negative integer, state index. The state represented by a particular index value can be derived from the “parent” basis 109  with half the number of spins, for example using the following loop structure. It is apparent that truncation can be performed at each stage of mixing by simply tuning the values of min2S, max2S, minSub2S, and maxSub2S. s i z e t s t a t e i n d e x =0; / / i n g e n e r a l t o t a l s p i n can be odd m u l t i p l e s o f 1 / 2 , so I worked w i t h 2S c o n s t i n t min2S= m i n i m u m 2 t i m e s t o t a l s p i n ; c o n s t i n t max2S= n u m b e r o f s p i n s ; c o n s t i n t minSub2S : = m i n i m u m 2 t i m e s s u b l a t t i c e s p i n ; c o n s t i n t maxSub2S : = n u m b e r o f s u b l a t t i c e s p i n s ; / / s t a t e s w i t h h i g h e r t o t a l s p i n have h i g h e r i n d i c e s f o r ( i n t ss=min2S ; ss<=max2S ; ss +=2) { / / l o o p over t h e l e f t k e t f o r ( i n t ssa =0; ssa<=maxSub2S ; ssa +=2) { / / l o o p over t h e r i g h t k e t f o r ( i n t ssb =0; ssb<=maxSub2S ; ssb +=2) { i f ( C l e b s c h G o r d a n c o e f f i c i e n t s ==0) c o n t i n u e ; c o n s t s i z e t na= number of ( ssa +1) b l o c k s w i t h t o t a l s p i n s s a ; c o n s t s i z e t nb= number of ( ssb +1) b l o c k s w i t h t o t a l s p i n s s b ; f o r ( s i z e t aa =0; aa<na ;++ aa ) { f o r ( s i z e t bb =0; bb<nb ;++ bb ) { / / l o o p over t h e ss+1 components o f z−p r o j e c t i o n f o r ( i n t zz=−ss ; zz<=ss ; zz +=2 ,++ s t a t e i n d e x ) { f p r i n t f ( stderr , ” S t a t e %l l u has s=%4 f and sz=%+4f , p r o d u c t o f ” , s t a t e i n d e x , 0 . 5 ∗ ss , 0 . 5 ∗ zz ) ; f p r i n t f ( stderr , ” b l o c k %l l u /% l l u o f t h e s=%4 f s e c t o r o f l e f t k e t and ” , aa , na , 0 . 5 ∗ ssa ) ; f p r i n t f ( stderr , ” b l o c k %l l u /% l l u o f t h e s=%4 f s e c t o r o f r i g h t k e t ” , bb , nb , 0 . 5 ∗ ssb ) ;  } / / zz } / / bb } / / aa } / / ssb } / / ssa } / / ss  Each state indexed by state index is associated with two blocks of states in the smaller parent basis. Each block contains the 2S + 1 projections related by S+/− . This is done to accomodate the state’s Clebsch-Gordan series, which has a delta function, δszT =sza +szb in the summation over z-projections. This blocking scheme immediately reveals the “parent” basis states which are releveant to that particular 110  Clebsch-Gordan series. The benefit of this enumeration is an instantaneous “reverse lookup”. When performing Hamiltonian operations, one is really interested in the non-zero overlap between the outgoing states and those in the orthonormal basis. The naive way is to compute the dot product againsts all basis states, but this adds a O(N) layer on top everything else and is detrimental in the case of large systems. Under the above looping scheme, an increasing state index is associated with increasing values of ssa, aa, ssb, and bb. Any arbitrary ket, |σa , σaz ⟩|σb , σbz ⟩, is trivially associated with these four indices so states with non-zero overlap is known immediately, with an O(1) reverse lookup operation (in analogy to Algorithm 151[122]). This hierarchy allows an efficient computation of the overall matrix in terms of z . The initial information needed is some small matrices O K∆L (α , β ), Tl−l ′ , and S∆L  within the smallest  N 8  spin basis. The matrices for the larger children basis required  z , one needs to evaluate two general operations. For K∆L (α , β ) and S∆L  ⟨σa , σaz |OL |σa , σaz ⟩ × ⟨σb , σbz |σb , σbz ⟩ ⟨σa , σaz |σa , σaz ⟩ × ⟨σb , σbz |OR |σb , σbz ⟩  (B.8) (B.9)  for the N8 -spin basis, use the result for the same procedure for the N4 -spin basis, and eventually propagated to the N-spin basis. The operation needed for each stage of the heirachy is the same and the computation is essentially the aforementioned trivial “revesere lookup” after having the indices and value of one ket modified by O of the smaller “parent” basis. For Jdd in the largest N-spin basis, the indices and values of both left and right kets are modified: ⟨σa , σaz |OL |σa , σaz ⟩⟨σb , σbz |OR |σb , σbz ⟩,  (B.10)  and again these are just trivial lookup operations. The octapartite division groups the sites such that sites from a particular group are neighbored by sites of only one other group in a particular direction (Fig.3.1). If one expand the spin basis in terms of the eight sublattice (Eq. 3.13), one can immediately see that the special neighboring pattern leads to simple computation which is the above simultaneous  111  OL/R lookup followed by a possible swap of left and right ket. ⟨σa , σaz |OR |σb , σbz ⟩⟨σb , σbz |OL |σa , σaz ⟩,  (B.11)  Due to the fact that the reverse lookups are performed according to the indices ssa, aa, ssb, and bb, the need of swaping the left and right ket does not introduce any complication. z Finally, I note that the computation for K∆L (α , β ), Tl−l ′ , S∆L , and each bond  of the Jdd terms are indepedent of one another and thus trivially parallelizable. The bottleneck of the overall computation is actually the vector size that can be accomodate by the eigen-solvers. The coarse-grain truncation discussed here leads to vector size of < 10GB for the undoped case with 64 spins and the two-hole case with 32 CuO2 unit cells.  112  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0071682/manifest

Comment

Related Items