2., 2x231; .,
.,.......s......€..
~ ‘ ; , . . , .T!
ﬁﬁéaaﬁﬁ s.
a , a
i.-
n;

Ismsmg’m

. ma
3... ;
.15.“...

: .

t .1 .f.
:9 2....

23:. .-
1.5.5.1.
p a

a5 .2
3 3.3.1991.
. (5.....1. .
2... 12:.
.:J‘
..

(I;

be" a; _. 2:55...

._ .. . :.
..Z.§.,.i...t

1.6.: . .1.
Litrf... .t?

r 4
.7

£5. 9%,.

. ‘Z
.

 

, aﬁwﬂgsﬁ

.7. . a ((7-17,

,. ..>..,

L.
$4.}!

 

 

LlBRARY

Michigan State
Universuty

 

 

 

This is to certify that the

dissertation entitled

Study of Electronic Structures of Solids
With Strongly Interacting Electrons

presented by

YEN—SEEING SU

has been accepted towards fulﬁllment
of the requirements for

Ph . D - degree in PHYSICS

 

 

 

 

\l/iﬂmttg A; iCOMaM

Major professor V
T .A . Kaplan

Date 06/06/2000

MSU is an Afﬁrmative Action/Equal Opportunity Institution 0-12771

   

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJClRC/DaieDuepGS-m 5

STUDY OF ELECTRONIC
STRUCTURES OF SOLIDS WITH
STRONGLY INTERACTING
ELECTRON S

By

Yen-Sheng Su

A DISSERTATION

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

2000

ABSTRACT

Study Of Electronic Structures Of Solids With Strongly
Interacting Electrons

By
Yen-Sheng Su

This work contains studies of two classes of perovskite transition metal oxides.
The ﬁrst class is the layered perovskite cuprates and the related nickelate. The second
class is the three dimensional perovskite manganites. Both model and ab initio
calculations are carried out for the two classes of systems. The dissertation is therefore
divided into the following four parts.

The ﬁrst part is about the 3-band Hubbard model. The model is commonly used
for describing the electronic properties of the important Cqu layers in the crystals of
high-Tc superconducting cuprates, such as doped LazCuO4 and YBazCu307. The
straightforward perturbation expansion on the model taking tpd/epd (~0.36 for the
cuprates) as the small parameter does not converge. In this work, I Show that there exist
canonical transformations on the model Hamiltonian such that the perturbation expansion
based on the transformed Hamiltonians converges.

In the second part, crystal Hartree—Fock calculations are carried out for LazNiO4
and LaZCuO4. The results predict correctly that these two materials are antiferromagnetic
insulators, in contrast to the wrong predictions made by the density functional
calculations using the local spin density approximation (LSDA). The spin form factors of

the materials are also calculated. The results agree with previous theoretical works using

an embedded cluster model. The calculated spin form factor of LazCuO4 is consistent
with the few experimental data currently available, while the results for LazNiO4 show a
large discrepancy between theory and experiment. We question the accuracy of the
experimental results of LazNiO4 and call for more experiments to settle the issue.

In the third part, crystal Hartree-Fock calculations are carried out for LaMnO3.
Our main focus is on the magnetic and orbital orderings, the effect of the crystal
distortion from the cubic perovskite structure, and the analysis of the projected density of
states. In addition, we also ﬁnd unexpected charge associated with the Mn ions. The
effective spin Hamiltonian obtained by mapping the energies of different magnetically
ordered states is consistent with experimental data of the spin wave dispersion, while it is
qualitatively different from that of a LSDA calculation which was also claimed to be
consistent with the experiment. We show that the current experimental accuracy is not
enough to distinguish these two theories.

In the fourth part, I report a study on possible extensions of the double-exchange
(DE) model for the doped manganites in the hope of understanding the recently observed
softening of the spin wave dispersion near the Brillouin zone edge. I also argue that most
models in the literature are ﬂawed by assuming uncoupled normal mode vibrations of the
MnO6 octahedral clusters in the crystal. As a ﬁrst trial, I consider several possible
electron-phonon couplings based on the single-band DE model. The results show that the
spin wave states are robust and the spin wave dispersion gets very slightly hardened in
the presence of the electron-phonon couplings. This suggests that the observed softening

of the spin wave dispersion may be beyond the scope of the single-band DE model.

 

 

 

T 0 My Parents and My Wife

iv

T 0 My Parents and My Wife

iv

ACKNOWLEDGEMENTS

It has been six years since I came to Michigan State University, East Lansing, in
1994. This is a bit longer than I would have liked to get my Ph.D. study done. But the
delay is mostly caused by a pleasant reason: The life here in East Lansing is so cheerful
that I simply lost the motivation to rush for the degree. ©

I was very fortunate to join the group of Profs. T. A. Kaplan, S. D. Mahanti in
Physics and J. F. Harrison in Chemistry departments. All the three professors are
outstanding scholars and excellent teachers. They are all very knowledgeable in Physics
in general but they each also have their own specialties. The discussion in our weekly
group meeting is often full of sparks of intelligence. That I can learn from the best the
various subjects of physics is truly the most wonderful thing about my study here. I
would like to express my deep gratitude to them. They are essentially all my advisors,
although my formal advisor is nominally only Prof. Kaplan.

In addition, I want to thank Prof. Kaplan for being so kind to me and giving me
the maximal degree of freedom. He never pushed me. He is just there to give advice and
to compliment what I have achieved. I also want to thank him and his wife, Pat, for their
hospitality. They are a large part of the reason why my life here is so cheerful.

I am grateful to Stephanie Holland, Debbie Simmons, and Janet King, who are
secretaries in our department. They are lovable for their professional assistance as well

as their enthusiasm in helping students.

Finally, I thank the Center for Fundamental Materials Research of MSU for
supporting the research these years. Without the crucial ﬁnancial support, this work

would not be possible.

vi

Contents

List of tables .. ..... .. ...................................................................... x
List of ﬁgures ............................................................................ xi
Introduction 1

Small-bandwidth perturbation theory for highly-covalent

Mott Insulators 4
2. 1 Introduction ...................................................................................... 4
2.2 The model, effective-Hamiltonian perturbation theory, and

2.3

2.4

2.5

TSDA ................................................................................................ 8
2.2.1 Model Hamiltonian ................................................................................ 8
2.2.2 Formal perturbation theory .................................................................... 9
2.2.3 TSDA ................................................................................................... 11
Achievement of convergence via l-body transformations ............ 15
2.3.1 Site localization .................................................................................... 17
2.3.2 Cell localization ................................................................................... 22
2.3.3 No localization ..................................................................................... 25

2.3.4 Comments on the TSDA solutions for the three localization choices .28

Further improvement via 2-body transformation .......................... 29

Summary and discussion ................................................................ 34

Crystal Hartree-Fock calculations for LazNiO4 and
L32Cll04 36

vii

3.1

3.2

3.3

Introduction and methods ............................................................... 36
Results ............................................................................................ 40
3.2.1 Heisenberg exchange parameter I ........................................................ 40
3.2.2 Spin density and spin form factor ........................................................ 41
3.2.3 Band structure ...................................................................................... 52
Summary and discussion ................................................................ 56

Electronic structure of LaMnO3 in the ab initio crystal

Hartree-Fock approximation 58
4.1 Introduction .................................................................................... 58
4.2 Methods .......................................................................................... 60
4.3 Results ............................................................................................ 62
4.3.1 Magnetic properties ............................................................................. 62
4.3.2 Spinless charge backﬂow and density of states ................................... 68
4.3.3 Orbital ordering .................................................................................... 73
4.4 Discussion and summary ................................................................ 76

Effect of oxygen vibration on spin waves in the double-

exchange model 78
5. 1 Introduction .................................................................................... 78
5.2 A ﬂaw of popular models in the literature ..................................... 81
5.3 Model studied ................................................................................. 84
5.4 Results ............................................................................................ 86

5.4.1 Pseudo-Bom-Oppenheimer approximation ......................................... 86

viii

5.4.2 QM treatment of the oxygen vibrations with cut-off in the number

of phonons ............................................................................................ 89

5.5 Remarks and conclusions ............................................................... 90

6 Concluding Remarks 92

Appendices 96

A A physical interpretation of the TSDA equations .................................... 97

B Solving the TSDA equations .................................................................. 100
C On the determination of the Heisenberg J from electronic structure

calculations ............................................................................................. 1 09

References 1 1 2

ix

3.1

4.1

4.2

4.3

4.4

List of Tables

Comparison of the crystal HF, the cluster HF and the experiment
results of the spin form factor of KNiF3 .............................................. 52

HFA and LSDA energies of LaMnO3 with various magnetic
orderings .............................................................................................. 63

HFA J’s in the effective spin Hamiltonian of LaMnO3 ...................... 65

Comparison of the HFA, LSDA and experimental J’s in the spin
Hamiltonian of LaMnO3 ...................................................................... 66

HFA energies of different magnetic and orbital ordering states for
the cubic structure ................................................................................ 75

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

2.9

2.10

3.1

3.2

3.3

3.4

3.5

3.6

4.1

List of Figures

Cu02 plane ............................................................................................. 6
The cluster Cu207 ................................................................................ 15
The nth approximation to J vs n according to straightforward

perturbation theory .............................................................................. 16
TSDA orbitals for site localization ..................................................... 20
Perturbation result for site localization ............................................... 21
TSDA orbitals for cell localization ..................................................... 23
Perturbation result for cell localization ............................................... 24
TSDA orbitals for no localization ....................................................... 26
Perturbation result for no localization. ................................................ 27

Perturbation result for cell localization l-body + minimal 2-body

transformation ...................................................................................... 33
The centered-tetragonal structure of LaNiO4 and LaCuO4 ................. 38
The spin density of LaNiO4 and LaCuO4 in the (001) plane .............. 43
The spin density of LaNiO4 and LaCuO4 in the (100) plane .............. 44
The crystal HF spin form factor of LaNiO4 and LaCuO4 ................... 47
Comparison of the crystal HF and cluster calculations with

experiments on the spin form factor of LaNiO4 and LaCuO4 ............. 49
The band structure of LaNiO4 and LaCuO4 ........................................ 54
The cubic perovskite crystal structure ................................................. 61

xi

4.2

4.3

4.4

5.1

5.2

5.3

The spin wave dispersions of the 3-J and its effective 2-J spin

Hamiltonians ........................................................................................ 67
The (projected) density of states of LaMnO3, with AAF ordering ..... 71
The spin density on the xy-, yz-, and zx-planes at a Mn with the

d(3x2-r2) and three tzg orbitals occupied .............................................. 74
The Mn-O linear chain ........................................................................ 82

Density of states of the DE model + electron-phonon interactions
for a 4-site, 3-electron ring .................................................................. 88

The spin wave dispersion of the DE model + electron-phonon
interactions for a 4-site, 3-electron ring .............................................. 89

xii

Chapter 1

Introduction

The fundamental physical law at the atomic length scale and above has long been
known; namely, the Schrodinger Equation with Coulomb interaction between charged
particles plus secondary relativistic corrections like spin-orbit coupling. For electronic
systems, a wonderful feature of the theory is that often the electrons can be treated as
non-interacting particles moving in an effective potential to a good approximation,
provided the proper statistics are obeyed. This approximation is a kind of mean-ﬁeld
theory and is often called a single-particle picture, where the electron-electron interaction
is treated in an average way. In this picture, one ﬁrst calculates the energy levels of an
electron in the presence of an effective potential. Then the states of a system with N
electrons will be just all possible conﬁgurations obtained by ﬁlling the N electrons in the
set of energy levels in different ways obeying the Pauli exclusion principle. For example,
the band theory of solids with periodic structures and the interpretation of the atomic
emission/absorption spectra by transitions of electrons between energy levels are two of
the most well known and successful single-particle theories.

Corrections to the single-particle picture are often called correlation or many-body
effects and cause only quantitative changes in many physical systems. However, in some
circumstances electron correlation can surprise us by leading systems to unexpected new
states or phases. The Mott insulator [1] and the fractional Quantum Hall effect [2] are

two of the famous examples. Simply put, correlation is probably the very source of most

of the radically new physics at the atomic level that we are still discovering everyday
even after the fundamental physical law was discovered so long ago in 1926 by E.
Schrodinger.

In this dissertation, I study the perovskite La-cuprate (and isostructural La-
nickelate) and La—manganite. Both have been the focus of intensive studies by many
researchers [3-4]. When properly doped with some elements, the cuprate reveals high-Tc
superconductivity and the manganite reveals the so-called colossal magnetoresistance.
The interesting physics of the two systems is generally believed to be dictated by the
strong interaction among the transition metal (Cu, Ni, Mn) 3d electrons and the oxygen
2p electrons. To understand the observed physical properties of these systems,
correlation effects must be taken into account to some extent. For both systems, I study
model Hamiltonians commonly used in the literature and also carry out ab initio
calculations using the crystal Hartree-Fock (HF) method. For the latter we use the
unrestricted HF method, which allows a sensible treatment of the ground state even when
correlation effects are large. The detailed motivation for the selection of different models
and the calculation methods will be given in the respective chapters.

In Chapter 2, I study the 3-band Hubbard model for cuprate and nickelate using a
perturbation approach called the effective Hamiltonian method. In Chapter 3, I carry out
crystal HF calculations for the cuprate and nickelate and compare various physical
quantities obtained with those from the Local Spin Density Approximation (LSDA) [5]
and those from an embedded cluster approximation [6]. In Chapter 4, I carry out the
crystal HF calculation for the manganite and compare the results with those of LSDA

calculations. Some quantities, like the order of the oxygen p band and the manganese d

band, and the charge of the manganese ions in the crystal, are found to be dramatically
different from those of LSDA calculations or those of the commonly conceived picture.
The signiﬁcance of the discrepancy will be discussed. In Chapter 5, I report a study on
modiﬁed Double-Exchange models for the doped manganite motivated by the recent
intriguing experimental observation of phonon—induced softening of the spin wave
dispersion near Brillouin zone boundaries. In Chapter 6, some concluding remarks of the

work of this dissertation are given.

Chapter 2

Small-Bandwidth Perturbation Theory
for Highly-Covalent Mott Insulators"

2.1 Introduction

The question of how to set up a perturbation theory for Mott insulators like NiO,
KNiF3, F e203, was considered in 1959 by P. W. Anderson [7]. The object was to take
advantage of the obviously localized property of the magnetic electrons to construct the
theory, where the interatomic overlap of localized orbitals is treated as the small
parameter. Such a theory leads to low-lying states that are governed by the Heisenberg
spin Hamiltonian, in leading order. At the time, it was known that most of the low-
temperature properties of these materials could be understood on the basis of such a
Hamiltonian.

Anderson, on considering how to determine the localized orbitals, stated that there
exist "exact" localized orbitals (orthogonal, and thus called Wannier functions) which
will make the perturbation theory converge well. These orbitals were required to be
nonmagnetic (spin-up and -down orbitals the same), as is natural, since the magnetic
properties were to be derived from the spin Hamiltonian resulting from the perturbation

theory (p.t.). Anderson gave an explicit prescription, namely Wannier functions derived

 

‘ This work has been published partially in “T. A. Kaplan, S. D. Mahanti, Y—S Su, and K. Kubo, J. Appl.
Phys. 79, 6433 (1996)” and partially in “Y-S Su, T. A. Kaplan, and S. D. Mahanti, Phys. Rev. B 56,
15596 (1997)”.

from Hartree-Fock theory in which the spins of the magnetic ions were all parallel. It
was noted. some time later [8] that there was an inconsistency with this, since Hartree-
Fock eigenstates with up and down spins differ in such a situation. [9] An alternative
approach was studied [8] based on the so-called Thermal Single Determinant
Approximation (TSDA) [10], a variational generalization of thermal Hartree Fock theory.
Solving this TSDA for non-magnetic localized solutions at sufﬁciently low temperature
amounted to minimizing the mean detenninantal energy, averaged over all spin
conﬁgurations of the magnetic ions, an approach adopted by Gondaira and Tanabe [11].
It was shown for a small cluster model [8,11] which we can call H-He-H (hydrogen-
helium—hydrogen) and for a 1-D crystal model [8] (the corresponding H-He chain) that
this TSDA choice indeed improved, quite dramatically, the convergence rate of
perturbation theory for the Heisenberg exchange parameter J, through 4th order in the H-
He overlap.

The overlap, or more precisely, the amplitude for cation-anion hopping, is small
in the highly ionic materials considered at the time, e.g. NiO, KNiF3, so the pi. was
expected to converge rapidly. For this reason higher order perturbations were not
investigated. However, materials that have recently become of interest in connection
with high-temperature superconductivity, namely the cuprates and the related LaQNiO4,
are thought to be much more covalent. Doubts about the validity of perturbation theory
in this connection, speciﬁcally for the 3-band Hubbard model of a Cqu plane (Fig. 2.1),
have been expressed [12-14]. And, in fact, later we have shown [15] that straightforward

perturbation theory does not converge, speciﬁcally for the version of the 3-band model

due to Hybertsen et al. [12], which is rather realistic for many physical observables [16],

and is quite similar to other versions of the 3-band model [17-20].

0 O O O O

O O O O O O O 0 O

[0 Copper 0 Oxygen]

Figure 2.1: CuO2 plane.

Clearly then it is of interest to try the TSDA approach to ﬁnd new localized
orbitals, and to test the implied new perturbation expansion for convergence. In fact, the
question of convergence of p.t. within the 3-band model has been addressed [13,14], and
rapid convergence of a new p.t. was claimed. But the transformation in these works
involved both 1- and 2-body transformations. (Finding new localized orbitals is an 1-
body transformation.) It is important to understand how far one can go with the l-body
transformation alone, because of its relative conceptual simplicity and its historical
interest. This is the content of this chapter.

We have been motivated to explore these questions by difﬁculties in connection
with calculating the ground state spin density in the cuprates and LaQNiO4 [21,22]. The
type of perturbation theory discussed above, which in effect separates the spin and space

(or charge) degrees of freedom, might be the only hope at present of making tractable ab

initio calculations of this observable. This is due to the fact that the reduction of the
ordered spin due to zero-point or quantum spin ﬂuctuations (QSF), which is large in these
cases, constitutes a macroscopic correlation effect. That is, (i) in principle it requires the
thermodynamic limit for its very existence, and therefore one needs very large systems to
estimate the QSF [23], and (ii) the wave function representing such a state requires a
linear combination of many Slater determinants. Thus, even with the remarkable
advances in computational many-body physics (e. g. Refs. 24 and 25), the present problem
is still beyond the reach of those approaches. Also, standard band-theoretic approaches
fail to capture essential physics (the local spin density approximation misses the
antiferromagnetism and insulating property of LaQCuO4 [26], and the Hartree-Fock
method (unrestricted) misses the QSF).

Since the work reported here involves calculations only for a small cluster, the
spin density of a crystal can't be addressed (the spin density in a ﬁnite cluster with an
even number of electrons is zero, since the ground state is a singlet). Instead we will be
investigating the Heisenberg exchange parameter J as an indicator of convergence of the
perturbation expansions. In Sec. 2.2 we deﬁne the Hamiltonian to be considered, review
effective-Hamiltonian perturbation theory, and the TSDA. In Sec. 2.3 we describe the
three types of single-particle transformations considered, which we call site localization,
cell localization, and no localization (the names are after the transformed ligand p
orbitals surrounding the magnetic cations), and ﬁnd the best in each case according to the
TSDA. Also in this section, J is calculated to high order in p.t. to examine its
convergence properties. In Sec. 2.4 we carry out the particular 2-body transformation

suggested by the work of Ref. 13, for the special model they considered, as a check; we

also carry this out for the model of Hybertsen et al. [16]. A summary and conclusions are
given in Sec. 2.5. A brief overview of the results obtained can be seen by glancing at
Figs. 2.3, 2.5, 2.7, 2.9, 2.10 which show the non-convergent series of straightforward p.t.,
and the rather dramatic improvement obtained in each of the four modiﬁcations,

respectively.

2.2 The model, effective-Hamiltonian
perturbation theory, and TSDA

2.2.1 Model Hamiltonian

We consider the 3-band Hubbard Hamiltonian (sometimes referred to as the

Anderson lattice model) as parametrized by Hybertsen et al. [12,16]:

H=H,+HU+HK (2.1)
where
HI = a}: n,” +de Eager; +h.c.)+tpp Eager, +h.c.), (2.1a)
i <il>a <ll'>0’
d I d
HU =UdZniTnE]+UpZnﬁnlﬂ+UPdZni n1”, (2.1b)
i I <iI>
and
HK = KP, chgcfocgc; + K”, Zc,’;+c,’;.c,’.’;.c,’.’a . (2.1c)
<il>aa’ <Il'>00"

V+

cm creates a hole in a Wannier function W”. of type v at site i with spin 0,

10' i0 ia‘

