Effects of the lattice distortion on themagnetic order in rare-earth nickelatesbyStepan Olegovich FomichevB.Sc., The University of Toronto, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2019c© Stepan Olegovich Fomichev 2019The following individuals certify that they have read, and recommend tothe Faculty of Graduate and Postdoctoral Studies for acceptance, the thesisentitled:Effects of the lattice distortion on the magnetic order in rare-earth nickelatessubmitted by Stepan Olegovich Fomichev in partial fulfillment of therequirements for the degree of Master of Science in Physics.Examining Committee:Mona Berciu, Physics and AstronomySupervisorGeorge Sawatzky, Physics and AstronomyExamining Committee MemberiiAbstractThe low-temperature magnetic order in the rare-earth nickelates is a subjectof vigorous debate in the literature. Recent work emphasized the primaryrole of the electron-phonon coupling for the metal-insulator transition inthe nickelates, and suggested that lattice distortions are the driver of thetransition, leading to the observed charge order. However, to our knowl-edge there has been little work on the impact of lattice distortions on themagnetic order, in particular whether distortions favour some orders overothers. In this thesis, we study the magnetic order in the nickelates at zerotemperature, and investigate whether the breathing-mode lattice distortionsselect a preferred ground state. An effective two-band Hubbard model forthe nickelates is constructed and coupled to the lattice distortions with anon-site Holstein-like term. The distortions are treated semiclassically. Us-ing the Hartree-Fock approximation, we obtain the magnetic phase diagram,then turn on the coupling to the lattice to observe its impact on the var-ious phases. Our model reproduces the earlier work showing the strongercharge disproportionation and insulating behaviour in the phase space dueto increased coupling to the lattice. Furthermore, we find numerous 4-sitemagnetic orders that are self-consistent within the model, including all of themain suggestions in the literature (states such as ↑↑↓↓, ↑→↓← and ⇑ 0 ⇓ 0).However, in this model a magnetic order can only couple to the lattice dis-tortions if there is nonzero charge disproportionation. As a result, we findthat coupling to the lattice distortions broadly favours the ⇑ 0 ⇓ 0 order inlarge sectors of the parameter space. Finally, we considered the impact oflonger range hopping on the magnetic order: we find that the shape of thedensity of states, rather than overall bandwidth, primarily determines themagnetic ground state. A van Hove singularity arises even for small 2nd-nearest hopping amplitudes, which results in robust ferromagnetism acrossmost of the phase diagram in a Stoner-like fashion. On the contrary, evensmall 4th-nearest amplitudes decrease the Fermi level density of states, re-sulting in ballooning of the metallic phase despite a barely renormalizedbandwidth.iiiLay SummaryWhile the study of magnetism began as a hobby for philosophers, it hassince led to deep insights about our physical world. Aside from the obviousutility of fridge magnets, understanding the atomic origin of magnetismhas given humanity many useful tools: MRI machines, loudspeakers, harddrives. Magnetism originates with an electronic property – spin: electronsare essentially tiny bar magnets. In most magnets, electrons align in onedirection, producing strong magnetism on the macroscopic scale. But rare-earth nickelates – at low temperature – become a peculiar kind of a magnet.Instead of aligning, the electrons exhibit an unknown pattern repeating everyfour lattice sites, potentially ↑↑↓↓ or ↑→↓←. In this thesis, we constructa quantum-mechanical model for the nickelates to study their magnetism.Using the magnetic pattern of this material to encode bits of information isjust one potential technological byproduct that could become available asour understanding of the material grows.ivPreface• I wrote this thesis in its entirety. The literature review and summaryin Chapter 1 arise from my own reading and analysis. The model forthe nickelates suggested in Chapter 2 is due to my supervisor for theproject, Mona Berciu, and our collaborator Giniyat Khaliullin fromthe Max Planck Institute at Stuttgart. The analytical calculationspertaining to the model, including the Hartree-Fock calculations, inChapter 3 have been carried out originally by Mona Berciu and sub-sequently repeated by me. I envisioned and constructed the numericalimplementation of the Hartree-Fock equations as described in Chapter4, including all of the additional techniques beyond the textbook ap-proach, such as Pulay mixing. All the results in Chapter 5 have beenobtained by me using the program described in Chapter 4. Finally,the analysis in Chapters 5 and 6 was conducted by me, in consultationwith Mona Berciu.• A version of the work discussed in Chapters 2, 3, 4, 5 is to be submit-ted for publication under the authorship of Stepan Fomichev, GiniyatKhaliullin, and Mona Berciu, with the manuscript presently in prepa-ration.• UBC Research Ethics Board approval was not required for this work.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction: Experiments and Calculations on the Nicke-lates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Crystal structure . . . . . . . . . . . . . . . . . . . . . . . . . 41.2 Structural Transition . . . . . . . . . . . . . . . . . . . . . . 61.3 Metal-Insulator Transition: Introduction . . . . . . . . . . . 71.3.1 Formal Valence rules . . . . . . . . . . . . . . . . . . 91.3.2 Crystal field splitting of RNiO3 . . . . . . . . . . . . 101.3.3 Zaanen-Sawatzky-Allen scheme . . . . . . . . . . . . 131.4 Metal-Insulator Transition: Mechanisms . . . . . . . . . . . . 151.4.1 Mechanism: Charge-transfer gap . . . . . . . . . . . . 151.4.2 Mechanism: Jahn-Teller effect . . . . . . . . . . . . . 151.4.3 Mechanism: Charge disproportionation . . . . . . . . 171.4.4 Mechanism: Negative charge-transfer theory . . . . . 181.5 Magnetic Transition: Introduction . . . . . . . . . . . . . . . 191.6 Magnetic Transition: Mechanisms . . . . . . . . . . . . . . . 222 Effective Two-Band Model for the Rare-Earth Nickelates 242.1 Hopping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.1.1 1st nearest neighbour hopping . . . . . . . . . . . . . 27viTable of Contents2.1.2 4st nearest neighbour hopping . . . . . . . . . . . . . 282.1.3 2nd nearest neighbour hopping . . . . . . . . . . . . . 282.2 Lattice contributions . . . . . . . . . . . . . . . . . . . . . . 302.3 Electron-electron interactions . . . . . . . . . . . . . . . . . . 302.4 Electron-lattice interaction . . . . . . . . . . . . . . . . . . . 312.5 Final form of the Hamiltonian Hˆ . . . . . . . . . . . . . . . . 323 Hartree-Fock Approximation for Rare-Earth Nickelates . 333.1 Minimizing lattice contributions . . . . . . . . . . . . . . . . 333.2 Minimizing electronic contributions . . . . . . . . . . . . . . 363.3 Mean-Field Parameters . . . . . . . . . . . . . . . . . . . . . 403.4 Energy of the Hartree-Fock ground state . . . . . . . . . . . 443.5 Hartree-Fock Hamiltonian . . . . . . . . . . . . . . . . . . . 474 Numerical implementation of the Hartree-Fock approxima-tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.1 The Hartree-Fock iterative algorithm . . . . . . . . . . . . . 544.2 Strategies for improving computation . . . . . . . . . . . . . 564.2.1 Pulay mixing . . . . . . . . . . . . . . . . . . . . . . . 564.2.2 Convergence in iteration count and momentum sam-pling . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.2.3 Choosing initial guesses . . . . . . . . . . . . . . . . . 594.2.4 Parallelization . . . . . . . . . . . . . . . . . . . . . . 604.2.5 Boot-strapping with neighbouring points . . . . . . . 605 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 615.1 Charge order and the Metal-Insulator Transition . . . . . . . 625.2 Magnetic order: types and the transition . . . . . . . . . . . 675.3 Impact of the hopping integrals . . . . . . . . . . . . . . . . 735.4 The anharmonicity parameter α . . . . . . . . . . . . . . . . 785.5 Convergence and stability . . . . . . . . . . . . . . . . . . . . 786 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86AppendicesA Calculating the hopping coefficients tab(k) . . . . . . . . . . 94A.1 1st and 4th neighbour hopping . . . . . . . . . . . . . . . . . 94viiTable of ContentsA.2 2nd neighbour hopping . . . . . . . . . . . . . . . . . . . . . 98A.3 Full hopping operator . . . . . . . . . . . . . . . . . . . . . . 101B Minimizing lattice contributions: solving the cubic equation 102C Calculation of 〈Hˆe〉 . . . . . . . . . . . . . . . . . . . . . . . . . 103D Fourier Transforming HˆHF . . . . . . . . . . . . . . . . . . . . 107E Re-calculating the Hartree-Fock order parameters . . . . . 111viiiList of Figures1.1 A temperature vs. type of rare-earth ion phase diagram ofthe rare-earth nickelates, adapted from Ref. [1]. . . . . . . . 21.2 A graphic of the rocksalt breathing-mode distortion in therare-earth nickelates. Adapted from Ref. [2]. . . . . . . . . . 21.3 From Medarde’s excellent 1997 review, Ref. [3]. (a) Ideal per-ovskite structure; (b) Rhombohedral distortion twists aboutthe [111] direction for t ∼ 1 (R3c point group symmetry class)(c) Orthorhombic distortion twists about the [110] direction –the case for most of the nickelates (dominant mode at smallert). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.4 Ni-O-Ni bond angle in the rare-earth nickelates is a measureof the perosvkite distortion, equivalent to the tolerance factor.Figure from Ref. [4]. . . . . . . . . . . . . . . . . . . . . . . 51.5 Various distortion modes available to the nickelates. Figurefrom Ref. [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . 81.6 Electrical resistivity of PrNiO3 as a function of temperature:the metal-insulator transition is clearly visible. Data and fig-ure from Ref. [6]. . . . . . . . . . . . . . . . . . . . . . . . . 91.7 The eg and t2g manifolds resulting from crystal field splitting.Figure from Ref. [7]. . . . . . . . . . . . . . . . . . . . . . . 121.8 A graphic depicting the Zaanen-Sawatzky-Allen classificationscheme, using featureless densities of electronic states withouthybridization. Figure from Ref. [8]. . . . . . . . . . . . . . . 141.9 The lowest order excitations within the Zaanen-Sawatzky-Allen classification scheme. Figure from Ref. [8]. . . . . . . . 151.10 The unit cell of KCuF3 and the associated JahnTeller distor-tion modes. Figure from Ref. [7]. . . . . . . . . . . . . . . . 161.11 The Ni magnetic moment in PrNiO3 as a function of tem-perature, with the onset of magnetism clearly visible. Figurefrom Ref. [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . 20ixList of Figures1.12 The appearance of a new ordering peak in neutron diffrac-tion data for DyNiO3 as the temperature is lowered past themagnetic transition temperature TN. Figure from Ref. [9]. . 212.1 A schematic diagram of the various hoppings included in thenickelates two-band model. The green arrow corresponds tothe first neighrest neighbour hopping; the yellow to secondnearest neighbour hopping; third nearest neighbour hopping,represented by the red arrow, is not included due to mini-mal orbital overlap; and fourth nearest neighbour hopping isshown by the purple arrow. . . . . . . . . . . . . . . . . . . . 265.1 Charge modulation δ (the magnitude is reflected in the colourbar) of the ground state in the U–J plane. Other parametersare t1 = 1, t2 = 0.15, t4 = 0, b = 0. Resolution is 40× 40. . . 635.2 Representative density of states diagrams for the various phasesin Fig. 5.1: (a) metallic phase, U = 1.077, J = 0.538; (b) itin-erant magnetic phase (see the magnetic phase analysis belowfor more details about this phase), U = 3.231, J = 0.385;(c) charge modulated phase, with a clear gap at the Fermilevel, U = 2.0, J = 1.538; (d) Mott insulating phase, U =5.846, J = 0.231. Other parameters as in Fig. 5.1, with theaddition of α = 1. Notice the increased weight (called thevan Hove singularity) at the lower band edge: its presence isdue to the nonzero t2 parameter, which introduces a strongasymmetry to the DOS. More on this below. . . . . . . . . . 645.3 Charge modulation δ (see colour scale) in the HF ground-state, as a function of U and b, for J/U = 0.2, 0.3 (top leftand right, respectively) and 0.4, 0.5, (bottom left and right,respectively). Other parameters are t1 = 1, t2 = 0.15, t4 =0.25, α = 1. Resolution is 20× 20. . . . . . . . . . . . . . . . 655.4 The HF ground-state is metallic (deep blue) or insulating(yellow). The results are shown in the U -b space, for J/U =0.2, 0.3, 0.4 and 0.5, respectively (panels arranged as in Fig.5.3). All other parameters are as in Fig. 5.3. . . . . . . . . . 66xList of Figures5.5 Top: U -J phase diagram of magnetic order, in the absence ofcoupling to the lattice (b = 0). Bottom: Same when couplingto the lattice is turned on (b = 0.8). Other parameters aret1 = 1, t2 = 0.15, t4 = 0, α = 1 for both. The black-line con-tours indicate the value of charge modulation δ. Resolutionis 40× 40. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.6 Magnetic order in the HF ground state, as a function of Uand b, for J/U = 0.2, 0.3 ((a) and (b), respectively) and0.4, 0.5, ((c) and (d), respectively). Other parameters aret1 = 1, t2 = 0.15, t4 = 0.25, α = 1. . . . . . . . . . . . . . . . 715.7 Energies of the various converged states (y axis) as a functionof b (x axis). The parameters are U = 3.158, J = 0.2U, t1 =1, t2 = 0.15, t4 = 0.25. Resolution is 20× 20. . . . . . . . . . 725.8 Magnetic order in the HF ground state, as a function of U andJ , for t2 = 0, 0.15, 0.25 (top to bottom). Other parametersare t1 = 1, t4 = 0, b = 0. . . . . . . . . . . . . . . . . . . . . 745.9 The non-interacting (U = J = 0) DOS, corresponding tothe systems in Fig.5.8. Notice that increasing the bandwidth(even if just modestly by at most %15), paradoxically, leadsto more robust ferromagnetism at lower U — a consequenceof the van Hove singularity at the lower band edge. Alsonotice how when the next-nearest neighbour frustration ismaximally reduced (t2 = 0), the non-collinear 4-site magneticorder dominates the collinear one. Currently it is not clearto us why this would be the case, given how introducing t2, t4seems to affect them equally based on pure lattice frustrationarguments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 755.10 Magnetic order in the HF ground state, as a function of Uand J , for t4 = 0, 0.1 ((a) and (b), respectively), and t4 =0.25, 0.35 ((c) and (d)). The bandwidth is W = 6.4 in allcases. Other parameters are t1 = 1, t2 = 0.15, b = 0. Thegrowth of the metallic region is clearly not the effect of arenormalized bandwidth, but rather due to the changes ofthe shape of the DOS. . . . . . . . . . . . . . . . . . . . . . . 76xiList of Figures5.11 Energies of several converged self-consistent HF states as afunction of t4, relative to the energy of the metallic state.The different colours correspond to different magnetic orders(see the legend). The parameters are U = 4, J = 0.316, t1 =1, t4 = 0, b = 0 for (a) and U = 5, J = 0.316, t1 = 1, t2 =0.15, b = 0 for (b). Notice how in (b) the relative ener-gies of the chief magnetic ground state contenders are notaffected by the change in the hopping rate t4 – except forferromagnetism, which gets strongly frustrated with the ad-ditional t4 hopping and, paradoxically, “unfrustrated” withthe introduction of t2 hopping due to DOS effects (see thetext for details). Meanwhile, in (a) with tuning the t2 rateaway from 0 the non-collinear ↑→↓← gets briefly favoured,but then quickly loses out to the collinear ↑↑↓↓, before ferro-magnetism begins to reign supreme. . . . . . . . . . . . . . . 795.12 Magnetic order in the HF ground state, as a function of U andb, for J/U = 0.2, and two different values of the anharmonic-ity α: α = 0 (top) and α = 1 (bottom). Other parametersare t1 = 1, t2 = 0.15, t4 = 0.25. . . . . . . . . . . . . . . . . . 805.13 A phase portrait of the iterative sequences for various start-ing parameters in the δ, S1z, S2z parameter subspace. On thex axis we plot the difference between the spin parametersof the sublattices, S1z − S2z: hence the perfect symmetrypoint is at the origin. On the y axis is the disproportion-ation parameter δ: the stable point ⇑ 0 ⇓ 0 is thus at thetop right and equivalently bottom left. The parameters areU = 1.5, J = 3, t1 = 1, t2 = 0.15, t4 = 0.55, b = 0. De-spite the relatively dense sampling of the phase space, the⇑ 0 ⇓ attractors dominate the dynamics, with no other states,specifically of the partially disproportionated kind ⇑↑⇓↓, be-ing convergent. . . . . . . . . . . . . . . . . . . . . . . . . . . 82xiiAcknowledgmentsConducting research and writing a thesis has been a massive undertaking,and I am extremely grateful to Mona Berciu for guiding me through it.From suggesting the problem (together with Giniyat Khaliullin), to slowlyshowing me the ropes of being a theoretical physicist, and to supporting meat conferences, seminars, and group meetings, Mona has been exceptionallykind, attentive, and encouraging.And the tightly-knit, warm and welcoming environment she nurturedwithin her group was indispensable to all of our academic pursuits, and toour sense of well-being. I am really grateful for all the discussions, hikes,ice-cream breaks, and escape room pursuits that helped me to get to knowbetter so many wonderful physicists, including Nathan Cheng, Mirko Mo¨ller,and John Sous. Without their peer-to-peer guidance, sharp questions andhelpful answers, I would be hard pressed to learn as much as I did: andwithout the coffee breaks and the jokes, it would have been hard to keepgoing.I would especially like to thank my long-time friend Kyle Wamer, whohas valiantly endured two years of sitting across from me and debating thevery foundations of condensed matter, by asking some silly questions (“so,what are phonons, anyway?”); and who has single-handedly introduced meto Vancouver – the city, the people, and the mountains.There are many wonderful people at UBC that have made this an un-forgettable journey, both in and outside the classroom. I am thankful toAli, Jordan, Christian, Connor, Adam, Tarun, Victor, Ian and Khoi, for theconferences, hikes, classes, ski hills, exams, lake-side cabins, boiler rooms,and coffee shops we enjoyed together.Lastly, I would like to thank Evgenia, for her loving support and cheerfulcompany, which made every moment of this journey exciting and memorable.xiiiTo my mother, without whose selfless sacrifices I would never be where I amtoday. And to my father, whose endless patience for my endless questionsnurtured my love of physics.xivChapter 1Introduction: Experimentsand Calculations on theNickelatesThe ultimate court of appeal isobservation and experiment...not authority.Thomas Huxley (1875)Rare-earth nickelates (chemical formula RNiO3, with R any rare-earthelement from La to Lu) are a fascinating perovskite transition-metal oxide[10] material series that have not only generated much excitement in thecondensed matter community as a matter of fundamental interest [3, 4], butare also now inspiring more and more potential applications, specificallyin the memory domain [11–13]. Obtained for the first time in 1971 byDemazeau et al. [14], they were largely ignored until the general surge ofinterest in perovskite materials in the 90s, generated by the discovery in oneof them – the cuprates – of high-temperature superconductivity in the late80s [15].In terms of basic science, the nickelates exhibit a variety of phase tran-sitions and emergent orders. All nickelates except LaNiO3 undergo – in factare considered the prime examples of – a sharp metal-to-insulator (MIT)transition [16] (see Fig. 1.1 for a phase diagram), which is generally under-stood to be a consequence of charge order (although the exact mechanism isstill disputed [17–22]). The transition is concurrent with a breathing-mode,rock-salt pattern lattice distortion [23], as shown in Fig. 1.2. In addition,for all nickelates except LaNiO31, at still lower temperatures (and for Pr andNd, at the same temperature as the MIT) there is an additional, magneticphase transition, resulting in an as-yet poorly-understood magnetic order,with various competing proposals [4, 25–28]. While the phases themselves1although this has recently been challenged [24].1Chapter 1. Introduction: Experiments and Calculations on the NickelatesFigure 1.1: A temperature vs. type of rare-earth ion phase diagram of therare-earth nickelates, adapted from Ref. [1].Figure 1.2: A graphic of the rocksalt breathing-mode distortion in the rare-earth nickelates. Adapted from Ref. [2].2Chapter 1. Introduction: Experiments and Calculations on the Nickelatesare reasonably well-understood and characterized (with the exception of thelow-temperature magnetism, which is largely due to the difficulty of synthe-sizing large enough single crystal samples), the precise mechanisms of thecharge and magnetic phase transitions are not yet entirely clear. Moreover,almost all of these transitions can be tuned with the type of rare-earth ion(as in Fig. 1.1), pressure [29], in heterostructures [30], and even swappingout the oxygen for its isotopes [31].One of the thorniest questions in the field in recent years has been thenature and origin of the magnetic order that obtains in this material seriesat low temperatures. Specifically, while neutron diffraction measurementshave firmly established the ordering wavevector of the magnetic phase asQm =2pia [14 ,14 ,14 ]pseudocubic2, the community is still debating the origin andprecise alignment of the local magnetic moments. Some of the magneticorders potentially consistent with the above ordering wavevector include⇑ 0 ⇓ 0 [25], as well as ⇑↑⇓↓ [26] and even the non-collinear option ⇑→⇓←(throughout this thesis we use fat arrows to indicate larger spin magnitudes).In light of recent work demonstrating that the lattice breathing-mode distor-tion helps stabilize charge order in the nickelates [32], it seemed reasonableto suggest that the lattice distortion could also preferentially favour one orthe other kind of magnetic orders suggested in the literature. In this thesis,we address the magnetic order question by constructing an effective two-band Hubbard Hamiltonian for the nickelates, and coupling it through anon-site interaction to the lattice distortion, which we treat semiclassically.We then search for the magnetic ground state in the Hartree-Fock approxi-mation, and consider how the presence of coupling to the lattice distortionaffects the magnetic phases that arise.This thesis is structured as follows: in the rest of this chapter we beginby describing, at the appropriate level of detail, the phenomenology of thisNi material series, including structural, charge and magnetic properties – asgleaned from experiment. In Chapter 2, we construct in detail the micro-scopic model used for the nickelates in our study. In Chapters 3 and 4, weapply and numerically implement the Hartree-Fock approximation for themodel. In Chapter 5, we present our results and answer the main researchquestions of the thesis, as well as comment on any unexpected discoveries.Our conclusions and outlook are summarized in Chapter 6.2 2pia[ 12, 0, 12] in perovskite notation.31.1. Crystal structureFigure 1.3: From Medarde’s excellent 1997 review, Ref. [3]. (a) Ideal per-ovskite structure; (b) Rhombohedral distortion twists about the [111] direc-tion for t ∼ 1 (R3c point group symmetry class) (c) Orthorhombic distortiontwists about the [110] direction – the case for most of the nickelates (domi-nant mode at smaller t).1.1 Crystal structureAt room temperature, most nickelates – whether found in nature or synthe-sized in the laboratory – consist of the usual columns of oxygen octahedraengulfing the nickel atoms, with the rare-earth ions sitting in-between them(see Fig. 1.3a for a diagram of the crystal structure). For all the nickelates,the structure actually is strictly speaking not an ideal perovskite lattice: ithas a slight deviation from perfect cubic lattice symmetry due to “twists” ofthe octahedra, caused by the too-small -to-fit rare-earth ions (many of thelattice parameters are tabulated in [33]; for more general information aboutthe crystal structure, see [3]). In practice, these deviations – characterizedby the Goldschmidt tolerance factor t [34], defined ast =dR-O√2dNi-O, (1.1)where d are the bond lengths of R-O and Ni-O, respectively, – are rathersmall and often either entirely ignored [26], or incorporated phenomenolog-ically (e.g. through modifying the hopping rates in a tight-binding model).The range of variation of t is fairly small – from t = 0.932 for LuNiO3 tot = 0.986 for LaNiO3. As an aside, sometimes also the Ni-O-Ni bond angleis used in place of the tolerance factor to characterize the distortion: it ismarked in Fig. 1.4.The most stark phase change in the nickelate phase diagram is betweenthe high-temperature metallic phase, and the low temperature insulating41.1. Crystal structureFigure 1.4: Ni-O-Ni bond angle in the rare-earth nickelates is a measure ofthe perosvkite distortion, equivalent to the tolerance factor. Figure fromRef. [4].51.2. Structural Transitionphase. (The only element in the series that stays metallic at all tempera-tures is the La based compound.) Before we discuss this metal-to-insulatortransition (MIT), which is still at times a source of controversy in the liter-ature, we first consider the structural changes that occur in the nickelatesacross that same phase boundary, as the former is widely thought to beheavily influenced, if not outright triggered, by the latter.1.2 Structural TransitionAs the temperature is lowered, all nickelates (except LaNiO3) undergo asymmetry reducing lattice distortion, which changes their point group classfrom the orthorhombic Pbnm to monoclinic P21/n [5, 9, 35–37]. In addition,there is a slight volume reduction [23]. As was just mentioned, the struc-tural deformations may be rather complex, combining twists and tilts aboutvarious crystal axes [3, 4, 23, 27, 38]: however, it is of prime importanceto determine exactly the kind of distortion type (or symmetry class) thatdominates the overall structural evolution, as different types would supportqualitatively different explanations for the MIT (see Secs. 1.3,1.4). This canbe accomplished through the use of symmetry-mode based decompositionanalysis (in this section we follow Refs. [5, 39]; for a general review of themethod in crystallography, see Ref. [40]). In short, the atomic positions ina distorted crystal, labelled rdisti , are written in terms of the coordinates ofthe undistorted, high-symmetry crystal r0i and a linear combination of par-ticular distortion modes with unit displacement vectors dim and distortionmode amplitudes Amrdisti = r0i +∑mAmdim. (1.2)The modes dim are defined based on symmetry considerations, whichmakes them independent from each other and allows for the decomposition.This approach helps process experimental data, which often comes as variousNi-O bond lengths, Ni-O-Ni angles, and other crystallographic parameters,in a transparent way and quickly identify the dominant distortion modesfrom change in the amplitudes Am across a phase line. Balachandran andRondinelli [5] have carried out a careful symmetry-mode based decomposi-tion analysis for the nickelates: they noted that, while there is significantdistortion of types R+4 (an out-of-phase rotation, see Fig. 1.5f), M+3 (an in-phase twist, see Fig. 1.5c), and X+5 (in-phase tilting, see Fig. 1.5h) already61.3. Metal-Insulator Transition: Introductionpresent in the nickelates at high temperatures, the biggest change acrossthe MIT phase line is in the sudden appearance from zero of the breathing-mode type distortion R+1 (already depicted above in Fig. 1.2). Meanwhilethe various Jahn-Teller like modes (see Sec. 1.4.2 for a discussion of whythey were thought to be important) are negligibly small.Some of the earlier structural analysis [41] seemed to suggest that Jahn-Teller distortions were present, but the conclusions were drawn based ona more rudimentary technique – that of calculating the distortion param-eter ∆d = (1/6)∑n=1,...6[(dn − 〈d〉)/〈d〉2]2, where 〈d〉 is the average Ni-Obond length and dn are the various Ni-O bond lengths measured in neutron-scattering studies. Moreover, in the analysis the change of this parameterwith temperature for a fixed R member of the nickelate series was neverconsidered. Rather, the distortion parameter was calculated for the nick-elate series and correlated with the transition temperatures TMIT : but itis the changes across the MIT phase line that are most significant, anda cross-sectional view of the nickelate series cannot conclusively establishwhat distortions are relevant across it. In light of the more sophisticatedtreatment via eq. 1.2, it is clear that the earlier approach (a) could noteasily differentiate between the various contributions to the distortions inthe Ni-O bonds and thus risked mis-attributing the distortion to the wrongmode, and (b) it is the changes across the MIT phase line that determinethe dominant distortion mode, which were not considered.The conclusion to take away from this discussion is that while the nick-elates may be octahedrally twisted and otherwise distorted in a number ofways, the primary change across the MIT phase boundary results from arocksalt, breathing-mode lattice distortion of type R+1 , wherein alternat-ing Ni octahedra expand and collapse isotropically and roughly by an equalmeasure, and that furthermore there are no Jahn-Teller distortions observed(once again, see Sec. 1.4.2 to see why the Jahn-Teller modes generated somuch controversy).1.3 Metal-Insulator Transition: IntroductionThe “smoking gun” kind of evidence for an MIT in the nickelates is thestandard measurement of electrical resistance: this first-order phase tran-sition is clearly indicated by the resistance curves, which shoot straightup once critical temperature is reached (see the measurements for severalof the compounds in Fig. 1.6, from [6]). While the presence of the MITis firmly established by measurements from several groups over the past71.3. Metal-Insulator Transition: IntroductionFigure 1.5: Various distortion modes available to the nickelates. Figure fromRef. [5].81.3. Metal-Insulator Transition: IntroductionFigure 1.6: Electrical resistivity of PrNiO3 as a function of temperature:the metal-insulator transition is clearly visible. Data and figure from Ref.[6].decades [35, 42, 43], its nature and origin are not immediately obvious fromthese measurements alone. Before we launch into a discuss of the variouspossibilities, let us consider the basic electronic structure of the nickelates.1.3.1 Formal Valence rulesWhen it comes to transition metal oxides and determining the formal3 ox-idation state4 of the various components, for instance in a perovskite withthe chemical formula ABX3, it is customary to do the following. Given thatthe unit cell ought to be electrically neutral – otherwise the material wouldnot be electrically stable – the valences of all the components must add upto zero. Furthermore, due to extremely high electronegativity of chalcogens(O, S, Se) and halides (F, Cl, Br, I), assume their outer shells are closed:they have stolen enough electrons to satisfy their outer shells. One can makea similar assumption about the alkali metals like Li, alkaline earth elementsin the Be column, and those in the column under Sc: their low electroneg-ativity leads them to give up their electrons in most situations. Once these3i.e. nominal, not physically true but a useful starting point.4that is, how many electrons an outer electronic subshell is gaining/losing.91.3. Metal-Insulator Transition: Introductionare fixed, the remaining valences are determined by setting the sum total tozero: this condition can typically be satisfied because the transition metalsoften have a variety of diverse stable oxidation states (e.g. for vanadiumthe oxidation state can range from +2 all the way to +5). In the case ofthe rare-earth nickelates RNiO3 we can assign formal valences as follows: R,being the in Sc column, has +3. In a formula unit of RNiO3 there is 1 Niand 3 O: if we insist on O−2, and R is +3, that forces Ni+3, or, in terms ofshell structure, Ni = [Ar] 3d7.Knowing the number of electrons residing in the outer shell of the Niions, one can, in order to describe the nickelate behaviour, attempt to con-struct a Hubbard-like model for the nickelates. The simplest approach isto forget about the rare-earth ions and the O bands, assuming they are faraway from the Fermi level and the Ni bands, and to focus purely on the Niconduction electrons. A Ni atom in a spherically symmetric vacuum pos-sesses the standard hydrogen-like electronic orbitals: in fact, according tothe hydrogen atom solution, all the orbitals in the Ni n = 3 subshell (wheren is the principal quantum number5) have the same energy (assuming theelectron-electron interactions are neglected!) and thus would all have to beincluded in a Hubbard model. However, the Ni atoms are not really in avacuum, nor are they in a spherically-symmetric environment: they are sur-rounded by a periodic array of O and R and other Ni. In particular, eachNi atom is encased in an octahedral cage of oxygen atoms: this breakingof full rotational symmetry down to the octahedral group lifts the orbitaldegeneracy and significantly simplifies the model, reducing the number ofstates that need to be included. The effects of the surrounding O on theelectronic structure of Ni are obtained within crystal field theory.1.3.2 Crystal field splitting of RNiO3There are excellent books available detailing the crystal field calculations forvarious geometries and point group symmetries, as well as the group theoryinvolved [7, 44]. Here we will only provide a qualitative description of thecalculation, omitting some of the calculation details.Consider the total electric potential V (r) on the electron in the outershell of an ion of interest – say, Ni (positioned at R0). Treating all other5which is a consequence of conservation of energy in the hydrogen atom problem.Other quantum numbers to remember are the angular quantum number l from rotationalsymmetry (l = 0 is the s orbital, l = 1 the p orbital, l = 3 the d orbital, l = 4 the forbital, ...), the angular projection number mz, and the spin projection quantum numbersz.101.3. Metal-Insulator Transition: Introductionions (oxygens) as point electric charges of magnitude qa for the moment, onecan writeV (r) =∑aqa|Ra − r| =q0|R0 − r| +∑anearestneighboursqa|Ra − r| +∑a nonnearestneighboursqa|Ra − r| = (1.3)= V0(r) + Vlig(r) + Vcr/lig(r). (1.4)The first term is simply the spherically symmetric potential from theelectron’s own core ion. The second is the so-called ligand field, producedby the ligands or atoms that are immediately surrounding (bonding to) theoriginal Ni ion. The second and the last term together comprise the crystalfield.It is reasonable to assume that the nearest neighbours exert the strongestinfluence on the electron of the ion at R0. Specializing to the case of a Ni ionin an octahedral cage of O, we have point charges qO at positions (±a, 0, 0),(0,±a, 0) and (0, 0,±a) (the assumption of perfect octahedra with all equalprimary axes can be relaxed, but it does not change the results qualitatively).The ligand field on the electron isVlig(r) =qO√(x− a)2 + y2 + z2 +qO√(x+ a)2 + y2 + z2++qO√x2 + (y − a)2 + z2 +qO√x2 + (y + a)2 + z2++qO√x2 + y2 + (z − a)2 +qO√x2 + y2 + (z + a)2. (1.5)Taylor expanding about r = 0 to first order givesVlig(r) ≈ 354qOa5(x4 + y4 + z4 − 35r4)= D(x4 + y4 + z4 − 35r4). (1.6)This term manifestly breaks spherical symmetry by virtue of the x4, y4, z4factors: the hydrogen atom problem thus acquires a perturbation to thespherical Coulomb potential 1/r and the large degeneracy between the dif-ferent subshell orbitals is lifted. One clever way to figure out how the degen-eracy is lifted is using representation theory: however, the exact derivationis lengthy and interested readers are directed to Ref. [7].To calculate the actual energy splittings and the make-up of the statesin terms of the original atomic orbitals |n, l,mz, sz〉 it is necessary to do de-generate perturbation theory on the ligand potential Vlig(r) and diagonalize111.3. Metal-Insulator Transition: IntroductionFigure 1.7: The eg and t2g manifolds resulting from crystal field splitting.Figure from Ref. [7].the resulting crystal-field matrix to obtain the splittings. Without going inthe details of the diagonalization (see Ref. [7]), we simply quote the results.After taking into account the crystal field, the degenerate d orbitals of the Niion split into two manifolds: a higher energy manifold called eg (for group-theoretic reasons), which encompasses the orbitals d3z2−r2 , dx2−y2 ; and t2gwith the orbitals dxy, dyz, dzx. Intuitively, the eg orbitals end up being higherin energy because they are oriented towards the negatively charged O ions,whereas the t2g are largely in interstitial space (see Fig. 1.7).Combining this knowledge of the new d orbitals with the previous insightabout the Ni valency of 3d7, we conclude that the lower-lying t2g triplet iscompletely filled, while the eg doublet remains degenerate, with a singleelectron splitting its time between the two orbitals – and two possible spinprojections sz = ±1/2 – resulting in quarter-filling at zero temperature.This is often written as d7 = t62ge1g. Thus the tight-binding model, if con-structed, ought to include two bands, |z〉 ≡ d3z2−r2 , |z¯〉 ≡ dx2−y2 .While the eg − t2g splitting is an adequate starting place for under-standing the impact of the ligand O on the electronic properties of Ni, thesimple point charge approximation leaves out important band effects. Totake those into account, it is necessary to take stock of the charge-transferenergy gap between the Ni and O bands, which we discuss next within the121.3. Metal-Insulator Transition: IntroductionZaanen-Sawatzky-Allen framework [45].1.3.3 Zaanen-Sawatzky-Allen schemeIn this discussion we draw on an excellent reference by George Sawatzky andRobert Green [8]. Standard band theory, motivated by the tight-bindingmodel, assigns only one criterion for judging whether a material is a metalor an insulator: the parity of the number of electrons per unit cell. If thenumber is odd, then the outermost band is only partially occupied and thematerial is a metal. Conversely, an even number of electrons means the bandis fully occupied, and in the presence of a finite band gap between the occu-pied (valence) band and the next higher, unoccupied (conduction) band, thematerial is insulating (this is well-described in many solid-state textbooks).The success of this model was astounding, but brief: it was quickly realizedthat there were a variety of materials, most notably transition-metal oxides,which were insulating despite formally having an odd number of electronsin the unit cell.To address this discrepancy, Mott and Hubbard pointed out that inter-actions between electrons must be taken into account [46, 47]. The Mott-Hubbard model, which incorporates the electron-electron interactions onlybetween electrons at a given site, was already a significant improvement: itwas able to show exactly how electron-electron interactions were forcing thelocalization of electrons, thus rendering the materials insulating. Withoutinteractions the Fermi level would sit in the middle of the outer electronband: the system would have many states which could potentially respondto an applied current (which is what conduction is). Instead, the electrons’mutual Coulomb repulsion (whose strength is characterized in the model bythe value U) opened a gap in the density of states at the Fermi level, split-ting the band in two, and eliminating any states that were once accessiblefor conduction, effectively confining the electrons to their “home cell” (seeFig. 1.8 for a useful graphic).However, soon this was also revealed not to be the whole story. Transi-tion metal oxides nominally incorporate both transition metal 3d bands andoxygen 2p bands in their bandstructure. Typically it is assumed that theoxygen bands are valence and are way below the transition-metallic bands:however, in many oxides that is not the case. In fact the bandstructurecan be such that the lowest energy excitations actually involve the an elec-tron hopping from the oxygen band to the transition-metal band, or yetmore exotic possibilities. This was the chief insight of Zaanen, Sawatzkyand Allen [45], who in 1985 classified the possibilities for such compounds131.3. Metal-Insulator Transition: IntroductionFigure 1.8: A graphic depicting the Zaanen-Sawatzky-Allen classificationscheme, using featureless densities of electronic states without hybridization.Figure from Ref. [8].in terms of the Hubbard U and the charge-transfer energy ∆: this latterenergy scale characterizes the energy cost of the electron hopping from thetransition-metal atom to the ligand oxygen.The different possibilities coming from the interrelationship of U and∆ are depicted in Fig. 1.8. They range from a full-on Mott-Hubbard sit-uation, where the oxygen band indeed plays no role; to a positive chargetransfer case where the lowest-order excitations are of dnL0 → dn+1L1 type(L indicates a hole state on the ligand oxygen); to a mixed valence statethat is somewhat difficult to understand and we will not discuss much here;and, finally, negative charge transfer insulators, which is what the nickelatesRNiO3 are believed to be in the ZSA scheme. The idea behind the negativecharge transfer insulator is that ∆ < 0 and it is thus energetically favourablefor the oxygen 2p electron to migrate toward the Ni : this situation is some-times called self-doping. The structure of the lowest order excitations inthe various scenarios is depicted in Fig. 1.9.The crucial takeaway from the ZSA analysis is that, unlike in the formalvalence counting in Sec. 1.3.1, the correct starting point for constructing aHubbard model for the Ni should be the d8L1 state, not d7. This has seri-ous implications for the picture of the charge order and the metal-insulatortransition, as we will see shortly.141.4. Metal-Insulator Transition: MechanismsFigure 1.9: The lowest order excitations within the Zaanen-Sawatzky-Allenclassification scheme. Figure from Ref. [8].1.4 Metal-Insulator Transition: MechanismsArmed with the understanding of the basic facts about the Ni valency, or-bitals and bandstructure, we can now start addressing the question of themetal-insulator transition and the various mechanisms that have been pro-posed.1.4.1 Mechanism: Charge-transfer gapIn the simplest and one of the earliest pictures, as the temperature is lowered,a charge-transfer gap appears in the band structure of the nickelates, leadingto insulating behaviour [16]. However, this does not readily account in anyway for the structural lattice distortions which are present in this materialand thus is not considered an adequate theory at this point.1.4.2 Mechanism: Jahn-Teller effectA good rule of thumb when it comes to electronic systems is that “natureabhors an orbital degeneracy”[48]. Whenever there is a possibility of liftinga degeneracy, it is generally done, unless the degeneracy is protected by sym-metry. So it should come as no surprise that the remaining degenerate egdoublet could have a natural way of resolving its degeneracy: in the case of aperovskite lattice, one such possibility is known as the Jahn-Teller effect. Inshort, it comes from the interaction between the electronic orbitals and lat-tice vibrational modes. The interaction produces a “cooperative distortion”of the lattice, elongating the O octahedra preferentially along one axis. This151.4. Metal-Insulator Transition: MechanismsFigure 1.10: The unit cell of KCuF3 and the associated JahnTeller distortionmodes. Figure from Ref. [7].splits the eg doublet and produces a kind of orbital order. A derivation ofthe Jahn-Teller effect from a Theory-of-Everything Hamitlonian is beyondthe scope of this thesis: an excellent introduction is given in Ref. [7].A typical example of a Jahn-Teller system is KCuF3, another perovskite:see the two possible distortion modes in Fig. 1.10. Intuitively, the Q2 modebrings the negatively charged O closer to the branches of dx2−y2 and so raisesits energy, while pushing other O away from d3z2−r2 and lowering its energy– and in so doing, lifts the degeneracy.In light of this well-studied phenomenon, and given that with the openouter subshell structure t62ge1g the nickelates are expected to be Jahn-Telleractive [38] (due to the unpaired electron in a degenerate pair of orbitals), itwas reasonable to suggest that a Jahn-Teller distortion is what resolves theelectronic degeneracy and results in insulating behaviour below a certaintemperature [31, 49]. At high temperatures, itinerant6 electrons possesslarge kinetic energy (the bandwidth W is large) that dominates the Jahn-Teller distortion and suppresses it, so the (electric) conduction is largelyunaffected.Many in the literature assert that, despite decades of research, no dis-cernible Jahn-Teller like (orthorhombic) distortion was found [27, 35, 39,42, 50, 51] – even though they might cite papers that explicitly say thatJahn-Teller distortions are present (such as in [41]). Moreover, some re-searchers suggest that neutron diffraction studies show that the change ofthe point group across the MIT phase line would not even be consistentwith a Jahn-Teller distortion: the change is such that the unit cell appearsto double, resulting in two inequivalent Ni sites – whereas a Jahn-Teller dis-6i.e. not localized – typical Bloch states in metals that are distributed over almost theentire sample.161.4. Metal-Insulator Transition: Mechanismstortion would operate in the same way within all octahedral O cages [52].However, others are suggesting that a Jahn-Teller distortion is indeed ob-served [4, 41, 49, 53]: they point to an increase in the distortion parameter∆d = (1/6)∑n=1,...6[(dn−〈d〉)/〈d〉2]2 as one follows the nickelate series fromLa to Lu, at room temperature. However, as we already discussed abovein Sec. 1.2 on the structural transition in the nickelates, a variety of dif-ferent types of distortion could potentially occur in the nickelates. Any ofthem could affect this very general parameter ∆d, without necessarily serv-ing as evidence for Jahn-Teller distortions. The more recent studies alreadymentioned that apply a sophisticated symmetry based mode decompositionanalysis for the distortions seem to concur that a breathing-mode distor-tion, rather than a Jahn-Teller one, is the main contribution to the latticedistortion of the nickelates across the MIT. Ref. [54] explains nicely how thesame data (e.g. Ni-O bonds lengths) could have been interpreted in favourof different theories by different researchers.In summary, with some degree or certainty one could conclude that Jahn-Teller distortions are not responsible for charge order, as there are no Jahn-Teller distortions observed: although a possibility of a dynamic Jahn-Tellereffect, which would not be readily picked up by the various X-ray and neu-tron diffraction experiments (as they time-average the measurements), hasnot yet been conclusively ruled out [41].1.4.3 Mechanism: Charge disproportionationAn alternative proposed microscopic explanation accounting for the MIT isthe effect of charge disproportionation (CD). Neighbouring Ni ions, nom-inally each being of d7 valency, spontaneously form pairs with alternatingvalence d7−δd7+δ (equivalently this is often written Ni3− δNi3+δ), where δis the magnitude of the charge disproportionation [4, 35–37, 42, 48, 50, 55].This happens to lift the eg degeneracy without requiring Jahn-Teller distor-tion, and would preferentially couple to breathing-mode type distortions7,as opposed to Jahn-Teller like distortions [54]. The breathing mode distor-tion makes sense in light of the changes of ionic radii with the introductionof a non-zero δ. This charge density wave then impedes conduction in theusual Mott-Hubbard insulator sense, opening a gap at the Fermi surface.However, there are some issues with this approach. Firstly, a charge densitywave can only occur for Hund’s coupling J8 large enough to compensate7i.e. ones where alternating O octahedra expand and contract.8Hund’s coupling J is the strength of the exchange interaction in an electronic system:typically smaller than U , it prefers to arrange the electrons into a state that maximizes171.4. Metal-Insulator Transition: Mechanismsfor the prohibitive cost of double occupancy U [48]: this indeed seems tobe the case in the nickelates [32]. But that kind of ratio between U and Jwould make it energetically favourable to transfer the electron away from Oand onto Ni (alternatively, a hole from Ni to O), thus making a 3d7 config-uration an inappropriate starting point (as discussed in terms of the ZSAframework in Sec. 1.3.3). Second, on more experimental grounds, it appearsto directly contradict resonant soft-x-ray diffraction measurements probingthe Ni d orbital occupation [56].1.4.4 Mechanism: Negative charge-transfer theoryAn alternative view, motivated by the Zaanen-Sawatzky-Allen framework[45], holds that the nickelates are a negative-charge transfer insulator: thatis, the energy cost ∆ (the charge-transfer energy) of occupying a Ni d or-bital versus an oxygen p orbital is such that ∆ < U and in addition ∆ < 0[18–22]. This forces a hole into the oxygen band (effectively self-doping thematerial), resulting in a d8L Ni valency state (L denotes a ligand (oxy-gen) hole). Subsequently, electron-phonon coupling results in a breathing-mode distortion of the lattice, which couples to a charge order of the form(d8L)i(d8L)j → (d8L2)S=0(d8)S=1. Since each hole on the oxygens is sharedequally by both Ni that are bonded to it, there is no actual overall move-ment of charge: both O cages have the same average hole concentration[32]. But when two holes are present on a Ni (collapsed) octahedron, thephase of their wavefunctions is modified to acquire the symmetry of the Nieg orbitals. The pair of holes then locks into a singlet with the Ni electrons,resulting in an effectively 3d6 Ni configuration, which is not susceptible toJahn-Teller distortions. Meanwhile the expanded Ni octahedron, left withno holes, orders as a triplet with S = 1, which allows the material to ordermagnetically (and acting in effect just like charge disproportionation wouldduring most experimental probes).These different mechanisms of how the MIT actually proceeds are notmerely questions of interpretation of data: different modeling approaches,treating the electronic structure to varying degrees of sophistication, have tobe implemented depending on which of the mechanisms one believes to becorrect. We use many of the considerations laid out here (Ni valence count-ing, crystal-field splitting, the Zaanen-Sawatzky-Allen framework, absence(maybe?) of Jahn-Teller distortions and orbital order) when we constructour model in Ch. 2.the local spin – justifying that old adage from chemistry, Hund’s rule.181.5. Magnetic Transition: Introduction1.5 Magnetic Transition: IntroductionThere are a few key experimental signatures that signal the onset of mag-netic ordering, such as a divergence in magnetic susceptibility, or the sud-den appearance of a nonzero atomic magnetic moment (picked up usuallythrough neutron scattering). The first indications of magnetic behaviour inthe nickelates were observed by the team Demazeau et al. that originallysynthesized the material series in the early 70s [14]. It was not until thelate 80s and early 90s that magnetic phenomena in the nickelates started toattract more widespread attention [38, 57–59]. An excellent plot of the Nimagnetic moment with temperature is depicted in Fig. 1.11: the sponta-neous transition to a magnetically ordered state with temperature is clearlyvisible.After a great deal of work it emerged that all of the members of the Nifamily possess nontrivial, magnetically ordered ground states at low enoughtemperatures (except for LaNiO3 which is assumed to stay paramagneticdown to the lowest temperatures – although there is very recent work thatis calling this into question [24]). Two of the members of the family, Ndand Pr, experience the magnetic transition simultaneously with the MIT:the transition is first order [3]. Meanwhile the rest of the family up toLu undergoes a second-order transition, with the transition temperature TNtrending upward as one approaches Lu. This is encapsulated in the canonicalphase diagram showed in an earlier section 1.1.Much work since those days has been devoted to characterizing this low-temperature magnetically ordered state, primarily though neutron diffrac-tion measurements [9, 28, 38, 49, 56, 60, 61]. Due to the difficulty of produc-ing large enough single crystal samples, the alignment of the local magneticmoments (both from the Ni and the rare-earth sublattices) remains a mys-tery: however, it has been conclusively established that whatever the align-ment may be, it obeys the ordering wavevector Qm =2pia [14 ,14 ,14 ] (2pia [12 , 0,12 ]in the perovskite lattice notation). This effectively means that the unit cellof the nickelate ought to include 4 separate, linearly coordinated Ni octahe-dral sites: only after the fourth site does the magnetic order repeat itself.The emergence of this ordering peak can be clearly seen in Fig. 1.12, whereit is marked in the perovskite notation.191.5. Magnetic Transition: IntroductionFigure 1.11: The Ni magnetic moment in PrNiO3 as a function of tempera-ture, with the onset of magnetism clearly visible. Figure from Ref. [4].201.5. Magnetic Transition: IntroductionFigure 1.12: The appearance of a new ordering peak in neutron diffractiondata for DyNiO3 as the temperature is lowered past the magnetic transitiontemperature TN. Figure from Ref. [9].211.6. Magnetic Transition: Mechanisms1.6 Magnetic Transition: MechanismsAt present it is not known what the mechanism of magnetization is in thenickelates: nor is it known what the alignment of local magnetic momentsis. This state of affairs fuels much of the current interest in the nickelates:it ranges from attempts to grow larger crystals9 [24] to be able to get atthe order experimentally – say through a scanning-tunneling microscopy(STM) experiment or inelastic neutron diffraction – to theoretical attemptsto model the nickelates and thus make a prediction for the alignment usingeither effective band models, like in this thesis or as Lee et al. for in Ref. [26],or ab initio methods [1, 2, 19, 25, 39]. There are a variety of propositions asto how the magnetic order arises, which are separate from the propositionsregarding what the order actually is.One idea for the appearance of the order is that it is due to an orbitalsuperlattice, triggered, in turn, by the orthorhombic distortions of the nick-elates away from perfect perovskite lattice structure [62]. The presence ofthis distortion results, potentially, in (very weak) orbital order: but if theenergy splitting introduced by it is comparable to the exchange couplingbetween the different orbitals on the same Ni site, a new energy-favourableconfiguration becomes available wherein two sublattices with different or-bital occupancy emerge. Then, given a Ni site, if its neighbour has thesame orbital occupancy, they are antiferromagnetically coupled through theoxygen 2p bands (superexchange); if it has the opposite orbital occupancy,then it is preferential for the spins to align ferromagnetically. This scenariowould reproduce the experimentally observable wavevector Qm: locally, itwould result in the magnetic order of type ↑↑↓↓.Another possibility is, once again, cooperative Jahn-Teller distortions.The mechanism here would be very similar: orbital order would couple withmagnetic exchange interactions to spontaneously generate two inequivalentsublattices with AF coupling between same-orbital occupancy neighboursand FM coupling between opposite-orbital ones. The crucial distinctionhere is the origin of the orbital order: it would be assumed due to theelongation of the Ni octahedra. However, we again reiterate here that noJahn-Teller like distortion has been observed for the nickelates.It would appear that one could disentangle the order generated by theorbital versus the spin degrees of freedom by studying the neutron scatteringpeaks of Ni family members with TMIT > TN: one might hope that the9The problem is that a variety of different orders, such as ⇑ 0 ⇓ 0 or ⇑⇑⇓⇓ cannot bedistinguished from one another during a neutron diffraction study on a powder.221.6. Magnetic Transition: Mechanismsprecursors of orbital order exist already above TN. However, a study withSm placed very stringent bounds on the magnitude of these peaks, if theyexist (< 10−4 of the largest peak typically observed in the study).Clearly, much more work needs to be done to fully understand the char-acter of the magnetic order and the magnetic transition at low temperature:that is the main purpose of this thesis.23Chapter 2Effective Two-Band Modelfor the Rare-Earth NickelatesAll models are wrong, but someare useful.George Box (1987)Now that we have extensively reviewed the status of experimental andtheoretical work to date, especially when it comes to the metal-insulatortransition and the magnetic order, the stage is set for constructing a phe-nomenological, microscopic tight-binding model for the nickelates. The con-siderations in the previous chapter should make clear as to what aspectsof the physics ought to be included, and what can be neglected or treatedphenomenologically as a renormalization to the model parameters.Given that the octahedral twists do not appear to change drasticallyacross the MIT or the magnetic phase line and thus apparently do notcouple much to the electronic/magnetic degrees of freedom, for simplicitywe take RNiO3 as having a perfect (pseudo)cubic perovskite lattice, withlattice constant a. In light of the crystal field and formal valence discussions,we will concentrate on two eg “effective” Ni orbitals for the tight-bindingmodel, |z〉 ≡ d3z2−r2 , |z〉 ≡ dx2−y2 .We call the orbitals effective, because we believe that in reality the point-charge model for the Ni orbital crystal field splitting is incomplete: it omitsany mention of the O 2p bands, or of the holes on the O due to negativecharge transfer. If one believes in the charge disproportionation mechanismof the MIT, the O bands are assumed inert and thus the simplistic crystalfield treatment described above is essentially accurate: so the eg bands inthe model would really be the eg bands of the Ni.However, that is not quite so in the negative charge transfer picture.When considering a single Ni octahedron, the holes on the O are found to bethose linear combinations of the pσ O orbitals that exhibit eg character [32].This makes sense, as eg-like O orbital combinations are the ones with max-24Chapter 2. Effective Two-Band Model for the Rare-Earth Nickelatesimum overlap with the Ni eg doublet: they are most likely to hybridize anddonate an electron to the Ni, leaving behind the ligand hole and producingthe 3d8L state in the first place. Hence the ligand hole states can be thoughtof as “effective” eg orbitals: the analogy does not fully bear out because suchorbitals would not actually be orthogonal between nearest-neighbour sitesand in fact need to be built up in the Wannier states manner from a varietyof more extended Ni and O states. This could potentially strongly renor-malize the hopping and interaction parameters as compared to their barevalues. More than that, some recent results suggest that even the make-upof the bands – i.e. what eg, t2g and other bands are active in the model –could dynamically respond to things like site occupation and crystal latticedistortion [63]. We deal with this by letting the model parameters vary andseeing what states obtain within the larger parameter space.Now that we settled what the lattice is and what orbitals to include,we can readily write down all the ingredients of a microscopic tight-bindingmodel. We will work at zero temperature for simplicity, as we are inter-ested mostly in the magnetic order that obtains at very low temperatures inthe nickelates. At the same time, we can represent sweeping the nickelateseries by adjustments in the model parameters. Roughly speaking, whileswapping out the rare-earth ion affects the bond angles in the crystal andso significantly modulates the bandwidth (which is determined by variousnearest-neighbour hopping amplitudes Ti in our model), it does little tochange the on-site Coulomb repulsion. Thus changing the hopping rate canrepresent the effect of continuously sweeping through the entire nickelateseries.As the importance of the lattice has already been established [31, 32, 39,51], we introduce the lattice semiclassically and couple it to the tight-bindingmodel via an on-site coupling, with the goal of demonstrating explicitly theforcing of charge disproportionation by the lattice and to investigate the im-pact of the lattice on the ground state magnetic order. Finally, we treat theelectron-electron interactions by employing the standard spherically sym-metric Kanamori Hamiltonian for the case of two bands. Our Hamiltonianwill thus have the general formHˆ = Tˆ + Hˆe−e + Hˆlatt + Hˆe−latt. (2.1)We use electronic creation and annihilation operators d†iaσ, diaσ to createand destroy electrons in an orbital a with spin σ on site i. In what follows,we build up the Hamiltonian 2.1 piece by piece.252.1. HoppingFigure 2.1: A schematic diagram of the various hoppings included inthe nickelates two-band model. The green arrow corresponds to the firstneighrest neighbour hopping; the yellow to second nearest neighbour hop-ping; third nearest neighbour hopping, represented by the red arrow, is notincluded due to minimal orbital overlap; and fourth nearest neighbour hop-ping is shown by the purple arrow.2.1 HoppingWe know from neutron scattering measurements to expect a linearly coor-dinated 4-site magnetic order in the nickelates It is necessary then to in-clude first nearest neighbour (with amplitude t1), second-nearest neighbour(t2) and fourth-nearest neighbour (t4) hopping, where the “nearness” of theneighbours is indicated in Fig. 2.1 and is purely spatial. The hopping partof the Hamiltonian readsTˆ = Tˆ1 + Tˆ2 + Tˆ4. (2.2)Let us work through these sequentially. We will explicitly write out262.1. Hoppingall the hoppings in the various directions: we will use notation like di+x,xσwhere the first x is a vector translating us one lattice constant length alongx to the next Ni site (with respect to the original site, for which we use acondensed 3-dimensional label i), and the second x is the orbital∣∣3x2 − r2〉,which is the equivalent of |z〉 = ∣∣3z2 − r2〉 “polarized” along the x direction,which allows hopping along x. The reason we have states like∣∣3x2 − r2〉 isbecause while we have only two true orbitals |z〉 and |z〉, the choice of thez axis is arbitrary: we could have equally well have chosen x as the z axisand still would have had an eg member “polarized” along that axis. We canalways write all these orbitals in terms of |z〉 = ∣∣3z2 − r2〉 , |z¯〉 = ∣∣x2 − y2〉by remembering what the decomposition of these orbitals is in terms of (real)spherical harmonics|x〉 ≡ ∣∣3x2 − r2〉 = −12|z〉+√32|z¯〉 , |y〉 ≡ ∣∣3y2 − r2〉 = −12|z〉 −√32|z¯〉 .(2.3)2.1.1 1st nearest neighbour hoppingWith these orbitals, the nearest neighbour hopping contribution isTˆ1 = −t1∑iσ[d†izσ (di+z,zσ + di−z,zσ) +d†ixσ (di+x,zσ + di−x,zσ) ++ d†iyσ (di+y,zσ + di−y,zσ)]. (2.4)Of course, nearest neighbour hopping is not possible for the |z〉 orbital alongthe z direction, because these orbitals have no overlap with the ligand O inthat direction. The hopping is possible along the x and y directions, but wechoose to consider the hopping through |z〉 (and |z〉-like orbitals polarizedalong different axes) as dominant and so neglect this fact.This quadratic Hamiltonian can be readily brought to almost diagonalform (save for the orbital subspace) by Fourier transforming to k basis, usingthe prescriptiond†kaσ =∑ieikRi√Nd†iaσ, (2.5)where k ranges over the full cubic Brillouin zone −pia < kη < pia , η = x, y, z.272.1. HoppingUsing this gives (the details of the calculation are in Appendix A):Tˆ1 = −t1∑kσ[t1zz(k)d†kzσdkzσ + t1z¯z¯(k)d†kzσdkzσ++ t1zz¯(k)(d†kzσdkz¯σ + d†kz¯σdkzσ)],(2.6)witht1zz(k) = 2 cos(kza) +12[cos(kxa) + cos(kya)],t1z¯z¯(k) =32[cos(kxa) + cos(kya)], (2.7)t1zz¯(k) = −√32[cos(kxa)− cos(kya)].2.1.2 4st nearest neighbour hoppingWith the Tˆ1 prescription in our pocket, we can immediately write down theTˆ4 operator. Since the 4th neighbour hopping is through the same orbitalsas the 1st neighbour hopping, the form of the operators is still the same, savefor the displacement vector being 2a · xˆ, yˆ, zˆ. Hence the result ought to bethe same as Tˆ1, except the cosines are now over twice the lattice constant,cos kxa→ cos 2kxa and so on. Explicitly,Tˆ4 = −t4∑kσ[t4zz(k)d†kzσdkzσ+t4z¯z¯(k)d†kzσdkzσ++ t4zz¯(k)(d†kzσdkz¯σ + d†kz¯σdkzσ)], (2.8)witht4zz(k) = 2 cos(2kza) +12[cos(2kxa) + cos(2kya)],t4z¯z¯(k) =32[cos(2kxa) + cos(2kya)], (2.9)t4zz¯(k) = −√32[cos(2kxa)− cos(2kya)].2.1.3 2nd nearest neighbour hoppingThis leaves the slightly more complicated 2nd neighbour hopping, where thehopping proceeds through mixed orbitals along multiple directions at once,282.1. Hoppingin one of the three planes xy, yz, zx. Only the origin and destination matterfor the hopping: we catalog the hopping contributions in terms of them.There are 12 possible destinations, 4 in each of the principal planes. Ingeneral the hopping term will look like thisTˆ2 = −t2∑iabσ,τˆ 6=µˆd†i+τˆ+µˆ,aσdibσ + h.c., (2.10)with τˆ , µˆ ranging over ±a · xˆ, yˆ, zˆ. However, notice that during such a hopthe electron always ends up in an orbital that is different from the one itstarted in. For instance, hopping xˆ+ yˆ forces the orbital to go either x→ yor y → x, depending on what orbital the electron started in (the orbitalwill change “virtually”, while the electron hops through a virtual site, butthe absence of any overlap between |z〉 and |z〉 does not allow orbitals tochange on real sites). Let us adopt the convention that in the site indices,order matters: in other words, if we have an index i + xˆ + yˆ, that meanswe first hop along xˆ and then along yˆ. Then it is clear that the startingorbital index b will be the same as the first hop displacement τˆ – after all,if an electron starts in the y orbital, it can only (initially!) hop along y.This allows us to set b = τ and a = µ. After that, it is a matter of carefullyFourier transforming all the 12 hopping destinations. The final result is (thedetails are in Appendix A)Tˆ2 = −2t2∑kσ[t2zz(k)d†kzσdkzσ+t2z¯z¯(k)d†kz¯σdkz¯σ++ t2zz¯(k)(d†kzσdkz¯σ + d†kz¯σdkzσ)], (2.11)witht2zz(k) = 2 cos(kxa) cos(kya)− 2 cos(kza)(cos(kya) + cos(kxa)),t2zz¯(k) =√3 cos(kza)(cos(kxa)− cos(kya)), (2.12)t2z¯z¯(k) = −3 cos(kxa) cos(kya).All of the different hopping contributions can be combined to give theusual form for TˆTˆ =∑kabσtab(k)d†kaσdkbσ, (2.13)The coefficients can be found at the end of Appendix A.292.2. Lattice contributions2.2 Lattice contributionsLet Uj be a distortion parameter characterizing the structural distortion atsite j. At the semiclassical level, where in addition the phonon modes of thematerial are assumed to be incredibly heavy, the dynamical (momentum)contribution to the phonon energy can be neglected, and the lattice energyis primarily determined by the distortion parameter Uj . In accordance withexperimental data described in an earlier section, it appears entirely justifiedto treat the main contribution to the lattice energy as that coming from anisotropic breathing-mode distortion mode R+1 (for a depiction of the modesee again Fig. 1.5). Specifically, the Ni octahedra are assumed to expand-contract in a checkerboard-like (also called breathing-mode like, also calledrocksalt like) pattern, with the expansion/contraction isotropic and roughlyequal in magnitude (opposite in sign). Thus the distortion parameter Uj isa simple scalar, and its energy can be writtenHlatt =∑j(k2U2j +A4U4j). (2.14)where for greater generality the anharmonic terms U4j have been included.For convenience, re-scale the distortion magnitude to extract the lattice cou-pling energy scale by non-dimensionalizing the lattice distortion parameteruj = kUj/g: this results inHlatt = b∑j(12u2j +α4u4j), b = g2/2k, α = Ag2/2k3. (2.15)So long as the value of the distortion is known, the contribution to theenergy can be readily evaluated with this expression.2.3 Electron-electron interactionsAs there are assumed to be two degenerate bands in the eg manifold, thereare markedly more options for electron-electron processes spurred on by theCoulomb interaction than there are in a single-band Hubbard model. Inter-action models of this kind have been studied previously: the generalizationto two bands is called the Kanamori [64–66] (sometimes Slater-Kanamori302.4. Electron-lattice interaction[19]) Hamiltonian, with the general formHU = U∑iαnˆiα↑nˆiα↓ + U ′∑iσnˆizσnˆiz¯σ¯+ (2.16)+ (U ′ − J)∑iσnˆizσnˆiz¯σ − J∑iσd†izσdizσ¯d†iz¯σ¯diz¯σ+ (2.17)+ J∑iad†ia↑dia¯↑d†ia↓dia¯↓. (2.18)The derivation of this interaction is beyond the scope of this thesis: inter-ested readers may be directed to [67]. Here we merely outline an intuitiveunderstanding of its behaviour. The first two terms should be familiar: theyare merely the usual Hubbard-like density-density interaction terms. Thepresence of J signifies Hund’s exchange interaction that arises in the pres-ence of more than one orbital, reflecting the varying cost of occupying thetwo orbitals with various spin orientations. A useful heuristic is Hund’s rule,which states that in the absence of other energy splitting, Hund’s J prefersto order the electrons by maximizing the total spin on the ion: thus in thecase of two degenerate orbitals, the preferred configuration for two electronswould be the triplet, not the singlet. Hund’s rule is a heuristic: if it wasfollowed to the letter in this situation, it would actually “prefer” the config-uration with truly maximum spin, namely t52ge2g, which has three unpairedspin up electrons and thus spin S = 32 . The magnitude of J essentially signi-fies how “strongly” the rule is followed: thus in this paper we do not considervery large values of J , as that would invalidate the two-band premise andwe would have to include the full t2g manifold into the calculation.The last two terms are somewhat more exotic: they represent the si-multaneous spin-flip and the pair hopping processes. Notice that we willmake the spherically symmetric choice U ′ = U − 2J and neglect the sym-metry reduction due to crystal fields: this is standard procedure [60]. Inthis construction the values of U and J should still be thought of as phe-nomenological parameters, instead of the proper Coulomb atomic matrixelements.2.4 Electron-lattice interactionNow comes the crucial piece: an interaction between the electrons and thecrystal lattice distortions. The collapsed or expanded octahedra should af-fect an electron’s affinity to occupy a given site: coupling the electronic312.5. Final form of the Hamiltonian Hˆdensity to the on-site distortion, we obtain the electron-phonon interactiontermHˆe−latt = −g∑iUi (nˆi − 1) . (2.19)Again applying the prescription ui = kUi/g, we obtain the final formHˆe−latt = −2b∑iui(nˆi − 1). (2.20)Notice that in this picture the electron-lattice interactions are assumed thesame for the two orbitals |z〉 , |z¯〉: this is because we assume that the onlyactive lattice distortion mode is the breathing-mode, which interacts in thesame way with both the orbitals. We thus implicitly assume that none ofthe other modes (like Jahn-Teller distortion) are active. We feel justified inthis assumption given that there appears to be no experimental evidence ofJahn-Teller distortions, as discussed in the introduction in Ch. 1.2.5 Final form of the Hamiltonian HˆNow that the stage is set and the model is fully defined, we have the followingHamiltonian for the rare-earth nickelatesHˆ =∑kabσtab(k)d†kaσdkbσ + b∑j(12u2j +α4u4j)− 2b∑iui(nˆi − 1)++ U∑iαnˆiα↑nˆiα↓ + U ′∑iσnˆizσnˆiz¯σ¯ + (U′ − J)∑iσnˆizσnˆiz¯σ− (2.21)− J∑iσd†izσdizσ¯d†iz¯σ¯diz¯σ + γJ∑iad†ia↑dia¯↑d†ia↓dia¯↓.All in all there are 7 parameters characterizing this Hamiltonian: the threehoppings t1, t2, t4, buried inside the coefficients tab(k); the on-site Coulombrepulsion U and Hund’s exchange J ; the breathing-mode distortion energyb and the dimensionless anharmonicity α.Next, we perform a Hartree-Fock calculation on this Hamiltonian. Oneobvious difficulty we will face is that it is easiest to treat the hopping part inFourier space, but it is simplest to treat the electron-electron (and electron-lattice) interactions in real space. Fortunately, there is a trick that can helpus make use of both simplifications at once.32Chapter 3Hartree-Fock Approximationfor Rare-Earth NickelatesThe most important part ofdoing physics is the skill toneglect.Lev LandauWe attack the Hamiltonian in 2.1 with the unrestricted Hartree-Fockapproximation (i.e. the mean-field approach – in this thesis we will oftenuse the terms interchangeably). It is unrestricted because we are allowingthe spin occupancy to be independent of the orbital occupancy, and notdemanding that for every occupied single-particle state that its spin-reversedpartner is also occupied. The ground state wavefunction is assumed to split|Ψ〉 = |Ψe〉 ⊗ |Φlatt{u}〉 , (3.1)with a Slater determinant |Ψe〉 for the electronic part and a semiclassicaldistortion |Φlatt{u}〉 for the lattice, which is a function of all the distortions{uj}. The electronic part is complicated and will produce many mean-fieldparameters (e.g. charge disproportionation δ =∑aσ〈niaσ−ni−1,aσ〉, orbitalorder∑σ〈nizσ − niz¯σ〉, spin density 〈Si〉 and so on), while the lattice partonly depends on the distortions uj . To find the Hartree-Fock ground state,we will derive the self-consistency equations, obtained by minimizing thetotal energyE = 〈Ψ| Hˆ |Ψ〉 , (3.2)with respect to the mean-field parameters.3.1 Minimizing lattice contributionsFirst, we minimize the energy with respect to the lattice distortions ui. Inprinciple, a straightforward and foolproof approach is to compute the energy333.1. Minimizing lattice contributionsas a function of ui and do a simple (or not so simple) parameter sweep tofind the lowest energy (and thus determine the ground state). However,this turns out to be an arduous and time-consuming process: luckily itis simplified significantly with the aid of the Hellmann-Feynman theorem[68, 69], which allows us to pass the derivative inside the expectation valueduring minimization∂∂uj〈Φlatt(ui)| 〈Ψe|H |Ψe〉 |Φlatt(ui)〉 = (3.3)= 〈Φlatt(ui)| 〈Ψe| ∂H∂uj|Ψe〉 |Φlatt(ui)〉 , (3.4)despite the fact that |Ψe〉 |Φlatt(ui)〉 is generally not an eigenstate of Hˆ (asit needs to be for the usual proof [70] of the theorem to apply). In fact,there is a generalization of the theorem to any kind of variational state, notmerely an eigenstate [71]. Consider the energy functionalET [Ψ{u}, {u}] = 〈Ψ{u}| Hˆ{u} |Ψ{u}〉 .By definition, |Ψe〉 |Φlatt{u}〉 is a stationary point of the functional withrespect to a variation δΨ{u}, i.e. |Ψe〉 |Φlatt{u}〉 is the solution toδETδΨ{u} = 0.On the other hand, taking a total u derivative of ET ,dduET (Ψ{u}, {u}) = ∂ET∂{u} +δETδΨ{u}∂Ψ{u}∂{u} =∂ET∂{u} ,we see that the Hellman-Feynman result still holds, thanks to the station-arity of the variational state. (In our calculation, we also checked this nu-merically, by comparing the results of Hellmann-Feynman calculation withexplicit minimization with respect to u.)Utilizing the Hellmann-Feynman theorem, we minimize with respect toui∑i〈Ψe| 〈Φ| ∂∂uj2b(12u2i +a4u4i − ui (nˆi − 1))|Φ〉 |Ψe〉 == 2b[uj + au3j − 〈nˆj〉Ψe − 1] = 0,343.1. Minimizing lattice contributionswhere we defined 〈Ψe| Oˆ |Ψe〉 ≡ 〈Oˆ〉Ψe . This gives us our first self-consistencyequationuj + au3j + 1 = 〈nj〉Ψe . (3.5)Careful neutron scattering measurements from multiple groups suggestthat the lattice distortions are isotropic, and alternate between the octahe-dra along all three crystallographic axes. Moreover, the collapse magnitudeu of one octahedra is exactly equal to the expansion of another. Thereforewe adopt the ansatz uj = ueiQc·Rj , where Qc = (2pi/a)(12 ,12 ,12) is the order-ing wavevector corresponding to a breathing-mode distortion. Let us alsochoose to represent the charge disproportionation behaviour with a param-eter δ, so that 〈nˆj〉 = 1 + δeiQc·Rj (charge disproportionation follows thesame ordering wavevector as the lattice distortion — this will be discussedin much more detail in Sec. 3.3). The equation above then reduces to adepressed cubic equationu+ au3 = δ, (3.6)which, interestingly, admits a closed-form solution (see Appendix B for de-tails):u = δ32β13[(1 +√1 +1β)13 + (1−√1 +1β)13], (3.7)with β = 274 aδ2. In other words, in our model the lattice distortion is inher-ently linked to charge disproportionation, just as observed in experiments (ofcourse, this was already clear in the general Eq. 3.5, before any assumptionsabout the specific type of lattice/charge order). This is helpful: it allows usto eliminate any u dependence in the Hartree-Fock equations, leaving onlydependence on the electronic mean-fields and thus reducing the problem toa purely electronic one.It is important to note that whilst the argument above establishes thevalidity of a Hellman-Feynman like result for any variational stationary stateof the energy functional, in reality there can be small discrepancies arisingdue to round-off and other errors inherent to any numerical implementationof the calculation [71]. The “stationary state” thus obtained will not beexactly the true stationary state and thus in principle the theorem does notapply. However, assuming the errors can be made sufficiently small, thetheorem should still apply within the error margins: and our calculationsconfirm this, showing that there is no difference (other than speed) betweenrelying on this extended Hellmann-Feynman result and doing the foolproofparameter sweep.353.2. Minimizing electronic contributions3.2 Minimizing electronic contributionsWe now derive the self-consistency equations for the electronic degrees offreedom. Assume, as usual in the unrestricted Hartree-Fock approximation[72], that the true ground state is a Slater determinant of some (so farunknown) states with creation operators a†p|Ψe〉 =∏pa†p |0〉 , (3.8)with the new states being related to the original d†iaσ by a unitary transfor-mationd†iaσ =∑nφ∗n(iaσ)a†n. (3.9)This is the first assumption of mean-field (or Hartree-Fock) theory. The goalis to determine this unitary transformation, by calculating the energy in thisstate and subsequently minimizing it with respect to the variational (mean-field) parameters φn(iaσ) (equivalently, expectation values 〈d†iaσdjbτ 〉Ψe).The Hlatt term has no electron operators and vanishes during minimiza-tion. The electron-phonon term Hˆe−latt immediately yields〈Ψe| Hˆe−latt |Ψe〉 = 〈Ψe| − 2b∑iui(nˆi − 1) |Ψe〉 == −2b∑iui[∑paσφ∗p(iaσ)φp(iaσ)− 1]. (3.10)The hopping term Tˆ produces, in real space〈Tˆ 〉Ψe =∑p∑ijabσtabij φ∗p(iaσ)φp(jbσ). (3.11)Finally, consider the Kanamori interaction operator Hˆe. As this is a four-operator interaction and the state |Ψe〉 is a Slater determinant, to evaluatethe expectation value we can use Wick’s theorem [73]. In short, it says thatwhenever we have an expectation value with respect to a Slater determi-nant of an even number of creation/annihilation operators (say 〈d†↑d↑d†↓d↓〉),the result is a sum of products of all possible pairwise expectation val-ues (called contractions), with the sign reflecting the number of permu-tations of the operators that are required to attain that given order (so363.2. Minimizing electronic contributions〈d†↑d↑d†↓d↓〉 = 〈d†↑d↑〉〈d†↓d↓〉 − 〈d†↑d↓〉〈d†↓d↑〉). In principle even contractionslike 〈d†↑d†↓〉 could be allowed for a state more complex than a Slater deter-minant: this expectation value would be nonzero if the ground state againstwhich this expectation value was taken would allow two-electron excitations(for example, this expectation value would be nonzero for the BCS state,if we were trying to construct a model for superconductivity [74]). In ourcase, since we assume a normal state for the electrons, we automatically dropsuch terms from the Wick expansion. Using this technique on the KanamoriHamiltonian immediately leads to〈Hˆe〉Ψe =U∑ia[〈d†ia↑dia↑〉〈d†ia↓dia↓〉 − 〈d†ia↑dia↓〉〈d†ia↓dia↑〉〉]++ U ′∑iσ[〈d†izσdizσ〉〈d†iz¯σ¯diz¯σ¯〉 − 〈d†izσdiz¯σ¯〉〈d†iz¯σ¯dizσ〉]++ (U ′ − J)∑iσ[〈d†izσdizσ〉〈d†iz¯σdiz¯σ〉 − 〈d†izσdiz¯σ〉〈d†iz¯σdizσ〉]−− J∑iσ[〈d†izσdizσ¯〉〈d†iz¯σ¯diz¯σ〉 − 〈d†izσdiz¯σ〉〈d†iz¯σ¯dizσ¯〉]++ γJ∑ia[〈d†ia↑dia¯↑〉〈d†ia↓dia¯↓〉 − 〈d†ia↑dia¯↓〉〈d†ia↓dia¯↑〉]. (3.12)Since all interactions are on-site, there will be no matrix elements 〈d†iaσdjbτ 〉for i 6= j. With this, the full set of unrestricted Hartree-Fock equationsmay be obtained by varying the energy E = 〈Tˆ + Hˆe + Hlatt + Hˆe−latt〉with respect to the amplitudes φ∗p(iaσ), together with the constraint thatthe set φp(iaσ) are normalized amplitudes. Mathematically, that last partmeans 1 =∑iaσ φ∗p(iaσ)φp(iaσ) for all p. This constraint can be accom-plished by what is essentially a Lagrange multiplier technique: simply addthe (constant) term∑liaσ Elφ∗l (iaσ)φl(iaσ) to the Hamiltonian,H˜ = Hˆ −∑liaσElφ∗l (iaσ)φl(iaσ), (3.13)where El is the Lagrange multiplier, one for every set of amplitudes φl (and,simultaneously, the Hartree-Fock energy of that state). The normalizationcondition of φ can be recovered from the total expression in the usual way– by differentiating with respect to the Lagrange multiplier and setting theresult equal to zero∂∂EpH˜ = 1−∑iaσφ∗p(iaσ)φp(iaσ) = 0.373.2. Minimizing electronic contributionsAt the same time, varying H˜ with respect to the amplitudes φ gives theequations which determine these amplitudes, subject to the constraint thatthey are normalized:Epφp(iaσ) =δδφ∗p(iaσ)〈Hˆ〉. (3.14)These are the celebrated Hartree-Fock equations. We derive them step bystep.First, noticing that 〈d†iaσdjbσ′〉 =∑p φ∗p(iaσ)φp(jbσ′), a key identity maybe obtained that simplifies the variation considerably (especially for theKanamori term)δδφ∗p(lcτ)∑ijabσσ′〈d†iaσdjbσ′〉 =∑ijabσσ′δlcτ,iaσφp(jbσ′), (3.15)where δlcτ,iaσ is the Kronecker delta function which is zero unless l = i, c =a, τ = σ, in which case it is 1. Using this identity, we calculate the variationterm by term. Varying the hopping term we immediately obtainδδφ∗p(iaσ)(∑s∑ljcbτtcbljφ∗s(lcτ)φs(jbτ))=∑jbtabij φp(jbσ). (3.16)Varying the Ee−latt,δδφ∗p(iaσ)(− 2b∑juj[∑lbτφ∗l (jbτ)φl(jbτ)− 1])= −2buiφp(iaσ).(3.17)Finally, tackle varying the Kanamori contribution, term by term, using theidentity 3.15:first term of 〈Hˆ〉 −→ δδφ∗p(iaσ)(U∑jb[〈d†jb↑djb↑〉〈d†jb↓djb↓〉 − 〈d†jb↑djb↓〉×× 〈d†jb↓djb↑〉])= U(φp(ia ↑)δσ,↑〈d†ia↓dia↓〉+ φp(ia ↓)δσ,↓〈d†ia↑dia↑〉−− φp(ia ↓)δσ,↑〈d†ia↓dia↑〉 − φp(ia ↑)δσ,↓〈d†ia↑dia↓〉) == U(φp(iaσ)〈d†iaσ¯diaσ¯〉 − φp(iaσ¯)〈d†iaσ¯diaσ〉). (3.18)In the last line we made the replacements of the type φp(ia ↑)δσ,↑〈d†ia↓dia↓〉 →φp(iaσ)〈d†iaσ¯diaσ¯〉. These are valid because only one of every pair of terms383.2. Minimizing electronic contributionsin the second line is nonzero for a given σ: and for that nonzero term, thespins can be deduced relative to the original spin σ.Using the same tricks, the variations of the other terms can be evaluated:second term of 〈Hˆ〉 −→ δδφ∗p(iaσ)(U ′∑jσ′[〈d†jzσ′djzσ′〉〈d†jz¯σ¯′djz¯σ¯′〉−− 〈d†jzσ′djz¯σ¯′〉〈d†jz¯σ¯′djzσ′〉])= U ′(φp(iaσ)〈d†ia¯σ¯dia¯σ¯〉 − φp(ia¯σ¯)〈d†ia¯σ¯diaσ〉),(3.19)third term of 〈Hˆ〉 −→ δδφ∗p(iaσ)((U ′ − J)∑jσ′[〈d†jzσ′djzσ′〉〈d†jz¯σ′diz¯σ〉−− 〈d†izσdiz¯σ〉〈d†jz¯σ′djzσ′〉])= (U ′ − J)(φp(iaσ)〈d†ia¯σdia¯σ〉−− φp(ia¯σ)〈d†ia¯σdiaσ〉). (3.20)The simultaneous spin-flip term results infourth term of 〈Hˆ〉 −→ δδφ∗p(iaσ)(− J∑jσ′[〈d†jzσ′djzσ¯′〉〈d†jz¯σ¯′djz¯σ′〉−− 〈d†jzσ′djz¯σ′〉〈d†jz¯σ¯′djzσ¯′〉])= −J(φp(iaσ¯)〈d†ia¯σ¯dia¯σ〉 − φp(ia¯σ)〈d†ia¯σ¯diaσ¯〉).(3.21)Finally, the pair-hopping process givesfifth term of 〈Hˆ〉 −→ δδφ∗p(iaσ)(γJ∑jb[〈d†jb↑djb¯↑〉〈d†jb↓djb¯↓〉−− 〈d†jb↑djb¯↓〉〈d†jb↓djb¯↑〉])= γJ(φp(ia¯σ)〈d†iaσ¯dia¯σ¯〉 − φp(ia¯σ¯)〈d†iaσ¯dia¯σ〉).(3.22)393.3. Mean-Field ParametersCombining these variations, we obtain the Hartree-Fock equationsEpφp(iaσ) =∑jbtabij φp(jbσ) +[− 2bui + U〈d†iaσ¯diaσ¯〉+ U ′〈d†ia¯σ¯dia¯σ¯〉++ (U ′ − J)〈d†ia¯σdia¯σ〉]φp(iaσ) +[− U〈d†iaσ¯diaσ〉−− J〈d†ia¯σ¯dia¯σ〉]φp(iaσ¯) +[(U ′ − J)〈d†ia¯σdiaσ〉++ J〈d†ia¯σ¯diaσ¯〉+ γJ〈d†iaσ¯dia¯σ¯〉]φp(ia¯σ)++[U ′〈d†ia¯σ¯diaσ〉+ γJ〈d†iaσ¯dia¯σ〉]φp(ia¯σ¯). (3.23)This is clearly a nonlinear eigenvalue problem, as we have terms of orderφ and of order φ3 on the right-hand side of the equation. The easiest (andoften only) way to solve it is by iteration. Given an initial guess for the mean-field parameters w(0) = {〈d†iaσdibτ 〉}, this set of equations can be solved forthe Hartree-Fock energies Ep and the amplitudes (eigenstates) φp(iaσ). Inturn, the energies can be used to determine which “states” φp(iaσ) areoccupied (remember that since our calculation is done at zero energy, thereis a well-defined Fermi level). At zero energy this is determined by filling:in our case, there is one electron per Ni site, which corresponds to quarter-filling (2 orbitals and 2 spins per orbital = 4 states, 1 electron). The bottomquarter of all the states is then considered occupied: from these occupiedstates, the mean-field parameters v(0) = {〈d†iaσdibτ 〉′} can be re-calculatedand compared to the initial guesses for these parameters. If the two do notagree (to within some pre-specified margin of error), the newly calculatedparameters become guesses for the next iteration step, w(1) = v0. Thisprocess is then repeated until convergence.3.3 Mean-Field ParametersOf course, the (infinite-dimensional) set of mean-field parameters 〈d†iaσdibτ 〉is too general and needs to be specialized to the problem at hand. Thisis where the second mean-field assumption comes in: we will assume thatthe set of relevant parameters can be truncated past a particular value ofi. Specifically, we assume that only the parameters within a periodicallyrepeating unit cell are non-trivial. The particular choice of this periodicstructure is motivated by phase phenomena experimentally found in rare-earth nickelates: charge disproportionation, (potential for) orbital order,403.3. Mean-Field Parametersand the various magnetic possibilities, from ferromagnetism to exotic non-collinear 4-site orders. Given the ordering wavevectors for these behaviours(Qc and Qm, respectively), we conclude that (one choice for) the unit cellis a 4-site linearly-coordinated one. Thus only mean-field parameters fromi = (i, j, k) to i+ 4xˆ = (i+ 4, j, k) should matter: the rest of the mean-fieldparameters can be obtained by translation of the unit cell by the primitivelattice vectors.However, it is not clear a priori what starting guesses for the mean-fieldparameters 〈d†iaσdibτ 〉 are appropriate: we do not have any intuition as towhat values of these matrix elements would correspond to particular kindsof actual charge/lattice/magnetic orders. Instead, it would be useful if wehad some combinations of these mean-fields (let us call the combinationsorder parameters), which would readily signal the presence or absence ofa particular order.Consider, for example, the phenomenon of charge order, with a 3Dcheckerboard pattern represented by the ordering wavevectorQc =2pia(12,12,12).The average charge at site i is 〈nˆi〉 =∑aσ〈d†iaσdiaσ〉: suppose there is a dif-ference δ between neighbouring sites in all three crystallographic directions,namely〈nˆ(i,j,k)〉 − 〈nˆ(i,j,k)+η〉 = δ, η ∈ {xˆ, yˆ, zˆ}. (3.24)A shorthand way of writing this is by using the ordering wavevector: set〈nˆi〉 = 1 + δeiQc·Ri = 1 + δeipi(i+j+k), Qc = 2pia(12,12,12). (3.25)It is clear from this expression that on each alternate site along the threecrystallographic axes, the magnitude of the on-site charge will be eitherenhanced (if the integer in the exponent is even) or reduced (if it is odd)by the same factor of δ. Thus a δ 6= 0 readily indicates the presence of(rocksalt-type) charge order. On the other hand, as we just said above, theon-site average electron density is the sum of the densities in the respectiveorbital and spin states〈nˆi〉 =∑aσ〈d†iaσdiaσ〉. (3.26)413.3. Mean-Field ParametersSo one order parameter δ wraps into itself a variety of different matrixelements 〈d†iaσdibτ 〉, while simultaneously being a signature of a given phase:it then appears more prudent to search for such order parameters than torandomly guess values for the[×2(orbitals)×2(spins)](states) choose 2 (each matrix elementincludes two states) −→ 6 times 4, as there are 4 unique sites−−−−−−−−−−−−−−−−−−−−−→ 24different allowed matrix elements.Given that no experimental signatures of orbital order have been ob-served to date in the nickelates, a reasonable assumption to make is that atleast the average occupancy of the original orbitals |z〉 , |z¯〉 is the same,〈nˆiz〉 = 〈nˆiz¯〉. (3.27)This does not entirely preclude the possibility of orbital order, which couldstill occur due to nonzero cross-orbital matrix elements 〈d†iaσdia¯τ 〉. (Muchlike the cross-spin matrix elements 〈d†iaσdibσ¯〉 correspond to spin order off-the-axis, namely along the x or y axis instead of the z, cross-orbital matrixelements would result in orbital order across some combination of the orbitalorbitals |z〉 , |z¯〉.) Combined with the expression that defined δ in eq: 3.26,the assumption of equal a priori orbital occupancy means〈d†ia↑dia↑〉+ 〈d†ia↓dia↓〉 =12[1 + δeiQc·Ri]. (3.28)To proceed further, we need to consider the average spin on site i, 〈Si〉 =∑aττ ′〈d†iaτ σττ ′2 diaτ ′〉 (σττ ′ is the vector of Pauli matrices). Applying thesame sort of thinking, we write down an ansatz for the average spin per sitebased on the possible magnetic orders〈Si〉 = SFMzˆ + SAFMeiQc·Ri zˆ + S1 cos(Qm ·Ri) + S2 sin(Qm ·Ri). (3.29)The ferromagnetic and antiferromagnetic orders should be self-explanatory.The last two terms are the ones that produce the possibility of a 4-sitemagnetic order: this is because the “period” of the cosine or a sine with thewavevector Qm is going to be 4 sites. For instance, a cosine will evaluateto 0, then 1, then 0, then -1: supposing that S1 = (0, 0, S1z), the resultingmagnetic pattern would be ↑ 0 ↓ 0, precisely one of the contenders for theground state suggested in the literature. The sine term, on the other hand,is offset by one lattice site: an equivalent setup for the sine term would423.3. Mean-Field Parametersresult in the magnetic pattern being 0 ↑ 0 ↓. Thus when put together, theycan produce another candidate, ↑↑↓↓. Finally, picking S1 = (0, 0, S1z) andS2 = (S2x, 0, 0) would produce the non-collinear order ↑→↓←.Once again we will resort to experiment and assume equal orbital occupa-tion of the original orbitals |z〉 , |z¯〉. We also assume that the magnetic ordervectors S1/2 = (S1/2x, 0, S1/2z) can only have nonzero x and z components,so that there are only two independent axes involved in the noncollinearorder. Unpacking the expression for the average spin density in terms ofPauli matrices 〈Si〉 =∑aττ ′〈d†iaτ σττ ′2 diaτ ′〉 for the particular axes we find12[〈d†ia↑dia↑〉 − 〈d†ia↓dia↓〉]=12[SFM + SAFMeiQc·Ri++ S1z cos(Qm ·Ri) + S2z sin(Qm ·Ri)],(3.30)12[〈d†ia↑dia↓〉+ 〈d†ia↓dia↑〉]=12[S1x cos(Qm ·Ri) + S2x sin(Qm ·Ri)] .(3.31)The y axis calculation yields an identity− i2[〈d†ia↑dia↓〉 − 〈d†ia↓dia↑〉]= 0, (3.32)due to our assumption of only 2 independent non-collinear axes. Thus〈d†ia↑dia↓〉 = 〈d†ia↓dia↑〉. Combining these results with the previous insightsabout δ, we can write a general form for the direct and cross-spin matrix ele-ments 〈d†iaσdiaσ〉, 〈d†iaσdiaσ¯〉, with both spin-dependent and spin-independentparts, as follows〈d†iaσd†iaσ〉 =14[1 + δeiQc·Ri]+σ2[SFM + SAFMeiQc·Ri++ S1z cos(Qm ·Ri) + S2z sin(Qm ·Ri)],(3.33)〈d†iaσd†iaσ¯〉 =12[S1x cos(Qm ·Ri) + S2x sin(Qm ·Ri)] . (3.34)In a very similar manner, directly by analogy we establish the mean-fieldorder parameters for the orbital order, that arises from the cross-orbital433.4. Energy of the Hartree-Fock ground statematrix elements 〈d†iaσd†ia¯σ〉. Altogether, we have the following definitions〈d†iaσd†iaσ〉 =14[1 + δeiQc·Ri]+σ2[SFM + SAFMeiQc·Ri++ S1z cos(Qm ·Ri) + S2z sin(Qm ·Ri)],〈d†iaσd†iaσ¯〉 =12[S1x cos(Qm ·Ri) + S2x sin(Qm ·Ri)] ,〈d†iaσd†ia¯σ〉 =O1 +O2eiQc·Ri + σ[Z1 + Z2eiQc·Ri++ Z3 cos(Qm ·Ri) + Z4 sin(Qm ·Ri)],〈d†iaσd†ia¯σ¯〉 =X1 cos(Qm ·Ri) +X2 sin(Qm ·Ri). (3.35)The parameters thus defined transparently allow for the presence of differentkinds of ordering. For instance, if SAFM = 1 and the other S are identicallyzero, that corresponds directly to antiferromagnetic ordering〈Szi 〉 =∑a〈d†ia↑dia↑ − d†ia↓dia↓〉 = SFM + SAFMeiQc·Ri++ S1z cos(Qm ·Ri) + S2z sin(Qm ·Ri) = SAFMeiQc·Ri ∼ SAFM(↑↓↑↓).In the last line, and throughout this thesis we frequently use this arrownotation to refer to the magnetic order that obtains in the lattice, by simplywriting down the arrows for the spin alignment within a single 4-site linearly-coordinated unit cell. Notice that there are only 15 free parameters here(16 when counting the lattice distortion u, which however is fixed by δ),as opposed to the 24 matrix elements we started with: this is due to theassumptions we made when defining the possible order parameters based onwhat could be expected to be seen from experiment.The various suggested guesses for the magnetic order can be obtained asfollows: the ⇑ 0 ⇓ 0 obtains if δ > 0 and S1z > 0, while all others are zero.The sublattice-symmetric collinear order ↑↑↓↓ obtains if S1z = S2z 6= 0.Finally, a non-collinear state ↑→↓← can be obtained if S1z = S2x 6= 0.3.4 Energy of the Hartree-Fock ground stateIn general, multiple solutions can be found to 3.23 for any given set of pa-rameters t1, t2, t4, U, J etc.. To identify the true ground state, it is necessaryto calculate the total energy of each of the solutions, then pick the statewith the lowest energy. As the Hartree-Fock approximation is a variational443.4. Energy of the Hartree-Fock ground stateapproximation, such a solution will be an upper-bound on the true groundstate energy of the model. The ground state energy per site of a trial state|Ψe〉 is given byEGSN=1N〈Ψe| Tˆ + Hˆe + Hˆe−latt +Hlatt |Ψe〉 . (3.36)It is immediately clear that there is a problem10. The definitions of theorder parameters w = δ, SFM, ... rely on the single-site (localized) matrixelements 〈d†iaσdibτ 〉: thus it is best to evaluate the interaction contribution tothe energy 〈Hˆe + Hˆe−latt〉 in the localized, site-centered single-particle basis|iaσ〉. Meanwhile the hopping energy 〈Tˆ 〉 is best expressed in the extendedBloch state basis |kaσ〉 = ∑i eik·Ri√N |iaσ〉. To circumvent the need to converteither of the energy contributions to a cumbersome representation, we resortto a trick that allows us to re-express the hopping average 〈Tˆ 〉 in terms of theHartree-Fock energies Ep and thus eliminate it from the energy expressionentirely. Multiplying on the left by φ∗p(iaσ) in the Hartree-Fock equations3.23 and summing over i, we find∑pφ∗p(iaσ)δδφ∗p(iaσ)〈Hˆ〉 =∑pφ∗p(iaσ)Epφp(iaσ) =∑pEp. (3.37)Generally speaking, in terms of the amplitudes φp(iaσ) any energy functionalE[φp] with at most quartic interactions may be writtenE = C+∑pq∑ijabσtabij φ∗p(iaσ)φq(jbσ)++∑pqrs∑ijkl,abcd,σσ′gabcdijkl (σ, σ′)φ∗p(iaσ)φ∗q(jbσ)φr(kcσ′)φs(ldσ′).(3.38)Plugging this ansatz into Eq. 3.37, we notice that whatever the variationalderivative extracts in terms of the φp amplitude factors, the subsequentmultiplication puts right back in. So the left-hand side essentially yieldsthe energy, except that each term in the average is multiplied by a prefactorbased on the power of φ∗p in that term (remember that in this variation φ andφ∗ are treated as separate dynamical variables, as is standard in variationalcalculations involving complex variables). Thus the hopping term 〈Tˆ 〉 does10one inherent to any Hubbard model type Hamiltonian, where a competition betweenthe localized and itinerant behaviours is present.453.4. Energy of the Hartree-Fock ground statenot acquire a prefactor, but the electron-electron interaction 〈Hˆe〉 picks upa factor of 2: ∑pEp = 〈Tˆ 〉+ 2〈Hˆe〉+ 〈Hˆe−latt〉.Re-expressing the hopping contribution,〈Tˆ 〉 =∑pEp − 2〈Hˆe〉 − 〈Hˆe−latt〉,it becomes possible to write the expression for the total energy per sitewithout the Tˆ contributionEGSN=1N(∑pEp + Elatt − 〈Hˆe〉). (3.39)Using the definitions of the Hartree-Fock parameters 3.35 and the expressionfor 〈Hˆe〉 from 3.12, we obtain a closed form expression for the mean-fieldenergy in terms of the mean-field parameters onlyEGSN=1N∑pEp + 2b(u22+au44)− 1N〈Hˆe〉, (3.40)where1N〈Hˆe〉 =3U − 5J8(1 + δ2)− U + J2(S2FM + S2AFM++S21x + S21z + S22x + S22z2)− 2(U − J(4 + γ))(O21 +O22)−− 2(U + (γ − 2)J)(Z21 + Z22 +Z23 + Z242)−− (U + (γ − 2)J)(X21 +X22 ). (3.41)(The N comes from the fact that the sum over the sites i yields N identicalterms, as per the mean-field assumption. See Appendix C for the completecalculation of 〈Hˆe〉.) Once a self-consistent solution w is generated, itsenergy can be immediately evaluated using this expression.While the order parameters w can be guessed initially, it is still notentirely clear how to re-calculate them – or how to find the energies Ep in thefirst place. In the next section we construct the Hartree-Fock Hamiltonian,which will address both those points.463.5. Hartree-Fock Hamiltonian3.5 Hartree-Fock HamiltonianGiven the set of Hartree-Fock equationsEpφp(iaσ) =∑jbtabij φp(jbσ) +[− 2bui + U〈d†iaσ¯diaσ¯〉+ U ′〈d†ia¯σ¯dia¯σ¯〉++ (U ′ − J)〈d†ia¯σdia¯σ〉]φp(iaσ) +[− U〈d†iaσ¯diaσ〉−− J〈d†ia¯σ¯dia¯σ〉]φp(iaσ¯) +[(U ′ − J)〈d†ia¯σdiaσ〉++ J〈d†ia¯σ¯diaσ¯〉+ γJ〈d†iaσ¯dia¯σ¯〉]φp(ia¯σ)++[U ′〈d†ia¯σ¯diaσ〉+ γJ〈d†iaσ¯dia¯σ〉]φp(ia¯σ¯). (3.42)or, having in mind the definitions of the various order parameters 3.35 andthe spherical symmetry U ′ = U − 2JEpφp(iaσ) =∑jbtabij φp(jbσ) +{3U − 5J4+[3U − 5J4δ − 2bu]eiQc·Ri−−σ2(U+J)[SFM+SAFMeiQc·Ri+S1z cos Qm ·Ri+S2z sin Qm ·Ri]}φp(iaσ)−− U + J2[S1x cos Qm ·Ri + S2x sin Qm ·Ri]φp(iaσ¯)−− (U + J(γ − 2))[X1 cos Qm ·Ri +X2 sin Qm ·Ri]φp(ia¯σ¯)++{(J(4 + γ)− U)[O1 +O2eiQc·Ri]− σ[U + J(γ − 2)][Z1 + Z2eiQc·Ri++ Z3 cos Qm ·Ri + Z4 sin Qm ·Ri]}φp(ia¯σ),it is always reasonable to ask: what quadratic Hamiltonian could they haveoriginated from? In other words, what is the Hamiltonian that results af-ter one applies the Hartree-Fock approximation to the original, interactingHamiltonian in Eq. 2.21? Such a Hamiltonian – now quadratic – and thuseasily diagonalizable – is termed factorized : it is relatively easy to recon-struct it from the equations. Recalling that 〈d†iaσdjbτ 〉 =∑p φ∗p(iaσ)φp(jbτ),we can visualize what the quadratic Hamiltonian “should have been” in or-der to produce the equations 3.23 after the expectation value with the trialstate |Ψe〉 =∏p a†p |0〉 and the variation with respect to φ∗p(iaσ). For eachterm, we need to recover the factor φ∗p(iaσ): moreover, their sum over pought to be replaced by the appropriate matrix element 〈d†iaσdjbτ 〉 and thenthe expectation values stripped away, leaving the operators d. Following473.5. Hartree-Fock Hamiltonianthis procedure, it is straightforward to reconstruct the appropriate quadraticHamiltonianHˆHF =∑ijabσtabij d†iaσdjbσ +{3U − 5J4+[3U − 5J4δ − 2bu]eiQc·Ri−−σ2(U+J)[SFM+SAFMeiQc·Ri+S1z cos Qm ·Ri+S2z sin Qm ·Ri]}d†iaσdiaσ−− U + J2[S1x cos Qm ·Ri + S2x sin Qm ·Ri] d†iaσdiaσ¯−− (U + J(γ − 2))[X1 cos Qm ·Ri +X2 sin Qm ·Ri]d†iaσdia¯σ¯++{(J(4 + γ)− U)[O1 +O2eiQc·Ri]−σ[U+J(γ−2)][Z1+Z2eiQc·Ri+Z3 cos Qm·Ri+Z4 sin Qm·Ri]}d†iaσdia¯σ.This quadratic Hamiltonian is equivalent, in its ground state and eigen-functions/eigenvalues, to the original Hamiltonian under the Hartree-Fockapproximation: it is the Hartree-Fock Hamiltonian. It has the advantageof being quadratic and thus can be diagonalized, obtaining said eigenval-ues/eigenfunctions. Making use of Bloch’s theorem is very advantageoushere, as the Hamiltonian admits significant block-diagonalization with trans-forming to Fourier (k) space (which was not the case for the original Hamil-tonian because of the 4-operator terms). The details of this are carried outin Appendix D: the result, after some mundane substitutions of the usualFourier transformd†kaσ =∑ieik·Ri√Nd†iaσ (3.43)becomes483.5. Hartree-Fock HamiltonianHˆHF =∑kabσtab(k)d†kaσdkbσ +∑kaσ(3U − 5J4− σ2(U + J)SFM)d†kaσdkaσ++(3U − 5J4δ − 2bu− σ2(U + J)SAFM)d†k+Qc,aσdkaσ−− σ4(U + J)((S1z − iS2z)d†k+Qm,aσdkaσ + (S1z + iS2z)d†k−Qm,aσdkaσ)−− U + J4((S1x − iS2x)d†k+Qm,aσdkaσ¯ + (S1x + iS2x)d†k−Qm,aσdkaσ¯)− U + J(γ − 2)2((X1 − iX2)d†k+Qm,aσdka¯σ¯ + (S1x + iS2x)d†k−Qm,aσdka¯σ¯)++ [(J(4 + γ)− U)O1 − σ(U + J(γ − 2))Z1] d†kaσdka¯σ++ [(J(4 + γ)− U)O2 − σ(U + J(γ − 2))Z2] d†k+Qc,aσdka¯σ− σ2(U +J(γ−2))[(Z3 − iZ4)d†k+Qm,aσdka¯σ + (Z3 + iZ4)d†k−Qm,aσdka¯σ].Notice that the different Q prefactors have resulted in the mixing ofcertain momenta within the Brillouin zone. Instead of folding the Brillouinzone and labeling the resulting states (a so-called reduced zone scheme treat-ment of the Brillouin zone), we continue to work in the original Brillouinzone (so-called extended scheme [75]).At this point, the Hamiltonian is ready to be diagonalized: thanks totranslational symmetry, it is already block-diagonal in k and can be writtenHˆHF =∑kψ†kh(k)ψkwhereψ†k = (ψ†kz↑, ψ†kz¯↑, ψ†kz↓, ψ†kz¯↓)andψ†kaσ = (d†kaσ, d†k+Qm,aσ, d†k+Qc,aσ, d†k−Qm,aσ).Hence only separate 16× 16 subblocks need to be explicitly diagonalized –a process we implemented in Python and discuss in a subsequent chapter.Explicit expressions for the matrices are available in Appendix D. (Notethat this matrix is, due to the properties of the Fourier transform, almost aToeplitz matrix.)Once the Hamiltonian is diagonalized, to complete the iteration it isnecessary to re-calculate the order parameters w from the eigenvalues and493.5. Hartree-Fock Hamiltonianeigenfunctions of HˆHF. First off, the eigenfunctions are ordered by theircorresponding eigenvalues: according to our quarter-filling scheme, only thebottom quarter are selected to participate in the calculation of the orderparameters. The Slater determinant (Hartree-Fock trial wavefunction) is|Ψe〉 =′∏kaσd†kaσ |0〉where the prime indicates the product over the occupied single-particlestates. The matrix elements 〈Ψe| d†kaσdkbτ |Ψe〉 need to be calculated andcompared to the values that were guessed for them initially in order tocomplete the iterative process. The easiest way to proceed is by using thedensity matrix approach [72]. Define the density matrix associated to theHartree-Fock ground state byρk(ζbτ, ηaσ) = 〈Ψe| d†k+ηQm,aσdk+ζQm,bτ |Ψe〉 . (3.44)Insert the identity decomposition 1 =∑kaσ |kaσ〉 〈kaσ| in between the cre-ation/annihilation operators,ρk(ζbτ, ηaσ) =∑pbσ′〈Ψe| d†k+ηQm,aσ∣∣pbσ′〉 〈pbσ′∣∣ dk+ζQm,bτ |Ψe〉 . (3.45)Each of the matrix elements can only be nonzero if the state |pbσ′〉 is occu-pied, i.e. it is part of the ground state. So the sum is now primed – only overthe occupied states. Moreover, since the expectation value is taken with re-spect to the ground state of the system, 〈Ψe| d†k+ηQm,aσ and dk+ζQm,bτ |Ψe〉should be considered as single-state excitations. This leads toρk(ζbτ, ηaσ) = 〈k + ηQmaσ| ′∑pbσ′∣∣pbσ′〉 〈pbσ′∣∣ |k + ζQm, bτ〉 . (3.46)This equation shows that the matrix elements 〈Ψe| d†k+ηQm,aσdk+ζQm,bτ |Ψe〉are the (reversed) matrix elements of the operator ρ =∑′pbσ′ |pbσ′〉 〈pbσ′|,constructed out of occupied single-particle Hartree-Fock states.The single-particle states themselves are 16-tuples |kaσ〉 of numbers.They can be used to construct the density matrix, which can then be usedto find the order parameters δ, SFM, ... via the Fourier transform rule. Thisis most easily understood when demonstrated on an explicit example. Recall503.5. Hartree-Fock Hamiltonianthe definitions of the order parameters〈d†iaσdiaσ〉 =14[1 + δeiQc·Ri]+σ2[SFM + SAFMeiQc·Ri++ S1z cos(Qm ·Ri) + S2z sin(Qm ·Ri)].Re-arranging into a more convenient form, sorting by the order of the expo-nential〈d†iaσdiaσ〉 =[14+σ2SFM]e0 +[δ4+σ2SAFM]eiQc·Ri++σ2[S1z − iS2z2]eiQm·Ri +σ2[S1z + iS2z2]e−iQm·Ri .Notice that by hitting the matrix element with the appropriate exponentialand summing over the sites, we can “project out” the quantities we areinterested in: the terms which do not have the same exponential will sumto zero. For instance, simply summing over sites i and orbitals a leavesonly the non-varying terms (other terms oscillate in sign and cancel out, atypical trick in condensed matter calculations – see [75] appendix for thisshown explicitly) – the sum over a only contributes an overall factor of 2∑ia〈d†iaσdiaσ〉 = N[12+ σSFM].Using the fact that the sign is different for the up and down spins, we canisolate for SFM by multiplying by σ and summing over itSFM =12N∑iaσσ〈d†iaσdiaσ〉.Finally, to make connection with the |k + ηQm, aσ〉 = d†k+ηQm,aσ |0〉 , η ∈{0, 1, 2,−1} eigenfunctions that are actually found during diagonalizationof the Hartree-Fock Hamiltonian, Fourier-transform the operators inside theexpectation values. Since in this particular case only on-site operators areinvolved and there are no exponential prefactors, the Fourier transform isimmediate and yieldsSFM =12N∑kaσσ〈d†kaσdkaσ〉 =12N∑kηaσσρk(ηaσ, ηaσ). (3.47)The matrix element 〈d†kaσdkaσ〉 is an expectation value of the number oper-ator nˆkaσ in the Hartree-Fock ground state |Ψe〉. It is also simply an entry513.5. Hartree-Fock Hamiltonianof the density matrix 3.46. Let us explicitly illustrate how to evaluate such asum relative to the 16×16 block density matrix at each value of momentum.In this case, the sum is over diagonal elements only: yet even for a singlevalue of kaσ, there are still four entries that need summing over, namelySFM =12N∑k in [−pia,pia]3and aσσ(〈d†kaσdkaσ〉+ 〈d†k+Qm,aσdk+Qm,aσ〉++ 〈d†k+2Qm,aσdk+2Qm,aσ〉+ 〈d†k−Qm,aσdk−Qm,aσ〉)==12N∑k in [−pia,pia]3and ηaσσ(ρk(0aσ, 0aσ) + ρk(1aσ, 1aσ)++ ρk(2aσ, 2aσ) + ρk(−1aσ,−1aσ)). (3.48)Similar expressions can be obtained for the other order parameters. Notethat for certain sums with displaced indices, it is necessary to sum over“looped” entries. For instance, consider obtaining the charge dispropor-tionation order parameter: after the usual exponential multiplication andsummation trick,∑ie2iQm·Ri〈d†iaσdiaσ〉 = N[δ4+σ2SAFM].Fourier transforming on the left and summing over the spin and orbitaldegrees of freedom,δ =1N∑kaσ〈d†k+Qc,aσdkaσ〉. (3.49)Summing over all allowed values of k, we inevitably stumble upon entrieslike 〈d†k+2Qc,aσdk+Qc,aσ〉 = 〈d†k,aσdk+Qc,aσ〉: such entries are from the uppertriangle of the matrix subblock and need to be summed together with thosefrom the lower triangle. Thus, the matrix 16× 16 subblock can be thoughtof as having “bands” that correspond to particular momentum space sepa-ration between the column and row. In the example above, the separation isQc: the matrix bands that correspond to that are of (i, j + 2) and (i+ 2, j)type, that is, the second off-diagonals. In the case of SFM there was zeroseparation (this corresponds to the true diagonal of the matrix subblock).For the curious case of Qm separation, the bands consists of the nearest523.5. Hartree-Fock Hamiltonianoff-diagonal and together with the opposite corner of the matrix. Inter-estingly, for the case of both δ and SAFM, this sort of summation ensuresthey are both real without requiring the matrix entries themselves to be real(because the sum is over the two diagonals that are each other’s Hermitianconjugates). Naturally, the reality of SFM is ensured by the fact that it onlyincludes diagonal entries. All of the order parameters are re-calculated andthe expressions are cataloged in Appendix E.One final point: since it is clear from the above that there is no mo-mentum mixing (this is a feature of the Hartree-Fock approximation), froma computational perspective it is useful to first quickly calculate the totaldensity matrix ρ =∑′ |kaσ〉 〈kaσ|, then merely extract the components onerequires from it.53Chapter 4Numerical implementation ofthe Hartree-FockapproximationPart of the inhumanity of thecomputer is that, once it iscompetently programmed andworking smoothly, it iscompletely honest.Isaac Asimov (1981)The equations have been derived, the expression for the total energyobtained: we have reached the extent of what can be done by hand. Thenext step is to implement the solution procedure (and, much more impor-tantly, automated sweeping/analysis) on the computer. This is necessarynot only because it is infeasible to diagonalize 16 × 16 matrices by hand(especially ones that involve complex expressions, such as the ansatzes for〈d†iaσdibτ 〉 and the hopping coefficients tab(k)), but also because the numberof mean-field order parameters w and adjustable parameters U, J, ... is ex-tremely large when compared to the typical mean-field model. Its thoroughanalysis demands a computerized approach, such that the true ground statein the multi-parameter w space is surely identified, and all the relationshipsbetween the different adjustable parameters and how they affect the groundstate is duly cataloged and categorized. In this chapter we describe thisprocess in detail.4.1 The Hartree-Fock iterative algorithmThe task at hand is to find the ground state of the Hamiltonian Hˆ 2.21. Asalready mentioned, in this thesis this is accomplished within the Hartree-Fock approximation: this produces a properly factorized companion Hamil-544.1. The Hartree-Fock iterative algorithmtonian HˆHF. Even though the Hamiltonian thus obtained is quadratic andso nominally exactly solvable, its coefficients depend on its own eigenfunc-tions: the Hartree-Fock equations are nonlinear and thus must be solvediteratively. The Hamiltonian to diagonalize isHˆHF =∑kψ†kh(k)ψk (4.1)whereψ†k = (ψ†kz↑, ψ†kz¯↑, ψ†kz↓, ψ†kz¯↓), ψ†kaσ = (d†kaσ, d†k+Qm,aσ, d†k+Qc,aσ, d†k−Qm,aσ),and the matrix h(k) is given in Appendix D. The energy can be evaluatedwith the help of the expression from Chapter 3EGSN=1N∑pEp + 2bN(u22+au44− 〈Hˆe〉), (4.2)where1N〈Hˆe〉 =3U − 5J8(1 + δ2)− U + J2(S2FM + S2AFM++S21x + S21z + S22x + S22z2)− 2(U − J(4 + γ))(O21 +O22)−− 2(U + (γ − 2)J)(Z21 + Z22 +Z23 + Z242)−− (U + (γ − 2)J)(X21 +X22 ). (4.3)Finally, the matrix elements that the Hartree-Fock Hamiltonian relieson can be re-calculated from the density matrix using the expressions in theAppendix E.Let us sketch out the iterative solution process by which the ground stateis found.1. Choose parameter values U, J, ... from the parameter regime of interest.2. Create a pool of initial guess values for the mean-field parametersw(0) = (δ(0), S(0)FM, S(0)AFM, ...). Think physically about what they meanand what states/symmetries could be present in this parameter region.Without good guesswork, it is likely the iterative solution will notwork.554.2. Strategies for improving computation(a) Pick one of the initial guesses and start the iterative solutionsequence:(b) i. Initialize the Hartree-Fock matrix HˆHF using the guessed or-der parameter values.ii. Diagonalize the Hamiltonian, then build the ground statedensity matrix ρ =∑′kaσ |kaσ〉 〈kaσ| at quarter-filing. Thatis, order the eigenvectors obtained during minimization bytheir eigenvalues, then pick the bottom quarter.iii. Re-compute the order parameters v(0) using the expressionsin the Appendix E: evaluate the residual (0) =∣∣w(0) − v(0)∣∣2.iv. Check if (0) < 0, the pre-established convergence condition.If yes, output the eigenvectors |kaσ〉 and eigenvalues Ekaσ.If not, repeat this loop with the initial guess being the newlyfound solution, w(1) = v(0) (or choose it in an more involvedway – more on such techniques later). If the maximum al-lowed number of iterations nmaxiter is exceeded, quit.(c) If the sequence converged, evaluate the total state energy EHF.3. When all the guesses from the pool have been explored, pick the lowestenergy solution as the best approximation of ground state.In principle, this algorithm, together with the Eqs. 4.1–4.2, should besufficient to fully solve the problem. However, in principle the energy land-scape can be complex and have multiple local minima, manifesting as differ-ent unique solutions to the self-consistency equations. We found that thereare many possibilities for improving convergence for difficult types of guessstates, ensuring that the true ground state is identified, as well as simplyspeeding up the convergence process.4.2 Strategies for improving computation4.2.1 Pulay mixingAlthough the self-consistent algorithm (4.1) suggests to always re-start thealgorithm on step n + 1 with the updated guess w(n) being the solutionv(n−1) from the previous step, in fact one can do much better. The simplestgeneralization is to use, as the guess for the next step w(n), an interpolationbetween the guess w(n−1) and the solution v(n−1) obtained with this guess,w(n) = v(n−1)α+ w(n−1)(1− α). (4.4)564.2. Strategies for improving computationThe literature on nonlinear iterative equation solution techniques suggeststhe optimal choice of α that ensures fastest convergence is highly problemspecific [76]. Our experience appears to be that α = 0.3 produces thesmoothest and most reliable iterations. The key advantage of this “cautiousconvergence” technique that only admixes part of the newly-found solutioninto the guess is being able to avoid the iteration getting stuck in a loop, orhelplessly hopping on either edge of a flat minimum.However, this approach works best in conjunction with another ideacalled Pulay mixing, or direct inversion in the iterative subspace (DIIS). Theidea behind DIIS, originally developed for high-dimensional11 Hartree-Fockquantum chemistry calculations by Pulay [77], is the following: suppose a se-quence of N solutions v(1),v(2), ...v(N) to the nonlinear system has been gen-erated. Together, they span a linear subspace VN = span{v(1),w(2), ...w(N)}within the higher-dimensional (possibly nonlinear) solution manifold V . Itis possible to pick the “best possible” vector w(N+1) within this linear sub-space, by minimizing the total error E.To construct this best possible vector, consider a linear combination ofthe vectors spanning the solution subspace with some a priori unknowncoefficients cmw(N+1) =N∑m=1cmv(m). (4.5)It turns out that a good way to pick the coefficients cm is so as to minimizethe residuum vector ∆w(N+1), defined as∆w(N+1) =N∑m=1cm(v(m) −w(m)), (4.6)subject to the constraint that the coefficients sum to one (i.e. the linearcombination in 4.5 is actually a convex combination, picking a vector frominside the space VN )N∑m=1cm = 1. (4.7)Define the total error E as the magnitude squared of the residuum vector,E = ||∑Nm cm(v(m)−w(m))||2, and minimize it, together with the Lagrange11in terms of the number of order parameters involved.574.2. Strategies for improving computationmultiplier factor λ(1−∑Nm cm) (which ensures the condition 4.7 holds) withrespect to the coefficients cm. Differentiating with respect to cj gives a seriesof equations∂∂cj(||N∑mcm(v(m) −w(m))||2 + λ(1−N∑mcm))=∂∂cj( N∑mc2m(v(m)−−w(m))2 + 2N∑m<ncmcn(v(n) −w(n)) · (v(m) −w(m)) + λ(1−N∑mcm))== 2cj(v(j) −w(j))2 + 2∑m 6=jcm(v(j) −w(j)) · (v(m) −w(m))− λ ==∑mcm(v(j) −w(j)) · (v(m) −w(m))− λ. (4.8)In the last step a re-scaling of the (arbitrary) Lagrange multiplier λ allowedus to get rid of the factor of 2. Adopting the notation Bjm = (v(j) −w(j)) ·(v(m)−w(m)), the system of equations, together with the constraint 4.7 maybe written B11 B12 .. B1N −1B21 B22 .. B2N −1. ... .. ... ..BN1 BN2 .. BNN −1−1 −1 .. −1 0c1c2...cNλ =00..0−1 . (4.9)So the search for the best possible vector within the subspace VN amountsto solving an N + 1×N + 1 system of linear equations (hence the “inversionof the iterative subspace” – for more details, see Ref. [77]). This approachalready yields an improvement on the iterative process.But one can do even better than this. Of course, the Pulay approach canquickly maximize the potential of the vectors within the subspace VN : butwhat if the true solution is outside the subspace by more than the allowedresidual 0 to begin with? Then no matter how many times Pulay mixing iscarried out, it will not result in improved convergence – instead, it will slowdown the search. The fix to this, as shown by Banerjee and co-workers [76],is to intersperse Pulay mixing with regular updates of the form 4.4, witha given periodicity k (typically k = 3). The algorithm thus alternates be-tween expanding its iterative subspace, and finding the lowest residual vectorwithin it, resulting in optimal convergence for most parameter values. Infact, our calculations show that this approach can even arrest divergences584.2. Strategies for improving computationfrom round-off error accumulation. Small errors (e.g. some order param-eters being nonzero while they should be zero for this particular type ofstate) often cause rapid divergence of the Hartree-Fock iterative trajectoryaway from a solution. Yet the Pulay approach appears to counteract thattendency, gently guiding the trajectory towards the self-consistent point byresetting the creeping small errors in irrelevant order parameters back tozero.4.2.2 Convergence in iteration count and momentumsamplingFrom the standpoint of standard iterative solution theory, there are two keyparameters that need to be tested for convergence: the cutoff iteration countnmaxiter and the number of points sampled in the Brillouin zone nnum-k-pts.It is clear that too fine of a sampling of momentum space can slow downthe calculation time as O(N3): as for the cutoff iteration count, it needs tobe large enough to allow convergence for the interesting solutions, while notso large so as to waste computation time on “dead-end” iterative trajecto-ries that never converge. Long-time convergence studies, studying both thenumber of identified self-consistent solutions, as well as the behaviour of theper-site ground state energy for those solutions, show that a k-space sam-pling rate of 25 points per dimension, or nnum-k-pts = 253 ∼ 15, 000 pointstotal, along with the cutoff iteration count nmax-iter = 500 are optimal forour purposes, taming the residual error (n) to less that 0 = 10−4 and thedifferences in energy between subsequent iterations to less than 10−3.4.2.3 Choosing initial guessesA good update choice strategy for picking the next guess for the itera-tion, such as Pulay mixing, is paramount to reaching converged solutions.However, equally important is the choice of the initial guesses w(0) for theiterative trajectory: our early calculations quickly showed that a lion’s shareof the initial guesses chosen at random never lead to interesting convergedsolutions even if they are known to exist, no matter how high the cutoffnmatiter is. Through a combination of physical thinking and grueling calcu-lations with extremely large nmax-iter and fine random sampling of the initialguesses, we identified a multitude of key converged states, which include theferro/antiferromagnetic solutions (↑↑↑↑ and ↑↓↑↓), the non-interacting zero-spin solution 0000, fully charge-disproportionated 4-site antiferromagneticsolution (⇑ 0 ⇓ 0), and the collinear and non-collinear 4-site solutions ↑↑↓↓,594.2. Strategies for improving computation↑→↓←, along with a few others. For subsequent calculations these high-symmetry solutions are always included as initial guesses: this enhancedboth the reliability of converging one of these states, as well as the speed ofsweeping the parameter space.4.2.4 ParallelizationAs can be seen from the self-consistent algorithm (4.1), the calculation ofan iterative trajectory is not an easily parallelalizable process. The diag-onalization of an individual 16 × 16 block submatrix for each momentumcan be done in parallel, but the additional difficulties of structuring thestoring and writing of the matrix and its eigenvalues/eigenvectors and thefact that each diagonalization process is very fast on its own makes thisoption unattractive. On the other hand, there are more high-level paral-lelization opportunities possible: (a) multiple initial guesses can be madeand followed through on for the same set of parameter values U, J, ..., and(b) ground states can be found at multiple parameter values at once. In ourcalculation we opted for option (b) due to ease of implementation and datamanagement and because it resulted in immediate linear speed-up of thecalculation (of course, ground states corresponding to different parametervalues are independent from each other, hence there is no overhead to theparallelization).4.2.5 Boot-strapping with neighbouring pointsWhile the calculations at neighbouring points in the parameter space maybe independent from each other, they can benefit from using each other’sconverged solutions as initial guesses for starting new iterative trajectories.On the grounds of adiabaticity (which should apply at least far away fromany phase boundaries), we expect small changes in the values of various pa-rameters U, J, t1, ... to lead to small changes in the nature of the ground state(characterized by its order parameter vector w): as such, an already con-verged solution from a neighbouring point can result in fast convergence toa similar kind of converged solution at the current point, further improvingthe reliability of the calculation and reducing computation time.Armed with the algorithm 4.1, the analytical equations 4.1 – 4.2 and thenumerical tricks from this chapter, we are ready to tackle the problem andlook at the results.60Chapter 5ResultsAn approximate answer to theright problem is worth a gooddeal more than an exact answerto an approximate problem.John Tukey (1962)Before discussing the results of applying the Hartree-Fock approximationto the nickelates, let us once again reiterate our inspiration for and goalsfor our research. Our study of the nickelates is largely motivated by therecent observation that the coupling to the lattice distortion has a signifi-cant impact on electronic behaviour and might ultimately be responsible fortriggering the MIT, with no need for actual charge disproportionation [32].In light of this, and the peculiar fact that the magnetic transition only coin-cides with the MIT for sufficiently weak lattice coupling (manifested by theEu to Lu part of the nickelate series), we wanted to investigate what impact,if any, the lattice coupling can have on the magnetic transition. So we builtan effective two-band Hubbard model for the nickelates, in the spirit of thecalculation of Lee et al. [26], and coupled it to the lattice semiclassically,with three goals in mind:1. Confirm quantitatively the earlier findings that the lattice couplingforces charge disproportionation;2. See what magnetic orders are self-consistent, in a mean-field calcula-tion, within our model;3. Observe whether one of the orders suggested in the literature, or an-other order that obtains in the calculation, is stabilized by the couplingto the lattice.615.1. Charge order and the Metal-Insulator Transition5.1 Charge order and the Metal-InsulatorTransitionThe experimental observation that the temperature of the metal-insulatortransition TMIT is higher than the temperature of the magnetic transition TNfor the heavier nickelates, along with the strong effect that higher neutroncount isotope substitutions of O18 have on this temperature [31], suggest theprimacy of electron-phonon interactions in the charge ordering process overthat of the electron-electron interactions. A recent study [32] clearly demon-strated the possibility of charge disproportionation without overt chargetransfer in the nickelates (see Sec. 1.4), being driven mostly by the cou-pling of ligand holes to the (static) lattice distortions. Here we reproducethis result, showing that even a Holstein-like on-site coupling between thelattice and the electronic degrees of freedom can drive the system through ametal-insulator transition. This is because the on-site coupling mimics theeffect of encouraging double occupancy of the shrunk site, as is obtainedmore naturally from the ligand hole picture.To first get an idea of what the phase space is like and how chargedisproportionation δ depends on the various adjustable parameters, we plotthe value of δ for the ground state in the U–J plane. Other parameters areset to the values used in the Lee et al. study [26]: t1 = 1, t2 = 0.15, t4 = 0and there is no coupling to the lattice (b = 0) for the moment. The data inFig. 5.1 consists of a 40×40 grid of solutions, where at each point in the gridwe carried out the full self-consistent calculation and obtained the groundstate. In total, to obtain the figure, we carried out 1,600 calculations, each ofwhich started from a pool of 20-30 initial guesses (or more, depending on itsneighbours) and proceeded through anywhere from 5 to 250 iterative steps(where each iterative step involves the diagonalization of roughly 15,00016 × 16 subblocks of the Hartree-Fock matrix, as well as the calculation ofthe density matrix and the order parameters w).All in all, Fig. 5.1 shows qualitative agreement with earlier studies (suchas those of Lee el al. [26] or Peters [78]), as well as with basic physicalintuition. The small U , small J region exhibits no charge modulation, as itis largely metallic (the hopping dominates the interactions), with a robustweight at the Fermi level (see density of state plots for various representa-tive points in the U–J phase diagram in Fig. 5.2). Upon increasing U nomodulation arises either: strong on-site Coulomb repulsion opposes chargedensity waves, instead favouring on-site localization. Finally, regardless ofthe magnitude of U , for J large enough it becomes possible to stabilize625.1. Charge order and the Metal-Insulator TransitionFigure 5.1: Charge modulation δ (the magnitude is reflected in the colourbar) of the ground state in the U–J plane. Other parameters are t1 = 1, t2 =0.15, t4 = 0, b = 0. Resolution is 40× 40.635.1. Charge order and the Metal-Insulator TransitionFigure 5.2: Representative density of states diagrams for the various phasesin Fig. 5.1: (a) metallic phase, U = 1.077, J = 0.538; (b) itinerant magneticphase (see the magnetic phase analysis below for more details about thisphase), U = 3.231, J = 0.385; (c) charge modulated phase, with a cleargap at the Fermi level, U = 2.0, J = 1.538; (d) Mott insulating phase,U = 5.846, J = 0.231. Other parameters as in Fig. 5.1, with the additionof α = 1. Notice the increased weight (called the van Hove singularity) atthe lower band edge: its presence is due to the nonzero t2 parameter, whichintroduces a strong asymmetry to the DOS. More on this below.charge density wave order thanks to Hund’s rule. It minimizes the energyby “poaching” an electron from a neighbouring site, then organizing thetwo electrons into a triplet state across the two orbitals |z〉 , |z¯〉. This isfavourable even despite the increased on-site inter-orbital repulsion from U(although it does take a larger J to accomplish this if U is large).What of the lattice? To investigate the impact from turning on thelattice coupling b we turn to Figs. 5.3 and 5.4, where the U–b diagrams ofcharge disproportionation and metal-insulator behaviour are presented forvarious ratios of J/U , from a modest 0.2 to a more experimentally reasonable0.5. These diagrams essentially represent line slices of the U–J diagramemanating from the origin in the lower left corner, with the new dimension645.1. Charge order and the Metal-Insulator TransitionFigure 5.3: Charge modulation δ (see colour scale) in the HF ground-state,as a function of U and b, for J/U = 0.2, 0.3 (top left and right, respectively)and 0.4, 0.5, (bottom left and right, respectively). Other parameters aret1 = 1, t2 = 0.15, t4 = 0.25, α = 1. Resolution is 20× 20.showing the impact of lattice coupling b. Here we let t4 = 0.25: the impactof this on the diagram is not extensive and will be reviewed later. The reasonwe include it here is to ensure that all of the pieces of our model (includingthe 4th neighbour hopping) are present during the analysis.Overall the behaviour again conforms to expectations from earlier work.For small U (and thus an even smaller J) the system is fully metallic (ascan be seem from the metal-insulator diagrams), and δ = u = 0 (recall thatin our model u is essentially fixed by δ as per the self-consistency equationfor u). In the absence of coupling to the lattice, just as was seen in Fig. 5.1,there is already a transition to a finite δ 6= 0 state with increasing J , whichoccurs sooner for larger slopes J/U . This can be readily explained, as forlarger J/U the propensity for two-orbital triplet occupancy (Hund’s order)grows faster than the push towards single-site localization (Mott-Hubbard655.1. Charge order and the Metal-Insulator TransitionFigure 5.4: The HF ground-state is metallic (deep blue) or insulating (yel-low). The results are shown in the U -b space, for J/U = 0.2, 0.3, 0.4 and0.5, respectively (panels arranged as in Fig. 5.3). All other parameters areas in Fig. 5.3.665.2. Magnetic order: types and the transitionlocalization). And once the coupling to the lattice is introduced, no matterwhat U and slope J/U we are considering, at large enough coupling strengthsthe system will gravitate towards the charge modulated state.To summarize, we find that our model supports the existence of the MIT,triggered by modulating material parameters (in particular, electron inter-action strength U, J). It manifests as a charge density wave with nonzerocharge modulation δ 6= 0. The introduction of the lattice further enhancesthe favourability of the charge density wave order across the parameterspace, the onset being more dramatic for larger values of J/U : the resultinglattice distortion u 6= 0 is readily obtained from the charge disproportiona-tion δ via the self-consistency condition derived earlier in Eq. 3.5.5.2 Magnetic order: types and the transitionHaving confirmed the existence of the MIT within the model and “soft-benchmarked” it against earlier works by confirming that charge modulationarises for large J , we move on the primary subject of this thesis: the issue ofthe magnetic order and its interaction with the lattice degrees of freedom.Once again, we begin by analyzing the large-scale picture of phases in theU–J plane. Before we discuss the insights it provides, note that in general atmost points in the diagram there are many different types of self-consistentsolutions from which to pick the ground state. Below is a list of most of thetypes of magnetic order states that we identified during our searches:1. The “no magnetic order”, fully metallic solution 0000;2. The ferromagnetic solution ↑↑↑↑;3. The antiferromagnetic solution ↑↓↑↓;4. The fully disproportionated antiferromagnetic solution ⇑ 0 ⇓ 0 order,favoured by some in the literature [25];5. The non-disproportionated collinear solution ↑↑↓↓, the focus of thework of Lee et al. [26];6. The non-disproportionated non-collinear solution ↑→↓←, suggestedby several recent studies [27, 28];7. A variety of unexpected orbital-order states such as spin ordered ⇑↓⇑↓with ferroorbital order of type XxXx (an x indicates a preference forthe orbital |x〉 = (1/2)(|z〉+ |z¯〉), with magnitude again reflecting the675.2. Magnetic order: types and the transitionsize of the preference): we will refer to this state as the disproportion-ated ferroorbital order;8. Moreover, we see the states ⇑↑⇑↑ (the aligned ferrimagnetic state), aswell as various odd asymmetric ones such as ⇑↓↑↓ (but they are alwayshigher energy states and involve orbital order).9. There are even more possibilities, which sometimes even manifest asground states, that are associated with complex combinations of thespin and orbital order: states such as ↑→↓←, 0x0x¯, or→→←← xxx¯x¯.In this thesis we do not spend much time studying such states, partlybecause they rarely occupy much space in the phase diagram. Theymanifest mostly when the hopping integrals are adjusted. All suchstates are marked with “misc” label in the phase diagrams and theirstudy is relegated to future work.Only a select few of these states actually end up having a robust presencein the magnetic phase diagrams. Consider Figs. 5.5, the magnetic phasediagrams in the U–J plane at the reference parameter values of Lee et al.,with the coupling to the lattice b = 0 in the first and b = 0.8 in the second.The black contours represent the (interpolated) value of δ. There are severalitems of note here: first of all, the magnetic diagram with no lattice cou-pling again exhibits qualitative agreement with prior work: a non-magnetic,metallic phase is found for small U/small J ; large values of Hubbard U re-sult in the presence of a ferromagnetic state ↑↑↑↑ due to localization of theelectrons; and large values of J lead to charge modulation and as a con-sequence two inequivalent spin sublattices, one of which is almost entirelydevoid of magnetic moment and the other is antiferromagnetically ordered,⇑ 0 ⇓ 0.One crucial distinction we point out at this stage: Lee et al. find a thinribbon-like region of intermediate-magnet phase ⇑↑⇓↓ between the ⇑ 0 ⇓ 0large J phase and the ↑↑↓↓ medium-to-large U phase. We do not find suchan intermediate phase at all: in fact, in all our calculations we have neverbeen able to successfully converge a magnetic state with (a) nonzero chargedisproportionation and (b) both magnetic sublattices having nonvanishingmagnetic moments. It appears that solutions of this type are simply notself-consistent within our model: in fact, they seem to be among the leaststable points in the phase space (see Sec. 5.5 about stability below). Whatinstead obtains at the boundary between the two phases is a typically narrowregion where the relative energies of the two states are the same to within thenumerical precision: such regions are designated “degenerate” in the phase685.2. Magnetic order: types and the transitionFigure 5.5: Top: U -J phase diagram of magnetic order, in the absence ofcoupling to the lattice (b = 0). Bottom: Same when coupling to the latticeis turned on (b = 0.8). Other parameters are t1 = 1, t2 = 0.15, t4 = 0, α = 1for both. The black-line contours indicate the value of charge modulationδ. Resolution is 40× 40.695.2. Magnetic order: types and the transitiondiagram legend. It is possible that the only reason Lee et al. saw such aregion in the first place was because of their specific technique for carryingout the mean-field approximation by explicitly minimizing the energy as afunction of the mean-field parameters (that is, via a parameter search), asopposed to solving the Hartree-Fock equations and looking for self-consistentpoints. These two techniques, while formally ought to yield the same results,can sometimes involve hidden assumptions and thus not quite produce thesame outcomes (as was discussed in Chapter 3).Aside from the intermediate phase issue, the qualitative agreement isclear (although there are sizable quantitative differences, which we attributeto the somewhat different form of the on-site interaction Hamiltonian Hˆe).As the coupling to the lattice is turned up, it is immediately obvious thatthe charge modulated ⇑ 0 ⇓ 0 order benefits extraordinarily from it: itexpands in reach down to J = 0 in some regions, suppressing the equal-sublattice order ↑↑↓↓ entirely. It also expands across the board, taking overabout half of what used to be the metallic phase, and making large gainsat the expense of the ferromagnetic phase as well. The charge modulation,however, expands even independently of its associated magnetic order, as isevidenced by the relative growth of the aligned-ferrimagnetic phase ⇑↑⇑↑.A curious feature emerges at the leftmost boundary between the metallicand the ⇑ 0 ⇓ 0 order directly on the U = 0 axis: a small phase pocketopens up, with the ground state there exhibiting charge order without theusually associated magnetic order ⇑ 0 ⇓ 0. In effect, the lattice distortionhelps decouple the charge and orbital degrees of freedom, which is stronglyreminiscent of the behaviour of the heavier rare-earth ions (Sm to Lu) inthe traditional nickelate phase diagram. This could be evidence for thedominant role of the lattice coupling during the MIT, and suggests thatthe strength of the lattice interaction might be responsible for determiningwhether the MIT and magnetic transitions are concurrent or not. Lookingback at the metal-insulator diagrams in Fig. 5.4, we notice that despiteacquiring charge density wave character, the ground state for small U is stillmetallic: despite the strong gapping out of the spectrum and flattening ofthe bands by the lattice distortion, the Fermi level stays well within theoccupied band, so there is no electron localization and magnetic order doesnot arise.If we fix charge modulation to the experimentally relevant range δ ∼0.3-0.4, we notice that most of the time these contours follow the phaselines between the ⇑ 0 ⇓ 0 phase and some other phase (metallic for Uup to 1.5, ferromagnetic after U ∼ 4). They, of course, stay on the ⇑0 ⇓ 0 side. However, right in the middle region (the more experimentally705.2. Magnetic order: types and the transitionFigure 5.6: Magnetic order in the HF ground state, as a function of U andb, for J/U = 0.2, 0.3 ((a) and (b), respectively) and 0.4, 0.5, ((c) and (d),respectively). Other parameters are t1 = 1, t2 = 0.15, t4 = 0.25, α = 1.reasonable values of U ∼ 1.5-3) these contours fall squarely into the bulk ofthe ⇑ 0 ⇓ 0 phase, suggesting that this indeed could be the dominant 4-sitelocal magnetic alignment for the real system. And this dominance is onlyenhanced by the coupling to the lattice, as the region of the overall phasediagram occupied by this state grows.Zooming in once again to specific values for the ratio J/U , we can in-terrogate the impact of the lattice more closely. Just as we argued fromthe large-scale diagrams of Fig. 5.5, it is clear from Fig. 5.6 that with theincrease of b the phase line for the ⇑ 0 ⇓ 0 order creeps further and fur-ther to the left towards smaller U and J values: again the effect is morepronounced for larger ratios J/U . The suppression of the other nontrivial4-site alignments is also evident with increasing lattice coupling: eventuallythe only possibilities for the ground state become a charge density wave,with or without magnetic order.An even more explicit picture of the preference of the lattice for the⇑ 0 ⇓ 0 type order can be seen in Fig. 5.7. Here we plot the energies of thevarious states at a particular set of parameter values relative to the metallicsolution 0000, as we sweep b: ⇑ 0 ⇓ 0 emerges a clear winner (and also715.2. Magnetic order: types and the transitionFigure 5.7: Energies of the various converged states (y axis) as a function ofb (x axis). The parameters are U = 3.158, J = 0.2U, t1 = 1, t2 = 0.15, t4 =0.25. Resolution is 20× 20.the only one responsive to the presence of the lattice). While this plot ismade at a specific set of parameter values, this pattern seems to hold acrossvarious parameter regimes of U, J, ti, as is evidenced by the magnetic phasediagrams where ⇑ 0 ⇓ 0 phase conquers an ever-larger area with increasingb. What is curious about these plots at a first glance is that relative tothe 0000 state, none of the states’ energies are affected by changing b,aside from ⇑ 0 ⇓ 0. But there is a straightforward explanation to this:any state that includes no charge disproportionation also includes no latticedistortion, as δ ∼ u through the lattice self-consistency equation 3.5. Thuschanging b has no effect neither in the total energy nor in the Hartree-FockHamiltonian, as it always enters in the combination bu. The real questionis: why are the states with partial charge disproportionation, e.g. ⇑↑⇓↓not convergent within this mean-field model? Such states could potentiallycompete with ⇑ 0 ⇓ 0 order for ground state status, as the coupling to thelattice is cranked higher. It is a fact, however, that such states are not onlynot convergent during the iterations no matter the starting point, but infact they are among the least stable (this is discussed more in Sec. 5.5).725.3. Impact of the hopping integrals5.3 Impact of the hopping integralsOne clear impact of adjusting the hopping integrals ti is to affect the overallbandwidth W , defined in the usual way as the difference between the largestand the smallest Hartree-Fock eigenvalues at zero interaction (U = J = b =0). We know from numerous studies of the single-band Hubbard model thatthe parameter that results in the onset of Mott insulator phase is the ratioof the interaction strength to the bandwidth, U/t (U/W with W the totalbandwidth if more than one hopping pathway is present): thus one couldanticipate changes to the phase diagram when any of the hopping rates tiincrease enough so as to change W significantly.On experimental grounds, however, given the spatial extent of the nick-elate orbitals, it is unrealistic to expect the farther-neighbour hopping t2, t4to be larger or even same magnitude as the nearest neighbour t1. The anal-ysis of Lee et al. put the rate t2 somewhere in the range ≈ 0.05-0.15t1for best agreement with LDA calculations and Fermi surface measurements[26]: thus we can reasonably restrict t2, t4 < 0.3. However, the bandwidthchanges by at most %10 for this range of values, which should lead to onlymodest shifts of the phases in the diagram.Yet we observe much more dramatic effects than could be expected on thebasis of bandwidth renormalization. Looking to Fig. 5.8, where a magneticphase diagram in the U–J plane is produced for varying levels of the second-nearest neighbour hopping (t2 = 0, 0.1, 0.25), we notice a profound changeof the phase boundaries. Moreover, this change favours ferromagnetism, inspite of the increasing bandwidth, which would be expected to push the insu-lating boundary to larger U . Looking to Fig. 5.10 with t4 = 0, 0.1, 0.25, 0.35,the opposite effect is observed: a change for the metallic region that is waybeyond the meager %10-%15 one could expect based on the U/W Mott pic-ture. Instead, it appears that the shape, rather than the width, of the DOSis to blame: both t2 and t4 hopping pathways affect the lop-sidedness of thedistribution strongly even when the overall width is barely renormalized, aphenomenon noted in the literature on Hubbard models [79].Both the t1 and t4 hopping pathways appear to produce a symmetric(about the midpoint) DOS: this is not entirely surprising, given that t4 hasexactly the same hopping coefficients, just with the lattice constant being2a. The t2 hopping, on the other hand, is different, producing an asymmetricDOS: even a small positive t2 value leads to a pile-up of DOS weight towardsthe lower band edge (see Fig. 5.9). This extra weight, known as a van Hovesingularity, increases the DOS at the Fermi level (given that at quarter-filling the Fermi level is in the lower quartile of the distribution). One735.3. Impact of the hopping integralsFigure 5.8: Magnetic order in the HF ground state, as a function of U andJ , for t2 = 0, 0.15, 0.25 (top to bottom). Other parameters are t1 = 1, t4 =0, b = 0.745.3. Impact of the hopping integralsFigure 5.9: The non-interacting (U = J = 0) DOS, corresponding to thesystems in Fig.5.8. Notice that increasing the bandwidth (even if just mod-estly by at most %15), paradoxically, leads to more robust ferromagnetismat lower U — a consequence of the van Hove singularity at the lower bandedge. Also notice how when the next-nearest neighbour frustration is max-imally reduced (t2 = 0), the non-collinear 4-site magnetic order dominatesthe collinear one. Currently it is not clear to us why this would be the case,given how introducing t2, t4 seems to affect them equally based on purelattice frustration arguments.755.3. Impact of the hopping integralsFigure 5.10: Magnetic order in the HF ground state, as a function of Uand J , for t4 = 0, 0.1 ((a) and (b), respectively), and t4 = 0.25, 0.35 ((c)and (d)). The bandwidth is W = 6.4 in all cases. Other parameters aret1 = 1, t2 = 0.15, b = 0. The growth of the metallic region is clearly notthe effect of a renormalized bandwidth, but rather due to the changes of theshape of the DOS.765.3. Impact of the hopping integralswell-known pathway to ferromagnetism is the Stoner criterion [80]: derivedwithin the Hubbard model, it says that ferromagnetism can spontaneouslyarise in an electronic system when the condition UD(EF) > 1 is met (D(EF)is the density of states at the Fermi level). This condition is derived for adifferent Hubbard model than the one employed here: however, it appearsto be consistent with that we observe. The increase/decrease in the densityof states at the Fermi level with the changes of t2, t4 is on the order of%30-40, much more commensurate with the observed shifts of the phaseboundaries in Figs. 5.8,5.10. In particular, t2 acts to increase the DOS atthe Fermi level: this means U (or some effective measure of the electron-electron interactions) at which ferromagnetism sets in can be significantlylower. Conversely, increasing t4 has the effect of boosting the DOS abovethe Fermi level, thus depleting it at the Fermi level and increasing the Urequired to set off ferromagnetism. This kind of behaviour is consistentwith what has been seen before with regards to the next-nearest neighbourhopping [79].Despite the seemingly strong impact of the hopping rates on the magneticphase diagram, it appears that it almost does not affect the energy of the4-site states. Fig. 5.11 shows how the energies of the three main 4-sitecontenders ⇑ 0 ⇓ 0, ↑↑↓↓ and ↑→↓← scale with varying t2 and t4. Theseplots represent the trend that is observed across most of the parameterregimes: that the hopping rates generally do not affect the energy orderof the various 4-site magnetic states, rather they modulate their energiesmore or less equally. This can be generally understood from frustration costarguments: the added frustration introduced by a nonzero t4 is identicalfor all the three candidates. For t2 the effect is slightly more complicated,as increasing it actually brings down the energy for all of the magneticstates: we believe this should be seen in light of the Stoner criterion whichbenefits all kinds of magnetic order, but prefers ferromagnetism above allelse. Unlike t4, increasing t2 actually does show a preference for magneticorders: namely it prefers collinear orders like ⇑ 0 ⇓ 0 and ↑↑↓↓ above thenon-collinear version ↑→↓← – however, its preference for ferromagnetismreigns above all else. Conversely, the frustration costs from t4 are strongestin ferromagnets, as can be seen by the steep increase in energy for theferromagnetic state in Fig. 5.11. Once again, we would like to point outthat these effects are all purely a consequence of the shape of the DOS, asthe bandwidth is barely renormalized (see Fig. 5.8).Another small caveat to the notion that the relative standing of thevarious 4-site candidates is unaffected by the various hopping rates is theemergence of the non-collinear ↑→↓← state as the ground state in what775.4. The anharmonicity parameter αused to be the ↑↑↓↓ domain at t2 = 0, visible in the top panel in Fig. 5.8.This phenomenon is still not entirely clear to us and is to be studied furtherin a later work: however, the fact that it only happens when t2 = 0 andis washed out at all other values of the hopping rates suggests that it isan isolated phenomenon and should not affect the general discussion of theground state in rare-earth nickelates.5.4 The anharmonicity parameter αThroughout this work the value α = 1 was mostly used during the calcu-lations: however, it is certainly important to reflect on the impact of thisadditional parameter in our model. Looking at the form of the lattice energyfor the typical value of the distortion u = 0.5,Elatt = b(u2i + αu4i /2) ∼ b(0.25 + 0.03α),we observe that the anharmonic contribution for α = 1 is merely %10 of thetotal value of the lattice energy. Hence its main impact is to renormalizethe ground state distortion value u away from the identity u = δ, as per theself-consistency equation 3.5,uj + au3j + 1 = 〈nj〉Ψe ,without affecting the physics of the problem. An example calculation fora fixed ratio J/U = 0.2 in Fig. 5.12 supports this intuition: there is onlya moderate displacement of the phase boundaries with the anharmonicityparameter going from α = 0 to α = 1, with the main change being theshift of the charge modulation δ contours. And it is mostly the contours forlargest δ (i.e. largest u) that are significantly displaced, as could be expectedfrom the energy scale comparisons above. In the interest of controlling themagnitude of the lattice distortion u, we employ α = 1 throughout most ofthe calculation in this thesis.5.5 Convergence and stabilityIt should be noted that many of these magnetic orders, excluding perhapsthe ⇑ 0 ⇓ 0 state, are very unstable numerically. Consider the no dispropor-tionation collinear state ↑↑↓↓. In our model it will have (for example) thefollowing mean-field description,δ, SFM, SAFM, S1x, S2x = 0, S1z = S2z 6= 0.785.5. Convergence and stabilityFigure 5.11: Energies of several converged self-consistent HF states as afunction of t4, relative to the energy of the metallic state. The differ-ent colours correspond to different magnetic orders (see the legend). Theparameters are U = 4, J = 0.316, t1 = 1, t4 = 0, b = 0 for (a) andU = 5, J = 0.316, t1 = 1, t2 = 0.15, b = 0 for (b). Notice how in (b)the relative energies of the chief magnetic ground state contenders are notaffected by the change in the hopping rate t4 – except for ferromagnetism,which gets strongly frustrated with the additional t4 hopping and, paradoxi-cally, “unfrustrated” with the introduction of t2 hopping due to DOS effects(see the text for details). Meanwhile, in (a) with tuning the t2 rate awayfrom 0 the non-collinear ↑→↓← gets briefly favoured, but then quickly losesout to the collinear ↑↑↓↓, before ferromagnetism begins to reign supreme.795.5. Convergence and stabilityFigure 5.12: Magnetic order in the HF ground state, as a function of U andb, for J/U = 0.2, and two different values of the anharmonicity α: α = 0(top) and α = 1 (bottom). Other parameters are t1 = 1, t2 = 0.15, t4 = 0.25.805.5. Convergence and stabilityTo make a reasonable guess for its convergence, we might guess S1z = S2z =0.5, and all other parameters zero. Then we usually get the desired state.However, if there is even a slight deviation from either perfect symmetrybetween the sublattices (S1z 6= S2z), or a small non-zero δ, or a differentparameter (O1, Z4, SFM, ...) being nonzero, most likely the iterations wouldrun away from the solution of the desired type. Sometimes there are linesof stability, along which deviations generally do not result in catastrophicrunaway, but there never are any sizable convergence basins for most of thestates on the list 5.2. In fact, iteration divergence occured (before the in-troduction of the Pulay mixing method) even if the iterations simply tooktoo long: the accumulated numerical errors after a sufficiently large numberof iterations automatically result in some of the other parameters becom-ing non-zero, which causes runaway effects and the desired solution is notobtained. This can be seen most clearly from cuts in the 15-dimensionalparameter space that show iterative trajectories of the Hartree-Fock cal-culation (such diagrams are called “Poincare sections” in the dynamicalsystems literature). In Fig. 5.13, we plot the difference S1z−S2x versus thecharge disproportionation δ during the iterative process, as the algorithmsearches for the solution to the Hartree-Fock equations. Both ↑→↓← and⇑ 0 ⇓ 0 are well-defined points in this plane, corresponding to S1z = S2z andone of the S being zero, respectively. The intermediate state ⇑→⇓← wouldthen lie somewhere on the diagonal, between the two. From the iterativetrajectories, we can see that any solution not already on the S1z − S2z = 0line (the y-axis) converges to ⇑ 0 ⇓ 0 instead: and the intermediate typestates on the diagonal have some of the fastest decay toward the ⇑ 0 ⇓ 0states out there.815.5. Convergence and stabilityFigure 5.13: A phase portrait of the iterative sequences for various startingparameters in the δ, S1z, S2z parameter subspace. On the x axis we plotthe difference between the spin parameters of the sublattices, S1z − S2z:hence the perfect symmetry point is at the origin. On the y axis is thedisproportionation parameter δ: the stable point ⇑ 0 ⇓ 0 is thus at thetop right and equivalently bottom left. The parameters are U = 1.5, J =3, t1 = 1, t2 = 0.15, t4 = 0.55, b = 0. Despite the relatively dense samplingof the phase space, the ⇑ 0 ⇓ attractors dominate the dynamics, with noother states, specifically of the partially disproportionated kind ⇑↑⇓↓, beingconvergent.82Chapter 6ConclusionsWhat we call the beginning isoften the end. And to make anend is to make a beginning. Theend is where we start from.T. S. Eliot (1942)Rare-earth nickelates, ever since being synthesized in the early 70s, havepresented us with numerous fascinating puzzles that speak to various areasof condensed matter physics. While the questions about its crystal struc-ture and the lattice distortion across the MIT phase boundary are consideredlargely settled, the conversations surrounding the MIT itself – including itsprecise mechanism and physical manifestation vis-a´-vis the orbital structurein the Ni-O cages – as well as questions about the nature of the magneticorder at low temperatures and its interrelation to charge and lattice degreesof freedom, are still ongoing. In this thesis we focussed our efforts on in-vestigating the relationship between the lattice distortions in the materialand the charge/magnetic order found at low temperatures. Our approachwas to construct an effective two-band Hubbard model, based on the phe-nomenology of the material that is known from experiments, that focusseson the nominally degenerate eg doublet of the Ni outer electron subshell.We included first, second and fourth neighbour hopping, and incorporatedthe electron-electron interactions in the two-band context in the usual way– by employing the Kanamori Hamiltonian. The lattice was introducedsemiclassically, characterized by a single scalar u representing the largelyisotropic distortion of the Ni-O octahedral cage: it was coupled to the elec-trons through a Holstein-like on-site term. The final model ends up hav-ing 7 independent adjustable parameters: the electron-electron interactionstrengths U, J , the hopping amplitudes t1, t2, t4, the lattice energy scale band anharmonicity α.Upon constructing the model, we applied the Hartree-Fock approxima-tion and solved the resulting Hartree-Fock equations to obtain the expectedground state of the model at zero temperature. Studying the ground state83Chapter 6. Conclusionsin a variety of slices of the high-dimensional parameter space, we were ableto confirm that the net effect of the electron-phonon coupling on the elec-tronic behaviour – even in this simple semiclassical model for the lattice –is to encourage the formation of insulating charge order, just as found byearlier investigators [32]. With the increase of the lattice coupling strengthb, we found that the area of phase diagrams occupied by the phase withcharge order increases dramatically, as the order becomes more energeticallyfavourable: in this simple model the charge order is most closely identifiedwith that of a charge density wave type Ni3+δNi3−δ, as described in Sec.1.4.3, although effective parallels can be drawn with the negative-chargetransfer picture as well, as discussed briefly during the construction of themodel in Ch. 2.Next we proceeded to interrogate the magnetic behaviour of the nicke-lates within our model. It turns out there is a great variety of states that areself-consistent within the model: the relevant ones are those that are lowestin total energy. We identified several 4-site self-consistent candidate statesthat could fit the experimentally observed wavevector Qm = (2pia )(14 ,14 ,14):these include ⇑ 0 ⇓ 0, ↑↑↓↓, and ↑→↓←. Even without the coupling to thelattice, for most hoppings ti and a wide range of the electron-electron inter-action parameters U, J , the ⇑ 0 ⇓ 0 state emerges as the dominant groundstate. The two other candidates only arise at intermediate U and small Jand zero lattice coupling b, and do not fare well when the lattice couplingis increased, quickly ceding their entire phase space to the ⇑ 0 ⇓ 0. Thissort of behaviour makes sense in light of the fact that the states ↑↑↓↓, and↑→↓← cannot couple to the lattice distortion, given that they have perfectsublattice symmetry and no charge disproportionation: and if δ = 0, thenthe self-consistency condition for the lattice degree of freedom demands thelattice distortion also vanish u = 0. So while these two states are unaffectedby the presence of the lattice, the state ⇑ 0 ⇓ 0 can couple to the distortionsand lower its energy, thus becoming the favoured ground state.In principle one could imagine states with only partial charge dispropor-tionation, say of type ⇑→⇓←, emerging within the model to compete ener-getically with the fully disproportionated state ⇑ 0 ⇓ 0. However, that doesnot seem to happen within our model: we never find such states among theself-consistent solutions. In fact, if we follow the iterative trajectories dur-ing the solution of the Hartree-Fock equations and start with initial guessesclose to such a partially disproportionated order, such points appear amongthe least stable in the phase space, the solution quickly diverging away eitherto the metallic 0000 state or to the ⇑ 0 ⇓ 0 one.During this study, an unexpected new feature emerged in the phase dia-84Chapter 6. Conclusionsgram when the lattice coupling b was turned on: a novel phase wherein thecharge density wave (δ 6= 0) appeared without the accompanying magneticorder ⇑ 0 ⇓ 0. This de-coupling of the charge and magnetic transitions inthe phase diagram appears to be mediated by the electron-phonon coupling:the energy decrease afforded by the charge disproportionation is reinforcedby the lattice distortion so far so as to make the additional appearance ofthe magnetic order (and simultaneously the onset of insulating behaviour)unfavourable. This feature is strongly reminiscent of the well-known phe-nomenon in the nickelates where the MIT and the magnetic transition oc-cur at different temperatures, TMIT > TN, for all the nickelates from Euto Lu, leaving out only the lightest rare-earth ions Pr and Nd, for whichTMIT = TN. The observation that this decoupling can be brought on by theelectron-phonon coupling even in this simple model further underscores theimportance of the lattice degrees of freedom and could potentially be seen asan explanation of the phenomenon, given that the lattice coupling strengthcan be expected to increase for the heavier rare-earth ions.However, plenty of open questions remain. For one, the negative chargetransfer picture of the MIT in the nickelates calls for the involvement of theoxygen 2p bands in the model: it is unclear to what extent the variationof the hopping and interaction parameters can capture the effects of intro-ducing the oxygen bands. Moreover, the Hartree-Fock approximation doesnot take into account any correlation effects (aside from Pauli’s exclusionprinciple): while it is generally believed that the effects of correlation inthe nickelates are not very important, it would still be useful to check thisexplicitly, e.g. via a dynamical mean-field theory (DMFT) approach.Aside from the questions about the applicability of the model and theaccuracy of the Hartree-Fock approximation, even within our study there area few ideas that can be pursued further. For instance, a clear understand-ing of the orbitally ordered state (the antiferromagnetic disproportionated-ferroorbital state) that obtains at intermediate to high values of U in themagnetic order phase diagrams would be very valuable, along with otherorbitally ordered states, especially considering that there have been no signsof any orbital order in most experimental studies.Overall, while we have made significant strides in addressing the questionof the impact of the electron-phonon coupling on the charge and magneticorder in the nickelates, further work, using more detailed Hamiltonians andmore accurate approximations, is required before the matter is fully put torest.85Bibliography[1] J. Varignon, M. N. Grisolia, J. I´n˜iguez, A. Barthe´le´my, and M. Bibes,“Complete phase diagram of rare-earth nickelates from first-principles,”npj Quantum Materials, vol. 2, p. 21, December 2017.[2] A. Hampel, P. Liu, C. Franchini, and C. Ederer, “Energetics of thecoupled electronic–structural transition in the rare-earth nickelates,”npj Quantum Materials, vol. 4, p. 5, December 2019.[3] M. L. Medarde, “Structural, magnetic and electronic properties of per-ovskites (R = rare earth),” Journal of Physics: Condensed Matter,vol. 9, pp. 1679–1707, February 1997.[4] G. Catalan, “Progress in perovskite nickelate research,” Phase Transi-tions, vol. 81, pp. 729–749, July 2008.[5] P. V. Balachandran and J. M. Rondinelli, “Interplay of octahedral ro-tations and breathing distortions in charge-ordering perovskite oxides,”Physical Review B, vol. 88, p. 054101, August 2013.[6] X. Granados, J. Fontcuberta, X. Obradors, and J. B. Torrance,“Metastable metallic state and hysteresis below the metal-insulatortransition in PrNiO3,” Physical Review B, vol. 46, pp. 15683–15688,December 1992.[7] E. Pavarini, “Crystal-Field Theory, Tight-Binding Method, and Jahn-Teller Effect,” in Correlated Electrons: From Models to Materials(E. Pavarini, E. Koch, F. Anders, and M. Jarrell, eds.), ch. 6, pp. 147–187, Forschungszentrum Ju¨lich: Forschungszentrum Ju¨lich GmbH In-stitute for Advanced Simulation, 2012.[8] G. Sawatzky and R. Green, “The Explicit Role of Anion States in High-Valence Metal Oxides,” in Quantum Materials: Experiments and The-ory (E. Pavarini, E. Koch, J. Van Den Brink, and G. Sawatzky, eds.),vol. 6, ch. 1, pp. 1–36, Forschungszentrum Ju¨lich: ForschungszentrumJu¨lich GmbH Institute for Advanced Simulation, 2016.86Bibliography[9] A. Mun˜oz, J. Alonso, M. Mart´ınez-Lope, and M. Ferna´ndez-Dı´az, “Onthe magnetic structure of DyNiO3,” Journal of Solid State Chemistry,vol. 182, pp. 1982–1989, July 2009.[10] I. B. Bersuker, Electronic Structure and Properties of Transition MetalCompounds. Hoboken, NJ, USA: John Wiley & Sons, Inc., April 2010.[11] N. A. Spaldin, “Multiferroics: Past, present, and future,” MRS Bul-letin, vol. 42, pp. 385–390, May 2017.[12] G. Giovannetti, S. Kumar, D. Khomskii, S. Picozzi, and J. van denBrink, “Multiferroicity in Rare-Earth Nickelates RNiO3,” Physical Re-view Letters, vol. 103, p. 156401, October 2009.[13] J. F. Scott, “Multiferroic memories,” Nature Materials, vol. 6, pp. 256–257, April 2007.[14] G. Demazeau, A. Marbeuf, M. Pouchard, and P. Hagenmuller, “Sur unese´rie de compose´s oxyge`nes du nickel trivalent derive´s de la perovskite,”Journal of Solid State Chemistry, vol. 3, pp. 582–589, November 1971.[15] J. G. Bednorz and K. A. Mu¨ller, “Possible high Tc superconductivity inthe Ba-La-Cu-O system,” Zeitschrift fu¨r Physik B Condensed Matter,vol. 64, pp. 189–193, June 1986.[16] J. Torrance, P. Lacorre, A. Nazzal, E. Ansaldo, and C. Niedermayer,“Systematic study of insulator-metal transitions in perovskites RNiO3(R=Pr,Nd,Sm,Eu) due to closing of charge-transfer gap,” Physical Re-view B, vol. 45, pp. 8209–8212, April 1992.[17] J. B. Torrance, P. Lacorro, C. Asavaroengchai, and R. M. Metzger,“Simple and perovskite oxides of transition-metals: Why some aremetallic, while most are insulating,” Journal of Solid State Chemistry,vol. 90, pp. 168–172, January 1991.[18] T. Mizokawa, D. I. Khomskii, and G. A. Sawatzky, “Spin and chargeordering in self-doped Mott insulators,” Physical Review B, vol. 61,pp. 11263–11266, May 2000.[19] H. Park, A. J. Millis, and C. A. Marianetti, “Site-Selective MottTransition in Rare-Earth-Element Nickelates,” Physical Review Letters,vol. 109, p. 156402, October 2012.87Bibliography[20] B. Lau and A. J. Millis, “Theory of the Magnetic and Metal-InsulatorTransitions in RNiO3 Bulk and Layered Structures,” Physical ReviewLetters, vol. 110, p. 126404, March 2013.[21] D. Puggioni, A. Filippetti, and V. Fiorentini, “Ordering and multiplephase transitions in ultrathin nickelate superlattices,” Physical ReviewB, vol. 86, p. 195132, November 2012.[22] A. D. Caviglia, R. Scherwitzl, P. Popovich, W. Hu, H. Bromberger,R. Singla, M. Mitrano, M. C. Hoffmann, S. Kaiser, P. Zubko,S. Gariglio, J.-M. Triscone, M. Fo¨rst, and A. Cavalleri, “Ultrafast StrainEngineering in Complex Oxide Heterostructures,” Physical Review Let-ters, vol. 108, p. 136801, March 2012.[23] J. L. Garc´ıa-Mun˜oz, J. Rodr´ıguez-Carvajal, P. Lacorre, and J. B. Tor-rance, “Neutron-diffraction study of RNiO3 (R = La,Pr,Nd,Sm): Elec-tronically induced structural changes across the metal-insulator transi-tion,” Physical Review B, vol. 46, pp. 4414–4425, August 1992.[24] H. Guo, Z. W. Li, L. Zhao, Z. Hu, C. F. Chang, C.-Y. Kuo, W. Schmidt,A. Piovano, T. W. Pi, O. Sobolev, D. I. Khomskii, L. H. Tjeng, andA. C. Komarek, “Antiferromagnetic correlations in the metallic stronglycorrelated transition metal oxide LaNiO3,” Nature Communications,vol. 9, p. 43, December 2018.[25] K. Haule and G. L. Pascut, “Mott Transition and Magnetism in RareEarth Nickelates and its Fingerprint on the X-ray Scattering,” ScientificReports, vol. 7, p. 10375, December 2017.[26] S. Lee, R. Chen, and L. Balents, “Metal-insulator transition in a two-band model for the perovskite nickelates,” Physical Review B, vol. 84,2011.[27] V. Scagnoli, U. Staub, A. M. Mulders, M. Janousch, G. I. Meijer,G. Hammerl, J. M. Tonnerre, and N. Stojic, “Role of magnetic andorbital ordering at the metal-insulator transition in NdNiO3,” PhysicalReview B, vol. 73, p. 100409, March 2006.[28] Y. Lu, D. Betto, K. Fu¨rsich, H. Suzuki, H.-H. Kim, G. Cristiani,G. Logvenov, N. B. Brookes, E. Benckiser, M. W. Haverkort, G. Khal-iullin, M. Le Tacon, M. Minola, and B. Keimer, “Site-Selective Probeof Magnetic Excitations in Rare-Earth Nickelates Using Resonant In-elastic X-ray Scattering,” Physical Review X, vol. 8, 2018.88Bibliography[29] X. Obradors, L. M. Paulius, M. B. Maple, J. B. Torrance, A. I. Nazzal,J. Fontcuberta, and X. Granados, “Pressure dependence of the metal-insulator transition in the charge-transfer oxides RNiO3 (R = Pr,Nd,Nd0.7 La0.3 ),” Physical Review B, vol. 47, pp. 12353–12356, May 1993.[30] M. Hepting, Ordering Phenomena in Rare-Earth Nickelate Heterostruc-tures. Springer Theses, Cham: Springer International Publishing, 2017.[31] M. Medarde, P. Lacorre, K. Conder, F. Fauth, and A. Furrer, “Giant16O-18O Isotope Effect on the Metal-Insulator Transition of RNiO3 Per-ovskites (R = Rare Earth),” Physical Review Letters, vol. 80, pp. 2397–2400, March 1998.[32] S. Johnston, A. Mukherjee, I. Elfimov, M. Berciu, and G. A. Sawatzky,“Charge Disproportionation without Charge Transfer in the Rare-Earth-Element Nickelates as a Possible Mechanism for the Metal-Insulator Transition,” Physical Review Letters, vol. 112, p. 106404,March 2014.[33] A. Wold, B. Post, and E. Banks, “Rare Earth Nickel Oxides,” Journalof the American Chemical Society, vol. 79, pp. 4911–4913, September1957.[34] V. M. Goldschmidt, “Die Gesetze der Krystallochemie,” Die Naturwis-senschaften, vol. 14, pp. 477–485, May 1926.[35] J. A. Alonso, M. J. Mart´ınez-Lope, M. T. Casais, J. L. Garc´ıa-Mun˜oz,and M. T. Ferna´ndez-Dı´az, “Room-temperature monoclinic distortiondue to charge disproportionation in RNiO3 perovskites with small rare-earth cations (R = Ho, Y, Er, Tm, Yb and Lu): A neutron diffractionstudy,” Physical Review B, vol. 61, pp. 1756–1763, January 2000.[36] M. Medarde, M. T. Ferna´ndez-Dı´az, and P. Lacorre, “Long-rangecharge order in the low-temperature insulating phase of PrNiO3,” Phys-ical Review B, vol. 78, p. 212101, December 2008.[37] J. L. Garc´ıa-Mun˜oz, M. A. G. Aranda, J. A. Alonso, and M. J.Mart´ınez-Lope, “Structure and charge order in the antiferromag-netic band-insulating phase of NdNiO3,” Physical Review B, vol. 79,p. 134432, April 2009.[38] J. Rodr´ıguez-Carvajal, S. Rosenkranz, M. Medarde, P. Lacorre, M. T.Fernandez-Dı´az, F. Fauth, and V. Trounov, “Neutron-diffraction study89Bibliographyof the magnetic and orbital ordering in SmNiO3,” Physical Review B,vol. 57, pp. 456–464, January 1998.[39] A. Hampel and C. Ederer, “Interplay between breathing mode dis-tortion and magnetic order in rare-earth nickelates RNiO3 withinDFT+U ,” Physical Review B, vol. 96, p. 165130, October 2017.[40] J. M. Perez-Mato, D. Orobengoa, and M. I. Aroyo, “Mode crystallog-raphy of distorted structures,” Acta Crystallographica Section A Foun-dations of Crystallography, vol. 66, pp. 558–590, September 2010.[41] J.-S. Zhou and J. B. Goodenough, “Chemical bonding and electronicstructure of RNiO3 (R = rare earth),” Physical Review B, vol. 69,p. 153105, April 2004.[42] J. A. Alonso, J. L. Garc´ıa-Mun˜oz, M. T. Ferna´ndez-Dı´az, M. A. G.Aranda, M. J. Mart´ınez-Lope, and M. T. Casais, “Charge Dispro-portionation in RNiO3 Perovskites: Simultaneous Metal-Insulator andStructural Transition in YNiO3,” Physical Review Letters, vol. 82,pp. 3871–3874, May 1999.[43] U. Staub, G. I. Meijer, F. Fauth, R. Allenspach, J. G. Bednorz,J. Karpinski, S. M. Kazakov, L. Paolasini, and F. D’Acapito, “DirectObservation of Charge Order in an Epitaxial NdNiO3 Film,” PhysicalReview Letters, vol. 88, p. 126402, March 2002.[44] C. Johan, Introduction to ligand field theory. New York: McGraw-Hill,1962.[45] J. Zaanen, G. A. Sawatzky, and J. W. Allen, “Band gaps and elec-tronic structure of transition-metal compounds,” Physical Review Let-ters, vol. 55, pp. 418–421, July 1985.[46] N. F. Mott, “The Basis of the Electron Theory of Metals, with Spe-cial Reference to the Transition Metals,” Proceedings of the PhysicalSociety. Section A, vol. 62, pp. 416–422, July 1949.[47] J. Hubbard, “Electron correlations in narrow energy bands,” Proceed-ings of the Royal Society of London. Series A. Mathematical and Phys-ical Sciences, vol. 276, pp. 238–257, November 1963.[48] I. I. Mazin, D. I. Khomskii, R. Lengsdorf, J. A. Alonso, W. G. Marshall,R. M. Ibberson, A. Podlesnyak, M. J. Mart´ınez-Lope, and M. M. Abd-90BibliographyElmeguid, “Charge Ordering as Alternative to Jahn-Teller Distortion,”Physical Review Letters, vol. 98, p. 176406, April 2007.[49] J. L. Garc´ıa-Mun˜oz, J. Rodr´ıguez-Carvajal, and P. Lacorre, “Neutron-diffraction study of the magnetic ordering in the insulating regime ofthe perovskites RNiO3 (R = Pr and Nd),” Physical Review B, vol. 50,pp. 978–992, July 1994.[50] J. A. Alonso, M. J. Mart´ınez-Lope, I. A. Presniakov, A. V. Sobolev,V. S. Rusakov, A. M. Gapochka, G. Demazeau, and M. T. Ferna´ndez-Dı´az, “Charge disproportionation in RNiO3 (R = Tm, Yb) perovskitesobserved in situ by neutron diffraction and 57Fe probe Mo¨ssbauer spec-troscopy,” Physical Review B, vol. 87, p. 184111, 2013.[51] A. Subedi, O. E. Peil, and A. Georges, “Low-energy description of themetal-insulator transition in the rare-earth nickelates,” Physical ReviewB, vol. 91, p. 075128, February 2015.[52] D. J. Gawryluk, J. Rodr´ıguez-Carvajal, P. Lacorre, M. T. Ferna´ndez-Dı´az, and M. Medarde, “Distortion mode anomalies in bulk PrNiO3,”arXiv e-prints, p. arXiv:1809.10914, Sep 2018.[53] J.-G. Cheng, J.-S. Zhou, J. B. Goodenough, J. A. Alonso, and M. J.Martinez-Lope, “Pressure dependence of metal-insulator transition inperovskites RNiO3 ( R = Eu, Y, Lu),” Physical Review B, vol. 82,p. 085107, August 2010.[54] N. Wagner, D. Puggioni, and J. M. Rondinelli, “Learning from Cor-relations Based on Local Structure: Rare-Earth Nickelates Revisited,”Journal of Chemical Information and Modeling, vol. 58, pp. 2491–2501,December 2018.[55] M. Medarde, C. Dallera, M. Grioni, B. Delley, F. Vernay, J. Mesot,M. Sikora, J. A. Alonso, and M. J. Mart´ınez-Lope, “Charge dispropor-tionation in RNiO3 perovskites (R = rare earth),” Physical Review B,vol. 80, p. 245105, December 2009.[56] Y. Bodenthin, U. Staub, C. Piamonteze, M. Garc´ıa-Ferna´ndez, M. J.Mart´ınez-Lope, and J. A. Alonso, “Magnetic and electronic propertiesof RNiO3 (R = Pr, Nd, Eu, Ho and Y) perovskites studied by resonantsoft x-ray magnetic powder diffraction,” Journal of Physics: CondensedMatter, vol. 23, p. 036002, January 2011.91Bibliography[57] J. Garc´ıa, J. Blasco, M. G. Proietti, and M. Benfatto, “Analysis ofthe x-ray-absorption near-edge-structure spectra of La1−xNdxNiO3 andLaNi1−xFexO3 at the nickel K edge,” Physical Review B, vol. 52,pp. 15823–15828, December 1995.[58] X. Q. Xu, J. L. Peng, Z. Y. Li, H. L. Ju, and R. L. Greene, “Resisitiv-ity, thermopower, and susceptibility of RNiO3 (R = La,Pr),” PhysicalReview B, vol. 48, pp. 1112–1118, July 1993.[59] J. L. Garc´ıa-Mun˜oz, J. Rodr´ıguez-Carvajal, and P. Lacorre, “SuddenAppearance of an Unusual Spin Density Wave At the Metal-InsulatorTransition in the Perovskites RNiO3 (R = Pr, Nd),” Europhysics Let-ters (EPL), vol. 20, pp. 241–247, October 1992.[60] C. Castellani, C. R. Natoli, and J. Ranninger, “Magnetic structure ofV2O3,” Physical Review B, vol. 18, pp. 4945–4966, November 1978.[61] M. T. Ferna´ndez-Dı´az, J. A. Alonso, M. J. Mart´ınez-Lope, M. T. Ca-sais, and J. L. Garc´ıa-Mun˜oz, “Magnetic structure of the HoNiO3 per-ovskite,” Physical Review B, vol. 64, p. 144417, September 2001.[62] M. Cyrot and C. Lyon-Caen, “Orbital superlattice in the degenerateHubbard model,” Journal de Physique, vol. 36, no. 3, pp. 253–266,1975.[63] K. Foyevtsova and G. A. Sawatzky, “A Band Theory Perspective onMolecular Orbitals in Complex Oxides,” Journal of Modern Physics,vol. 10, no. 08, pp. 953–965, 2019.[64] J. Kanamori, “Electron Correlation and Ferromagnetism of TransitionMetals,” Progress of Theoretical Physics, vol. 30, pp. 275–289, 09 1963.[65] A. M. Oles´, “Antiferromagnetism and correlation of electrons in tran-sition metals,” Physical Review B, vol. 28, pp. 327–339, July 1983.[66] A. B. Georgescu and S. Ismail-Beigi, “Generalized slave-particle methodfor extended Hubbard models,” Physical Review B, vol. 92, p. 235117,December 2015.[67] J. B. Goodenough, Magnetism and the Chemical Bond. John Wiley &Sons Inc., 1963.[68] H. Hellmann, “Einfu¨hrung in die Quantenchemie,” in Hans Hellmann:Einfu¨hrung in die Quantenchemie, pp. 19–376, Berlin, Heidelberg:Springer Berlin Heidelberg, 2015.92[69] R. P. Feynman, “Forces in Molecules,” Physical Review, vol. 56,pp. 340–343, August 1939.[70] D. J. Griffiths, Introduction to Quantum Mechanics. Toronto: PearsonPrentice Hall, 2nd ed., 2005.[71] F. Jensen, Introduction to Computational Chemistry. Mississauga:John Wiley & Sons Inc., 2nd ed., 2007.[72] J.-P. Blaizot and G. Ripka, Quantum Theory of Finite Systems. Lon-don: MIT Press, 1st ed., 1986.[73] M. E. Peskin and D. V. Schroeder, An Introduction to Quantum FieldTheory. Boca Raton, FL: CRC Press, 1995.[74] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods ofQuantum Field Theory in Statistical Physics. Dover, 1975.[75] M. P. Marder, Condensed Matter Physics. Hoboken, NJ, USA: JohnWiley & Sons, Inc., October 2010.[76] A. S. Banerjee, P. Suryanarayana, and J. E. Pask, “Periodic Pu-lay method for robust and efficient convergence acceleration of self-consistent field iterations,” Chemical Physics Letters, vol. 647, pp. 31–35, March 2016.[77] P. Pulay, “Convergence acceleration of iterative sequences. the case ofscf iteration,” Chemical Physics Letters, vol. 73, pp. 393–398, July 1980.[78] R. Peters, Magnetic Phases in the Hubbard Model. PhD thesis, derGeorg-August-Universitat Gottingen, 2009.[79] R. Peters and T. Pruschke, “Half-filled Hubbard model on a Bethelattice with next-nearest-neighbor hopping,” Physical Review B, vol. 79,no. 4, pp. 1–7, 2009.[80] E. C. Stoner, “LXXX. Atomic moments in ferromagnetic metals andalloys with non-ferromagnetic elements,” The London, Edinburgh,and Dublin Philosophical Magazine and Journal of Science, vol. 15,pp. 1018–1034, May 1933.93Appendix ACalculating the hoppingcoefficients tab(k)In this appendix we start with the basic expressions for the hopping withinthe tight-binding model that we are constructing for the nickelates, and wecalculate the hopping matrix elements to put them eventually into a simpleform in Fourier space.A.1 1st and 4th neighbour hoppingFor 1st order hopping, there will not be any mixed orbital hopping betweendifferent directions, as discussed in the text, so we haveTˆ1 = −t1∑iσ[d†izσ(di−z,zσ + di+z,zσ) + d†ixσ(di−x,xσ + di+x,xσ)++ d†iyσ(di−y,yσ + di+y,yσ)]. (A.1)The first step it to put everything in terms of the |z〉 and |z¯〉 states, usingthe expressionsd†ixσ = −12d†izσ +√32d†iz¯σd†iyσ = −12d†izσ −√32d†iz¯σ. (A.2)Let us massage the terms separately. The y hopping givesd†iyσ(di−y,yσ+di+y,yσ) =(12d†izσ +√32d†iz¯σ)××(12di−y,zσ +√32di−y,z¯σ +12di+y,zσ +√32di+y,z¯σ)94A.1. 1st and 4th neighbour hoppingwhile the x terms ared†ixσ(di−x,xσ+di+x,xσ) =(12d†izσ −√32d†iz¯σ)××(12di−x,zσ −√32di−x,z¯σ +12di+x,zσ −√32di+x,z¯σ)Now all the operators are expressed in terms of the same orbitals, but hop-ping still occurs along different directions. Fourier transforming into mo-mentum space should fix this. Notice that what orbital the operator haswill have no influence on this transformation: the only thing that matters isthe site index. Letting τ = x, y, z denote the displacement index for a site(and simultaneously a physical unit displacement vector axˆ, ayˆ, azˆ), we seethat a standard term will transform in momentum space as∑iσd†iaσdi+τ,bσ =1N∑ikqσe−ik·Rid†kaσeiq·Ri+τdqbσ ==1N∑ikqσe−i(k−q)·Rid†kaσeiq·τdqbσ =∑kσeik·τd†kaσdkbσSince the Hamiltonian is Hermitian, the hermitian conjugate of this termwill also be present, which is of course hopping in the opposite direction— in our Hamiltonian it is represented by the hopping from the other sideonto the same site i (plus orbital swap), as can be seen by shifting the indexi→ i+ τ∑iσd†ibσdi−τ,aσ =∑iσd†i+τ,bσdiaσsame as above−−−−−−−−→∑kσe−ik·τd†kaσdkbσ.Together, of course, these terms produce the usual band structure∑iσ(d†iaσdi+τ,bσ + h.c.) = 2∑kσcos(k · τ)d†kaσdkbσ.To utilize these insights, let us group the appropriate terms from the x, y, zdirections into nice groups based on orbitals. First, these are all the zz95A.1. 1st and 4th neighbour hoppingterms from Tˆ1:zz terms hopping along z =∑iσd†izσ(di−z,zσ + di+z,zσ) =∑iσd†izσdi−z,zσ++∑iσd†izσdi+z,zσ︸ ︷︷ ︸re-sum to h.c. form=∑iσd†izσdi−z,zσ +∑iσd†i−z,zσdizσ ==∑iσ(d†izσdi−z,zσ + h.c. )applying the prescription−−−−−−−−−−−−−−−→2∑kσcos(k · τ︸︷︷︸=azˆ)d†kzσdkzσ = 2∑kσcos(kza)d†kzσdkzσ.Now we collect zz terms hopping along x. We havezz terms hopping along x =14∑iσ(d†izσdi−x,zσ + d†izσdi+x,zσ)apply prescription−−−−−−−−−−−→re-sum as before12∑kσcos(kxa)d†kzσdkzσ.Now for zz terms hopping along y. We have exactly the same outcome asabovezz terms hopping along y =14∑iσ(d†izσdi−y,zσ + d†izσdi+y,zσ) −→12∑kσcos(kya)d†kzσdkzσ.Now let us do the z¯z¯ terms. There are none with z hopping (for z these arethe in-plane orbitals), so we go to x:z¯z¯ terms hopping along x =34∑iσ(d†iz¯σdi+x,z¯σ + d†iz¯σdi+x,z¯σ) −→32∑kσcos(kxa)d†kz¯σdkz¯σ.Exactly the same for yz¯z¯ terms hopping along y =34∑iσ(d†iz¯σdi+y,z¯σ + d†iz¯σdi+y,z¯σ) −→32∑kσcos(kya)d†kz¯σdkz¯σ.96A.1. 1st and 4th neighbour hoppingNow for the cross terms. Again, there are no cross terms with z hopping, sowe go directly to x. There will be four such terms, in conjugate pairszz¯ terms hopping along x = −√34∑iσ((d†izσdi−x,z¯σ + d†iz¯σdi+x,zσ)︸ ︷︷ ︸h.c.’s++ (d†iz¯σdi−x,zσ + d†izσdi+x,z¯σ)︸ ︷︷ ︸h.c.’s)apply prescription−−−−−−−−−−−→same resummation−√32∑kσcos(kxa)(d†kzσdkz¯σ + d†kz¯σdkzσ).Since the terms in the expansion of y are exactly the same save for thedisplacement vector yˆ and that there are no minus signs, we get exactly thesame outcome, with a + signzz¯ terms hopping along y =√34∑iσ((d†izσdi−y,z¯σ + d†iz¯σdi+y,zσ)++ (d†iz¯σdi−y,zσ + d†izσdi+y,z¯σ))−→√32∑kσcos(kya)(d†kzσdkz¯σ + d†kz¯σdkzσ).And that’s it! We have put our Tˆ1 Hamiltonian into momentum space, allin terms of z and z¯ orbitals. We writeTˆ1 = −t1∑kσ[t1zz(k)d†kzσdkzσ+t1z¯z¯(k)d†kzσdkzσ++ t1zz¯(k)(d†kzσdkz¯σ + d†kz¯σdkzσ)], (A.3)witht1zz(k) = 2 cos(kza) +12[cos(kxa) + cos(kya)],t1z¯z¯(k) =32[cos(kxa) + cos(kya)], (A.4)t1zz¯(k) = −√32[cos(kxa)− cos(kya)].As mentioned in the text, we can immediately write down the corre-sponding result for the 4th neighbour hopping by merely doubling the hop-97A.2. 2nd neighbour hoppingping distance a→ 2a. Explicitly, we haveTˆ4 = −t4∑kσ[t4zz(k)d†kzσdkzσ+t4z¯z¯(k)d†kzσdkzσ++ t4zz¯(k)(d†kzσdkz¯σ + d†kz¯σdkzσ)], (A.5)witht4zz(k) = 2 cos(2kza) +12[cos(2kxa) + cos(2kya)],t4z¯z¯(k) =32[cos(2kxa) + cos(2kya)], (A.6)t4zz¯(k) = −√32[cos(2kxa)− cos(2kya)].A.2 2nd neighbour hoppingStart withTˆ2 = −t2∑iabσ,τˆ 6=µˆd†i+τˆ+µˆ,aσdibσ. (A.7)As discussed in the text, we can take b = τ and a = µ, and then take theFourier transform just as before, obtaining∑iabσ,τˆ 6=µˆd†i+τˆ+µˆ,aσdibσ =∑iσ,a6=bd†i+bˆ+aˆ,aσdibσ =∑kσ,a6=beik·(bˆ+aˆ)d†kaσdkbσ.There are multiple possible summands arising from the above, but they willcome in conjugate pairs. There are a total of 24 terms (6 ·6 = 36 overall, butcannot have a and b be equal. There are 6 identical terms, but we sum overboth a and b, so in fact we leave out 12 terms, so 36− 12 = 24), which areconjugate pairs based on the overall sign, like +x+ y and −x− y, or +x− yand −x + y. Grouping them into conjugate pairs produces the following98A.2. 2nd neighbour hopping(note the standard cos dispersion)∑kσ,a6=beik·(bˆ+aˆ)d†kaσdkbσ =∑kσ2 cos(kxa+ kya)[d†kxσdkyσ + d†kyσdkxσ]++ 2 cos(kya+ kza)[d†kyσdkzσ + d†kzσdkyσ]++ 2 cos(kza+ kxa)[d†kzσdkxσ + d†kxσdkzσ]++ 2 cos(kxa− kya)[d†kxσdkyσ + d†kyσdkxσ]++ 2 cos(kya− kza)[d†kyσdkzσ + d†kzσdkyσ]++ 2 cos(kza− kxa)[d†kzσdkxσ + d†kxσdkzσ].Using the trigonometric identity cos(x+ y) + cos(x− y) = 2 cos(x) cos(y),we can reduce this even further to obtain∑kσ,a6=beik·(bˆ+aˆ)d†kaσdkbσ =∑kσ(4 cos(kxa) cos(kya)××[d†kxσdkyσ + d†kyσdkxσ]+ 4 cos(kya) cos(kza)[d†kyσdkzσ + d†kzσdkyσ]++ 4 cos(kza) cos(kxa)[d†kzσdkxσ + d†kxσdkzσ] ).To make further progress, we need to express everything in terms of ourpreferred orbital basis |z〉 and |z¯〉. Notice that, simply inverting the matrixin favour of |x〉 , |y〉 from before, we find[|x〉|y〉]=[−12√32−12 −√32]︸ ︷︷ ︸=U[|z〉|z¯〉]→[|z〉|z¯〉]=[−1 −11√3− 1√3] [|x〉|y〉]. (A.8)We can represent the three terms in our Hamiltonian in the following way(using the shorthand dx ≡ dkxσ and so on)[d†kxσdkyσ + d†kyσdkxσ]→[d†x, d†y] [dydx](A.9)Symbolically changing basis, we have[d†x, d†y] [dydx]=[d†z, d†z¯](U−1)tU s[dz¯dz]==[d†z, d†z¯] [−12 −12√32 −√32][−√32 −12√32 −12] [dz¯dz]=[d†z, d†z¯] [ 0 1−32 0] [dz¯dz],99A.2. 2nd neighbour hoppingwhere the notation U s means we swapped the two rows. Instead of pluggingin expressions for each term, we can just perform this matrix product. Theother terms looks similar, except the matrix will have one simple row forthe terms that already include z:[d†y, d†z] [dzdy]=[d†z¯, d†z] [−√32 0−12 1][1 0−12 −√32] [dzdz¯]==[d†z¯, d†z] [−√32 0−1 −√32] [dzdz¯]=[d†z, d†z¯] [−√32 −10 −√32] [dz¯dz]just to straighten it out in accordance with the rest of the terms. Finally,for the last term we find (Adopt the shorthand cx ≡ cosx)[d†z, d†x] [dxdz]=[d†z, d†z¯] [1 −120√32][√32 −120 1] [dz¯dz]==[d†z, d†z¯] [√32 −10√32] [dz¯dz].Adding them all together, we findcos(kxa) cos(kya)[d†kxσdkyσ + d†kyσdkxσ]+ cos(kya) cos(kza)××[d†kyσdkzσ + d†kzσdkyσ]+ cos(kza) cos(kxa)[d†kzσdkxσ + d†kxσdkzσ]== cos(kxa) cos(kya)[d†z, d†z¯] [ 0 1−32 0] [dz¯dz]++ cos(kya) cos(kza)[d†z, d†z¯] [−√32 −10 −√32] [dz¯dz]++ cos(kza) cos(kxa)[d†z, d†z¯] [√32 −10√32] [dz¯dz]==[d†z, d†z¯] [√32 cz(cx − cy) cxcy − cz(cy + cx)−32cxcy√32 cz(cx − cy)] [dz¯dz]== [cxcy − cz(cy + cx)] d†zdz +[−32cxcy]d†z¯dz¯++[√32cz(cx − cy)](d†zdz¯ + d†z¯dz).100A.3. Full hopping operatorDone! Putting this all together, the form of the Hamiltonian isTˆ2 = −2t2∑kσ[t2zz(k)d†kzσdkzσ + t2z¯z¯(k)d†kz¯σdkz¯σ++ t2zz¯(k)(d†kzσdkz¯σ + d†kz¯σdkzσ)], (A.10)witht2zz(k) = 2 cos(kxa) cos(kya)− 2 cos(kza)(cos(kya) + cos(kxa)),t2zz¯(k) =√3 cos(kza)(cos(kxa)− cos(kya)), (A.11)t2z¯z¯(k) = −3 cos(kxa) cos(kya).A.3 Full hopping operatorAltogether, then, the hopping operator has the formTˆ =∑kσ[tzz(k)d†kzσdkzσ + tz¯z¯(k)d†kz¯σdkz¯σ++ tzz¯(k)(d†kzσdkz¯σ + d†kz¯σdkzσ)], (A.12)with the definitionstzz(k) =− 2t1(cos(kza) +14[cos(kxa) + cos(kya)])−− 2t4(cos(2kza) +14[cos(2kxa) + cos(2kya)])−− 2t2 [cos(kxa) cos(kya)− 2 cos(kza)(cos(kya) + cos(kxa))] ,tz¯z¯(k) =− 3t12[cos(kxa) + cos(kya)]− 3t42[cos(2kxa) + cos(2kya)]+ (A.13)+ 6t2 cos(kxa) cos(kya),tzz¯(k) =√3t12[cos(kxa)− cos(kya)] +√3t42[cos(2kxa)− cos(2kya)]−− 2√3t2 cos(kza)(cos(kxa)− cos(kya)).101Appendix BMinimizing latticecontributions: solving thecubic equationIn the text we were faced with the need to solve the cubic equationu3 +1au− δa= 0 (B.1)which arose during lattice energy minimization in the Hartree-Fock process.It turns out that it admits a closed form solution: writeu3 = (s− t)3, 1a= 3st,δa= s3 − t3Combining the latter two equations into one fo t, we findt6 + t3(δa)−(13a)3= 0which is easily solved as a quadratict3 =12−δa±√(δa)2+427a3s3 = t3 +δa=12δa±√(δa)2+427a3from which we can recover the expression for uu = δ32β1/3[(1 +√1 +1β)1/3+(1−√1 +1β)1/3](B.2)where we set β = 274 aδ2. So we see that the charge disproportionation factorδ couples to the lattice distortions.102Appendix CCalculation of 〈Hˆe〉In this appendix we express the electron-electron interaction contribution tothe Hartree-Fock energy, 〈Hˆe〉, which we found in the text to be〈Hˆe〉Ψe = U∑ia[〈d†ia↑dia↑〉〈d†ia↓dia↓〉 − 〈d†ia↑dia↓〉〈d†ia↓dia↑〉]++ U ′∑iσ[〈d†izσdizσ〉〈d†iz¯σ¯diz¯σ¯〉 − 〈d†izσdiz¯σ¯〉〈d†iz¯σ¯dizσ〉]++ (U ′ − J)∑iσ[〈d†izσdizσ〉〈d†iz¯σdiz¯σ〉 − 〈d†izσdiz¯σ〉〈d†iz¯σdizσ〉]−− J∑iσ[〈d†izσdizσ¯〉〈d†iz¯σ¯diz¯σ〉 − 〈d†izσdiz¯σ〉〈d†iz¯σ¯dizσ¯〉]++ γJ∑ia[〈d†ia↑dia¯↑〉〈d†ia↓dia¯↓〉 − 〈d†ia↑dia¯↓〉〈d†ia↓dia¯↑〉], (C.1)in terms of the newly-defined order parameters δ, SFM, ..., which are set by〈d†iaσd†iaσ〉 =14[1 + δeiQc·Ri]+σ2[SFM + SAFMeiQc·Ri++ S1z cos(Qm ·Ri) + S2z sin(Qm ·Ri)],〈d†iaσd†iaσ¯〉 =12[S1x cos(Qm ·Ri) + S2x sin(Qm ·Ri)] ,〈d†iaσd†ia¯σ〉 =O1 +O2eiQc·Ri + σ[Z1 + Z2eiQc·Ri++ Z3 cos(Qm ·Ri) + Z4 sin(Qm ·Ri)],〈d†iaσd†ia¯σ¯〉 =X1 cos(Qm ·Ri) +X2 sin(Qm ·Ri). (C.2)We calculate the contributions term by term. The important trick toremember is that∑i eis·Ri(...) = 0 unless s ·Ri = 0 mod 2pi. Otherwise theexponential keeps alternating sign and the sum essentially cancels itself out.An exact derivation of this fact, as well as many important examples of itsapplication are given in the appendices of Marder’s textbook on condensed103Appendix C. Calculation of 〈Hˆe〉matter physics [75]. Below, whenever we encounter sums over terms thathave non-vanishing exponentials we will immediately drop them from theexpression. And whenever there is a cosine or a sine, expressing it in termsof the exponentials will allow us to apply that same rule. Let us derive thoserules: for a product of cosines∑icos2(Qm ·Ri) =∑i14(ei(Qm·Ri) + e−i(Qm·Ri))2 ==∑i14(eiQc·Ri + 2 + e−iQc·Ri) =142N =N2. (C.3)For sines this is identical, except there is a minus sign which is fixed by thepresence of the i in the denominator that gets squared to −1∑isin2(Qm ·Ri) =∑i1−4(ei(Qm·Ri) − e−i(Qm·Ri))2 = N2. (C.4)Finally, for a product of sines and cosines we find∑icos(Qm ·Ri) sin(Qm ·Ri) =∑i14(ei(Qm·Ri) + e−i(Qm·Ri))××(ei(Qm·Ri) − e−i(Qm·Ri)) =∑i14(eiQc·Ri + 0− e−iQc·Ri) = 0. (C.5)The first term givesU∑ia[〈d†ia↑dia↑〉〈d†ia↓dia↓〉 − 〈d†ia↑dia↓〉〈d†ia↓dia↑〉]== U∑ia[ 116(1 + δeiQc·Ri)2 − σ24[SFM + SAFMeiQc·Ri + S1z cos(Qm ·Ri)++ S2z sin(Qm ·Ri)]2 − 14(S1x cos(Qm ·Ri) + S2x sin(Qm ·Ri))2]== U∑i[18(1 + δ2)− 14[S2FM + S2AFM +12S21z +12S22z]−− 14(12S21x +12S22x)]== 2UN[116(1 + δ2)− 14(S2FM + S2AFM +S21z + S22z + S21x + S22x2)].(C.6)104Appendix C. Calculation of 〈Hˆe〉The second term gives, by analogy (it is almost exactly the same as thefirst, same for X instead of S1,2x in the last term)U ′∑iσ[〈d†izσdizσ〉〈d†iz¯σ¯diz¯σ¯〉 − 〈d†izσdiz¯σ¯〉〈d†iz¯σ¯dizσ〉]== U ′∑iσ[ 116(1 + δeiQc·Ri)2 − σ24[SFM + SAFMeiQc·Ri + S1z cos(Qm ·Ri)++ S2z sin(Qm ·Ri)]2 − 14(X1 cos(Qm ·Ri) +X2 sin(Qm ·Ri))2]== 2U ′N[116(1 + δ2)− 14(S2FM + S2AFM +S21z + S22z2)− 12(X21 +X22 )].(C.7)The third term gives(U ′ − J)∑iσ[〈d†izσdizσ〉〈d†iz¯σdiz¯σ〉 − 〈d†izσdiz¯σ〉〈d†iz¯σdizσ〉]== (U ′−J)∑iσ{14(1+δeiQc·Ri)+σ2[SFM+SAFMeiQc·Ri+S1z cos(Qm ·Ri)++ S2z sin(Qm ·Ri)]}2+{O1 +O2eiQc·Ri + σ[Z1 + Z2eiQc·Ri++ Z3 cos(Qm ·Ri) + Z4 sin(Qm ·Ri)]}2= (U ′ − J)∑iσ{ 116(1 + δ2)++14[S2FM + S2AFM +S21z + S22z2]+σ4[SFM + SAFMδ]}+{O21 +O22 + Z21++ Z22 +Z23 + Z242+ σ(O1Z1 +O2Z2)}= 2(U ′ − J)N[ 116(1 + δ2)++14[S2FM + S2AFM +S21z + S22z2]+O21 +O22 + Z21 + Z22 +Z23 + Z242]. (C.8)where towards the end the σ-dependent terms dropped out because ofthe sum over the spins.105Appendix C. Calculation of 〈Hˆe〉The fourth term gives− J∑iσ[〈d†izσdizσ¯〉〈d†iz¯σ¯diz¯σ〉 − 〈d†izσdiz¯σ〉〈d†iz¯σ¯dizσ¯〉]== −J∑iσ[14(S1x cos(Qm ·Ri) + S2x sin(Qm ·Ri))2 − (O1 +O2eiQc·Ri)2++(Z1 + Z2eiQc·Ri + Z3 cos(iQc ·Ri) + Z4 sin(iQc ·Ri))2]== −2JN[S21x + S22x8− (O21 +O22) +(Z21 + Z22 +Z23 + Z242)]. (C.9)And, finally, the fifth termγJ∑ia[〈d†ia↑dia¯↑〉〈d†ia↓dia¯↓〉 − 〈d†ia↑dia¯↓〉〈d†ia↓dia¯↑〉]== γJ∑ia[(O1 +O2eiQc·Ri)2 −(Z1 + Z2eiQc·Ri + Z3 cos(iQc ·Ri)++ Z4 sin(iQc ·Ri))2 − (X1 cos(Qm ·Ri) +X2 sin(Qm ·Ri))2] == 2γJN[O21 +O22 −(Z21 + Z22 +Z23 + Z24 +X21 +X222)]. (C.10)Combining all of these together, we find〈Hˆe〉N=3U − 5J8(1 + δ2)− U + J2(S2FM + S2AFM +S21 + S222)−− 2(U ′ − (2 + γ)J)(O21 +O22)− 2(U ′ + γJ)(Z21 + Z22++Z23 + Z242)− (U ′ + γJ)(X21 +X22 ). (C.11)In this thesis we make the assumption of spherical symmetry for thetwo-band Hubbard terms: this assumption means that U ′ = U − 2J . Withthis in mind, we obtain the final form〈Hˆe〉N=3U − 5J8(1 + δ2)− U + J2(S2FM + S2AFM +S21 + S222)−− 2(U − (4 + γ)J)(O21 +O22)− 2(U + (γ − 2)J)(Z21 + Z22++Z23 + Z242)− (U + (γ − 2)J)(X21 +X22 ). (C.12)106Appendix DFourier Transforming HˆHFThe task of this appendix is to make use of the Bloch symmetry of theHartree-Fock HamiltonianHˆHF =∑ijabσtabij d†iaσdjbσ +{3U − 5J4+[3U − 5J4δ − 2bu]eiQc·Ri−−σ2(U+J)[SFM+SAFMeiQc·Ri+S1z cos Qm ·Ri+S2z sin Qm ·Ri]}d†iaσdiaσ−− U + J2[S1x cos Qm ·Ri + S2x sin Qm ·Ri] d†iaσdiaσ¯−− (U + J(γ − 2))[X1 cos Qm ·Ri +X2 sin Qm ·Ri]d†iaσdia¯σ¯++{(J(4 + γ)− U)[O1 +O2eiQc·Ri]− σ[U + J(γ − 2)]××[Z1 + Z2eiQc·Ri + Z3 cos Qm ·Ri + Z4 sin Qm ·Ri]}d†iaσdia¯σ. (D.1)and Fourier transform it into k-space according to the ruled†kaσ =∑ieikRi√Nd†iaσ. (D.2)The Fourier transform of the hopping operator was already carried out inAppendix A, resulting inTˆ =∑kabσtab(k)d†kaσdkbσ,with the coefficients tab(k) tabulated in Appendix A. In order to deal withthe rest, we need a quick rule of thumb. Fourier transforming a term of theform∑i eis·Rid†iaσdibτ simply offsets the momentum label of one of the doperators by the vector s∑ieis·Rid†iaσdibτFourier transform−−−−−−−−−−−→with D.2∑kd†k+s,aσdkbτ . (D.3)107Appendix D. Fourier Transforming HˆHFArmed with this rule, we tackle the Hamiltonian term by term. The firstterm gives{3U − 5J4+[3U − 5J4δ − 2bu]eiQc·Ri− σ2(U +J)[SFM +SAFMeiQc·Ri++ S1z cos Qm ·Ri + S2z sin Qm ·Ri]}d†iaσdiaσ ==[3U − 5J4− σ2(U + J)SFM]d†kaσdkaσ +[3U − 5J4δ − 2bu−− σ2(U + J)SAFM]d†k+Qc,aσdkaσ −σ4(U + J)[(S1z − iS2z)d†k+Qm,aσdkaσ++ (S1z + iS2z)d†k−Qm,aσdkaσ].The second term gives− U + J2[S1x cos Qm ·Ri + S2x sin Qm ·Ri] d†iaσdiaσ¯ == −14(U + J)[(S1x − iS2x)d†k+Qm,aσdkaσ¯ + (S1x + iS2x)d†k−Qm,aσdkaσ¯].The third term gives(U + J(γ − 2))[X1 cos Qm ·Ri +X2 sin Qm ·Ri]d†iaσdia¯σ¯ ==12(U +J(γ−2))[(X1− iX2)d†k+Qm,aσdka¯σ¯ + (X1 + iX2)d†k−Qm,aσdka¯σ¯].Finally, the last term gives{(J(4 + γ)− U)[O1 +O2eiQc·Ri]− σ[U + J(γ − 2)][Z1 + Z2eiQc·Ri++ Z3 cos Qm ·Ri + Z4 sin Qm ·Ri]}d†iaσdia¯σ ==[(J(4 + γ)−U)O1−σ(U + J(γ− 2))Z1]d†k,aσdka¯σ +[(J(4 + γ)−U)O2−−σ(U+J(γ−2))Z2]d†k+Qc,aσdka¯σ−σ(U+J(γ−2))[(Z3−iZ4)d†k+Qm,aσdka¯σ++ (Z3 + iZ4)d†k−Qm,aσdka¯σ].Define some abbreviated interaction strengthsU¯ =3U − 5J4, U0 =U + J2, U1 =U + J(γ − 2)2, U2 = J(4 + γ)− U.(D.4)108Appendix D. Fourier Transforming HˆHFAlso defineΣz = S1z + iS2z, Σx = S1x + iS2x, χ = X1 + iX2, ζ = Z3 + iZ4. (D.5)Using these, we can write down the expressions for the various 4×4 subblocksof the 16× 16 matrices h(k) that we diagonalize at every k. Recall that thematrices h(k) are writtenHˆHF =∑kψ†kh(k)ψk,whereψ†k = (ψ†kz↑, ψ†kz¯↑, ψ†kz↓, ψ†kz¯↓)andψ†kaσ = (d†kaσ, d†k+Qm,aσ, d†k+Qc,aσ, d†k−Qm,aσ).The matrices for the same orbital, same spin subblock have the formhˆaσ,aσ =U¯−σU0SFM++taa(k)−σU02 Σz U¯δ−2bu−−σU0SAFM −σU02 Σ∗z−σU02 Σ∗z U¯−σU0SFM++taa(k+Qm) −σU02 Σz U¯δ−2bu−−σU0SAFMU¯δ−2bu−−σU0SAFM −σU02 Σ∗zU¯−σU0SFM++taa(k+Qc)−σU02 Σz−σU02 Σz U¯δ−2bu−−σU0SAFM −σU02 Σ∗zU¯−σU0SFM+taa(k−Qm). (D.6)The matrices for the same orbital, opposite spin subblocks arehˆaσ,aσ¯ =−U02 Σx −U02 Σ∗x−U02 Σ∗x −U02 Σx−U02 Σ∗x −U02 Σx−U02 Σx −U02 Σ∗x . (D.7)The matrices for the opposite orbital, opposite spin arehˆaσ,a¯σ¯ =−U1χ −U1χ∗−U1χ∗ −U1χ−U1χ∗ −U1χ−U1χ −U1χ∗ . (D.8)109Appendix D. Fourier Transforming HˆHFFinally, the opposite orbital, same spin matrices have the formhˆaσ,a¯σ =U2O1−2σU1Z1++taa¯(k)−σU1ζ U2O2−−2σU1Z2 −σU1ζ∗−σU1ζ∗ U2O1−2σU1Z1++taa¯(k+Qm) −σU1ζ U2O2−−2σU1Z2U2O2−−2σU1Z2 −σU1ζ∗ U2O1−2σU1Z1++taa¯(k+Qc) −σU1ζ−σU1ζ U2O2−−2σU1Z2 −σU1ζ∗ U2O1−2σU1Z1++taa¯(k−Qm).(D.9)110Appendix ERe-calculating theHartree-Fock orderparametersIn this appendix we demonstrate how to close the self-consistency loop whensolving the Hartree-Fock equations, by recalculating the order parametersδ, SFM, ... for which the initial values were simply guessed. The order pa-rameters are defined by the expressions〈d†iaσd†iaσ〉 =14[1 + δeiQc·Ri]+σ2[SFM + SAFMeiQc·Ri++ S1z cos(Qm ·Ri) + S2z sin(Qm ·Ri)],〈d†iaσd†iaσ¯〉 =12[S1x cos(Qm ·Ri) + S2x sin(Qm ·Ri)] ,〈d†iaσd†ia¯σ〉 =O1 +O2eiQc·Ri + σ[Z1 + Z2eiQc·Ri++ Z3 cos(Qm ·Ri) + Z4 sin(Qm ·Ri)],〈d†iaσd†ia¯σ¯〉 =X1 cos(Qm ·Ri) +X2 sin(Qm ·Ri). (E.1)When we diagonalize the Hamiltonian HˆHF from Appendix D, we are leftwith the 16-tuple momentum space eigenvectors arranged likeψ†k = (ψ†kz↑, ψ†kz¯↑, ψ†kz↓, ψ†kz¯↓)withψ†kaσ = (d†kaσ, d†k+Qm,aσ, d†k+Qc,aσ, d†k−Qm,aσ).As is discussed in the main text in Chapter 3, the density matrix,ρ(ζbτ, ηaσ) =∑k〈Ψe| d†k+ηQm,aσdk+ζQm,bτ |Ψe〉 , (E.2)111Appendix E. Re-calculating the Hartree-Fock order parametersis essentially the quantity we need to calculate, and it turns out that it canbe obtained from the occupied eigenvectors |k, n〉 (n runs from 1 to 16, asthere are 16 eigenvectors for each 16× 16 subblock) byρ(ζbτ, ηaσ) =′∑k|k, n〉 〈k, n| , (E.3)where the prime on the sum indicates that the summation is over only theoccupied eigenvectors (which are determined by ranking all the eigenvectorsin accordance with their corresponding eigenvalues, then taking the lowerquarter). Thus the crucial task is to connect the real-space matrix elementsin E.1 with the momentum space density matrix ρ(ζbτ, ηaσ). To do this,we use Fourier’s trick: this is the same trick we used when Fourier trans-forming the Hamiltonian in Appendix D. The trick is that any sum of theform∑i eis·Ri(...) is zero unless s ·Ri = 0 mod 2pi: essentially the sign ofthe summands keeps alternating and they effectively average out. We canuse this to invert the expressions in Eqs: E.1 for the various order parame-ters δ,O1, ... by multiplying by the inverse of the exponential factor of theparameter we want so as to cancel it out: the rest of the terms in that ex-pression will then sum to zero because of Fourier’s trick and we can isolatethe one we are interested in. This was already demonstrated in the text forSFM and for δ: we quote the expressions here for completenessδ =1N∑kaσ〈d†k+Qc,aσdkaσ〉 =∑aσ[ρ(0aσ, 2aσ) + ρ(1aσ,−1aσ)++ ρ(2aσ, 0aσ) + ρ(−1aσ, 1aσ)]. (E.4)SFM =12N∑kaσσ〈d†kaσdkaσ〉 =12N∑ηaσσρ(ηaσ, ηaσ). (E.5)We obtain the rest of the by analogy. The expressions for many areessentially immediateSAFM =12N∑kaσσ〈d†k+Qc,aσdkaσ〉 ==∑aσσ[ρ(0aσ, 2aσ) + ρ(1aσ,−1aσ) + ρ(2aσ, 0aσ) + ρ(−1aσ, 1aσ)],(E.6)112Appendix E. Re-calculating the Hartree-Fock order parametersO1 =14N∑kaσ〈d†kaσdkaσ〉 =14N∑ηaσρ(ηaσ, ηa¯σ), (E.7)Z1 =14N∑kaσσ〈d†kaσdkaσ〉 =14N∑ηaσσρ(ηaσ, ηa¯σ), (E.8)O2 =14N∑kaσσ〈d†k+Qc,aσdkaσ〉 ==∑aσσ[ρ(0aσ, 2aσ) + ρ(1aσ,−1aσ) + ρ(2aσ, 0aσ) + ρ(−1aσ, 1aσ)].(E.9)The other parameters require just a little more finesse. Looking to thesecond line in eqs. E.1, re-express the trigonometric functions in terms ofthe exponentials〈d†iaσd†iaσ¯〉 =14[(S1x − iS2x) eiQm·Ri + (S1x + iS2x) e−iQm·Ri]. (E.10)Multiplying by the two exponentials and carrying out the sum produces theexpressions∑iaσe−iQm·Ri〈d†iaσd†iaσ¯〉 = S1x − iS2x,∑iaσeiQm·Ri〈d†iaσd†iaσ¯〉 = S1x + iS2x.(E.11)Adding or subtracting these and carrying out the Fourier transform givesthe expressions for S1/2xS1x =12N∑kaσ[〈d†k+Qm,aσdkaσ〉+ 〈d†k−Qm,aσdkaσ〉]==12N∑aσ[(ρ(0aσ, 1aσ¯) + ρ(1aσ, 2aσ¯) + ρ(2aσ,−1aσ¯) + ρ(−1aσ, 0aσ¯))++(ρ(1aσ, 0aσ¯) + ρ(2aσ, 1aσ¯) + ρ(−1aσ, 2aσ¯) + ρ(0aσ,−1aσ¯))]. (E.12)S2x =12iN∑kaσ[〈d†k+Qm,aσdkaσ〉 − 〈d†k−Qm,aσdkaσ〉]==12iN∑aσ[(ρ(0aσ, 1aσ¯) + ρ(1aσ, 2aσ¯) + ρ(2aσ,−1aσ¯) + ρ(−1aσ, 0aσ¯))−−(ρ(1aσ, 0aσ¯) + ρ(2aσ, 1aσ¯) + ρ(−1aσ, 2aσ¯) + ρ(0aσ,−1aσ¯))]. (E.13)113Appendix E. Re-calculating the Hartree-Fock order parametersBy analogy with these we can immediately write down the rest of the ex-pressions, as the only ones remaining also only depend on the trigonometricfunctionsS1z =12N∑kaσσ[〈d†k+Qm,aσdkaσ〉+ 〈d†k−Qm,aσdkaσ〉]= (E.14)=12N∑aσσ[(ρ(0aσ, 1aσ) + ρ(1aσ, 2aσ) + ρ(2aσ,−1aσ) + ρ(−1aσ, 0aσ))+(E.15)+(ρ(1aσ, 0aσ) + ρ(2aσ, 1aσ) + ρ(−1aσ, 2aσ) + ρ(0aσ,−1aσ))]. (E.16)S2z =12iN∑kaσσ[〈d†k+Qm,aσdkaσ〉 − 〈d†k−Qm,aσdkaσ〉]==12iN∑aσσ[(ρ(0aσ, 1aσ) + ρ(1aσ, 2aσ) + ρ(2aσ,−1aσ) + ρ(−1aσ, 0aσ))−(ρ(1aσ, 0aσ) + ρ(2aσ, 1aσ) + ρ(−1aσ, 2aσ) + ρ(0aσ,−1aσ))]. (E.17)Z3 =14N∑kaσσ[〈d†k+Qm,aσdka¯σ〉+ 〈d†k−Qm,aσdka¯σ〉]==14N∑aσσ[(ρ(0aσ, 1a¯σ) + ρ(1aσ, 2a¯σ) + ρ(2aσ,−1a¯σ) + ρ(−1aσ, 0a¯σ))++(ρ(1aσ, 0a¯σ) + ρ(2aσ, 1a¯σ) + ρ(−1aσ, 2a¯σ) + ρ(0aσ,−1a¯σ))]. (E.18)Z4 =14iN∑kaσσ[〈d†k+Qm,aσdka¯σ〉 − 〈d†k−Qm,aσdka¯σ〉]==14iN∑aσσ[(ρ(0aσ, 1a¯σ) + ρ(1aσ, 2a¯σ) + ρ(2aσ,−1a¯σ) + ρ(−1aσ, 0a¯σ))−(ρ(1aσ, 0a¯σ) + ρ(2aσ, 1a¯σ) + ρ(−1aσ, 2a¯σ) + ρ(0aσ,−1a¯σ))]. (E.19)Finally,X1 =14N∑kaσσ[〈d†k+Qm,aσdka¯σ¯〉+ 〈d†k−Qm,aσdka¯σ¯〉]==14N∑aσσ[(ρ(0aσ, 1a¯σ¯) + ρ(1aσ, 2a¯σ¯) + ρ(2aσ,−1a¯σ¯) + ρ(−1aσ, 0a¯σ¯))++(ρ(1aσ, 0a¯σ¯) + ρ(2aσ, 1a¯σ¯) + ρ(−1aσ, 2a¯σ¯) + ρ(0aσ,−1a¯σ¯))]. (E.20)114Appendix E. Re-calculating the Hartree-Fock order parametersX2 =14iN∑kaσσ[〈d†k+Qm,aσdka¯σ¯〉 − 〈d†k−Qm,aσdka¯σ¯〉]==14iN∑aσσ[(ρ(0aσ, 1a¯σ¯) + ρ(1aσ, 2a¯σ¯) + ρ(2aσ,−1a¯σ¯) + ρ(−1aσ, 0a¯σ¯))−(ρ(1aσ, 0a¯σ¯) + ρ(2aσ, 1a¯σ¯) + ρ(−1aσ, 2a¯σ¯) + ρ(0aσ,−1a¯σ¯))]. (E.21)115
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Effects of the lattice distortion on the magnetic order...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Effects of the lattice distortion on the magnetic order in rare-earth nickelates Fomichev, Stepan Olegovich 2019
pdf
Page Metadata
Item Metadata
Title | Effects of the lattice distortion on the magnetic order in rare-earth nickelates |
Creator |
Fomichev, Stepan Olegovich |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | The low-temperature magnetic order in the rare-earth nickelates is a subject of vigorous debate in the literature. Recent work emphasized the primary role of the electron-phonon coupling for the metal-insulator transition in the nickelates, and suggested that lattice distortions are the driver of the transition, leading to the observed charge order. However, to our knowledge there has been little work on the impact of lattice distortions on the magnetic order, in particular whether distortions favour some orders over others. In this thesis, we study the magnetic order in the nickelates at zero temperature, and investigate whether the breathing-mode lattice distortions select a preferred ground state. An effective two-band Hubbard model for the nickelates is constructed and coupled to the lattice distortions with an on-site Holstein-like term. The distortions are treated semiclassically. Using the Hartree-Fock approximation, we obtain the magnetic phase diagram, then turn on the coupling to the lattice to observe its impact on the various phases. Our model reproduces the earlier work showing the stronger charge disproportionation and insulating behaviour in the phase space due to increased coupling to the lattice. Furthermore, we find numerous 4-site magnetic orders that are self-consistent within the model, including all of the main suggestions in the literature (states such as ↑↑↓↓,↑→↓← and ⇑0⇓0). However, in this model a magnetic order can only couple to the lattice distortions if there is nonzero charge disproportionation. As a result, we find that coupling to the lattice distortions broadly favours the ⇑0⇓0 order in large sectors of the parameter space. Finally, we considered the impact of longer range hopping on the magnetic order: we find that the shape of the density of states, rather than overall bandwidth, primarily determines the magnetic ground state. A van Hove singularity arises even for small 2nd-nearest hopping amplitudes, which results in robust ferromagnetism across most of the phase diagram in a Stoner-like fashion. On the contrary, even small 4th-nearest amplitudes decrease the Fermi level density of states, resulting in ballooning of the metallic phase despite a barely renormalized bandwidth. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-08-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0380543 |
URI | http://hdl.handle.net/2429/71390 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2019-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2019_september_fomichev_stepan.pdf [ 3.32MB ]
- Metadata
- JSON: 24-1.0380543.json
- JSON-LD: 24-1.0380543-ld.json
- RDF/XML (Pretty): 24-1.0380543-rdf.xml
- RDF/JSON: 24-1.0380543-rdf.json
- Turtle: 24-1.0380543-turtle.txt
- N-Triples: 24-1.0380543-rdf-ntriples.txt
- Original Record: 24-1.0380543-source.json
- Full Text
- 24-1.0380543-fulltext.txt
- Citation
- 24-1.0380543.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0380543/manifest