n,” = Z n}; , n.” = cwcv The orbital at a copper site is (1.2-,2 ; at each oxygen site there

is onep orbital. The parameter values are s= 3.6, tpd = -1.3, tpp = -0.65, Ud = 10.5, Up
= 4, Upd = 1.2, Kpd = -0.l8, Kpp = -0.04, all in eV [27]. Note that for the convenience
of having the d—p hopping parameters all the same, and similarly for the p-p hopping, as
in Eq. (2. la), one has to set up the underlying Wannier ﬁrnctions in the way that they may
change signs under the crystal or cluster symmetry. The signs corresponding to orbital

phases are chosen as follows: if die—)9 exists at a particular Cu, then at the O's

immediately to its right and left (along the x-axis) the orbitals are px and -px

respectively; similarly, the orbitals at the nearest 0's below and above are py and —py

respectively. The remaining phases are determined by having the nearest neighbor d—p
overlap always negative. We limit ourselves to the case where the number of holes

equals the number of Cu sites, i.e. the Mott insulator limit.

2.2.2 Formal perturbation theory

The form of the perturbation theory used involves an effective Hamiltonian.

Given
W = E‘P , (2.2)
then the effective Hamiltonian Heffsatisﬁes

Hcﬂ(E)P‘IJ = EP‘I’, (2.3)

namely

H,,,(E) = P{H+H[Q(E—H)Q]"H}P. (2.3a)

Here, P = l-Q is a projection operator which projects ‘P onto a given subspace. The

"inverse of the corner" [28] G(E) E [Q(E-H)Q]'1 has to be understood as the matrix with

zeroes everywhere except in the Q-subspace, where the matrix is the inverse of the Q-

projection of E - H. Similarly, Q E Q l is the matrix with zeroes everywhere except for
the Q-subspace, where it is the unit matrix. We choose the P-subspace as that deﬁned by
the ground states of some chosen unperturbed Hamiltonian H0, with eigenvalue E0.

Then PHQ = P(H0+V)Q = PVQ, and Eq. (2.3) can be rewritten as
P{V + VG(E)V} P? = (E PLP, (2.4)

where 6E = E-EO. Expanding G(E) in powers of (E -V, and substituting the full left—

hand side of Eq. (2.4) for oEP‘l’ everywhere it appears in the expansion, Eq. (2.4)
becomes [29]
P{H0 + V + VGOV + (VGOVGOV — VGOZVPV) + - - -}P\I’ = EP‘I’. (2.5)

where

G. = [Q(Eo — Hoar".

Straightforward perturbation theory takes the terms in H that involve hopping (tpd
and tpp) plus the exchange terms (HK) as the perturbation V (and H0 = H — V). We
studied [15] this p.t. in the case of a crystal (inﬁnite Cqu plane), where we obtained the
ﬁrst three terms (3rd, 4th and 5th order) for the nearest Cu—Cu exchange parameter in the

Heisenberg Hamiltonian [30]. We also studied this p.t. on the cluster CuzO7, using the

10

embedding scheme described in Ref. 16, and showed that the corrections to J increase in

magnitude and oscillate in sign through 14th order [15]. See Fig. 2.3 below.

2.2.3 TSDA

To help in choosing transformations which are to deﬁne new p.t.'s (i.e. new
choices of V), we will appeal to the TSDA [10] which, essentially, ﬁnds the “best” 1-
particle wave ﬁmctions according to the free energy variational principle (for the

canonical ensemble)

F(p) a tr(pH) + ﬂ“tr(plnp) Z — ,6" 1ntrexp(—,BH). (2.6)
Here p is any density operator and ,8" is Boltzmann's constant times absolute

temperature. Writing ,0: Zn” exp(—,BHa) with Z“ 2 tr exp(-ﬂHa), one can consider H,
as an approximate Hamiltonian. The TSDA is deﬁned by taking Ha as a general real

function of the occupation number operators for some complete orthonormal set of one-

particle wave functions luv and then varying the functional form of Ha and the \iIv’s to

minimize F. A complete set of eigenstates of Ha can be chosen as single Slater
determinants with all the various subsets of the wv’s occupied. It turns out [10] that the
best function Ha for a given set of wv is simply the determinantal energy of the exact
Hamiltonian H with respect to that set of WV. In other words, if H is expressed with

respect to an arbitrary complete orthonormal set of one-particle wave functions wv as

H = z hvﬂc:c# +% Z vV/LMCCC/jckc,I , (2.7)

141 milk“

11

where the subscripts v, u, .. each label both the spatial and spin quantum numbers and the

quantities kW and v are matrix elements (which, by properly choosing the phase

via/Ix
factors of the rpv’s, are taken to be real in this paper), respectively, of the one- and two-

body terms in the exact Hamiltonian H, then Ha will be

Ha = zhvvnv +%Z(vvp,vp — vv;1.;1v )nvny (28)

V}!
which contains terms linear and quadratic in the occupation-number operators n,

corresponding to the wv’s. (The distinction between this and thermal Hartree F ock theory
is that in the latter, H" is restricted to be linear in the occupation numbers, i.e. a one-

electron operator.) Ha turns out to be equal to the diagonal part of H with respect to the

single-determinantal states constructed from having the one-particle states u, v,
occupied in all possible ways. As an example, for the particular model Eq. (2.1) the
determinants with the original Wannier functions occupied in all possible ways have

energies given by

Ha = .92 nip + HU + KW, annf; + KW 2,2,13,12,30 , (2.9)

<il>0' <ll'>0

the values of the occupation numbers being deﬁned by the determinant being considered.

Furthermore, stationarity of F ( p) with respect to the \[Iv’s implies [10]

< n, —n# > hvﬂ +Z<(nv —nﬂ)n, >(vW1M — Van”): 0, (2.10)
,1

where the brackets < -- > mean average over the canonical ensemble with H" as

Hamiltonian, i.e. for any operator A, < A > means

12

”(e—W“ A)
< A >5. —-—-——. 2.11
”(e—3”“) ( )

A physical interpretation of Eq. (2.10) is presented in Appendix A.

It should be noted that TSDA allows localized solutions in crystals, or
analogously in molecules with high symmetry, in contrast to thermal HF approximation
[10,31]. This distinction is important in the Mott insulators we are considering. Further,
the form of the localized orbitals in TSDA is dictated by the Hamiltonian, H. As far as
we are aware TSDA is the only variational theory in the literature with this formal
property.

It would be natural to ask here why not simply use the TSDA as a variational
approximation to the Hamiltonian H of Eq. (2.1)? The answer is that it is known to give
a poor description of the low—energy physics when the hopping is very small (actually,
speciﬁcally in the case of the single-band Hubbard model), so that the low-lying energies
are accurately given by the Heisenberg model. E.g., it gives a Curie law susceptibility
rather than the correct Curie-Weiss law at temperatures above the antiferromagnetic
transition temperature, and it gives the latter about an order of magnitude too large. [32]
Thus in the present work we are not considering TSDA in this way. Rather we are using
it as a formal device to help discover a l-body transformation that (hopefully) will lead,
via a rapidly convergent perturbation theory, to the proper effective Hamiltonian (it is the
solution of the latter that will determine accurately the low-energy physics of the model).
The early success of this approach discussed in the Introduction has motivated our

present investigation.

13

Thus in our TSDA considerations here, [3 will be considered as a parameter whose
optimum value [30 would ideally give the “fastest” convergence of the perturbation
expansion Eq. (2.5). As is apparent from the above TSDA equations, the solution, i.e. the
guy’s, will depend on 00. If we take a ﬁnite number of terms in the perturbation expansion
(as one would do), then the resulting effective Hamiltonian will thus depend on [30. For
thermal properties, e.g., [50 would then enter the partition function

Z" exp[—En(ﬂo) / kBT] where En(,80) are the eigenvalues of the spin Hamiltonian and T

is the physical temperature. When working within the TSDA to ﬁnd the \pv’s, we will
think of B" intuitively as temperature to help ﬁnd a “good” TSDA solution, i.e. a solution

that gives rapid convergence of p.t.; in fact we will tend to consider large [3, since the

effective Hamiltonian is designed to give the low-energy states of H. (One could

probably not explore the whole range of [3 in practice.)

One should realize that even for the small-basis-set model of the small cluster
treated below, the 1-body Hilbert space is too large to explore completely. So one needs
some guide. There are alternate approaches to determining localized orbitals, widely used
in quantum chemistry; see Ref. 25 and references contained there, and Ref. 33. These
conceivably could be useful in the magnetism problem considered here, and should

probably be studied in this connection.

14

2.3 Achievement of convergence via 1-body
transformations

As mentioned above, we applied [15] pt to the Hamiltonian in the original
Wannier-function representation, i.e. Eq. (2.1), where the perturbation V incorporates all
the hopping and exchange terms and H0, the rest. For the crystal (an inﬁnite Cqu plane)
we obtained the ﬁrst 3 contributions to the nearest-neighbor exchange parameter J (they
occur in 3rd, 4th and 5th order), with no sign of convergence [34]. We also calculated J
to very high order for the embedded cluster CuzO7, shown in Fig. 2.2, with two holes.
The results of the latter work are given in Fig. 2.3. Plotted there is

J“) = ZJ", (2.12)

":1

where J” is the contribution to J (the splitting between the lowest singlet and triplet) from
n-th order p.t.. The non-convergence is apparent, the result oscillating with increasing

amplitude about the exact value (straight line) calculated by direct diagonalization.

O 0
P2 P5
0 O O O 0
P3 d1 P1 d2 P6
0 0
P4 P7
Figure 2.2: The cluster CuzO7. The symbols at each site serve both as labels and as

orbitals.

15

 

 

o 2 4 6 ' 8 1o 12m14
order

Figure 2.3: The nth approximation to J vs 12 according to straightforward perturbation
theory. The horizontal line shows the exact value of J. From Ref. 15.

16

To follow up on the idea [7] that an appropriate change of l-electron basis states,
with its consequent change in the partition of H into unperturbed and perturbed pieces,

might improve the convergence, we consider the unitary linear transformation,

qr, = Z Aye-1.0. (2.13)
.I

An essential assumption we make about this transformation is that it maintains
localization and symmetry of the d states and that all states should be nonmagnetic. The
reason for this is that we want to maintain the property of the straightforward p.t. that the
low-lying energies are governed by a Heisenberg Hamiltonian with the symmetry of the

nuclear structure. Av is a unitary matrix, which we take as real, since we assume the new
Wannier functions are real. The spin independence of AD. embodies the assumption that
the new Wannier functions are nonmagnetic (they are each a product of a spatial orbital

times a spin function, or or [3, the spatial orbital being the same for either spin). These

properties are assumed of course to hold for the original Wannier functions, i.e. those

created by the cg.

2.3.1 Site localization

We begin with our ﬁrst scheme which we call site localization, deﬁned by the
requirement that the new Wannier functions have the same symmetry properties as the
original ones. Referring to Fig. 2.2, let 6V and oh be reﬂections in the vertical and

horizontal symmetry axes, respectively. Then the d states are taken to satisfy

d" = _0.vd?t (2.14)

17

(a similar relation being assumed, of course, for the original, unprimed, orbitals; the sign
choice follows the discussion in Sec. 2.2.1). Similarly, we require

_ I I _ I I_ I
p2 —_O-vp5’ p3 _ —0-vp6, p] — _0-\’pl , etc.

I r r r t r (215)
p2 = 011194 , p5 : 01:1)7 a 1’3 : 0'11P3 2 etc'

Eqs. (2.14) and (2.15), satisﬁed also by the unprimed Wannier functions, deﬁne "site
localization".

One can see that these requirements plus orthonormality are not sufﬁcient to
completely determine the transformation. A simple example is the 3-site cluster
dl — pl — d2 or CuzO: Then the two symmetry requirements, Eq. (2.14) and pl’ 2 —a‘,p,',
give

d1, : ]V1(d1+ A1d2 + Azpi)
d,’ = N,(A,d1 +ai2 +A2p1).
pi = N2[B(d1+d2)+p1]

The Ni are normalization constants, determined directly in terms of the three coefﬁcients,
A 1, A 2, and B. Orthogonality clearly gives two equations, leaving one independent
coefﬁcient. For the cluster CuZO7, it turns out that there are 9 independent coefﬁcients. It
is useful to see this in detail. We can write the transformation as

‘11, = N1[d1+ A1d2 + A2p1+ A3(p2 + p4)+ A4P3 + A5(P5 + p7)+ A6p6]

pi = NziBr(d1+d2)+p1+Bz(P2 +p4 +P5 +p7)+Bs(p3 +P6)l
Pi = N3(C1d1 + Czdz + C3p1 + F2 + C4p3 + C5p4 + C6ps + C7176 + C8177)
p5 = NilDid. + Dzdz + D31). + D4092 +p4)+p3 +D5(p5 + p7)+D6p6]

(2.16)

The ﬁrst two and the last forms are dictated by the horizontal reﬂection symmetry; the

remaining ﬁve primed orbitals are obtained from those given here by operating with 0',

l8

and oh. Aside from the M there are 23 coefﬁcients A,, B,, Ci, Di. One can see that there

are 14 orthogonality conditions, leaving 9 free coefﬁcients.
To determine these coefﬁcients we turn to the TSDA. For every pair of states [1
and v, there is a corresponding TSDA equation (2.10) to be satisﬁed. Since both It and v

are spin independent, the states labeled v and u in the TSDA equation (2.10) must
correspond to the same spin, each such state also corresponding to an atomic site using
the site-localized orbitals deﬁned by Eq. (2.16). It is easy to see that if two sites are
equivalent by symmetry, e.g. the two Cu sites or the two O sites labeled p2 and p,, then
the corresponding equation vanishes identically, giving no information. Thus from Eq.
(2.10) the only equations with content are those where the sites in v and u are non-
equivalent, e.g. a Cu and an oxygen, or pz-p3. The number of such non-equivalent pairs
which are independent can be seen directly to be 9, the number of free variables to be
determined, i.e. there is the correct number of equations to determine these unknowns.
Because we are looking for basis functions on which a perturbation theory will be
based with the purpose of deriving a spin Hamiltonian which will yield the low-lying
energies and magnetic properties, we look for low-temperature nonmagnetic solutions
such that for every spin up state there is a spin down state with the same spatial orbital.
In the same spirit, we consider Kdd<<B" <<e, where Kdd is the splitting of the lowest
levels within the single-determinantal energies (It was 0 for the original Hamiltonian, Eq.
(2.1), but as we’ll see later it has value of about 0.03 eV after the 1-body

transformations). Unfortunately, the low temperatures of interest lead to large disparity

in the sizes of different contributions to the thermal averages <n V>, etc., this causing

19

serious numerical difﬁculties. An approximation has been introduced to get rid of the
exponential factors in, and thus simplify, the TSDA equations. Details of the calculation
and discussion about its validity are given in the Appendix B.

The results obtained after solving these approximate TSDA equations are shown
in Fig. 2.4 which gives an idea of the “shape” of the new orbitals relative to the original
ones. The Hubbard Hamiltonian, Eq. (2.1), is then rewritten in terms of the new orbitals
and thus the parameters in the model get renormalized. Note that the price to pay for the
transformation is mainly that the Hamiltonian, which was restricted up to nearest-
neighbor terms in Eq. (2.1), now resumes back to a general form like Eq. (2.7). We take
together all those terms which are associated only with number operators as unperturbed
Hamiltonian H0, and the rest of the transformed Hamiltonian as perturbation [35]. The
exchange parameter obtained from the new perturbation expansion is seen in Fig. 2.5 to

converge (to the exact value). The errors in 2nd and 4th order are 15.4% and 15.7%,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

respectively.
.24 400 .20 .20
[ .24 .88 .22 -.02 -.02 ] | —.15 -.28 .80 -.28 —.15 |
.24 -.00 .20 .20
d1, pi
.73 -00 -.55 .03
.57 -.26 -.08 .03 .01 | [ .52 .06 .34 -.08 -.04|
—.27 -.00 -.55 .03
05 p;

Figure 2.4: TSDA orbitals within the site localization assumption--the numbers indicate
the amplitudes of the transformed orbital at the original orbitals. (The positions of the
numbers correspond to the CILZO7 cluster in Fig. 2.2.)

20

 

 

 

 

 

 

0.15 ' T T 7
i— _____
0.1 i
t, 0.05;, [I
0]
-0.05
0 2 A 4 A6 8 101 ‘12
order

Figure 2.5: Perturbation result for site localization.

21

If

2.3.2 Cell localization

To see if we could do better, we considered another type of symmetry assumption,
named cell localization with respect to the transformed p orbitals. This was inspired by
the work of Zhang and Rice [36], which also played an important role in the work of Refs
13 and 14. Instead of demanding that each new p orbital is localized at a site, we require
each of them, with the exception of the central p, to be localized either within the cell on
the left (the d plus the 4 p’s) or the one on the right, maintaining a certain symmetry,
analogous to Eqs. (2.14) and (2.15). There is no distinction between site and cell

localization for the d orbitals, so they obey Eq. (2. 14):
dl’ = —o*vd2' (2.17)

The distinction between site and cell localization is seen in the following symmetry

requirements for the (cell localization) p orbitals:

I I I - .
p' =—O-vpir’ pi“) =O-lipi(r), WIthl:1’2

I

pi = -0.p§ . pi“ = ﬁat"), (218)

C

P = «7.00, p“ = 0119(-
The superscripts l and r stand for left and right, and c for central. Thus, for example,
rather than 07, taking orbital p; into an orthogonal orbital p], as in the site localization, it

takes p,’ into itself. The explicit transformation in this case is

d,’ = N1[dI +A,d2 + Azp, + A3(p2 +p4)+ A4!)3 +A5(p5 +p7)+ A6126]

1)” = B,(d. +6172)+Bzrl +B3(pz +04 +05 +p7)+Bi(p3 +126)

pi = Cidr + Czdz + C3171 + C4(P2 +P4)+ C5173 + C6(P5 +107) + C7136 (2-19)
p; = p,’ with C, replaced by D,

pi =7'7(pz -p.)

22

Following a counting similar to that done in the site localization case, one can show that,

for the cell localization, the number of free coefﬁcients after orthonormalization and the

number of independent nontrivial TSDA equations are equal, as in the previous case.

The analogous ﬁgures to the site localization case for the cell localization are

shown in Fig. 2.6 and Fig. 2.7. The errors in the exchange parameter in 2nd and 4th order

are 15.4% and 13.6%, respectively. It shows that the cell localization is a slightly better

choice for the perturbation calculation than the site localization.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.24 -00 .20 .20
.24 .88 .22 -.02 .02] [-15 -.28 .80 -.28 -.15
.24 —.00 .20 .20
d; pr
.32 -.01 —.55 .03
| .80 -.36 -.11 .04 .01 ] | .52 .06 .34 -.08 -.04]
.32 —.01 —.55 .03
pl pi

Figure 2.6: TSDA orbitals for cell localization.

23

 

 

 

 

.3 0-051
0
l
—0.05 g g
0 2 4 6 8 10 12 14

order

Figure 2.7: Perturbation result for cell localization.

24

2.3.3 No localization

The third choice of the symmetry of the new orbitals is named no localization.
The term means releasing any localization requirement on the new p orbitals (but still
keeping the new 0’ orbitals localized). This is analogous to the case in a crystal where the
p orbitals may form a p-band while the d orbitals remain localized. The new orbitals are

required to satisfy the following symmetry properties:
d,’ = —ovd2' (2.20)

which is the same as in the previous two cases, and

[717+ = —0-vpi—+ a P? = 07,1): Withi = 19 2’ 3
pl.++ = O'Vpl.” , i“ = 031p:+ withi = l, 2 (2.21)
p" = —0.p", p" = ﬂip", pi" = 0.19"“, 19““ = ﬂip“.

where the ﬁrst sign in the superscripts denotes the symmetry under 0“,, and the second

sign, the symmetry under 07,. The explicit transformations for some of the orbitals are

listed below

d1, : Nl[dl + Aidz + A2p1+ A3(P2 +174) + A4103 + A5(p5 +p7) + Aspol

p.” = 31(61, +d2)+sz. +3300: +194 +05 +07)+B4(03 +196)

p; (pf) = pf“ with B, replaced by Cl.(D,.)

pr++ : E1(d1 _d2)+E2P1 +E3(P2 +174 —p5 -p7)+E4(p3 _p6) (2'22)
p2++ = p,++ with E, replaced by E.

p" = 5002 -p4 +195 -p7)

+——

.0 —%(pz -p4 -ps+p7)

25

The transformed orbitals and the perturbation results for the no localization case is shown

in Fig. 2.8 and Fig. 2.9. The errors at 2nd and 4th order in the exchange parameter J are

15.5% and 14.8%, respectively, not as good as the cell localization.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.24 -.00
|.24 .88 .22 -02 -02]
.24 -00
d;
.31 .31
|.37 -.35 .31 —.35 .37]
.31 .31
p?
.19 —.19
|.60 -.27 0 .27 —.60]
.19 -.19
p1”

Figure 2.8:

TSDA orbitals for no localization.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.03 .03
]-.45 -.10 .76 -.10 —.45]
.03 .03
p?“
-.35 -.35
[.37 -.02 .48 -02 .37]
-.35 -.35
p?”
.43 -.43
]—.33 -.13 0 .13 .33]
.43 -.43
.02”

..__h___..__. s . t . A . ﬂ -ﬁ.——‘—____ - . . A » ,_ A__.——_. - . _.—...

0 2 4 6 8 10 12
order

Figure 2.9: Perturbation result for no localization.

27

14

2.3.4 Comments on the TSDA solutions for the three
localization choices

Because the symmetry requirements in these three choices are different for p
orbitals but the same for d orbitals, and since the lowest energy states in which only the d
orbitals are occupied dominate at low T, it turns out that the transformed d orbitals for the
three localization choices are the same. After TSDA transformation, the energies of the
lowest two single-determinantal states of the half-ﬁlled CuzO7 cluster are —2.7718 and
-2.7430 (all in eV) for the spins of the two holes, one at each transformed d orbital, being
parallel and antiparallel, respectively. This can be compared with the energies using the
original orbitals, 0 for both parallel and antiparallel spin conﬁgurations, and with the
exact eigen-energies from direct diagonalization of the Hamiltonian matrix, -2.7738 and
-2.9000 for the triplet and the singlet, respectively.

After transformations, all parameters hvp and v in the Hamiltonian get

vallx

renormalized. For example, some new parameters resulted from the TSDA solutions are:
1;, = —0.403, U; = 6.910, Uéd = 0.113, K], = —0.029, (2.23)

where the lst, 3rd and 4th ones were absent in the original form of the Hamiltonian, Eq.

(2.1).

28

2.4 Further improvement via 2-body
transformation

Although the TSDA guided l-body transformations did lead to convergent
perturbation expansions as we had hoped, we were somewhat disappointed about the
relatively slow convergence rate compared to the achievement of Refs. 13 and 14, in
which transformation is not restricted to be 1-body. An interesting question was then
asked: What is the main factor which limits the convergence rate of 1-body-
transformation-induced p.t.? The answer is given below and a minimal 2-body
transformation is added to improve the convergence rate.

Our attempt to incorporate both 1- and 2-body transformations is inspired by the
work of Zhang and Rice [36], as well as J. H. Jefferson et al. [13]. On the CuO2 planes, a
cell is deﬁned as any one copper and its four surrounding oxygens. Turning on the
perturbation (p-d hopping) will induce hybridization among the copper d orbital and the
oxygen p orbitals. If the perturbation is small, the new orbitals would remain localized
within a cell. The copper d orbital is expected, by symmetry, to equally hybridize with its

four surrounding oxygen p orbitals. Let us deﬁne the symmetric p orbital of each cell as
P: = 52 p, ,
le{i}

where the summation index indicates that l is taken over the (four) p sites in cell i. It is
intuitive to expect that the largest hybridization will occur between the d 1 and pf . The
symmetric pf orbital has its analogue in our cell localization choice of the l-body

transformation for the Cu2O7 cluster, namely p,’ and p,’ , for the left cell and the right cell

29

respectively. Therefore our following investigation is based on the cell localization
choice.

Because the Hamiltonian of Eq. (2.1) conserves the total spin S and its projection
on the z-axis $2, the Hamiltonian matrix elements form isolated blocks for different
values of S,. From here to the end of this section we only consider the SZ = 0 case.
Through the process of the l-body transformations, the hybridization between the d and p
orbitals are actually taken into account. Because of the expected large hybridization

between the d and the symmetric p5 orbitals, we are interested in those states which

consist of only d,’ (dz') and p,’ (Pil- Among the total 81 states for the half-ﬁlled CuzO7

cluster (i.e. only two holes on it) with SZ = 0, this means the states of interest are

diidit>’

 

dﬁph>t Phph>a (224)

 

 

idﬁphl
and four similar states for the right cell. Hereafter we consider only the left cell because

the expressions for the right cell are exactly the same. After the l-body transformation,

the Hamiltonian submatrix in the basis set of Eq. (2.24) is

’ 5.483 —1.149 —2.239 —0.148‘
- 1.149 5.483 2.239 0.148
—- 2.239 2.239 4.055 0.984
— 0.148 0.148 0.984 11.922

(2.25)

 

 

(row and column indices correspond to the states in Eq. (2.24) in order). Two things are
worth noting. First, the off-diagonal matrix elements are large. This is due to d 1' and p1'
being mainly strong mixtures of d, and pf . This hybridization leads to large

renorrnalized model parameters for hole-transfer terms (between “'1' and p,’) in the

30

transformed Hamiltonian, and thus contributes to the off-diagonal matrix elements. It
indicates that these four states are far from the eigenstates and may lead to slow
convergence if adopted as unperturbed states in the perturbation calculation. Second, the
energy of the state of double occupancy of the d orbital, 4.055, the third diagonal element
of Eq. (2.25), is dramatically reduced from its original value 10.5, also a result of strong
hybridization. Each new 61 orbital contains part of the original 61 orbital and part of the
original p orbitals. Therefore double occupancy in the new d orbital only gains energy
from a fraction of the large (original) Ud, and that causes the large change of this energy
level. Direct diagonalization of the Hamiltonian matrix shows that the energy levels
associated with double occupancy at the d orbitals remain high (~12.7). Thus, in spite of
giving good approximate low-lying energy levels, the l-body transformation also leads to
some high-lying levels which are not at all good approximations of corresponding
eigenstates, and these high-energy eigenstates do play am important role in the
perturbation theory. In conclusion, it is the strong hybridization (covalence effect) and
the large on-site Coulomb energy Ud (correlation effect) together that prevents a l-body
transformation from generating good approximation for the whole energy spectrum.

In order to take care of the above limitation, we diagonalize the matrix of Eq.
(2.25) and use the 4 resulting eigenstates (the 2-body transformation comes in here), plus
the 4 counterpart states for the right cell and the remaining 73 (unchanged) states, as the
basis set to perform the perturbation calculation. In some sense, by this we allow
relaxation in the intermediate states of the p.t., where one hole hops into the neighbor cell
which has already been occupied by another hole. We call this scheme the minimal 2-

body transformation where “minimal” means that, among the various p orbitals, only

31

p1’(p,’) are involved in the 2-body transformation. The result is shown in Fig. 2.10. We
have also tried including p§(p2') in the 2-body transformation [p§(p3’) are out of our
consideration because they possess different symmetry from d1’(d2')]. However, it only
made a tiny change from that of the minimal 2—body transformation and so the result is
omitted here.

In Ref. 13, J. H. Jefferson et al. described another transformation scheme, which
is, in spirit, similar to our l-body + minimal 2-body transformation, and demonstrated
rapid convergence to 2nd order p.t. on a model with a restricted parameter set: 8:3, (pd:-
I, tpp=-0. 5, Ud=oo, and all other parameters equal to zero. (The straightforward p.t. for
this parameter set is divergent, too.) For comparison, we applied our method to the same
parameter set. The result showed that the l-body transformation alone fails to lead to a
convergent perturbation expansion. This is not surprising since the renormalization in
Ud, the double-occupation energy at d orbitals, caused by the l-body transformation is
fatal to the p.t. in the case of Ud=oa Nevertheless, we achieved similar precision as

theirs after we added the minimal 2-body transformation to the p.t..

32

 

0 2W4 6 8- 10 12 14
order

 

 

Figure 2.10: Perturbation result for cell localization l-body + minimal 2-body
transformation.

33

2.5 Summary and discussion

For simplicity, we sought 1-body transformations to convert the pi. from
divergent to convergent in the context of strongly covalent Cqu planes. Taking TSDA
as a guide, several l-body transformations were found, in the CuzO7 case, to achieve this
goal, but the convergence rate in each case was rather slow. Besides those suggested by
TSDA, we have also tried some other l-body transformations, of which the best one gives
errors of about 13% and 6% at 2nd and 4th order p.t., respectively, better than the best
TSDA result (cell localization, where the errors at 2nd and 4th order p.t. are 15.4% and.
13.6%). However, the convergence is still slow, particularly in higher order where the
results oscillate with slowly decreasing amplitude, similarly to TSDA (see Fig. 2.5, 2.7,
2.9). The details of these other transformations are omitted here for the above reason and
the fact that they were obtained in a rather ad hoc way, not as systematically as in TSDA.
A striking difference between the best ad hoc trial and the cell localization TSDA is the
following. Two of the transformed p-orbitals of the ad hoc trial (the one with errors 13%
and 6%) were deliberately guided to be close to the Zhang-Rice-type orbital [36] (where a

transformed orbital localized in a cell consists of four oxygen p-orbitals in the cell with

equal weight). On the other hand, the p[ and p,’ in the cell localization choice of TSDA
are further from the Zhang-Rice-type orbital, which is a result of the low symmetry of the
CuzO7 cluster plus the free energy minimization. This suggests that being closer to the
Zhang-Rice-type orbital helps the convergence. However, one can see that this issue

becomes moot when the cell localization TSDA is applied to the crystal, where the 4-fold

34

symmetry at a Cu site holds and therefore the cell localization TSDA will give the
Zhang-Rice-type orbital.

After all our considerations of l-body transformations, we guess that we are close
to their limit in connection with improving convergence. If one wants to do better, 2- or
more-body transformation is probably necessary. It is the large covalence and strong
correlation in the problem that limits the effect of l-body transformations: a 1-body
transformation designed to give some zero-order energy levels accurately will likely spoil
other levels, all of which enter p.t. in high order.

We then found that a minimal 2-body transformation dramatically improved the
convergence rate. Only transformed di-like and pf—like orbitals are essential to be taken
into account in extracting the Hamiltonian submatrix to solve for “good” two-particle
basis states. In comparison, our cell-localization 1-body + minimal 2-body
transformation scheme achieves similar (high) precision by second order p.t. as obtained

in the work of Ref. 13.

35

Chapter 3

Crystal Hartree-Fock Calculations for
LazNiO4 and LaZCuO4*

3.1 Introduction and Methods

L212NiO4 and LaZCuO4, which have essentially the same crystal structure, have
attracted considerable theoretical and experimental research efforts, especially since the
latter was found to become a high-Tc superconductor when properly doped with Sr or Ba
[3 7]. Above their respective structural transition temperatures TS, both materials are in a
tetragonal structure [space group D4,,17 (I4/mmm), see Fig. 3.1 below]. Below TS, they are
slightly distorted to an orthorhombic structure. Each of them is an antiferromagnetic
(AFM) insulator below a Néel temperature TN ( < Ts ), and is a paramagnetic (PM)
insulator above TN.

The past theoretical studies of these two materials can be broadly classiﬁed into
two categories. The ﬁrst one is density functional calculations executed on the crystals.
This sort of approach only achieved partial success and failed to predict certain properties
of the materials. For instance, local spin density approximation (LSDA) calculations
incorrectly predicted LazCuO4 to be a PM metal [37] and LaQNiO4 to be an AFM metal

[38]. Since 1990, several modiﬁcations have been made to LSDA to produce AFM

 

' This work has been published in “Y-S Su, T. A. Kaplan, S. D. Mahanti, and J. F. Harrison, Phys. Rev. B
59, 10 521 (1999)”.

36

insulating solutions, like SIC-LSDA, LSDA+U [39,40]. However, since the
modiﬁcations are somewhat ad hoc, whether the modiﬁed LSDA approaches still fall into
the ab initio regime or not is arguable.

The second category is cluster calculations. This category itself can be further
divided into two subsets: single-magnetic-ion cluster (SMIC) calculations and multi-
magnetic—ion cluster (MMIC) calculations. The SMIC calculations were aimed at
obtaining the neutron scattering form factors of the crystals [41] and hyperﬁne properties
[42]. A more detailed description of this approach will be given in Sec. 3.2 when we
come to the comparison of results of the SMIC calculations with the results of our present
work. Note that, since the cluster in the SMIC calculations contains only one magnetic
ion, the AF M ground state of the crystals must be assumed a priori and quantities such as
the Heisenberg exchange parameter J cannot be extracted from the calculations. On the
other hand, MMIC calculations [43-46] have given information about J.

The present calculations were done using the program CR YSTAL95 [47], which is
designed to do Hartree-Fock (HF) calculations on inﬁnite crystalline materials, using
linear-combination-of-atomic-orbital (LCAO) basis sets. Clearly, the crystal Hartree-
Fock approximation (HFA) falls outside the two categories mentioned above and thus
may provide useful supplementary information. The crystal HFA has recently been
applied to a number of AFM Mott insulators including MnO, NiO [48], KNiF3 [49],
KCuF3 [50], and CaCuO2 [51]. These were successful in that they gave the materials to
be insulators (vs. metals), to be antiferromagnets (vs. ferro- or paramagnets); also they
gave energy differences between the various magnetic states in semiquantitative

agreement with experiment. In this paper, in addition to these considerations, we will

37

compare the spin density with neutron diffraction experiments. Our major emphasis is on
the ground state properties, although some discussion of band structure (l-electron

energies, nature of l-electron states) is presented.

 

 

 

 

 

 

Figure 3.1: The centered-tetragonal structure of LazNiO4. LaQCuO4 has the isosturcture
of LazNiO4 (Ni atoms are replaced by Cu atoms). The NiO6 octahedral cluster, an
important signature of perovskite crystal structure, is marked by the dotted lines.

In this work, we use the centered tetragonal crystal structures with lattice
constants taken from experimental data [52]: a=3.855, c=12.652 for LazNiO4, and
a=3.7793, 0:13.226 for LazCuO4 (in units of A). The basis sets of 0, Ni and Cu are taken
from Ref. 53, which are specially designed for crystal calculations. (Note that the usual

basis sets used in molecular and cluster calculations are often too diffuse for crystal

38

calculations; they may cause nonconvergence or a very slow convergence rate with not
much increment in the accuracy to the result one might get using tighter basis sets [47].)
Since La is a heavy atom (atomic number 57), relativistic effects are not negligible for the
inner shell electrons. Besides, an atomic basis set for La suitable for crystal calculation is
not available in the literature. We thus used the effective core potential (ECP) calculated
by Hay and Wadt [54] for the La+3 ion core [55], where the relativistic effect of the inner
shell electrons is supposed to have been taken into account. In Hay and Wadt’s papers,
basis sets for valence electrons designed for use with the ECPs were also given.
However, they are so diffuse that they cause a numerical problem in the crystal HF
calculation. Thus, La was treated in the present work as a bare La’”3 ion and represented
by the Hay and Wadt large core ECP, no valence orbitals (therefore no valence electrons)
being attached to it. We did test the signiﬁcance of the bare La“3 core approximation to
some extent by adding a d shell (consisting of a single, Optimized, Gaussian exponent) to
the La+3 and seeing how it affected the results. It turned out that about 0.45 electron per
La gathered into the added shell. Accompanying that, the total energy per formula unit
(including one Ni or Cu) changed by about 4.5 eV. However, the shape of the occupied
bands and the spin density changed negligibly, and the energy difference between
ferromagnetic (FM) and AFM solutions, which is around 36 meV, varied by only about
1.4 meV. The insensitivity of the FM-AFM energy difference to variation of the outer
shells of basis sets was also observed for other materials [49]. The properties that are

only slightly affected by the bare La+3 core approximation are the main focus in this

paper.

39

3.2 Results

3.2.1 Heisenberg exchange parameter J

Our results show that both LazNiO4 and LaQCuO4 are AFM insulators, in
agreement with experiments. The FM solutions also exist for both materials, with higher
energy. The energy difference per formula unit between FM and AF M states is 36.9 meV
for LaQNiO4 and 36.1 meV for LaQCuO4. It is a common and long-standing practice
[56,57,48-50,43] to map energies of FM and AFM electronic states to the mean ﬁeld
approximation for the nearest-neighbor Heisenberg Hamiltonian Hum-s. This is done for

materials such as these, where H is known to accurately describe the low-lying

Heis
excitations [56,58], and leads to an estimate of the Heisenberg exchange parameter J. In
Appendix C, we give a rationale for doing the mapping in HF calculations, while pointing
out a potential inconsistency if the same procedure is used in approximate density

functional theories.

To proceed then with the determination of J, the Heisenberg Hamiltonian is

HM = 1231.3, , (3.1)
<1 1)

the summation going over each nearest neighbor pair (1, j) of magnetic ions, taken once.

We equate the calculated FM-AFM energy difference in the HFA to the corresponding
FM-AFM energy difference in the mean ﬁeld approximation of Eq. (3.1), and thus deﬁne

the Heisenberg J in the HFA, J HF, as follows

NZ]. (3.2)

AEHF = (2J,.,Fs"~(7

40

AEHF is the FM-AFM energy difference for the whole crystal in the HFA, N is the
number of magnetic ions, s is the spin of each site (1 for Ni and 1/2 for Cu), and Z is the
coordination number. Eq. (3.2) (with J instead of JHF) appears elsewhere [43-46] for
dimers where N = 2, Z = 1, but the claim in some of the references that this gives the
correct or exact J is not warranted in general (see Appendix C). Since the interactions
within a Cqu plane are much stronger than the interplanar ones, we take Z = 4. It turns
out that J is 9.2 meV for LaQNiO4 and 36.1 meV for LaQCuO4. The LazCuO4 result agrees
with an embedded cluster calculation using HFA [43], which gives 37.8 meV,
strengthening the suggestion [43-46] that J appears to be a rather local property. These
values are to be compared with the experimental data, which give 18 meV for LaQNiO4
[59,60], and 134 (neutron) or 128 (Raman) meV for LaQCuO4 [61,62]. Although our
results are about a factor of 3 too small, it actually is remarkable to be so close, since we
are picking up the small difference between the FM and AFM energies with each on the
order of 4000 Hartrees. (So it is at an accuracy of about one out of a million.) Similar

underestimates for I using HFA were found for other materials [48-50].

3.2.2 Spin density and spin form factor

Spin density maps on the perovskite magnetic ion — oxygen layer [(001) plane] are
given in Fig. 3.2 for both materials. The x-y plane of the graphs matches the NiO2
(Cqu) plane and includes 3 x 3 magnetic ion sites. The large peaks on the graphs
coincide with the locations of the magnetic ions. The sign of the spin density alternates

from one site to another, a manifestation of the AFM character of the systems. The

41

magnitudes of spin density around oxygens, which are half way between each nearest

neighbor pair of magnetic ions, are extremely small as expected. Every peak in the

graphs has an approximate 4-fold symmetric shape, a result of the singly occupied die—)9

orbital. Although Fig. 3.2(a) (for LaQNiO4) and Fig. 3.2(b) (for LazCuO4) resemble each
other very much, a characteristic difference between the spin densities of these two
materials does exist and can be seen in a tomograph of other directions. For instance,
Fig. 3.3 shows the spin density tomographs of these two materials in the (100) plane
passing through a magnetic ion. For this plane, the peak in the graph of LaQNiO4 still has
four lobes, while the peak in the graph of LaQCuO4 has only two lobes aligned in the y

direction. This is expected because each copper in LazCuO4 has only a singly occupied

dx,_y, , while each nickel in LaQNiO4 has a singly occupied dx2_y2 and a singly occupied

d 322-12 whrch together grve a spin d1str1butlon wrth cubrc symmetry (small d1stortrons

from cubic symmetry are caused by the surroundings). Later on we will see that this
difference in spin densities leads to characteristically different spin form factors for these
two materials.

It has been shown [41] that, for the single-band Hubbard model assuming small

hopping t compared to the repulsion U, the Fourier transform of the spin density is

42>:<s.>...F<r>1<211+0100. 63>

42

 

La2NiO4

 

(a)

La2Cu04

 

(b)

Figure 3.2: The spin density in the magnetic ion - oxygen plane [(001) plane]: (a)
LazNiO4, AFM, (b) LazCuO4, AFM. The range of the basal plane in the graphics is 20 (a
is the lattice constant) in each direction and is centered at one magnetic ion. The vertical
axis is in atomic units (same as Fig. 3.3 in the following).

43

La2NiO4

    

o
o .0
- .-
O . -0 .o
. .~:9:'.:.
0
. .. -
ﬁat..- .
Tran-:0... .
m4". *0 "

   

.
O-
..

 
      
  

 

  
   
  
 
    

S 2 0 . 5
l
o . 2 5 \
0
‘5 1'
La2CuO4
l
i
l
ﬁéﬁafaz - f
~ i:::.:.. '1
~~ ‘ﬁkiﬁm ’
0.7 51 .32: ”’5"?
1‘ Q
0.51
82 '
0.25

(b)

Figure 3.3: The spin density in the (100) plane: (a) LazNiO4, AFM, (b) LazCuO4, AFM.
The range of the basal plane in the graphics is a in each direction and is centered on one
magnetic ion. The labels of the two axes of the basal plane refer to the y and 2 directions

in the conventional tetragonal lattice coordinates.

44

 

Here (52),“ is the expectation value of the spin at a site (an up site) in the
(spontaneously) symmetry-broken ground state of the corresponding Heisenberg model.
f (I; ), the form factor, is the Fourier transform of Us times the HF spin density over a
nonmagnetic unit cell centered at a magnetic ion, where s is the spin quantum number for

the free ion (1 for Ni+2 and 1/2 for Cu”). This implies f (0) = 1 . F (l; ) is given as

F02): get-W (3.4)

R

where R runs over all the Bravais lattice vectors and k A is an antiferromagnetic wave

_.

vector (e.g. for the square lattice with lattice constant a, k A is l i]). The correction

a’a

term 0((%,)2) is negligible for the cuprate [41] and, we expect, for the nickelate [63].

Although the derivation of Eq. (3.3) is somewhat lengthy, the result is physically

understandable. (sz>Heis accounts for the quantum spin ﬂuctuations (QSF) in the

Heisenberg-model ground state and is just an overall scale factor. F (I; ) accounts for the

occurrence of the Bragg peaks. f (If ) gives the relative sizes of the Bragg peaks, and

conveys detailed information about the spin density s(F) in the unit cell. Eq. (3.3) is
discussed further in Sec. 3.3.

In an earlier work [41] of our group, SMIC calculations on several AFM
insulating compounds including LaQNiO4 and LaQCuO4 were carried out based on this
formula. The results agreed extremely well for KNiF3 and NiO; for LaQCuO4 they gave
rough overall agreement for the shape of f (it) (for a sample well off stoichiometry).

Only in the case of LaZNiO4 did the results show a serious qualitative disagreement with

45

experiment [64]. This disagreement provided a strong part of our motivation for the
present study.

Thus we calculated the Fourier transform of the spin density that we had obtained
from the crystal HF calculations. The resulting form factors at the k-values for the
magnetic Bragg peaks are shown in Fig. 3.4.

Note that if the spin density is spherically symmetric, all the data will fall on a
smooth curve. In Fig. 3.4(b), the form factor of LazCuO4, the data separate into several

branches, indicating that the form factor varies dramatically in different directions of I; .

This is understood as due to the asphericity of the Cu dx2_y2 orbital of the unpaired

electron, as ﬁrst pointed out by S. Shamoto et al. [65]. On the other hand, in the Ni
compound both eg states are singly occupied, giving a nearly cubic spin density; and thus
s(F) is closer to spherical. Indeed, the corresponding form factor [Fig 3.4(a)] reveals no

obvious branches.

46

La2NiO4

1 lemﬁmmmmlllvm__m_l_”m__~_mm_mm.-.1l.1ﬂ.lq,pl

 

 

 

 

 

0.8 * .
‘ o
l * .~. 1
H ‘ .‘ i
o -t
“'k
80 6 m
g M
. I.
At
E i .i
(40.4} I
8 . n
I
0.2L
‘ 1
1mm -1-” _L ........ 1d
1 2 3 4 5 6
k
(a)
La2CuO4
1 V V a
0.8;
‘ *9 t o t i
O 1
i
i o ]
14 t
80 6 O
F o
*1
. .‘IA.
EO.4' ‘ . _-
o ‘ I
Lu ‘ A
0.2- i
3' i
l 2 3 4 5 6
k
(b)

Figure 3.4: The spin form factor: (a) LaQNiO4, AF M, (b) LaQCuO4, AFM. The
horizontal axis is the magnitude of it , in units of inverse Angstroms. The black stars are
for the family of (i —§,l) Bragg peaks, diamonds for (%,%,l), circles for (%,%,l), empty

2’

stars for (%,—%,l), triangles for (%,—%,l) and squares for (%,%,l), where the three

components of the triple numbers are in units of 27” , 27” , 27" , respectively. With our

choice of coordinate system, for the Bragg peaks the sum of the three components of each
triple number must be an even integer.

47

Our present results are compared to the previous cluster calculations and to

experiment in Fig. 3.5. Note that Fig. 3.5(a) for the nickelate shows the Bragg scattering

amplitude (in units of [33) and Fig. 3.5(b) for the cuprate shows the form factor. In Fig.

3.5(a), the Bragg scattering amplitude g(sz ) Hm f (I; l is calculated using g = 2.29 [66] and

(S, > Has = 0.8 from the spin wave theory [67]. The two theoretical results agree with each

other very well; however, both disagree qualitatively with experiment. In View of the
very good agreement between theory and experimental results for several other
compounds like KNiF3 and NiO [41], it seems odd that the theory would be so far off
from experiments for the case of LaZNiO4. The agreement between the cluster and the
crystal calculations ensures that the theoretical results are calculation-error free,
excluding one possible source of the discrepancy between theory and experiment. Then,
logically, more independent experiments are called for to cross-check with the single
existing experimental work on a presumably stoichiometric sample [64,68]. We should
note that the plot in Fig. 3.5(a) is on an absolute scale; thus the agreement between theory

and experiment at large k might be signiﬁcant [69].

48

La2NiO4

 

 

 

 

 

 

 

 

2! ’ r“ 1
1.75 i
1
. l
0
1.5 t l
t
O t
O
*

131.25 Q i
3 1 1. 1‘
0
.1
a 1 ll l ‘i in i l

9 ‘1
e 1 *. *1
‘H 0 ‘ J
10 75 l I; T]
‘ r:
0.5 ' 911*
i
0.25 i
_ . . l
l 2 3 _14 5 6
4nsine/A (A )
(a)
La2CuO4
l[ ll_l_ "T?
1 1
l .
1 ﬂ 1
0.8] g * 4 .
l 05 o 9 ]
| I ‘
0.6 l
E 31' ‘
W 99
0.4 ]
l
l 1
0.21 1
i .3
1 i
1 . __-_-_ . l
1 2 3 _14 5 6
4nsin6/A (A )
(b)

Figure 3.5: Comparison of the crystal and cluster calculations with experiments: (a)
Bragg scattering amplitude for LaZNiO4. (b) spin form factor for La2CuO4. The
diamonds are for the crystal, stars for the cluster, and circles with error bars for the
experimental data. Note that in (b) the experimental data has been scaled to make the
least-square ﬁt to the crystal results.

49

For the cuprate [Fig 3.5(b)], the experimental form factor data [70] have been
scaled to give the least-square ﬁt to the crystal HF result. It has been argued [71, 41] that
the moment (and Néel temperature) should be maximum at the stoichiometric compound;
the main experimental question here is the oxygen content [72,73]. Since the
experimental “moment” in Ref. 70 is appreciably smaller than that of other samples [72],
we expect the one used to measure the form factor shown is far from being
stoichiometric. Therefore, only comparison of the “shape” of the form factor might be

sensible. It is seen that the overall shape is in rough agreement with theory.

For a test of the absolute value of our calculated spin density, we consider the
sample with the highest “moment”, 0.6 i 0.05 Bohr magnetons (with Néel temperature of
298K) [72]. As discussed in Ref. 41, the “moment” (called moment in that reference)

deﬁned by the experimentalists is the amplitude g<s(lt)> at the Bragg peak with the

smallest k, divided by 0.835 (the form factor of KzCuF4, interpolated for the k—value

appropriate for the cuprate being measured [72]); g is the usual g-factor, E 2.2 for the

cuprate. Multiplying the above “moment” by 0.835 yields g<s(/;)> = 0.50 i 0.04.

exp

Using our calculated HF value of f (I; ) = 0.763, (sJHm = 0.3 for the square lattice
(appropriate to a Cqu plane) [67], and g = 2.2, Eq. (3.3) gives g<s(R)> theo = 0.504. This

exact agreement is surely fortuitous. In fact a sample with a slightly higher TN (325K)
has been reported [74] as being highly stoichiometric. Unfortunately an absolute
magnetic Bragg intensity measurement (which would yield 6. g. a value of the “moment”)

was not reported on this sample. Judging from the increase in “moment” with TN shown

50

in Ref. 72, the experimental amplitude will be larger than the 0.5 quoted above.
Nevertheless, this agreement between theory and experiment suggests that the HF
approach, corrected for quantum spin ﬂuctuations, gives a very good account of the
ground state spin density. It is clearly important to have the absolute intensities at several
Bragg peaks measured on an excellent sample like that of Ref. 74, for comparison with
our calculations.

It is noted that an appreciable discrepancy between the crystal and cluster HF
results exists. This is interesting in view of the very good agreement between the two
theoretical calculations for the nickelate as shown in Fig. 3.5(a). We believe that the
larger covalence in the cuprate is the cause of this discrepancy. Indeed, we performed the
integration of the HF spin density over a Wigner-Seitz cell on the magnetic ion sublattice,
as a measure of the ordered spin per magnetic ion (without QSF), on both materials. We
found 0.926 for the nickelate and 0.421 for the cuprate. In comparison with the free
magnetic ion spin (1 for Ni+2 and '/2 for Cu”) we see a 7.36% reduction for the nickelate
and 15.8% reduction for the cuprate [75]. This reduction is evidence of the covalence
between Ni (Cu) and Oxygen [76,77]. The ﬁnding here of appreciably less covalence in
the nickelate strongly suggests that the original explanation [64] of the extreme ﬂattening
of the form factor at small k as being due to strong covalence may not be correct [78].

As one further check on the accuracy of our crystal calculations, we considered

another Ni compound, KNiF3, where the cluster approach agreed excellently with

experiment [41,66]. Taking (s2) = 0.92 from the spin wave theory [67], we

Heis

summarize the results for the amplitude <s(lii)> at the three Bragg peaks measured in

51

Table 3.1. It is apparent that both the crystal and cluster results are slightly larger than

the experimental data. However, if we took (s2) .3 = 0.9, both of them predict the

Her

experimental data to within the error bars! It should be remembered that the spin wave

theory of Ref. 67 is an approximation to the ground state of the Heisenberg model [79].

Table 3.1: Comparison of the crystal HF, the cluster HF and the experiment results of
the spin form factor of KNiF3. The ﬁrst row shows the indices of the three Bragg peaks
measured in experiments up to date. For the convention of the indexing, see the caption
of Fig. 3.4.

 

A
NI—
N|--
NI—
v

 

 

 

Crystal HF 0.813 0.679 0.554
Cluster HF 0.807 0.683 0.564
Experiment 0.783 a 0.018 0.672 1 0.015 0.553 i 0.013

 

 

 

 

 

 

3.2.3 Band structure

The band structures of the AF M solutions are given in Fig. 3.6. The results show
a wide gap (~17 eV) for both materials; we also see a similar large gap for the FM
solutions. The experimental (optical) gap for La2CuO4 is about 2.0 eV [80]. That the gap
is much larger than the experimental gap has been found much more generally in the
HFA; see the papers of Refs. 48-50 for examples and Ref. 48 for discussion. The results

suggest that the insulating character of both materials is not strongly correlated with their

52

 

 

 

magnetic ordering, in agreement with experiment [81], and in contradiction to a band-
theory scheme proposed by Slater [51,82]. In the LaQNiO4 AF M solution, the F point of

the highest occupied band consists mainly of the Ni 3dx2—y2 and 0(1) 2 pm , and the F

point of the lowest unoccupied band consists mainly of the O(1) 3s and 0(2) 3s, plus a bit

of Ni 4s and 0(2) 3pz. In the LaQCuO4 AFM solution, the F point of the highest

occupied band consists mainly of the Cu 3dx2—yz and 0(1) 2 pm , and the F point of the

lowest unoccupied band consists mainly of the 0(1) 3s and 0(2) 3s, plus a bit of Cu 4s.
For both materials, the F point is the minimum in the lowest unoccupied band but not the
maximum in the highest occupied band, showing an indirect band gap. The width of the
highest occupied band is about 1.6 eV for LaQNiO4 and 0.6 eV for LaQCuO4. The
projected densities of occupied states (not shown) have many features in common with
those previously reported for other Mott insulators [48-50]. Two of these features are
that it is mainly anion p-states that exist near the Fermi energy, and that the bulk of the
magnetic ion eg states is lower in energy than the t2g states. These results need further

analysis and are suitable for a future publication.

53

 

 

0.4:
1 l , i l 1
’3? 0.2 l l l 1
. 1 1 j
(0 1
132 1
U 01

 

 

ll._.1.___ l

(a)

 

 

 

 

Figure 3.6: The band structure: (a) LazNiO4, AFM, (b) LaQCuO4, AFM (on next page).
The band structures of the spin-up and spin-down states are the same (because of the

antiferromagnetism), so the ﬁgure shown here is for either one. The (vertical) .13-axis is in
units of Hartrees. The (horizontal) k-axis goes through two contiguous paths: (%%0) —->
(000) —> (3,40) —> (440) —> (000) a (001) and (44,4) —) (000) —+ (404). For the
convention of the triple numbers, see the caption of Fig. 3.4. The horizontal line inside
each graph indicates the Fermi energy.

54

 

 

 

Figure 3.6 (continued from previous page)

55

3.3 Summary and Discussion

We have calculated the electronic structure of LaQNiO4 and LaQCuO4 using the
crystal HF approach. The major emphasis has been on ground state properties, although
some discussion of band structure has been given. The results provide supplementary
information to that obtained from crystal LSDA and embedded cluster calculations in
interpreting the experimental data. For both materials, our crystal HF results correctly
predict an AF M insulating ground state. The calculated Heisenberg exchange parameter J
is about a factor of 3 too small compared with experiments. This is consistent with the

interpretation of Martin and Illas [43], and Towler et al. [48] that the HFA overestimates

the on-site Coulomb interaction. The spin form factor f (I; ) basically agrees with earlier
cluster calculations. However, there is an appreciable discrepancy between the two
theoretical calculations in the cuprate case, which we believe is due to the larger

covalence in the cuprate. We found rough overall agreement with experiment for the

shape of f (It. ) vs. l]: | for a poor sample of the cuprate, and excellent agreement for the
absolute intensity of the one Bragg reﬂection measured on a good cuprate sample.
However in the case of LaQNiO4 the shape of the form factors in both the crystal and
cluster calculations disagree seriously with experiment. The fact that we found the
nickelate to be (appreciably) less covalent than the cuprate deepens the puzzling nature of
this disagreement with the nickelate together with agreement for the cuprate. For a
further check, we also calculated the form factor for KNiF3 and found excellent
agreement with experimental absolute intensities. These results should create strong

motivation for performing more experiments on a stoichiometric sample of the nickelate;

56

it is also important to measure the absolute Bragg intensities at more than one Bragg peak
for a stoichiometric sample of LaQCuO4 for a more stringent test of our theory.

We conclude with a discussion of the procedure used to calculate the spin density.
We used our earlier result calculated within the single-band Hubbard model as a guide.
Namely, following that result, we took the spin density to be that in the HFA reduced by
the ratio <sz>Heis / s, the reduction from 1 being due to the quantum spin ﬂuctuations in
the Heisenberg AF M ground state. In fact we have shown [83] that this procedure is not
exact for a model more general than the single-band Hubbard model. Thus we need to
understand, e.g., why our present and previous [41] results are in excellent agreement
with experiment for a number of cases using the HFA. Such studies are in progress. In
any case, we want to emphasize that the QSF (which give large effects) cannot be
ignored, as they were in Ref. 46 where a large moment reduction was attributed entirely
to covalence. We note also that in Ref. 43, the large QSF correction was made, showing

the confusion on this issue in the literature. The single-band Hubbard model result, Eq.

(3.3), is a good starting point for clariﬁcation, and the generalization is clearly important.

57

Chapter 4

Electronic Structure of LaMnO3 in the Ab
Initio Crystal Hartree-Fock
Approximation*

4.1 Introduction

The divalent-metal doped rare-earth-manganites have received a great deal of
attention in recent years largely due to the colossal magnetoresistance (CMR) in these
compounds. The phase diagrams of these manganites are very rich; the magnetic and
conducting properties, as well as crystallographic structure, can vary substantially with
temperature and doping concentration [84]. These properties alone, aside from the CMR,
constitute an interesting research area. LamCaanO3 is a typical family of these
manganites, with the CMR occurring in the region near x = 1/3 [84]. Here we focus on
LaMnO3 (the x=0 end member of the above family). Its properties, which involve
magnetic and orbital orderings plus strong Jahn-Teller distortion, are interesting in their
own right.

There has been a great deal of work on LaMnO3; the papers most directly relevant
here are Refs. 85 — 90. Saitoh et al. [85] obtained experimental photoemission results

and interpreted them through a cluster conﬁguration-interaction model. LSDA (local spin

 

* This work has been published in “Y-S Su, T. A. Kaplan, S. D. Mahanti, and J. F. Harrison, Phys. Rev. B
61, 1324 (2000)”.

58

density approximation) band calculations were reported in Refs. 86-88. In all three of the
band calculations, the experimentally observed ground state magnetic ordering was
found. The observed orbital ordering was obtained by Satpathy et al. [87] (this property
was not discussed in Refs. 86 and 88). The spin Hamiltonian, which governs the
magnetic properties including the low-lying excitations (spin waves or magnons), was
calculated by Solovyev et al. [89] within the LSDA. Hirota et al. [90] determined the
magnon dispersion via inelastic neutron scattering measurements, and claimed it to be
consistent with the theory of Ref. 89. Sarma et al. [86] found their LSDA density of
states and calculated photoemission intensity to agree with the measurements of Saitoh et
al. [85].

These facts would seem to have the theory of these properties of LaMnO3 in
satisfactory shape. However, the theories of Refs. 85 and 86, both of which apparently
explain the photoemission data, disagree with each other. In Ref. 86 (in agreement with
the other LSDA calculations of Refs. 87 and 88), the Mn d-band lies near the top of the
valence band, with the O p-band lower than the d-band, while in Ref. 85 the opposite
order occurs. Furthermore, an LSDA+U band calculation gave the Mn d-band below the
O p-band [87]. Thus there are major differences between existing pictures, and they are
of considerable importance, e.g. at stake is the nature of the band gap, O-p —> Mn-d or
Mn-d —> Mn-d (the two possibilities have been called charge-transfer insulators or Mott—
Hubbard insulators, respectively [91]). We were thereby motivated to carry out
calculations on LaMnO3 using the Hartree-Fock approximation (HFA). HFA is a well-

established theoretical approach, independent of the above methods, and thus can provide

59

valuable comparison with the other results. Also, it is known that for some 3d-transition-
metal oxides with perovskite-based structure, e. g. lanthanum cuprate and nickelate,
LSDA failed [92] and the HFA succeeded [93] to predict correctly the ground state
insulating property and the magnetic ordering.

Our HFA results show some surprising physical effects, and signiﬁcant
differences from LSDA. The correct magnetic and orbital orderings and insulating
character are found in both theories. However, we ﬁnd major differences in the occupied
densities of states (e. g. the O p bands lie close to the top of the valence band in the HFA,
more like the interpretation of Ref. 85 and the LSDA+U results of Ref. 87). There is also
a major difference in the effective spin Hamiltonians; yet the magnon dispersion curves
of both theories are consistent with the neutron scattering experiment [90], as we will
explain. Also, in apparent contrast to LSDA [87,94], a large spinless charge backﬂow,
0'2 2p ——> Mn+3 3d, is found in HFA. Finally, the fundamental type of insulator differs in
the two approximations: LSDA gives a band insulator (the gap doesn't exist for the cubic
structure), HFA yields a Mott insulator (the gap exists for both the cubic and distorted

structures).

4.2 Method

If there were no distortion, LaMnO3 would be a perfect cubic perovskite crystal as
shown in Fig. 4.1. Note that the MnO6 octahedron marked in the ﬁgure is connected with

every one of its neighboring octahedra by sharing a common oxygen, which is a key

60

signature of the perovskite structure. In reality, the MnO6 octahedra (denoted as [Mn06])
are strongly distorted due to Jahn-Teller effect and rotate from the crystal axes by an
appreciable amount, resulting in an orthorhombic crystal structure (space group ana)
with four Mn per unit cell [95]. The unit cell of the orthorhombic structure is close to a
J2 x J2 x 2 supercell of the undistorted cubic perovskite structure. The longest side of
the unit cell is chosen as c-axis [96] and the MnO2 planes perpendicular to it are called

basal planes.

 

 

 

 

Figure 4.1: The cubic perovskite crystal structure. This is the structure for compounds
like KNiF3 and CaMnO3. For LaMnO3, large distortion exists. The atom at the center of
the cube is Mn, the face centers are 0, and the corners are La. The important MnO6
octahedral cluster is marked by the dotted lines.

61

We carried out the calculations on both the observed orthorhombic structure and
the ﬁctitious iso-volume cubic structure in order to investigate the effect of the strong
lattice distortion. To our knowledge, ours is the ﬁrst ab initio HFA calculation for this
material [97]. The calculation makes use of the program CRYSTAL95 [98]. The basis
sets of Mn and O are those optimized for CaMnO3 [99]. La is treated as a bare La+3 ion
and represented by an effective core potential [100]; a test of this approximation as well

as other accuracy control parameters of the program will be discussed in Sec. 4.4.

4.3 Results

4.3.1 Magnetic properties

For the orthorhombic structure, it can be shown that there are ﬁve collinear
magnetic orderings that maintain the size of the unit cell. They are ferromagnetic (FM),
A-type, G-type, and C-type antiferromagnetic (AAF, GAF, and CAF), and ferrimagnetic
(FI) orderings, deﬁned as follows. AAF: the Mn spins are parallel in a basal plane and
antiparallel from plane to plane. GAF: each nearest neighbor (n.n.) pair of Mn are
antiparallel. CAF : each n.n. pair of Mn are antiparallel in a basal plane and parallel along
the c-axis. F I: one of the four Mn in a unit cell is antiparallel to the other three. HFA
results for these ordered states are listed in Table 4.1, together with LSDA results [88] for
comparison. It is seen that HFA predicts the ground state of LaMnO3 to be an AAF
insulator, in agreement with experiment [101]. LSDA also makes the same prediction.

In the cubic structure, both theories predict the FM state to have the lowest energy.

62

However, there are substantial differences between the theories. E.g., in HFA all states
are insulating (for both cubic and orthorhombic structures), while in LSDA all states of
the cubic structure are metallic [88]. The band gaps (not shown in the table) in HF A are
much larger than those in LSDA. From the results in the AAF column, it is seen that the
crystal distortion lowers the energy per Mn by ~1 eV in HF A, vs. 0.27 eV in LSDA. For
the cubic case, it is seen that LSDA predicts much larger energy differences among the

various magnetically ordered states.

Table 4.1: HFA and LSDA [88] energies of LaMnO3 with various magnetic orderings.
The energies shown are relative to the FM state of the cubic structure, in meV / Mn.
(ins.=insulator, met.=metal.) The HFA energies are for the states with the observed
orbital ordering (see text).

 

 

 

 

 

FM AAF GAF CAF FI
cubic 0, ins. 0.4, ins. 34, ins. 32, ins. 16, ins.
”FA orth. -1053, ins. -1055, ins. —1041, ins. 4039, ins. -1047, ins.
cubic 0, met. 110, met. 365, met. -- --
LSDA orth. -- -156, ins. -- —- —-

 

 

 

 

 

 

 

 

 

The HFA energies of the ﬁve magnetically ordered states for the orthorhombic
structure are used to map to an effective spin Hamiltonian (see Refs. 93 and 99 for
discussion of the mapping). We add 4-spin terms to the usual 2-spin terms in the spin

Hamiltonian [102]

H =constant+ZJij§i S]. + ZJiik,(4—spin tenns)+---; (4.1)

(131‘) (i,j,k,1)

63

where S, is the spin at (Mn) site i, 1,,- and J,,-,,, are exchange parameters, and each

combination of summation indices is summed once. Since there are ﬁve energies, the
mapping can determine four J’s plus the constant in H. We choose the J’s in the
following way. We consider the one-band Hubbard model with n.n. intra- and inter-
(basal) plane hopping t,, t2. Perturbation theory with the ti small implies the spin
Hamiltonian of Eq. (4.1). Keeping through fourth order terms, more than four J’s occur;
however, only a particular subset of four can be determined by the ﬁve ordered states
considered. This is the chosen set and is as follows: the intra-plane n.n. J ,, the inter-plane
n.n. 12, the inter-plane next n.n. J3, and the inter-plane 4-spin J4. Due to the distortion of
the crystal structure, there are actually two types of J3 and J 4, denoted as J 3“) and J 3(2) [89],
and J 4“) and J 4(2). In the mapping here, only the average values of the two types of J’s are
determined: J3 is the average of J3“) and J3“), and similarly for J4. Note that the same
mapping can also be done for the cubic case. This is possible because the existing orbital
ordering (discussed later on) breaks the cubic symmetry and make, for example, J, not
equal to J 2.

The results of the mapping are listed in Table 4.2. In both the cubic and
orthorhombic cases, J 3 and J 4 are negligible. This supports the neglect of further neighbor
terms, expected a priori to be small. The signs of J, and J2 reﬂect the ground state
magnetic properties (FM for the cubic and AAF for the orthorhombic). LSDA
calculations by Solovyev et al. for the orthorhombic case give a qualitatively different

picture [89]. As shown in Table 4.3, LSDA gives a much bigger J3; also J, and J2 are both

64

ferromagnetic and they alone would yield a wrong magnetic state. It is the J 3 that turns

the ground state from FM to AAF in the LSDA.

Table 4.2: HF A J’s in the effective spin Hamiltonian of LaMnO3, in meV.

 

 

 

J, J2 J, J,
cubic -2.1 -013 0.019 0.0049
orth. -0.88 0.21 0.0036 0.00051

 

 

 

 

 

 

 

The spin wave spectrum of the system has been measured in a recent inelastic
neutron scattering experiment and was ﬁtted very well using a 2-J (Jl and J2) spin
Hamiltonian by Hirota et al. [90]. The results are also included in Table 4.3 for
comparison. Our HFA results are smaller in magnitude than the experimental values by
about a factor of 2 for J, and 6 for J 2. At ﬁrst glance one would say that the experiment
shows a 2-J character of the system and thus favors the HFA picture. However, as we
analyze further, we ﬁnd that this conclusion is not solid enough yet. Hirota et al. state
that their experimental results are consistent with those of Solovyev et al. in the following

sense. If one maps the 3—J model by Solovyev et al. to the following effective 2-J model,

{me : JIM (4 2)

(21) _ (U) (U)

(the factor of 4 in the second equality comes from the coordination number involved with
J 3), then J ,(2’) and J2” have the right signs and are both about a factor of 2 too large in

magnitude compared with the experimental values. By Eq. (4.2), the 3—J model and its

65

effective 2-J model have identical spin wave dispersion along the c-axis. Dispersion
along several other directions is shown in Fig. 4.2. It is seen that the spin wave spectra of
the two models are quite close. The difference is estimated to be about the size of the
error bars of the experimental data cited. So at this point, the LSDA picture is not
deﬁnitely ruled out. The issue could be settled by somewhat more accurate

measurements.

Table 4.3: Comparison of the HFA, LSDA [89] and experimental [90] J’s in the spin
Hamiltonian of LaMnO3, in meV.

 

 

 

 

 

 

 

 

J, J2 J 3
Exp. -l.67 1.21 ~ 0
HFA -0.88 0.21 ~ 0
LSDA -2.28 -0.78 0.78

 

 

66

 

 

 

.//’ ' / /l a
t ”l ’ ' , ’
:1 l /// / Iii/ﬂ
° 10 f— ~ ,.
j / ’ ﬁ/i
, //g [1 10] j / , [1 12]
: /,//// //
O / ,1
, T ,, -
20 2 ,
3‘ / /
m , /
é , , / ,
c: 10 ” / I
= / A
g— / [114] 1 / [116]
/ ‘ _/'
0 1/ V

 

Figure 4.2: The spin wave dispersions of the 3-J (solid line) and its effective 2-J
(dashed) spin Hamiltonians, along four arbitrarily chosen directions. The values of J’s

are those given by Solovyev et al. [89]. The k-indexing follows the convention used in
Ref. 90.

67

4.3.2 Spinless Charge Backﬂow and Density of States

To our surprise, Mulliken population analysis (MPA) on the HFA results gives
charges that deviate substantially from the formal valence picture, although the spin
values do not show such deviation. We ﬁnd Mn+2'2, O,"'7 and O,,"'75, while the Mn spin is
1.98; the formal valence picture is Mn”, 02, with Mn spin of 2. (O, is the apical oxygen
and O,,, the basal-plane oxygen.) We note that, due to the bare core approximation for
La, these MPA results should only serve as an indication of deviation from the formal
valence picture, not taken as being very accurate. A similar but more severe departure
from the formal valence picture has been reported for an all-electron HFA calculation on
CaMnO3 [99], where Mn was found to be Mnﬂl7 with spin 3.25/2, compared with the
nominal value Mn+4 with spin 3/2. That is, there is a large nearly spinless backﬂow of
electrons from O'2 to Mni3 or Mn”. Such large changes in ionic charges would have
obviously important consequences, e.g., on Madelung energies, phonon spectra, dielectric
constant. How they would affect current simpliﬁed models is obscure, but certainly their
implications for such models would have to be considered.

It must be pointed out that the MPA is basis set dependent, and can be very
misleading when in the basis set there are diffuse basis functions that are appreciably
occupied [103,104]. However the following considerations suggest that the deviation
from the formal valence picture found here cannot be explained solely by basis set
dependence. (i) We repeated the calculation with the most diffuse Mn (1 and 0 sp basis
functions omitted from the basis set. This gave a very small change in the MPA, the

result actually being further from the formal valence picture. (ii) The way that the MPA

68

attributes an overlap charge to the two atoms involved is somewhat artiﬁcial and
therefore results in an inherent uncertainty in the meaning of the MPA results. Moreover,
when the MPA gives ridiculous results due to the presence of very diffuse basis ﬁmctions
in the basis set [103,104], there usually are large overlap charges. So smallness of the
overlap charges is an indication of the reliability of the MPA results. In the present case,
the overlap charges totally account for 0.06 electrons for Mn, considerably smaller than
the deviation of the MPA charge of Mn ﬁom the formal charge. (iii) In contrast to MPA,
the charge from the actual integration of the charge density over a reasonable volume
(vide infra) around an atom is much more basis-set insensitive, and gives a realistic value
of the HFA result (assuming a good basis set, of course). So this integration is a good
check for verifying the correctness of MPA results. We didn't do the integration for
LaMnO3 since the precise numbers will not be useful due to the bare core approximation
for the La, as mentioned in the previous paragraph. Instead, we redid the all-electron
calculation of Ref. 99 for CaMnO3, and then integrated the charge density over a cube
around a Mn. The faces of the cube are perpendicular to the Mn-O bonds and pass
through the minima of the charge density along the bonds, which gives the cube edge to
be 1.8 A (the cube is close in size to the Mn sphere of diameter 2.12 A used in Ref. 88).
The integrated charge for Mn was +2.52, close to the MPA result of +2.17. Further
discussion of the La bare core approximation is given in Sec. 4.4, where it is suggested
that a large deviation from formal valence, as found above, probably will hold true in
more accurate HFA calculations.

The projected density of occupied states of LaMnO3 is shown in Fig. 4.3. It is

seen that a small amount of Mn-d projected density of states (DOS), both up and down,

69

exists in the range of 0 to 6.9 eV below the top of the valence bands, E,, coinciding with
the range of the O p bands and accounting for the nearly spinless backﬂow. The large
peaks (spin up) in the Mn-d projected DOS in the 8.7-10.3 eV range below Et are
associated with the Mn (1 bands, which are spin polarized and give the moment of Mn
close to that predicted by the formal valence picture.

This HFA picture is rather different from previous LSDA calculations. First, to
our knowledge, no departure from the formal valence picture has been reported in LSDA
studies of the system. In fact, Satpathy et al. [87] state that the charge states are close to
nominal in their LSDA results. Second, the O p bands lie above the Mn d bands in HF A,
opposite to what LSDA predicts [86-88,94]. The disagreement between HFA and LSDA
in the order of O p and transition metal (1 bands has also been seen for other systems.
[92,93,88,99]

We also note that LSDA+U calculations [87] show ordering of the O p and Mn (1
bands much like our HFA results. Those calculations were disparaged [88] on the
grounds that the Mn spin turned out to be larger than the nominal value (3/2 for
CaMnO3). However, as far as we are aware, there is no theorem that the Mn spin must
not be greater than nominal. In fact intra-atomic exchange would tend to increase the Mn
spin by polarizing electron transfer from 0; also, we ﬁnd in HFA an increase in Mn spin
(as determined by integration of the spin density within a suitable cube centered at the

Mn) over the nominal value for CaMnO3.

70

 

 

 

 

 

 

 

6 . J ‘11 ' z *1

4 L ] ,i i

i" t] H‘ 1] Mn d ,

I 15 /J "111 1

0 L 1v -W A_,_\1 j

W ' 7

. 1

-2 g a , ; 1

i l

1 * On p , 11 i -

3“ V AH \V’quk‘x 1 J

t f 1

-1 r \fWML/A 1 41
—2 1 9 1 1 f . e ,

1 [' 0| p “

 

 

 

 

 

 

 

 

-2 it + 1 1 1 l l t 1
40 ~ I r
30 — 1 -
~ Total ' 1
20 ~ I -
10 F W J
r- I d
0 - : J
_10 1 1 1 1 1 1 1 1 . 1 1 1 1 . . 1 1 1 1 1 1 n 1 1 1
—12 -1O —8 —6 -4 _ —2 0

Energy (eV)

Figrue 4.3: The (projected) DOS of LaMnO3, with AAF ordering. Positive and negative
DOS are for up- and down-spin states, respectively. Energies are relative to the top of the
valence band. The projected Mn-d and O,,-p DOS are for Mn and O on an up-spin basal
plane. The projected O,-p and total DOS are symmetric for up- and down-spin. The
down-spin part of the total DOS is omitted from the ﬁgure.

71

The HFA predictions of the charges on ions and the characters of the valence
bands are consistent with the interpretation of a recent experiment by Saitoh et al., who
studied La,_,Sr,,MnO3 by photoemission and x-ray-absorption spectroscopy [85]. These
authors determined that the character of the band gap of LaMnO3 is of the p-to-d charge-
transfer type while that of SrMnO3 (corresponding to CaMnO3 in our discussion) has
considerable p-p character as well as p-d character. The HP A results (via the MPA) show
that for LaMnO3 and CaMnO3, the Mn-d electron population is roughly equal, while the
O-p population gets reduced for the latter. This implies that the valence bands of
CaMnO3 consist of a smaller amount of O-p character than LaMnO3, which in turn
suggests that there is a larger amount of O—p character in the conduction bands of
CaMnO3. Combining this with the other HFA prediction that the O p bands are the
highest occupied bands can explain the above experimental observation. Saitoh et al.
also reported that the Mn-d electron population is 4.5 for LaMnO, and 3.8 for SrMnO3, a
considerable deviation from the formal valence picture similar to (although not as severe
as) that predicted by the HFA, which we ﬁnd to be about 4.7 for both LaMnO3 and
CaMnO3. The most surprising aspect of the Saitoh et al. work [85] is that their
calculation showed a large reduction in photoemission intensity of the Mn d band as
compared to the density of states (due to matrix-element effects). We intend to calculate

that intensity using our HF A wave functions to check this.

72

4.3.3 Orbital Ordering

We now discuss orbital ordering — this is ordering of the single eg orbital occupied
at each Mn. In the observed structure, each [MnOé] is stretched substantially along one
axis and is rotated by 10-15° from its orientation in the cubic structure. Disregarding the
rotation for simplicity, the stretched axes lie in the basal plane, alternating in direction by
90°, and this pattern repeats along the c-axis. The stretching, being driven by the J ahn-
Teller splitting of the eg states, leads to an expected orbital ordering — the single occupied
eg orbital at each Mn is d(3zZ-r2)-type with its axial symmetry along the stretched axis of
the associated [MnO6]. This orbital ordering is indeed found in our HFA solutions, by
plotting the spin density. A combination of the three t2g and one eg electrons (nominally
the only unpaired electrons in the system) dictates a unique shape in the spin density
distribution, which therefore can be used to identify which eg orbital is actually occupied.
For example, Fig. 4.4 shows the spin density at a Mn with the d(3x2—r2) orbital occupied
(the z-direction is along the c-axis). This orbital ordering, which was predicted long ago
by Goodenough [105], has also been obtained in LSDA calculations [87], and was

recently conﬁrmed in experiments [106].

73

 

Spin Density at d(3x2-r2)-Mn

   
 
   
      
 
    

  
   

 

 
   
 

1.1
WWI“,
Y
X
l
s (r)
S (If)

 

Figure 4.4: The spin density on the xy-, yz— and zx-planes at a Mn with the d(3x2-r2) and
three t2g orbitals occupied. The plotting region is about (3.9 A)2, centered at the Mn. The
slightly off of the symmetry axes of the distribution from the x—, y- and z-axes is a result
of the rotations of the [Mn06] in the observed crystal structure.

74

In our study, we obtain the same orbital ordering for all trials with various
magnetic orderings as well as initial conditions used to start the Hartree-Fock
calculations, for the observed orthorhombic structure. However, for the cubic structure
we ﬁnd a variety of orbital orderings by starting with different initial conditions. This is
probably due to the absence of the Jahn-Teller distortion which would stabilize the
occupancy pattern of the eg orbitals [107]. Table 4.4 lists the energies of the various
states we have obtained in the cubic case. Roughly speaking, the results suggest that the
energy scales associated with the change in crystal structure (cubic to orthorhombic),
orbital ordering (in the cubic structure), and magnetic ordering are in the ranges of l, 0.1

and 0.01 eV / Mn, respectively.

Table 4.4: HFA energies of different magnetic and orbital ordering states for the cubic
structure. The symbol d(3z2-r2) means that the eg electron is of this type for all Mn,
similarly for d(x2-y2). The third one, d(3x2-r2)/d(3y2-r2), is the observed ordering of the
orthorhombic structure, as discussed in the text. Energies are relative to the AAF state of
the orthorhombic structure, in meV / Mn.

 

 

 

 

 

 

 

 

 

 

FM AAF GAF CAF
«322-12) - 1 147.4 - —
d(x2-y2) - 1165.9 1145.0 1157.5
d(3x2-r2)/d(3y2-r2) 1054.8 1055.2 1088.9 1087.2

 

75

4.4 Discussion and Summary

The extreme smallness of the various energy differences between magnetic states
requires a discussion of the numerical accuracy of our calculations. The sources of error
are [98] the tolerance in the direct-space summations of Coulomb and exchange series
(controlled in the program by ﬁve parameters called TOLINTEG), the number of
sampling points for the Brillouin zone integration (controlled by a shrinking factor IS),
and the ﬁnite basis sets used. The presented results of this work were calculated using
TOLINTEG of (7, 7, 7, 7, l4) and IS of 6. IS = 6 translates to 80 asymmetric k-points in
the Brillouin zone and is more than adequate for our need of accuracy; the total energy
obtained from using IS = 4 (30 asymmetric k-points) deviated by less than 0.003 meV /
Mn. The deviation of total energy due to change of TOLINTEG is much larger; by using
(8, 8, 8, 8, 16) for TOLINTEG, we observed that the total energy changed by about 50
meV / Mn. However, the energy differences between various magnetic states in the test
remained quite stable, typically only varying by ~ 0.1 meV / Mn. Concerning the errors
due to the ﬁnite basis sets, we expect the most severe should come from the La bare core
approximation, which we now discuss.

The test of the bare La+3 ion approximation consisted of adding to the La+3 core an
optimized d shell consisting of a single primitive Gaussian (0.32 Bohr‘2 is the optimized
exponent). The total energy decreased by about 2 eV / Mn; however, the energy
differences between the various magnetic states changed by only ~ 0.1 meV / Mn (a few
%), and the change in the occupied band structure (not shown) is small. Similar behavior

of the basis set dependence was also found in studies of other systems [93]. The charges

76

found in the test were La+2'58, Mnm", O,"'55 and OU'1'59. Interestingly, 0.42 electrons per
La occupied the added La (1 shell, but these electrons were taken from oxygen, making
the results even further from the formal valence picture.

To summarize, our I-[FA results on LaMnO3 give a dramatically different picture
from that of previous LSDA calculations, although for some properties, e. g. orbital and
magnetic ordering, the two theories agree. This agreement is surprising in view of, but
certainly not inconsistent with, the large differences. Some of the HFA predictions have
been supported by experiment. In particular, the HFA spin Hamiltonian is consistent
with a spin wave experiment, but the present accuracy of the experiment is not quite
sufﬁcient to distinguish conclusively between the HFA and LSDA results. The DOS in
HFA is consistent with one interpretation [85] of a photoemission experiment [85].
However it is not consistent with another interpretation [86] of the same experimental
results, based on LSDA. Further, the DOS in LSDA+U theory [87] is more like the HFA
result, and is deﬁnitely inconsistent with the LSDA result. We can conclude that our
HFA results have added weight to the picture where the top of the valence band consists
predominantly of O-p states, as found in the LSDA+U [87] and the cluster model
interpretation of photoemission results [85] (in contradiction to the LSDA results [86—

88]).

77

Chapter 5

Effect of Oxygen Vibration on Spin
Waves in the Double-Exchange Model

5.1 Introduction

In recent years, research activities on perovskite manganites such as A,_,,B,MnO3
(A represents a trivalent rare-earth element such as La, Pr or Nd; B represents a divalent

alkali or transition element such as Ca, Sr, or Pb) have increased dramatically due

primarily to the possibility of important technological applications of the colossal
magnetoresistance observed in about the 0.2 S x S 0.4 region. A,_,,B,,MnO3 has the typical
perovskite structure, where each Mn is surrounded by six oxygens forming an MnO6
octahedron, and each MnO6 octahedron shares one oxygen with each of its six nearest-
neighbor octahedra. Fig. 4.1 in Chap. 4 shows the basic cubic structure of A,_xBanO3,
while the actual crystal may exhibit distortion depending on the constituents and doping
concentration.

When x = 0, each Mn carries charge of +3 with electronic conﬁguration [Ar]3d4.
The substitution of a divalent element for the trivalent element introduces holes and
results in Mn+4 ([Ar]3d3) ions in the crystal. The electronic transport properties of these
manganites are usually described by the so called double-exchange (DE) model [108]. In

the model, only the d electrons of Mn’s are considered. Due to the cubic (or near cubic)

crystal ﬁeld, the Mn d shell splits into t2g (xy, yz, zx) and eg (xz-yz, 3zz-r2) levels with the

78

former lower in energy. In the DE model, three electrons are assumed to be in the tzg

orbitals, localized at each Mn, constituting a local spin 5‘, at i site, S, = 3/2. All

remaining (1 electrons, the number of them depending on the doping concentration, go
into the eg orbitals and can hop from Mn to Mn. The interactions between the t2g and eg
electrons are governed by the strong intra—atomic exchange J (Hund’s rule coupling) and
Coulomb repulsion. Historically, people considered the simpler single-band DE model

(assuming one eg orbital per Mn), which reads:

H=—z Z(c;;c,, +h.c.)—JZ§, -§,, (5.1)

<ij>,a'

where <ij> denotes nearest neighbors, c; (cw) creates (destroys) an eg electron with spin

0' at Mn of site i, §,(§,) is the spin of the localized t2g (conduction eg) electrons, and t (J)

is the hopping (exchange) parameter, both positive. In these materials, J/t >> 1. This

model can explain the concurrence of the metallic and ferromagnetic phases of the

materials in the following way: if all of the local spins 5;, are parallel, then the
conduction (eg) electrons with the same spin direction can move most freely because there
is no energy barrier (due to the exchange term) for them to hop between sites, thereby
minimizing the kinetic, as well as the exchange energy.

The DE model of Eq. (5.1) also predicts a shape of the spin wave dispersion very
much like that of the Heisenberg model with nearest neighbor spin-spin coupling in a
broad range of doping concentration, 0.1 S x S 0.6 [109]. This has been conﬁrmed by
experiments on high—Tc (Curie temperature) manganites, e.g. LamemMnO3 [110].

However, recent experiments on low-TC manganites such as LamCamMnO3 [111] showed

79

strong softening of the dispersion near zone boundaries compared with that of the
Heisenberg model with same spin wave stiffness. Moreover, results of Ref. 111 strongly
suggest that this softening is phonon-induced and the magnon-phonon band crossing
pattern is unusual and therefore intriguing. Under common circumstances, the spin wave
dispersion is a steadily increasing curve from Brillouin zone (BZ) center to the zone
boundary. (In one dimension, it is more or less close to an upside-down cosine curve.)
On the other hand, an optical phonon band remains quite ﬂat throughout the whole BZ.
If these two bands happen to cross and there is some inter-band interaction present,
usually we would expect that a gap opens up near the crossing and two new non-crossing
bands are formed. This expectation turns out not to be the scenario observed in the low-
Tc manganite experiment of Ref. 111. The experimental data show the following: The
optical phonon band remains roughly constant all the way to the BZ edge. The spin wave
band increases from zero at k = 0, but where it meets the optical phonon band, it ﬂattens
(“softens”) and merges with the phonon band, out to the El edge. Thus there is no band
above the phonon band near the B2 edge as would be expected in usual band-crossing
situations. Based on the DE model, we have tried several possible schemes for the
electron-phonon interaction to seek the mechanism which underlies the observed spin

wave softening.

80

5.2 A ﬂaw of popular models in the literature

In fact, even before the recently found spin wave softening many modiﬁcations of
the DE model had been studied in the literature for various reasons, e.g. isotope effect
[112]. Most of the works extended the DE model by including some form of the so
called J ahn-Teller (J T) coupling. According to the well known J ahn-Teller theorem, it is
energetically favorable to distort the MnO6 octahedra to lift the degeneracy of the two eg
orbitals. This interplay between the electronic energy and the geometry of a MnO6
octahedron is therefore called JT coupling and is introduced in many extensions of the

original DE model. A typical model including the JT coupling is the following:

H = H, +th +HJT, (5.2)
(lb + " .. U
He =— 21‘, (€11.00ij +h.c.) —JZ S, as,“ —?Zr,yr,, , (5.3)
<1_'/>,a.ab i,a i.7
M - K
th :—2_ZQ132+EZQIZCI’ (5'4)
HJT : g'Z Q1111 "' g2 (Q1271: + Qi3z-iz)’ (55)

This is a two-band model: it considers both the two eg orbitals of each Mn, which are
described by a pseudospin 1'. Due to some technical details, only x and 2 components of
r are present in the model. The ﬁrst two terms in Eq. (5.3) are just the direct
generalization of the single-band DE model. The third term in Eq. (5.3) takes into

account the on-site Coulomb repulsion between eg electrons, which is irrelevant in the

81

single-band DE model in the large J/t limit but important in the two-band DE model. Eq.
(5.4) is the Hamiltonian of the three relevant normal mode vibrations (with coordinates
Q,, Q2 and Q3, respectively) of each MnO6 octahedron [112]. Eq. (5.5) is the JT coupling
(an electron-phonon interaction), where the Q modes affect the on-site energies of the two
eg orbitals in a speciﬁc way deduced from symmetry considerations ( g' and g are two

coupling constants). Note that there is a deﬁnite proportionality among the components
(lb . . . o . .
of t,,. , but the actual relation rs not uniquely determined and its form may vary slightly

from model to model [113].

By Eq. (5.4), this model assumes that the vibrations of MnO6 octahedra of
different sites are uncoupled, which we consider as a rather unphysical assumption. The
fact that every pair of nearest-neighbor octahedra share a common oxygen implies that

their individual vibrations cannot be uncoupled. In other words, there should be some
cross terms (Q,Q_, and Q,Q, terms) in Eq. (5.4). Unfortunately, the signiﬁcance of this

assumption has not been discussed in most of the literature. [114]
One may speculate that maybe the cross terms are small in the ﬁrst place. This
turns out not to be the case. Furthermore, there are more complications than meet the

eyes. To illustrate the point, consider a 1D chain of Mn-O as shown in Fig. 5.1:

o O o O 0 O O O
i+2 u.

z u. 1+1 u. 1+2

1'" I [H
o manganese O oxygen

Figure 5.1: The Mn-O linear chain

82

The manganese atoms are assumed heavy and therefore still. The ui’s are displacements
of the oxygens from their equilibrium positions. Suppose there are N unit cells with
periodic boundary condition, i.e. uN+i = 2.1,. The Hamiltonian of these vibrations can be

written in terms of ui’s and their conjugate momenta:

H —_1_ﬁ 2+£ﬁu2 (56)
pk 2M i=1 pi 2 i=1 1' -
In this 1D model, the relevant normal mode coordinate of each unit cell would be

Q.- = $0.- —u.--1)- (5.7)

Now, can we transform the Hamiltonian to a form in terms of the Q variables? The
answer is NO. Because of the periodic boundary conditions, we have

2Q, = 0. (5.8)

1:]

This means that we only have N-l truly independent Q variables. Therefore the
transformation from the u variables to Q variables is impossible unless we add an
arbitrarily chosen new variable to the Q set. A natural choice for the Nth variable (let’s

call it Q,,) will be

Q0 172“; (59)

Now we are able to proceed with the transformation of Eq. (5.6) from {u} to {Q}. The
result is not in a simple form and therefore is omitted here. But the general properties of

the results can be described easily: After the transformation, terms Q,Q,- with any i and j

83

are present. For i and j nonzero, the coefﬁcients of Q,Q,- (i at j) are of the same order of
magnitude as those of the diagonal terms Q,2 when |i - j] are small, then their magnitudes
decrease gradually with Ii - j| and become near zero when [i - j] reach the maximum.
Furthermore, cross terms not only occur in the potential energy part (the Q,Q,- terms), but
also in the kinetic energy part (the P,P,- terms) with the same coefﬁcients. This means that
the neglect of cross terms in Eq. (5.4) is not justiﬁed.

There is actually a paper by G. Khaliullin and R. Kilian [115], which includes
cross terms of nearest-neighbor Q’s with the coefﬁcient as an adjustable parameter,
claiming to have solved the spin wave softening problem. However, there is a ﬂaw in

their model: cross terms of the kinetic energy part are missing.

5.3 Model Studied

In view of the question raised in the previous section, we suggest that the Hp, in

Eq. (5.4) should be replaced by
Hp, =%Zu§a +—I-2<—Zu;, (5.10)

where the u variables denote the displacements of individual oxygens. The H,T in Eq.
(5.5) remains unchanged but should be now written in terms of the u variables.

To seek the mechanism which underlies the spin wave softening near the zone
boundaries, we explored several possible schemes of inclusion of electron-phonon

interaction into the DE model. To get started, we considered the 1D single-band model

84

with only oxygen vibrations taken into account. The portion of the crystal structure of
interest to the model is the Mn-O linear chain, as shown in Fig. 5.1. The Hamiltonian

considered has the following form:

HZHe+th+He—ph’ (511)

H, =-tZ(c;;c-,,,a +h.c.)-J28, 5,, (5.12)
1 K

H =—— .2+— u.2. 5.13

ph 2MIZPJ 2 Z]: 1 ( )

The electronic part He is just the original single—band DE model. The phonon part Hp,
represents the simple-harmonic oxygen vibrations (u,- and pi are the displacement and
corresponding momentum of oxygen at site i). He_p,, stands for the electron-phonon
interaction part and hasn’t been speciﬁed above. Which form of the electron-phonon

interaction H would give rise to the spin wave softening is not known a priori. It

e-ph

could be that displacements of the neighboring oxygens would affect the local on—site

energy a of the eg orbital. Assuming linear response 8, oc (u,-u,_,), then H would have

e-ph

the form:
He—ph : ZEiCLCia = _g2(ui — ui—l )Citrcia ' (514)

This is analogous to Eq. (5.4) for the two-band model. In addition, we have also
considered another possibility where oxygen vibrations affect the local electron hopping
parameter t (in both magnitude and phase). The results of the ﬁrst scheme (changing on-

site energy) are presented in next section, while the results of the other scheme will only

85

be mentioned brieﬂy in Sec. 5.5, Remarks and Conclusions, since its results are not any

better than those of the ﬁrst scheme towards the explanation of the spin wave softening.

5.4 Results

We solved the Hamiltonian of Eq. (5.11) — (5.14) for a ring with 4 Mn’s and 4 0’s

in the S2101“, = Smaximum — 1 manifold in two different approximate ways.

5.4.1 Pseudo-Born-Oppenheimer approximation

In this approximation, the oxygen is assumed heavy and Hp, is neglected. Each u,-
in H,ph [Eq. (5.14)] is randomly taken from the interval [-6/2, 6/2] with a uniform
probability distribution. For a set of {u,} randomly taken, He + He_p,, is solved exactly,
giving a set of electronic eigenenergies. Then another set of {u,} is chosen, and the
procedure repeats. We ran many iterations to get a statistical distribution of the electronic
eigenenergies. A typical case is shown in Fig. 5.2 with ground state energy set to 0.

Fig. 5.2(a) shows the energy levels of the ordered system. The lowest three peaks
carry the wave vector quantum number k of 0, win/2 (two-fold degeneracy) and 7r,
respectively. These are excited sites due to creation of a single magnon, with wave vector
k, from the ground state. For an inﬁnite ring, these three levels will expand into a
complete band and become the familiar spin wave dispersion. Thus we refer to the states

associated with these three peaks as the spin wave states.

86

 

 

Fig. 5.2(b) shows clearly that each energy level of the ordered system is now
broadened by the random electron-phonon interaction. Compared with the width of the
random distribution of ui (0.5 in this example), the broadening of the spin wave levels is
surprisingly small. Other higher energy levels do have much larger broadening. The
results suggest that the spin wave dispersion remains robust in the presence of the
electron-phonon interaction of the form of Eq. (5.14). Moreover, we also calculated the
distribution centers of the spin wave states in Fig. 5.2(b) and found that the ratio of the
(averaged) k=7t magnon energy to the k=7r/2 magnon energy increases slightly compared
with that of the ordered system [Fig. 5.2(a)]. Therefore this model actually predicts a

small hardening in the spin wave dispersion, opposite to the experimental observation.

87

 

4 sﬂfs , 3 ehajzrns

-—‘— ._ —_ ,—. ~.~_. .. 14 -77 7V i—fﬁ

 

0
CD

O
ON

 

DOS (arb. unit)
0
:5

O
l\)

(b) disordered system

Figure 5.2: Density of states of the DE model for a 4-site, 3-electron ring.
(t=1,J=oo,g=1,6=0.5)

88

5.4.2 QM treatment of the oxygen vibrations with cut-off in
the number of phonons

In the second approximation, we treated the oxygen vibrations fully quantum-
mechanically but a cut-off in the number of phonons was imposed. This approximation
was checked for convergence by using different cut-offs. The energies of the spin wave
states are shown in Fig. 5.3, plotted as E versus k. The result is consistent with that of the
ﬁrst approximation; namely, the spin wave energies are very robust in the presence of the
electron-phonon interaction Eq. (5.14) and the spin wave dispersion is actually slightly

hardened in this model.

 

 

 

 

i _._ 9:0
0.15] * g=0 1
0.1: —I~ g=0 2

i 4- g=0.3
0.05

 

 

 

 

Figure 5.3: The spin wave dispersion of the DE model for a 4-site, 3-electron ring. (cut-
off ofphonon number = 4, t =1, J = 00, (up, = 0.15)

89

5.5 Remarks and Conclusions

We have criticized the unphysical assumption of uncoupled vibrations of the

MnO6 octahedra in studies of the DE model with JT coupling, which unfortunately is
prevalent in the literature. We suggest that the Hp, in question should be changed from
I Eq. (5.4) to Eq. (5.10).

In seeking the underlying mechanism of the recently observed spin wave
softening in low-Tc manganites, we considered the simple 1D single-band DE model with
electron-phonon coupling in the spirit of our suggestion mentioned above. Given the
Hamiltonian of Eq. (5.11)-(5.l4), the spin wave energies hardly change at all in the
presence of the electron-phonon interaction. A study of the ratio of the energies of the
k=7r and 7r/2 magnons shows that the model actually predicts slightly hardening in the
spin wave dispersion, contrast to the large softening found in experiments for low Tc
manganites.

Although omitted here, another possible scheme for the He_,,,,, in which oxygen
vibration can modify the local electron hopping, is also studied. Two cases in this regard

are considered: One is

HM, = a2 uf (4011“,, + h.c.), (5.15)

where a is an adjustable but positive parameter. Note that the sign of this term is
Opposite to that of the hopping term in Eq. (5.11); that is to say that the displacement of

oxygen, u,, can reduce the magnitude of the local electron hopping. The proportion to u,2

90

is due to symmetry consideration; the reduction of the hopping magnitude shouldn’t

depend on which way the oxygen moves. The other case is
Hm), = _tZ(e’“"‘ —1Xc,+ac,+,‘a + h.c.)z —iatZu,(c,;c,+,ﬂ + he). (5.16)

This assumes that ui can alter the phase of the local electron hopping. Despite the rather
different mechanisms deﬁned in Eqs. (5.14), (5.15) and Eq. (5.16), their effects are very
similar for the spin wave states alone (while they make larger difference for other states).
Namely, the spin wave states are robust and the dispersion gets very slightly hardened in
the presence of these electron-phonon interactions.

This strongly indicates that the observed softening of the spin wave dispersion
cannot be explained by the simple single-band model. A more realistic two-band model
which incorporates the Jahn-Teller effect such as Eqs. (5.2), (5.3), (5.5) and (5.10) may

be needed in this respect.

91

Chapter 6

Concluding Remarks

This work gives several new results which are likely to have a broader impact on
the research community on these subjects. First, although the three-band Hubbard model
with the model parameters ﬁtted to high-Tc superconducting cuprates [12] has a large
[pd/8pd (~0.36) ratio, the perturbation expansion using the effective Hamiltonian method
on this model for the undoped case is proved to converge well in the leading order. It
explains why the low-lying states of the undoped cuprates are found to be governed by
the nearest-neighbor Heisenberg spin Hamiltonian in a study of the model on small
clusters [16]. However, the rapid convergence of the expansion doesn’t come from the
straightforward perturbation. In fact, the straightforward perturbation does not converge
at all. To achieve convergence, one has to ﬁrst perform a proper canonical
transformation on the model Hamiltonian and then apply the perturbation method to the
transformed Hamiltonian. In this work, it is shown that transformations involving only
recombinations of the underlying single-particle orbitals (called l-body transformations)
lead only to a slow convergent expansion. The reason for the slow convergence is the
following: 1-body transformations are usually aimed to reduce the hopping between the
transformed orbitals. But at the same time, it also reduces the intra-atomic Coulomb
repulsion in the new (1 orbitals due to the p-d orbital mixing (the Coulomb repulsion in
the p orbitals is much less than that in the d orbitals). This repulsion energy occurs in the

energy denominators of the perturbation expansion, associated with the intermediate

92

states where two holes occupy the same copper site. Thus the reduction of this repulsion
energy has a negative impact on the convergence of the perturbation expansion. This is
an inherent limitation of l-body transformations. I have shown that an l-body + minimal
2-body transformation can simultaneously reduce the hopping and retain the large
repulsion energy, and leads to rapid convergence.

Second, the crystal Hartree-Fock (HF) calculations on the spin form factors of
LaQNiO4, LaQCuO4 and KNiF3 in this dissertation conﬁrm the results of Ref. 41 which
explored the embedded cluster approximation. In both the works, Eq. (3.3) is used as the
approximation to obtain the spin form factor. The authors of Ref. 41 actually have
calculated spin form factors for a much broader range of Mott insulators. Their results all
agree very well with experimental measurements except for the case of LaQNiO4 in which
a large discrepancy between the theory and an experiment [64] is observed. Based on the
agreement between this work and Ref. 41, and the fact that Eq. (3.3) has been
successfully applied to many other materials in Ref. 41, we cast doubt on the
experimental results of LazNiO4 and call for more experiments to settle the issue.

Third, the crystal Hartree-Fock approximation (HF A) is proved to be valuable for
the solids with strong interacting electrons. For example, the crystal HFA predicts
correctly that LazCuO4 and LazNiO4 are both antiferromagnetic insulators, in contrast to
the wrong prediction (paramagnetic metal for LaQCuO4 and antiferromagnetic metal for
LaQNiO4) made by the local spin density approximation (LSDA) [37,38]. It is commonly
perceived that, in general, the HFA overestimates band gap and the LSDA underestimates
or predicts no band gap for insulators. However, there are other important differences in

the predictions of the two theories which have not been tested by direct experimental

93

measurements yet. For example, the crystal HFA predicts the oxygen p band while the
LSDA predicts the magnetic ion (1 band to be the top valence band for the perovskite
Mott insulators. This is a fundamental problem in understanding the electronic structures
of these systems, but has not yet been resolved by direct experimental measurement.
Besides, it is found in the crystal HFA that the charge of the Mn ion is close to +2 for
both LaMnO3 and CaMnO3, in contrast to the formal valence picture where it is Mn"3 for
LaMnO3 and Mn’“4 for CaMnO3. What makes the matter more intriguing is that the spin
of the Mn ion predicted by the crystal HFA is very close to the value of the formal
valence picture (2 for LaMnO3 and 3/2 for CaMnO3). Based on these results, we
hypothesize that high charge ion states are not energetically favorable in reality and, for
systems containing ions with large formal charges, there is a spinless charge backﬂow
from the anion (such as O, F) p orbitals to the magnetic ion (such as Cr, Mn) d orbitals
(thus it only changes the charge but not the spin of the magnetic ion). Since the formal
valence picture has been a universal concept in general physics and chemistry, the
importance of examining this hypothesis, which is a departure from the formal valence
picture, is obvious.

Fourth, in this work it is argued that the common treatment of the double-
exchange (DE) model with Jahn-Teller (JT) coupling found in the literature [112] is
seriously ﬂawed. The DE model is commonly used in describing the electronic
properties of the perovskite manganites. The point of concern here is the vibrations of the
MnO6 clusters in the crystal, which affect the electronic properties through the JT
coupling. Obviously for the sake of convenience, most people treat the vibrations of the

MnO6 clusters as uncoupled. This leads to an approximation which essentially neglects

94

all the cross terms in both the kinetic and the potential energy parts of the Hamiltonian
for the MnO, vibrations. Using a simple example, I show that the coefﬁcients of the
cross terms are just as large as those of the diagonal terms. Thus the neglect of the cross
terms is not justiﬁed.

Fifth, the energies of the spin wave states in the single-band DE model are found
to be surprisingly robust in the presence of several trial electron-phonon couplings.
Based on this, it is believed that the mechanism responsible for the recently observed
phonon-induced softening of the spin wave dispersion near Brillouin zone boundaries in
several low-TC (Curie Temperature) manganites [111] lies beyond the simple single-band
DE model. It remains to be seen whether the two-band DE model with JT coupling can

explain this phenomenon satisfactorily.

95

APPENDICES

Appendix A

A physical interpretation of the TSDA
equations

In the one-body terms of the Hamiltonian, hvﬂc:c#(v¢ ,u) is the so called

hopping term which causes a particle initially in the one-particle state it to hop to the one-
particle state v, provided the one-particle state v is empty initially. Therefore, for an
arbitrary single-determinantal state |‘P>, the expectation value of the operator

hvﬂn #(1 — nv) gives the amplitude of the hopping u—w when [‘1’] is operated on by the

hopping terms of the Hamiltonian; i.e.,

<LIJ h n (1_ anLP> : <LP'IHf—body

w /1

 

‘1’), (A1)

 

where H,'_body E (H — Ha),_body = Zhaﬂckﬂ and I‘lj') is the single-determinantal state
(1:38

which differs from |‘IJ) by replacing, if there is any, the one-particle state it in [‘1’] with

the one-particle state v. For the same reason, the operator hv,,nv(1—n,,] (note that

hw=hw are assumed to be real here) corresponds to the reverse hopping v—>u. The

difference between these two operators:
hwnv(l—n,,)-hvﬂn,,(l—nv)=hv,,(nv—-n,,) (A2)

reﬂects, in some sense, the “net” hopping v—>p.. Therefore we deﬁne it as the net

hopping operator due to the one-body terms of the Hamiltonian.

97

Now, let’s turn to the two-body terms of the Hamiltonian. The analogous hopping
terms which cause a particle in the one-particle state u to hop to another one-particle state
v while other particles in the system don’t move are vv,‘,,,,c:c;c,c# and vv,‘,#c:c;c,,c,.
What’s different from that of the one-body terms is that this hopping must occur in the

presence of another particle in the orbital )1. Also keep in mind that the second term,
vMMC: c] cﬂc, , introduces an extra minus sign in the ﬁnal state from that of the ﬁrst term,
vv, 1.16:6: c ,, c ,1. Analogous to the case of the one-body terms of the Hamiltonian, for an
arbitrary single-determinantal state |LI’) the expectation value of the operator

2694.714 —vv,,,,, )n,n,,(l—nv) gives the amplitude of the hopping u—w when [‘11) is

A

operated on by the two-body terms of the Hamiltonian,

 

r
HZ-body

 

)1!) (A3)

<LIJ I; (Vt/2,112 _ V1722); )"anp (1 _ ”V X ‘1’) = (‘11,

where H2’_body denotes the two-body terms of H — Ha and H“) has the same meaning as

in Eq. (A1); and correspondingly, 2(1),,” —-v,,,’ 7111 )n ,nv (1 —n y) accounts for the reverse
1

hopping v—au. The difference,

2 (Vt/2,112 _ Vwmy >11n1(1 _ ”/1 )_ Z (Vt/71,111 _ Vv2,2)1)"2"11(1— ”v)
4 2
= 2644.14 “ Vvuy )"4 (”v " ”11)
2

(A4)

is deﬁned as the net hopping v—>u due to the two-body terms of the Hamiltonian. Then

it’s natural that the total net hopping operator is
hv;1(nv —ny)+Z(vt/A,M —va,/lp )nl (”v _n,u)’ (A5)
21

98

which is just the left-hand side of the TSDA equation (2. 10) without the thermal average
acting on it. Thus we get another view of TSDA; namely, the solutions of TSDA give
zero thermal average of the net hopping operator (again, note that here the thermal

averages are not with respect to the exact Hamiltonian H but to the approximation Ha).

99

Appendix B

Solving the TSDA equations

Since the detailed use of the TSDA in the present context is not commonly seen,
we present the details to aid in further pursuit of the matter. Because of the numerical
difﬁculty mentioned in the text and explained later in this appendix, what we have
actually solved are modiﬁed TSDA equations, where the modiﬁcation is based on
argument of the large [3 (or low temperature) region of interest. In the following, we
consider site localization as an example to explain the difﬁculty and show how we

circumvent it.

The TSDA equations are: [Eq. (2.10) of the text]

< n, -n,, > 11,, +Z< (n, —n,,)n,, > (11%,, —v,,,,,) = 0 (B1)
A

for any pair of orbitals labeled v, u. The thermal average bracket, < -- >, is deﬁned as, for

any operator A,

= tr(e"ﬂ”1A) _ Z < ”an'13-ﬂﬂ“ 41'",va >

< A >_ — , B2
tr(e‘ﬂ”" ) Z < .., NV,--|e_ﬂH“ I", NV," > ( )

 

 

where the summation is over all possible occupations. The subscripts in Eq. (B1) label
both the spatial and spin quantum numbers of the one-particle states. In the Hubbard

model, the one-body and two-body matrix elements, h,” and vm 7111-» preserve spin (i.e.

_ . ,. . . _ _ . , , .
h,“ h,0,ja.—01f0'¢0',51m11arlyv =v —01f0'¢0' ors¢s)andaresprn

v,u./1K i0js,kcrls'

100

independent (i.e. hm}, E 11,; similarly, v,,j,,,,,, a vW). So the TSDA equations in this

case can be written in the following form:

< n,, — 11,, > h, + Z< (n,, — n,,)n,,. > (v,,,,, — 5,,.v,,,,) = 0 (B3)
k.o'
for any pair of spin-orbitals ,, and 1,. In the large B regime in which we are

interested, the exponential factors in the thermal averages of the number operators in Eq.
(B3) cause numerical difﬁculty and need special treatment.
In viewing the diagram of the CuzO7 cluster (Fig. 2.2), it is apparent that all the

orbitals can be divided into four groups:

((1,311; ), (p1), (péaplapéapé), (19331012). (B4)

Each orbital is said to be equivalent by symmetry to any other orbital in the same group.
Because no original orbital (unprimed) will occur in this discussion, the prime on the new
orbitals will be suppressed hereafter for brevity. By symmetry requirements, it is easy to
see that the TSDA equation (B3) for any pair of orbitals is automatically satisﬁed if the
two orbitals of this pair are in the same group. These trivial equations give us no
information in determining the new orbitals. Further, careful counting shows that, among
those remaining non-trivial TSDA equations, there are precisely nine independent
equations; they are the TSDA equations for the pairs

(d. 11p. 0.0. tpz Md. 1.1230141 T1121 MA 1,1911),

(135)
(10. T1192 TM101 T1113 T),(p2 T’10:. T),(p2 T1106 T)-

(The corresponding equations for the down-spin pairs are identical to these.) If B" is far

below the transfer energy ap-g d but much larger than the parallel-antiparallel splitting of

101

the underlying states, i.e. K ,, << ,6" << 8,, — 8,, then the states which have one hole in

each of the Cu d orbitals will dominate the thermal averages and the weights they
contribute to the thermal average are essentially equal no matter what their spin

conﬁgurations are. For the CuZO7 cluster this means

— — _ ~ 1
<”d,t >—<”d.t >—<nd2, >—<nd2, >= 2
, . (B6)

< np, >='_- 0 for all z =1..7 and 0' =T,i

Because the TSDA equations for the last four pairs of Eq. (B5) only involve the thermal
averages < "M > and < n Mn, > , their numerical values for the left hand side of Eq. (B3)
are very small. This causes difﬁculty in solving numerically the nine simultaneous
TSDA equations resulting from Eq. (B5) since the computer program may only “see” the
ﬁrst ﬁve equations of Eq. (B5) and regard the last four as having been solved. To
circumvent this problem, we keep only the leading terms in the thermal averages,
substitute them into the TSDA equations, and then extract the limiting forms of the
equations in the large B regime, as shown below.

The lowest energy group of the half-ﬁlled CuZO7 cluster is the set of 4 single

determinants
‘IJ(d,0',dzo"), (B7)
where 0' and o" are either T or \L. The next three higher energy groups are
.‘I’(d,0',p,o"), LI’(d20',p,0")

2. ‘1’(d20', p20”), LI’(d20',p4o"), ‘I’(d,o,p50"), Ll’(a’,0',p70"). (B8)
3. L11(d20',p30"), lI’(d,0',p,o")

102

(We don’t know the energy order of the three groups in Eq. (B8) before we solve the

TSDA equations.) To circumvent the numerical difﬁculty mentioned above, we make

two assumptions about B: (i) B" is sufﬁciently low such that in the thermal average the

exponential weight due to a lower energy group numerically overwhelms the weights due

to other higher energy groups. (ii) But at the same time B'1 is sufﬁciently high such that
the parallel-antiparallel splittings are to be neglected and average energies (over 4
possible spin conﬁgurations) are to be used in substitution of real individual energies.
For example, deﬁne E (d, 0', p,o") as the energy of state LI’(d,o*,p,o") and E (d, p,) as

the average of E (d ,0', p,o") over the 4 possible spin conﬁgurations,

1 .
E<d1p1)= 2219811011210)

1

: 22(hdldl + hPIPI + vdlplsdlpl — aaa'vdlpl’pldl ) (B9)

_ __1.
— hdldl + hPIPI + vlervdIPI 2 vdlpl'pldl

Because there is essentially no energy resolution among the different spin conﬁgurations
by our assumption, we are to replace E (d,0', p,o") with E (d , p,) in our calculation. The

energies of other states are similarly treated. We will come back to our assumption about

B later. Some thermal averages of the number operators are listed below as examples:

103

l _ _ 1
< n >5 —(e ﬂw'a‘dm + e [’W'O'dz“) E 2 < n n . >5 —
d,0' Z (1,0' dzo 2

1 _ _ _ _
< n >5 —(e ﬂE<p.a.d.T) +e 415010.411) +e 450210.42?) +e 4402104120)E 4 < "plan“, >5 0

l _ _
< n >5 E02 Wham“ +e Who‘d”) 5 2 < n,,,nd . >5 0

20'

1 _ _
< n >2 —(e ”5(p30’d2T) +e Who‘d”) 5 2 < n, ,n, ,. >2 0
Z 3 2
(1310)

where Z denotes the partition function. The TSDA equation for the pair (d, T, p, 1)

therefore reduces approximately to a simple form in our assumed B regime:

k,0'

l _
I) (3 — O)hd;P1 + (< ndllndzl > —O)(vd1d21prd2 _ vdldzﬂzpl ) + (< nlendzl > ’0)(vdld2'pld2 ) _ O

_i =
3 hdlpl + (vdldzvpldz 2 vdldzsdzpl) 0

The TSDA equations for the 2nd to 5th pairs in Eq. (A5) are similar in form to that for

(d, T, p, T) . Now we examine the TSDA equation for the 6th pair (p, 1, p2 T):

(< ”911 > _ < ”I721 >)hPIP2 + Z“ ”mink“ > — < npzln’m >)(vp1k.p2k — 61avp1k1kpz ) = O
k,a

3 (< ”PIT > _ < "p21 >)hPIP2 + (< npllndll > _O)(vPid1~P2dl — vpldl’dlpz ) + (< nplTnle' > _0)(vp,d,,p2d,)

+(< npl,nd2, >—<np2,nd2, >)(v )+(< nPITndzl' >—<np2,nd2, >)(v )=0

Pldzspzdz _ vp.d2,d2p2 Pldz'Pzdz
_ _ L

:> (2 < "p1tnd1a > < np2,n,2, >)hp1p2 + < "p.tndp >(vl’1dhl’2dl 2 Vp1d.,d1pz)

) = 0

_ _ L
+ (< np,1nd,-0 > < npzl‘ndza >)(vpid2~P2dz 2 vp,d2,d2p2

where in the last expression d,, can be either (1,, , d,, , d2, , or d2,. If
E(p,d,) < E(p2d2) , then <npl,ndi,> >> <npz,n,2,> by the former assumption about B and

the previous equation reduces to

104

_ -1. _ 'l—V
2hPIP2 + (VPJdlapzdi 2 vpldl’dlpz ) + (vP1d2~P2d2 2 vprdz d2P2 )= 0’

on the other hand, if E(p,d,) > E(p,d,), then <"p.t”d.-o> << <np2,n,2,> and the equation

reduces to

)= O

hP1P2 + (VP1d2aP2d2 —2 VPrdz d2P2
The limiting forms of the TSDA equations for the 7th to 9th pairs in Eq. (B5) can be

obtained in a manner similar to that for (p, T,p2 1). All the nine equations are

summarized below:

(d, t, p, t) : 71,, + (v,,,2,,,2 — 4v,,,,,2,, ) = 0 for i = 1,2,3,5,6

(AT45TVWQA+WmMa‘2%am)+WMMA- 4v,,,m)=(nfmpwd)<mpd)
hp],2 +(Vp1d2.p2d2 —4vpld2,,2p2) =0 if E(p1d1)> E(p d)

(p, t, p, T) : 2’14)». +(v,,,,,,,l —4v,,,,,,,3)+(v,,,2 ,3, 412,, ,2”): 0 if E(p1d1)< E(p d )

mm+WMMﬂ- zvmgmg= OdEde)>ﬂpd)

)= 0

)= OﬁEwwd)<ﬂpd)

h
(p. Tm, t) : h
(p. 1.14.04

h

Iv
+ V
P2P3 ( Pzdszsz —2 vadz d2P3

__l_v
+ V
le’o ( Pzdszodz —2 vadz vdzPo

)= onEde)>Emd)

_ 1v
P2Po +(vP2d1rP6dl 2 vadl dIPo

(311)

There are 3! = 6 possible orderings of the energies E(p,d,), E(p,d,), E(p,d,). For each
possible order, we solve Eq. (B11) and then check the consistency of the solutions with
their respective presumed energy order. Those which are not self-consistent are, of
course, dropped. Although the forms of Eq. (B11) are pretty simple, the transformed

parameters h”, and v in Eq. (B11) are each a polynomial function of the

Vinita“

transformation coefﬁcients deﬁned in Eq. (2.16) of order 2 and 4, respectively. Therefore

105

there indeed exist many solutions. We solved Eq. (B11) numerically. By varying the
starting point and length of searching steps, we got several solutions of Eq. (B11). These
solutions were examined by comparing their resulting TSDA free energies and also by
cross checking with the solutions obtained in the cell localization and no localization
cases to conﬁrm that the real TSDA solution (which means the one which, within the
single-determinantal and the particular localization-choice framework, gives the lowest
free energy) has been found and identiﬁed. It turns out E(p,d,) < E(p,d,) < E(p3d2) for

the real TSDA solution.

Now we examine the consistency of the solution with our assumption about [3,
which was made to simplify the TSDA equations and therefore circumvent the numerical
difﬁculty it leads to. Based on the solution, we calculate the energies of various states.
The average energies (over 4 possible spin conﬁgurations) of the ﬁrst four lowest energy

groups, i.e. those of Eqs. (B7) and (B8), are

E(d.d2) E(p,d,) E(pzd.) E(p.d.)

. (312)
—2.75742 3.80976 3.40425 4.6475

The smallest energy gap among the above four groups is between the second and the third
ones, which is about 0.41. For our assumption to be right, the [3" must be well below
0.41 to make a lower energy group numerically overwhelm other higher energy groups in
contributing the weight in thermal average. However, the parallel-antiparallel splittings,
which we have neglected in order to replace the energy of every state of a group with the

average energy of that group, for the four groups in Eq. (B12) are

106

K(d1d2) K(p1d1) K(p2d2) K(p3d2)

0.0288 0.7605 0.0097 0.0823 ' (B 13)

(Note that the second one is big because the p1 and d1 are nearest neighbors to each other
so that they are from strong hybridization of the original p1 and d], and that strong
hybridization also leads to big exchange energy between the new p1 and d1 orbitals.) So
[3" must be much larger than 0.76 for the neglect of parallel-antiparallel splittings to be
legitimate. An inconsistency thus occurs.

We have two reasons to believe that this inconsistency does no harm to our work
here. First, in addition to what we presented so far (which we call the average case), we
have done the same calculation on two artiﬁcial cases: one is to assume the spins of the
two particles in the Cu207 cluster can only be parallel and the other, antiparallel. As a
matter of fact, the parallel case is not really artiﬁcial because it is exactly the limiting
situation when B"—>0. (Note that a state with a pair of spin-parallel particles is always
lower in energy than its antiparallel counterpart by an amount of the exchange energy K
between these two particles. When [5" is sufﬁciently low, the parallel states will
dominate the thermal average.) At ﬁrst glance, the inconsistency suggests that we
shouldn’t assume [3" is sufﬁciently high and neglect the parallel-antiparallel splittings. It
turns out that the TSDA solution of the average case is somewhere between the TSDA
solutions of the parallel and antiparallel cases, and, more importantly, they only slightly
differ from each other. This suggests that the solution we obtained is close to the solution
one might get by solving the exact (unsimpliﬁed) TSDA equations (for ,6" << 0.4 ).

Second, the average case itself can also be viewed as the exact calculation for another

107

artiﬁcial case, namely, a system which has no spin degrees of freedom and has the
average energies as the energies of its levels. Our purpose is to seek a deﬁnite guideline
which can lead us to a good transformation on the basis set such that the perturbation
converges rapidly. So the model which we use to derive the transformation may be
closely related to the system of interest, but it doesn’t have to be the same. Also, recall
that [3" here is rather a parameter than a physical temperature. Therefore, the

inconsistency is not really our concern as long as we can obtain a suitable transformation.

108

Appendix C

On the determination of the Heisenberg J
from electronic structure calculations

In this appendix, we discuss the theoretical justiﬁcation behind the commonly
used mapping scheme mentioned in Sec. 3.2.1 in estimating the Heisenberg exchange
parameter J from ab initio electronic calculations. We assume that the problems under
consideration have their low-lying states governed by a Heisenberg Hamiltonian. Of
course, if a problem is simple enough, one can calculate J exactly (using standard
quantum chemistry techniques such as conﬁguration interaction). For example, if one
were actually dealing with two magnetic electrons, with a limited number of non-
magnetic electrons, one could in practice calculate the energies of the true lowest energy
eigenstates, singlet and triplet. Equating their energy difference to the difference in
eigenvalues of the Heisenberg Hamiltonian, would then yield J. This approach has been
used [1 16-1 18,44-46].

However, for more complicated systems, we are forced to make approximations.
The mapping scheme to estimate J in Sec. 3.2.1 has been used for years, in both Hartree-
Fock and density functional calculations. We have shown [83] that for some special
cases, e.g. a single-band Hubbard model, JHF deﬁned by Eq. (3.2) equals J to the leading
order in perturbation with small hopping; but that does not hold in general (although one
may get the contrary impression from some papers [46]). Still, J HF serves as a reasonable

and well-deﬁned estimation of J.

109

We now discuss the determination of J within density functional theory (DF T)
approaches, which invoke the following procedure [119,120]. One obtains a symmetry-
broken (SB) solution for the Kohn-Sham determinant DKS for the “antiferromagnetic”
state and a triplet state (for simplicity we discuss the case of one magnetic electron per
magnetic site in a two-magnetic-site cluster). The triplet energy minus the SB AF M state
energy is equated to the triplet energy minus the (average) energy of the SB spin state
a(1),8(2) for the Heisenberg model, thereby obtaining a number for J [121]. Although
making the correspondence “SB electronic state H SB spin state” sounds reasonable, one
must realize that there is not a unique SB spin state, even in this simplest of cases.

The following spin state is SB for all values of ,u¢0 or 00.

Q. = —1——{a(l)ﬂ(2) — [301042) + #[a(1)ﬂ(2)+ 40104211} . (CI)

1/2]1+,u2]

The expectation value of the Heisenberg energy is

(C2)

 

If one used the exact density functional (DF), then the energies of the exact lowest
triplet and singlet states would be obtained. The difference in these energies would be
equated to the corresponding difference in Eq. (C2) with ,u = 0 and 00, and this would
give the exact J. However, given an approximate DF, as in almost all applications, one is
faced with the problem of what value of ,u to use in Eq. (C2) for the SB AFM state.
Many people [57,119,120] take [2 = 1 without giving a reason. But there is a potential

inconsistency in doing so: as one’s approximate DF gets closer to the exact DF, continued

110

 

use of ,u = 1 may lead to serious errors. Thus the appropriate value of ,u to use is

unknown.

A possible way out of this dilemma is as follows. From Eq. (C2) one gets

_ §2
u—\/< %2—«62 >), (C3)

where S = S, + 5,. One could demand that < S2 > be equal to the expectation value of

 

the square of the total electronic spin in DKS, thus determining ,u via Eq. (C3). One
should remember that there is no guarantee that DKS is an eigenfunction of S2 or that the

average value < S2 >in DKS is correct, even for the exact DKS, presenting a possible
inconsistency in this approach. Also, how to generalize this mapping idea for
determining the SB AFM spin state to more than two sites (two electrons) is not obvious
and remains to be investigated.

Interestingly, this idea raises a question about the standard mapping used for the

HFA where ,u is taken to be 1 for the AFM spin state. For this two-site two-electron
example, one can show easily that < S 2 > D”; = 1— 62 , where DHF is the HF determinant of

the AF M state and 6 is the overlap integral between the two localized orbitals. Since 6 is

usually quite small (< 0.1), it is seen that this probably would make little difference.

111

REFERENCES

[1]
[2]

[3]
[4]

[5]
[6]

[7]
[8]

[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]

[22]

[23]
[24]
[25]

[26]
[27]

[28]

N. F. Mott, Proc. R. Soc. London A 382, 1 (1982).

S. Das Sarrna, A. Pinczuk, eds., Perspectives in Quantum Hall Eﬂects, Wiley,
New York (1997).

W. Pickett, Rev. Mod. Phys. 61, 433 (1989).

T. A. Kaplan and S. D. Mahanti, eds., Physics of Manganites, Kluwer Academic /
Plenum Publishers, New York (1999).

J. Callaway and N. H. March, in Solid State Physics 38, 136 (1984).

T. A. Kaplan and S. D. Mahanti, eds., Electronic Properties of Solids using
Cluster Methods, Plenum Press, New York (1995).

P. W. Anderson, Phys. Rev. 115, 2 (1959), and in Solid State Physics 14, 99
(1963)

R. S. Tu and T. A. Kaplan, AIP Conf. Proc. No. 24, p. 218 (1974); R. S. Tu, Ph D.
dissertation, Michigan State University, 1975.

See B. Brandow, Advances in Physics 26, 651 (1977) for other criticisms and a
general history of the early work.

T. A. Kaplan and P. N. Argyres, Annals of Physics 92, 1 (1975).

K.-I. Gondaira and Y. Tanabe, J. Phys. Soc. Japan 21, 1527 (1966).

M. S. Hybertsen, M. Schluter, N. E. Christensen, Phys. Rev. 39, 9028 (1989).

J. H. Jefferson, H. Eskes, L. F. F einer, Phys. Rev. B 45, 7959 (1992).

V. I. Belinicher and A. L. Chemyshev, Phys. Rev. B 49, 9746 (1994).

T. A. Kaplan, S. D. Mahanti, Y-S Su, K. Kubo, J. Appl. Phys. 79, 6433 (1996).

M. S. Hybertsen, E. B. Stechel, M. Schluter, D. R. Jennison, Phys. Rev. B 41,
11068 (1990).

A. K. McMahan, R. M. Martin, and S. Satpathy, Phys. Rev. B 38, 6650 (1988).

H. Eskes, G. A. Sawatzky and L. F. Feiner, Physica C 160, 424 (1989).

V. Emery, Phys. Rev. Lett. 58, 2794 (1987).

C. M. Varma, S Schmitt-Rink, E. Abrahams, Solid State Comm. 62, 681 (1987).
Hyunju Chang, J. F. Harrison, T. A. Kaplan and S. D. Mahanti, Phys. Rev. B 49,
15753 (1994).

T. A. Kaplan, Hyunju Chang, S. D. Mahanti, J. F. Harrison, in Electronic
Properties of Solids using Cluster Methods, edited by T. A. Kaplan and S. D.
Mahanti, (Plenum Press, NY, 1995) p.73.

J. D. Reger, A. P. Young, Phys. Rev. B 37, 5978 (1989); M. Gross, E. Sanchez-
Velasco, E. Siggia, Phys. Rev. B 39, 2484 (1989).

L. Mitas, in Electronic Properties of Solids using Cluster Methods, edited by T.
A. Kaplan and S. D. Mahanti, (Plenum Press, NY, 1995), p.131.

R. B. Murphy, M. D. Beachy, R. A. Friesner, M. N. Ringnalda, J. Chem. Phys.
103, 1481 (1995).

W. Pickett, Rev. Mod. Phys. 61, 433 (1989).

These are the same parameter values as those given in Refs. 12 and 16 except that
the hopping parameters t are the negatives of those given there. We thank Dr.
Ellen Stechel (private communication) for verifying that these negative values are
the values actually used in Ref. 16.

P. O. Lowdin, J. Math. Phys. 3, 969 (1962).

113

[29]
[30]

[31]
[32]

[33]

[34]
[35]

[36]

[3 7]
[3 8]
[39]
[40]
[41]

I. Lindgren and J. Morrison, Atomic Many-Body Theory, Springer Series in
Chemical Physics, Vol. 13 (Springer-Verlag, NY 1982), p.207 ff.

A similar calculation was carried out for a Hamiltonian less general than Eq. (2.1),
by H. Eskes and J. H. Jefferson, Phys. Rev. B 48, 9788 (1993).

T. A. Kaplan and P. N. Argyres, Phys. Rev. B 11, 1775 (1975).

T. A. Kaplan and R. A. Bari, J. Appl. Phys. 41, 845 (1970); T. A. Kaplan, AIP
Conf. Proc. No.5, p.1305 (1972).

W. England, L. S. Salmon, K. Ruedenberg, Molecular Orbitals (Springer-Verlag,
New York, 1971) and H. Weinstein, R. Pauncz, M. Cohen, Advance in Atomic &
Molecular Physics, Vol. 7 (Academic Press, New York, 1971).

The same contributions were found in Ref. 13 for less general parameters in H.

There is an ambiguity for the exchange terms, Kﬁclgciockcja, because they can
also be expressed in terms of number operators when 0' = 0". For convenience,
we take all the exchange terms as a part of the perturbation. In fact, due to the
small values of the exchange parameters Ki]; this convention hardly affects the
convergence rate.

F. C. Zhang and T. M. Rice, Phys. Rev. B 37, 3757 (1988). These authors worked
within straightforward p.t..

W. E. Pickett, Rev. Mod. Phys. 61, 433 (1989).

G. Y. Guo and W. M. Temmerman, J. Phys. C: Solid State Phys. 21, L803 (1988).
A. Svane, Phys. Rev. Lett. 68, 1900 (1992).

P. Wei and Z. Q. Qi, Phys. Rev. B 49, 12159 (1994).

T. A. Kaplan, H. Chang, S. D. Mahanti, and J. F. Harrison, in Electronic
Properties of Solids Using Cluster Methods, edited by T. A. Kaplan and S. D.
Mahanti (Plenum, New York, 1995), pg. 73.

T. P. Das, in Electronic Properties of Solids Using Cluster Methods (see Ref. 41),
pg. 1.

R. L. Martin and F. Illas, Phys. Rev. Lett. 79, 1539 ( 1997).

F. Illas, I. de P. R. Moreira, C. de Graaf, O. Castell and J. Casanovas, Phys. Rev.
B 56, 5069 (1997).

I. de P. R. Moreira and F. Illas, Phys. Rev. B 55, 4129 (1997).

J. Casanovas, J. Rubio and F. Illas, Phys. Rev. B 53, 945 (1996).

R. Dovesi, V. R. Saunders, C. Roetti, M. Causa, N. M. Harrison, R. Orlando, E.
Apra, CR YSTAL95 User ’s Manual, University of Torino, Torino, 1996.

M. D. Towler, N. L. Allan, N. M. Harrison, V. R. Saunders, W. C. Mackrodt, and
E. Apra, Phys. Rev. B 50, 5041 (1994).

J. M. Ricart, R. Dovesi, C. Roetti, and V. R. Saunders, Phys. Rev. B 52, 2381
(1995)

M. D. Towler, R. Dovesi, and V. R. Saunders, Phys. Rev. B 52, 10150 (1995).

S. Massidda, M. Posternak, and A. Baldereschi, Phys. Rev. B 46, 11705 (1992).
R. W. G. Wyckoff, Crystal Strucuters, 2“d ed. (Interscience Publishers, New York,
1963)

The ofﬁcial web site of the CRYSTAL program is

http://www. dl. ac. uk/ T CS/Software/ CR Y S T AL,

114

[54]
[55]
[56]

[57]
[53]

[59]
[60]
[61]
[62]
[63]

[64]

[65]

[66]
[67]

[69]

[70]

[71]

which contains a collection of basis sets for various atoms suitable for crystal
calculations.

P. J. Hay and W. R. Wadt, J. Chem. Phys. 82, 270 (1985).

N. W. Winter and R. M. Pitzer, J. Chem. Phys. 89, 446 (1998).

P. W. Anderson, Solid State Physics 14, 99 (1963). The mapping scheme is
mentioned on pg. 166. Note also Anderson’s discussion of Kondo’s work on this
subject, on the same page.

L. Noodleman, J. Chem. Phys. 74, 5737 (1981).

S. Chakravaty, in High-Temperature Superconductivity, edited by K. S. Bedell, D.
Coffey, D. E. Meltzer, D. Pines, J. R. Schrieffer (Addison-Wesley, New York,
1990)

G. Aeppli, and D. J. Buttrey, Phys. Rev. Lett. 61, 203 (1988).

S. Itoh, K. Yamada, M. Arai, Y. Endoh, Y. Hidaka, and S. Hosoya, J. Phys. Soc.
Jpn. (Japan) 63, 4542 (1994).

G. Aeppli, S. M. Hayden, H. A. Mook, Z. Fisk, S-W. Cheong, D. Rytz, J. P.
Remeika, G. P. Espinosa, and A. S. Cooper, Phys. Rev. Lett. 62, 2052 (1989).

P. E. Sulewski, P. A. Fleury, K. B. Lyons, S-W. Cheong, and Z. Fisk, Phys. Rev.
B 41, 225 (1990).

The parameters in the effective l-band Hubbard models for La,CuO4 and LaQNiO4
are similar [A. K. McMahan, private communication].

X.-L. Wang, C. Stassis, D. C. Johnston, T. C. Leung, J. Ye, B. N. Harmon, G. H.
Lander, A. J. Schultz, C.—K. Loong, and J. M. Honig, Phys. Rev. B 45, 5645
(1992)

S. Shamoto, M. Sato, J. M. Tranquada, B. J. Sternlieb, and G. Shirane, Phys. Rev.
B 48, 13817 (1993).

M. T. Hutchings and H. J. Guggenheim, J. Phys. C 3, 1303 (1970).

P. W. Anderson, Phys. Rev. 86, 694 (1952). For more recent calculations
supporting Anderson’s original results see J. D. Deger and A. P. Young, Phys.
Rev. B 37, 5978 (1988); M. Gross, E. Sanchez-Velasco and E. Siggia, Phys. Rev.
B 39, 2484 (1989); G. Castilla and S. Chakravarty, Bull. Am. Phys. Soc. 37, 447
(1992)

A form factor measurement was made on LaQNiO,125 by J. M. Tranquada et al.,
Phys. Rev. B 52, 3581 (1995). The result showed a similar ﬂatness at small k,
being somewhat smaller than the stoichiometric case. This is large doping; the
spin ordering is radically changed.

In the face of this major overall disagreement, we investigated an exotic ground
state based on mixing of the high- and low-spin states on the Ni”, that could
occur if the tetragonal structure caused a sufﬁciently large splitting between the
two eg states [41]. Unfortunately, the state obtained did not correct the form factor
at the important, small-k, Bragg peaks [T. A. Kaplan and S. D. Mahanti,
unpublished]; this was actually consistent with the small eg splitting found
theoretically [A. K. McMahan, private communication].

T. Freltoﬁ, G. Shirane, S. Mitsuda, J. P. Remeika, and A. S. Cooper, Phys. Rev. B
37,137(1988)

T. A. Kaplan and S. D. Mahanti, J. Appl. Phys. 69(8), 5382 (1991).

115

[72]

[73]

[74]
[75]

[76]
[77]

[78]

[79]

[80]

[81]
[82]
[83]

K. Yamada, E. Kudo, Y. Endoh, Y. Hidaka, M. Oda, M. Suzuki, and T.
Murakami, Solid State Comm. 64, 753 (1987).

The arguments of Refs. 71 and 41 are empirical. A plausible rationale for the
conclusion is the following. If there is oxygen deﬁciency, there should be some
Cu“, which simply dilutes the localized Heisenberg spins with nonmagnetic sites,
leading to a reduction in the AFM moment and TN. On the other hand, excess
oxygen goes in interstitial positions located within the La,O2 bilayers as in the
nickelate [D. E. Rice and D. J. Buttrey, J. Solid State Chem. 105, 197 (1993)];
rather than causing Cu”, it causes 0'1 in the CuO2 planes, as is known for
YBa,Cu,Oé,5. But such an ion between two Cu"2 causes a ferromagnetic
interaction [T. A. Kaplan, Nuclear Phys. B (Proc. Suppl.) 5A, 151 (1988); A.
Aharony et al., Phys. Rev. Lett. 60, 1330 (1988)], so this tends again to make the
AF M moment and TN decrease on adding oxygen to the stoichiometric material.
B. Keimer, A. Aharony, A. Auerbach, B. J. Birgeneau, A. Cassanho, Y. Endoh, R.
W. Erwin, M. A. Kastner, and G. Shirane, Phys. Rev. B 45, 7430 (1992).

In Ref. 41, the spin reduction for La2CuO4 is larger than this 15.8%. We believe
the main source of the difference lies in the different deﬁnitions used.

J. Hubbard and W. Marshall, Proc. Phys. Soc. 86, 561 (1965).

This is qualitatively consistent with parameters in the multi-band models as
calculated by A. K. McMahan using constrained local density calculations. He
ﬁnds (see pg. 157 of the book in Ref. 41) the charge-transfer energy epd is
appreciably larger in the nickelate than the cuprate, while the hopping parameter
tpd is very similar in the two cases [private cormnunication].

Some of the present authors had accepted this: T. A. Kaplan, S. D. Mahanti and
H. Chang, Phys. Rev. B 45, 2565 (1992).

Comparison with experiment of the “magnetic moment” In calculated in the HFA
for MnO and NiO was made in the course of calculations similar to ours, by M. D.
Towler et al. [48]. We point out that their discrepancies (Inth > mcxp) would be
reduced if the quantum spin ﬂuctuations were taken into account (as they should
be, in principle, according to our argument above). Also, replacing their use of
the Mulliken spin population by the Hubbard-Marshall (HM) extrapolation [76]
would go in the direction of decreasing the theoretical moment. That the HM
theory is more appropriate for comparison with neutron difﬁaction results is
discussed in S. D. Mahanti et al., J. Appl. Phys. 73 (10), 6105 (1993). The best
procedure for comparing theory and experiment is to directly consider the
measured quantities, the intensities of the Bragg peaks (as we have done above),
avoiding perhaps controversial questions of how to deﬁne the antiferromagnetic
moment.

S. L. Cooper, G. A. Thomas, A. J. Millis, P. E. Sulewski, J. Orenstein, D. H.
Rapkine, S-W. Cheong, and P. L. Trevor, Phys. Rev. B 42, 10785 ( 1990).

D. Adler, Solid State Physics 21, 1 (1968).

J. C. Slater, Phys. Rev. 82, 538 (1951).

T. A. Kaplan, S. D. Mahanti, Y.-S. Su and J. F. Harrison, unpublished.

116

[84]

[85]
[86]
[87]
[88]
[89]
[90]
[91]

[92]
[93]

[94]

[95]

[96]

[97]
[98]
[99]
[100]

[101]

[102]
[103]
[104]

[105]
[106]

[107]

S-W. Cheong and C. H. Chen, in Colossal Magnetoresistance, Charge Ordering
and Related Properties of Manganese Oxides, edited by C. N. R. Rao and B.
Raveau (World Scientiﬁc, 1998).

T. Saitoh, A. E. Boucquet, T. Mizokawa, H. Namatame, A. Fujimori, M. Abbate,
Y. Takeda, M. Takano, Phys. Rev. B 51, 13942 (1995)

D. D. Sarrna, N. Shanthi, S. R. Barman, N. Hamada, H. Sawada, K. Terakura,
Phys. Rev. Lett. 75, 1126 (1995)

S. Satpathy, Z. S. Popovié, F. R. Vukajlovié , Phys. Rev. Lett. 76, 960 (1996).

W. E. Pickett and D. Singh, Phys. Rev. B 53, 1146 (1996).

I. Solovyev, N. Hamada, K. Terakura, Phys. Rev. Lett. 76, 4825 (1996).

K. Hirota, N. Kaneko, A. Nishizawa, Y. Endoh, J. Phys. Soc. Jpn. 65, 3736
( 1996).

J. Zaanen, G. Sawatzky, J. W. Allen, Phys. Rev. Lett. 55, 418 (1985); Phys. Rev.
Lett. 53 , 2339 (1984); J. Zaanen, G. Sawatzky, J. Solid State Chem. 88, 8 (1990).

W. E. Pickett, Rev. Mod. Phys. 61, 433 (1996)

Y-S. Su, T. A. Kaplan, S. D. Mahanti, J. F. Harrison, Phys. Rev. B 59, 10521
(1999)

W. E. Pickett and D. J. Singh, in Physics of Manganites, edited by T. A. Kaplan
and S. D. Mahanti (Kluwer Academic / Plenium Publishers, New York, 1999),
p.87.

J. B. A. A. Elemans, B. van Laar, K. R. van der Veen, and B. O. Loopstra, J. Solid
State Chem. 3, 238 (1971).

In the ana space group deﬁned in the International Tables For X-Ray
Crystallography, vol. I, edited by N. F. M. Henry and K. Lonsdale (The Kynoch
Press, Birmingham, England), the longest side is along the b-axis.

There have been HFA calculations on model systems, e.g., T. Mizokawa and A.
Fujimori, Phys. Rev. B 51, 12880 (1995).

R. Dovesi, V. R. Saunders, C. Roetti, M. Causa, N. M. Harrison, R. Orlando, E.
Apra, CR YSTAL95 User ’3 Manual, University of Torino, Torino, 1996.

F. F. Fava, Ph D’Arco, R. Orlando and R. Dovesi, J. Phys.: Condens. Matter 9,
489(1997)

P. J. Hay and W. R. Wadt, J. Chem. Phys. 82, 270 (1985).

G. H. J onker and J. H. van Santen, Physica (Amsterdam) 16, 337 (1950); J. H. van
Santen and G. H. Jonker, Physica (Amsterdam) 16, 599 (1950); E. O. Wollan and
W. C. Koehler, Phys. Rev. 100, 545 (1955).

M. Takahashi, J. Phys. C: Solid State Phys. 10, 1289 (1977).

T. A. Kaplan, H. Chang, S. D. Mahanti, J. F. Harrison, in Electronic Properties of
Solids Using Cluster Methods, edited by T. A. Kaplan and S. D. Mahanti
(Plenum, New York, 1995), pg. 73.

P. S. Bagus, F. Illas, C. Sousa, G. Pacchioni, in Electronic Properties of Solids
Using Cluster Methods (see Ref. 103), pg. 93.

J. B. Goodenough, Phys. Rev. 100, 564 (1955).

Y. Murakami, J. P. Hill, D. Gibbs, M. Blume, I. Koyama, M. Tanaka, H. Kawata,
T. Arima, Y. Tokura, K. Hirota, and Y. Endoh, Phys. Rev. Lett. 81, 582 (1998).

T. Mizokawa and A. Fujimori, see Ref. 97.

117

[108]
[109]
[110]
[111]

[112]
[113]

[114]

[115]
[116]
[117]
[118]

[119]
[120]

[121]

C. Zener, Phys. Rev. 82, 403 (1951).

T. A. Kaplan and S. D. Mahanti, in Physics of Manganites, edited by T. A. Kaplan
and S. D. Mahanti (Kluwer Academic / Plenum Publishers, New York, 1999), pg.
135.

T. G. Perring, G. Aeppli, S. M. Hayden, S. A. Carter, J. P. Remeika, and S-W.
Cheong, Phys. Rev. Lett. 77, 711 (1996).

P. Dai, H. Y. Hwang, J. Zhang, J. A. Femandez-Baca, S-W. Cheong, C. Kloc, Y.
Tomioka, and Y. Tokura, cond-mat / 9904372.

H. Meskine and S. Satpathy, J. of Appl. Phys. 85, 4346 (1999).

E. Dagotto, S. Yunoki, and A. Moreo, in Physics of Manganites (see Ref. 109),
pg. 39.

S. Yunoki et al., cond-mat / 9909254, have done calculations using both non-
cooperative and cooperative Jahn-Teller phonon models and found no essential
difference in the results for some static properties such as the charge, spin, and
orbital orderings. However, whether this is also true for dynamic properties like
the spin wave dispersion remains to be seen.

G. Khaliullin and R. Kilian, cond-mat / 9904316.

J. Hubbard, D. E. Rimmer, and F. R. A. Hopgood, Proc. Phys. Soc. 88, 13 (1966).
D. E. Rimmer, J. Phys. C: Solid State Phys. 2, 329 (1969).

A. B. van Oosten, R. Broer and W. C. Nieuwpoort, Chem. Phys. Lett. 257, 207
(1996)

F. Illas and R. L. Martin, J. Chem. Phys. 108 (6), 2519 ( 1998).

R. Caballol, O. Castell, F. Illas, I. de P. R. Moreira and J. P. Malrieu, J. Phys.
Chem. A 101, 7860 (1997).

In an earlier work, Heisenberg exchange constants (J) are also obtained by making
a correspondence between Kohn—Sham determinants and product spin states, T.
Oguchi, K. Terakura and A. R. Williams, Phys. Rev. B 28, 6443 (1983), with no
mention of symmetry-breaking.

118

 

 

 

.r. . . . . . ..:.:

1.2... .3 f

5;...

. 7 .
: 3.7.1.. .

1:13. .

.7. .

ii .2. .
.1. .r .

 

25

l
2
O
3
9
2

7 .
1:5... . . . ;
‘14:}. . . . .. . :

 

:
.. . . .1715
.C............. , 5.2.3.
. j: ..r.,:. 3...:
.. .... .
.s

._u...

3......
Z. .

,5"

'1! I"!

 

 

. . . . ... r.
. a. . 3 a... .;
z: .. . . . . . ......:...r.._....: