v... w.
t... :9...

vel- 11.1.1}.
‘L‘.t his!)

 

u

x:

.1:

.l

L

 

T155332?
”and: .d‘bnihfi

.1 and:

‘9:
’ .3.
:1; . an. .hw‘
anfiL} Q
..
‘33:"5. no I l

.1»: «9a..

:5
4 .117:
5»

I‘ll . I...»
gainflfilt‘.
1... ..

L.I .- k .

.>..‘ .ss,,.=
‘| v‘ ‘Ivniui‘ A
.qsggv

fi,. ‘

U3)!

{xiiin v.
I 83.12313. ..

.u

5“; $1 .1.
r?«n«:xf?.
.13.}...15; .S
.‘v)
:5.
VJ!

$1253.11 .
y k

:whfi

i, 3:...
:3? s

‘05.?!
ax~x

x. :1 l . 4
.t. 1 $13!.”

 

.{xz

 

LIBRARY
Michigan Stat;-
University

-Mr

 

 

This is to certify that the

dissertation entitled

Critical Behavior of a Class of Quantum Disordered
Systems at T = 0 and Finite Temperature Magnetization
Studies of Snall Magnetic Clusters

presented by

VIK'IOR Z . CEROVSKI

has been accepted towards fulfillment
of the requirements for

Ph. D . degree in .Bhysms_'

 

Major professor

S.D.. Mahanti

Date 08/21/01

MS U is an Affirmative 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 C'jClRC/DataDUOpGS-p. 15

 

III} I ‘I‘lllv. .I'\

 

OF CRITICAL BEHAVIOR OF A CLASS OF QUANTUM
DISCRDERED SYSTEMS AT T = 0 AND FINITE
TEMPERATURE MAGNETIZATION STUDIES OF SMALL
MAGNETIC CLUSTERS
By

Viktor Z. Cerovski

AN ABSTRACT OF A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy
2001

Prof. S.D. Mahanti

ABSTRACT

CRITICAL BEHAVIOR OF A CLASS OF QUANTUM
DISCRDERED SYSTEMS AT T = 0 AND FINITE
TEMPERATURE MAGNETIZATION STUDIES OF SMALL
MAGNETIC CLUSTERS

By

Viktor Z. Cerovski

Critical properties of several models from the class of quantum disordered systems
with chiral symmetry on a two-dimensional square lattice are studied using numerical
simulations. The first one, the Anderson bond disordered model (ABD), is known
to have a logarithmically divergent density of states (DOS) and the parameter deter-
mining this divergence is calculated for the first time in this dissertation. It is also
shown that the known asymptotic expression of DOS is inadequate for describing of
the calcualted part of the infinite system size DOS and a new asymptotic expression
is proposed. Multifractal analysis associated with the scaling of inverse participation
numbers is carried out for eigenstates near the band center and it is found that scaling
at the band center is described by a significantly different set of generalized dimen-
sions than of those at nearby energies. This discontinuity at the band center is shown
to be a pr0perty of the two central states of the spectrum. It is known from previous
studies that ABD model has a critical point at the band center, and this dissertation

studies the corresponding divergence of the correlation length, which is Shown to be

accurately described by a power-law, and the critical exponent is calculated for the
first time. This calculation is done using a standard transfer-matrix method (TMM)
as well as a new finite-size scaling analysis (FSS) related to the two central states of
the spectrum, and both methods give the same exponent within error bars. Another
model in the same class is the random flux model (RFM), describing an electron
moving on a square lattice in the presence of random magnetic flux through each
plaquette. RFM has attracted a lot of interest and previous studies have proposed
several mutually incompatible solutions to the problem of its scaling properties. The
model is studied using TMM, and the high-precision of calculations allowed us to con-
clude that scaling of the correlation length near the band center has similar behavior
to that of the ABD model. Therefore the FSS related to the two central states of the
spectrum is applied to obtain the critical exponent. The differences between ABD
and RFM are also discussed, the main one being a broken additional 0(2) symmetry
at the critical point of the RPM.

Another problem studied in this dissertation is the temperature (T) dependence
of the magnetic properties of a small classical magnetic cluster of Gd13. It was known
from previous T = 0 studies that the cluster has a hexagonal close-packed structure
and its magnetic properties could be described by a classical Heisenberg model with
competing ferromagnetic (among the nearest) and antiferromagnetic (among the next-
nearest neighbors) interactions. This model has been used here to probe finite T
properties of the cluster. Magnetization is calculated from the parallel tempering
Monte-Carlo numerical simulations, and it exhibits an anomalous increase with T
in a certain range of T and for particular ratios of fero- to antiferro—couplings. Our

simulation results are in good agreement with experimental findings.

ACKNOWLEDGMENTS

I am very pleased to acknowledge my gratefulness to, first and foremost, my ad-
visor prof. S. D. Mahanti for his guidance through vicissitudes of my graduate studies
—— from him I have learned what the physics is all about; to prof. T.A. K aplan for his
help — his enthusiasm and knowledge will be a lasting inspiration in my future physics
research; to prof. M. Thorpe for providing computer time; to prof. S. Billinge and his
group; to prof. V. Zelevinskii, prof. M.I. Dykman and prof. C. P. Yuan— their courses
I have enjoyed the most; to physicists who kindly found time to discuss various mat-
ters of importance to this dissertation: S.A. ngman, T.A. Kaplan, V. Zelevinskii,
A. Taraphder, I. Herbut, F.M. Izrailev, S. Khanna, R. Bhatt, D. Mulhall, N. Birge,
S. T essmer, B. Nikolic’, I. Smolyarenko, I. Mostovoy, E. Boiin, R. di Francesko,
J. Kondev, T. Nagy, M. Drndic’, V. Srdano'v, I. Dukovski; to always helpful people of
Department of Physics and Astronomy and National Superconducting Cyclotron Lab
here at Michigan State University for their support; to my good friend Declan for all
the good times; to one enormously cute mathematician, for things I would rather not
Specify; to J.L. Goddard, F. Fellini and A. Egoyan, who provided me with hundreds
of hours of quality television watching; and, at last but not at least, to F. Nietzsche,
M. Heidegger, and J. Baudrillard, who explained it all. This thesis is devoted to

whom it will mean the most, my mother.

iv

TABLE OF CONTENTS

LIST OF FIGURES vii
LIST OF TABLES ix
1 Overview 1
1.1 Strongly correlated systems ......................... 3
1.2 Disordered quantum systems with chiral symmetry ............ 4
1.3 Small magnetic clusters and mesoscopic physics .............. 6
2 Anderson localization and related phenomena 7
2.1 Models of disorder .............................. 8
2.2 Anderson localization ............................. 10
2.3 One-parameter scaling theory of localization ................ 12
2.3.1 Weak localization corrections ....................... 14
2.3.2 Scaling theory of localization: aftermath ................. 19
2.4 Universal conductance fluctuations and mesoscopic physics ........ 21
2.5 Random matrix theory and universality ................... 22
2.5.1 Quantum Chaos .............................. 25
3 Two dimensional chiral disordered systems on a square lattice 28
3.1 Chiral symmetry ............................... 29
3.2 Double-exchange model ........................... 31
3.2.1 Density of states .............................. 32
3.2.2 Scaling properties of eigenstates ...................... 34
3.2.3 Band center states ............................. 42
3.3 Anderson bond-disordered model ...................... 45
3.3.1 Density of states near the band center .................. 47
3.3.2 Distribution of the nearest-neighbor level spacings ............ 49
3.3.3 Multifractality of eigenstates ....................... 52
3.3.4 Localization length near the band center ................. 63
3.4 Random flux model .............................. 68
4 Monte-Carlo studies of cluster magnetization 78
4.1 Monte-Carlo method ............................. 82
4.1.1 Metropolis algorithm ............................ 84
4.1.2 Thermalization ............................... 85

4.1.3 Parallel tempering ............................. 86

4.2 Magnetization of Gd13 clusters ........................ 88
4.2.1 Results of the numerical simulations ................... 88
4.2.2 Results of experimental studies ...................... 90
5 Conclusions 94
BIBLIOGRAPHY 97
Bibliography 97

vi

LIST OF FIGURES

2.1 Schematic diagram showing the behavior of conductance, correlation
length and density of states near the mobility edge ..........
2.2 Scaling fl function of the dimensionless conductance ............
2.3 Quantum coherent paths in the semi—classical limit ............
2.4 A pair of paths contributing to the weak-localization corrections to the
conductivity ................................

14
16

18

2.5 Conductance fluctuations of a small two-dimensional disordered conductor 24

2.6 Crossover between random and chaotic spectra in d = 2 Anderson model.

3.1 Density of states for the two-dimensional DE, |DE|, and ABD models. . .
3.2 Participation ratio of ABD, DE and |DE| models .............
3.3 Dependence of the participation number of ABD, DE and |DE| models on
system size ................................
3.4 Enrgy depdendence of the exponent ,6 for the two-dimensional DE, |DE|,
and ABD models ..............................
3.5 Goodness of the power-law scaling of PN for the two-dimensional DE,
IDEI, and ABD models. .........................
3.6 Schematical representation of In PN vs ln N ................
3.7 Dependence of the participation number on systems size at the band center
of DE, |DE|, and ABD models. .....................
3.8 Density of states of the ABD model near the band center, for different
system sizes ................................

26

33
36

38

48

3.9 Extrapolated density of states of the ABD for an infinite size square lattice 50

3.10 Distribution of the nearest-neighbor level spacings of the ABD model . .
3.11 Dependence of the inverse participation numbers on the bin size for dif-
ferent values of q. .............................
3.12 Dependence of the inverse participation numbers on the bin size for dif-
ferent systems sizes L ...........................
3.13 Dependence of the average inverse participation numbers on system size .
3.14 Scaling exponent of the inverse participation numbers ...........
3.15 Multifractal spectrums D, at the band center and nearby energies . . . .
3.16 Dependence of the new length scale 5’ in the RPM model .........
3.17 Scaling of the largest inverse Lyapunov exponent of the ABD model . . .
3.18 Universal function A of the one-parameter scaling of the ABD model
3.19 Largest inverse Lyapunov exponent of the RPM model on a square lattice

vii

52

56

57
59
60
61
63
65
67
71

3.20 Scaling of the additional length scale é’ in the RFM model ........ 73
3.21 Scaling of the largest inverse Lyapunov exponent in the RF M model in

the presence of the on-site disorder ................... 75
3.22 Distribution of 8 individual eigenenergies closest to the band center in the

ABD and RFM models .......................... 76
4.1 Gd13 cluster with Dgh-Symmetric spin configuration ............ 81
4.2 Results of Monte-Carlo simulations for various numbers of spin updates . 87
4.3 Temperature dependence of magnetization of Gdlg cluster ........ 89

4.4 Distribution of magnetization of Gd13 cluster at different temperatures . 91

viii

CHAPTER 1

Overview

Field of condensed matter physics is concerned to a large extent with the pr0perties
of systems containing many interacting particles. Let us, for example, consider a
metallic bar, whose microscopic Hamiltonian is quite complicated and involves many
degrees of freedom. Yet the properties of the bar leading to its mechanical oscillations
are experimentally very well described by only one number, the Young modulus. From
a theorist’s viewpoint, at the microscopic level there is a crystalline lattice (with im-
purities and defects) of ions held together by the itinerant conduction electrons. Ions
and conduction electrons interact among themselves with the long-range Coulomb
interaction, and the very existence of the crystal is due to the quantum mechanical
nature of its constituents. And yet to describe the mechanical oscillations of a metal-
lic bar one can go quite far by distinguishing between important and unimportant
degrees of freedom. Assuming that the pipe is a uniformly distributed set of masses
(ions) connected by springs, and neglecting the degrees of freedom associated with
conduction electrons, one is able not only to get Young’s modulus expressed in terms
of the ion mass and the Spring constant, but also to understand the collective os-
cillations of all the ions in terms of simply behaving non-interacting quasi-particles,
phonons, that represent collective behavior of all the ions together.

As a system is cooled down to lower temperatures, the quantum mechanical nature

of the matter becomes more apparent. Let us therefore consider the ground state of
a many-particle quantum mechanical system, which fully describes the system at
zero temperature. In the absence of interactions, the ground state is a product of
single-particle eigenstates. AS the interaction is turned on and becomes stronger, the
ground state becomes more and more complicated superposition of states from the
Hilbert space of noninteracting particles. Since the number of states is exponentially
large in the number of particles in the system, it becomes generally impossible to
keep track of all of the required numbers describing exactly the ground-state. It
is nevertheless usually enough to calculate perturbative corrections to the ground
state of non-interacting system in order to quantitatively include the presence of the
interaction.

If the strength of the interaction among particles becomes very large, however, it
is natural to expect that the approximate description in terms of the initial particles
becomes less accurate. In particular, a phase transition may occur, when the ground
state changes dramatically its properties, and when it is often necessary to change
the description of the wave-function in terms of a different set of basis states, more
appropriate to the new phase. And yet, if one is interested only in the ground state
energy, the Kohn—Sham theorem [1] assures that its calculation requires only polyno-
mialy many operations in the number of particles, since one needs to know only the
electron density at various points in space rather than the whole wavefunction. This
therefore makes the calculation of the ground state energy tractable in a very general
case of quantum many-body systems.

At the level of understanding the transport properties of many-particle systems,
similar attempts have been made to give a general description, such as the Landau
theory of Fermi liquids [2], in which the cumulative effect of Coulomb interaction
among electrons in metals is quantitatively included by the renormalization of the

electron mass, after which one ends up with weakly-interacting quasi—particles with a

new, effective mass and a finite life-time. Landau’s conjecture was that the existence
of quasi-particles was a general feature of quantum many-body systems. Even in the
case when systems undergo phase transitions, for instance, one has to identify the
order parameter and construct Landau-Ginzburg functional [3], which then allows for
a quantitative study of the transition itself. The new quasi-particles in the ordered

phase then can be associated with the fluctuations of the order parameter.

1.1 Strongly correlated systems

Even though Landau theory already sets a very general framework for studying con-
densed matter systems, it has become clear in the last two decades that different ap-
proaches are needed in order to understand a variety of new and previously unexpected
phenomena, such as the fractional quantum Hall effect (FQHE) [4], high-temperature
superconductivity (HTc) [5], colossal magneto—resistance, or properties of quantum
systems in general whose size interpolates between microscopic and macroscopic. In
the case of FQHE, for instance, an approximate ground state wavefunction is an en-
tangled state involving very large number of electrons that cannot be adiabaticaly
mapped into a set of non-interacting quasi-particles [4].

These realizations led to the concept of “strongly correlated systems”. There is
no strict definition of strongly correlated systems, partly because the term is often
used to characterize systems (models) that are still often not very well understood.
Some important ingredients are nevertheless known: low (effective) dimensionality
of systems at low temperatures leads to strong effects due to quantum fluctuations
(such as quantum phase transitions); strong-coupling limits of various theories (where
perturbative and some non-perturbative approaches fail); strong coupling between
spin and charge degrees of freedom. The main realization underlying the necessity

of introducing the concept of strongly correlated systems is that we still don’t know

what are all the forms in which matter can appear in the nature.

The problem studied in this dissertation within the framework of “strongly cor-
related systems” is that of the ferromagnetic Kondo lattice model in the limit of
infinite ferromagnetic coupling (double exchange model). The model is of interest
in the context of colossal magneto-resistance. In this model there is a strong ferro-
magnetic coupling between local magnetic moments of localized electrons and spins
of the conduction electrons. Since the magnitude of local magnetic moments is large,
they are approximated by classical ones, which are then studied using the approach

described bellow.

1.2 Disordered quantum systems with chiral sym-
metry

The main part of this dissertation is concerned with the physics of disordered quantum
systems in two dimensions with chiral symmetry (to be explained below). Initial
interest in these systems came from a particular approach to the double exchange
model, when one approximates the dynamics of the itinerant electrons as one where
non-interacting electrons from the itinerant conduction band move in a given random
and frozen (quenched) configuration of local spins coming from the localized electrons.

One thus reduces the initial many-body problem to that of a single particle moving
on a square lattice with random hopping energies between sites that are nearest
neighbors. In this case there exists a symmetry such that all the eigenergies come in
opposite-energy pairs for every realization of disorder, and this symmetry is commonly
called the particle-hole or the chiral symmetry (Sec. 3).

The presence of disorder, on the other hand, significantly influences electronic

transport in low-dimensional systems through the phenomena of Anderson localiza-

tion [6], leading in certain cases to the complete absence of diffusion of electrons.
There are many excellent reviews of this subject, which are the basis for an overview
given in the Chapter 2. It turns out that symmetry plays very important role in dis-
ordered systems, much more important than, for instance, the actual distribution of
energies of impurities or whether these energies are uncorrelated or have a short-range
correlation. Symmetry is a basis for the identification of universality classes. Three
of these classes, corresponding to systems with time-reversal symmetry, without it,
and those with the random spin scattering, have been known for quite a while, but
only recently it has been recognized that other symmetry classes exist as well, chiral
symmetry being the defining feature for some of those.

This dissertation addresses properties of several models with chiral symmetry.
The first on is the double exchange model mentioned above. Next one the Anderson
model with random hopping energies, also called Anderson bond-disordered model
(ABD) on a two-dimensional (2d) square lattice. In particuler, I will be focussing on
the singularities of the density of states, the scaling of the eignestates and correlation
length in the neighbourhood of quantum critical points (which occur at the band
center) charactereistic of these disordered systems. A more realistic chiral disordered
model, appearing in certain effective theories of HTc and QHE of half-filled lowest
Landau level, is the random-flux model (RFM), describing a particle moving in a
random quenched magnetic flux. This problem has attracted a lot of interest and
has turned out to be a djflicult one to which numerous incompatible solutions have
been proposed. This dissertation addresses this problem in the limit of an infinite 2d
square lattice, partially using the insights gained from the study of the ABD model,
by proposing the solution of it and by pointing out sources of discrepancies among

previous numerical studies with contradicting conclusions.

1.3 Small magnetic clusters and mesoscopic
physics

Finally, the Chapter 5 of this dissertation is devoted to the study of small magnetic
clusters. Properties of matter change significantly when the system size is reduced.
Two regimes then become important: microscopic and mesoscopic. The first roughly
corresponds to the case of small magnetic clusters, up to several tens of atoms large.
Magnetic properties of a small cluster can be completely different than those exhibited
by the bulk of the same material, and this dissertation addresses the problem of
anomalous thermal behavior of 13 atom clusters of Gadolinium.

The problem is studied within the context of classical statistical physics and us-
ing the Monte—Carlo method. This method is one of the most important methods
in computational physics, because it allows accurate calculation of thermodynamic
properties of a given system. Since the only limitation is the size of the system stud-
ied, the method is very useful when all the other approaches fail, or as an independent
check of the validity of other approximate methods. It is often the method of choice
to probe the physical prOperties of strongly interacting systems.

Another important regime is the mesoscopic regime, when the number of atoms
can be thousands or more, and yet prOpertieS of this system are still quite different
from those of the bulk materials and those of small clusters. Electronic transport
properties of small disordered metallic grains and wires at low-temperatures is an
excellent example for a variety of novel phenomena due to the quantum interference
effects, such as magneto-fingerprints and the universal conductance fluctuations, as

discussed briefly in Sec. 2.4.

CHAPTER 2

Anderson localization and related

phenomena

When a material has crystalline structure, it is natural to use the translational sym-
metry to express the electron eigenstates in terms of Bloch waves, which simplify the
problem and lead to the band structure of electronic states. The crystalline structure
is, however, only an idealization and real materials have certain degree of disorder,
either in the form of impurities, or as a result of topological disorder due to the amor-
phous structure of a solid. The degree of disorder thus may vary from very small to
very large. Here we are primarily interested in the effects of static quenched disorder
due to, for instance, impurities.

In a metal, impurities act as scattering centers for electrons modifying the trans-
port properties of the metal. If the concentration of impurities is small, one may
expect that their effect can be treated perturbatively, through Boltzmann transport
formalism. The characteristic parameter is kpé, where kp is the Fermi wavelength and
8 the mean free path. So long as kpf >> 1 one expects the Boltzmann semi-classical
approach to be valid [8]. This, however, turns out to be incorrect, because of the ex-
istence of non-local quantum interference between two possible electron trajectories.

Such purely quantum-mechanical effects are not accounted for in a standard semi-

classical approach and lead to non-trivial corrections that must be included to obtain
correct transport properties. These are the well-known weak-localization corrections,
discussed in the Sec. 2.3.].

In semiconductors, on the other hand, energy levels of impurities may be in the
band gap. The corresponding electron state bound to an impurity is exponentially
decaying in space. Since impurities which give rise to states close in energy are ex-
pected to be spatially far apart, overlap of the corresponding states is exponentially
small for small concentration of impurities. But for sufliciently high concentrations,
one would expect strong overlap of different bound states and therefore metallic be-
havior. This reasoning is nevertheless invalid since the introduction of impurities
breaks the translational symmetry of the crystal, and with it the whole picture of
bands of plane-wave like quasi-particle states. The larger the impurity concentration,
more perturbed the band structure becomes, and one needs a different approach that

can take into account effects of these perturbations.

2.1 Models of disorder

One commonly used approach to probe the effects of disorder on Single particle prop-
erties is to treat electrons via the tight-binding approximation, in which one starts
from the expansion of electron wavefunctions in terms of the Wannier states corre-
sponding to orthogonalized one-electronic orbitals localized about different sites of a
crystal. In the approximation of one orbital per site, the quantum disordered system
is then described by the Hamiltonian

N

H = Zeiclc, + Z (ha-ck,- + tzjcl-ci) , (2.1)

i=1 (133')

where c,- is the electron annihilation operator for the state localized at site 1",- of a d-

dimensional lattice with N = Ld sites, Q’s are on-site and t,,j’s are hopping energies.

The symbol (...) denotes nearest neighbor hoppings only. Disorder is introduced
through randomly chosen on-site and/ or hOpping energies. An alternative approach

to treat disorder is to use the continuum description:
H=E+U26(F—R,)EE+V(F), (2.2)
m m

where Nimp is the number of impurities which are located at positions R4, distributed
as d-correlated white noise, i.e. their disorder average (V(7"'1 — 172)) = 76(7‘1 -— 7'3),
with U chosen such that '7 oc NimplUl2 —» const when Nimp —+ 00, U —> O.

The tight-binding model (2.1) has the advantage of being amenable to numerical
simulations and certain analytical studies. On the other hand, the continuum model
(2.2) is somewhat more realistic, since it directly describes an electron near the Fermi
level interacting with impurities, and is usually suitable for analytical studies in the
limit kpf >> 1.

To connect the two models, let us consider the Anderson model (AM), described
by (2.1) with t,-, ,- = t and e,- uniformly and independently distributed between —W/ 2
and W/ 2 on a 3D lattice [6]. The limit of weak disorder is t > W, when one can
approximately treat electrons as belonging to a band of free electrons scattering off
the impurities. The mean free path can be estimated as t’ = (SimpNimp)_la where
Simp is the scattering cross-section, estimated to be Simp = (b 6E/E)2, where (SE is
the fluctuation in energy due to disorder and b is the concentration of scatterers. For

b = N7”3 E ~ t, (SE ~ W, we have 6 ~ b(t/W)2. Assuming furthermore kp ~ 1/b,

Imp’
t 2
k€~ — . 2.
F (W) (3)

Finally, several models studied in this dissertation are all of the kind commonly

we get

referred to as the Anderson model with off-diagonal (or bond) disorder (ABD model),

described by Eq. (2.1) with e,- = 0. Since the on—site disorder is now absent, the

Hamiltonian takes the following form:

H = Z (ti’jCICj + tf,jc]c,-) . (2.4)
(131')

2.2 Anderson localization

When the kinetic energy is much smaller than the strength of the on-site disorder
measured by W (t < W), approximate eigenstates of the Anderson model to the
lowest order correspond to the electron being localized at individual sites of the lattice.
Kinetic energy then tends to delocalize electron states across the lattice. Anderson,
however, found that for sufficiently strong on-site disorder all eigenstates of the AM
remain localized in the infinite 3D system [6]. The localized eigenstate is described
by a wavefunction that decays exponentially fast at infinity, and the corresponding
length scale é describing the decay is called the localization length.

The main consequence of the Anderson localization phenomena is that, for sufli-
ciently strong disorder, non-interacting electron disordered systems become insulators
at temperature T = 0. With decreasing W/ t there is a transition to a metallic sys-
tem. The change from the insulator to the metal takes place through a quantum
phase transition. Quantum phase transition is a continuous phase transition that
occurs when some parameter of the Hamiltonian reaches a certain (critical) value.
Since it happens at temperature T = 0, there are no thermal fluctuations governing
the usual finite T phase transition in statistical physics. Instead, their role is taken
over by quantum fluctuations [14].

In the case of Anderson localization, according to the scaling theory of localiza-
tion [15] discussed in the next section, as the disorder strength is reduced, the AM
undergoes localization-delocalization quantum phase transition at W = Wc in d = 3
(energy scale is set by = 1). For W < We, a band of extended states appears around

the band center (for [E] < E), separated from the localized states by the mobility

10

 

 

 

\““~§
>
E

Figure 2.1: Schematic diagram showing the behavior of: conductance a, correlation
length 5 and density of states near the mobility edge of an infinite 3D system for
W < WC m 13.6.

edge EC(W), Exactly at the mobility edge (E 2 EC) there exists a critical state. The
transition is characterized by a divergent correlation length near the mobility edge,
where {(E) ~ |E — Eel—V. This is sumarised in Fig. 2.1.

Since the localization length is not a directly measurable quantity, one can alter-
natively describe the transition in terms of the transport properties of the system,
and such a scaling theory of this metal-insulator transition is described in Sec. 2.3,
the main consequences of which are that the dimensionless conductance is the only
relevant quantity describing the transition. In d = 1, 2 all states become localized in

the AM for any disorder strength W.

11

2.3 One-parameter scaling theory of localization

Abrahams, Anderson, Licciardello and Ramakrishnan [15] addressed the issue of the
Anderson transition. They studied change of the dimensionless conductance g as a
function of L, based on the Thouless approach [7, 16] to the quantum transport.

The Drude formula for conductivity 0 (equivalent to the Einstein relation) of a
diffusive electron with mass m and charge e is

Neezr

o = e2pD, (2.5)

where r is the mean scattering time, D is the diffusion constant, and p is the density
of states, p(E) = Z, 6(E — E.)
The conductance G E 00’“2 can be expressed using (2.5) as

33132:

G hA’

(2.6)

where A = (pLd)‘1 is the mean level spacing and ET}, = hD/L2 is the Thouless
energy. The latter can be understood as the inverse time that takes for a particle to
diffuse throughout the sample. Due to the relation G = (e2/h)g, one then arrives at

the relation

_ ETh
g— A . (2.7)

Dimensionless conductance g(E) is thus equal to the number of energy levels inside
the energy interval of width En, around a given energy E. Abrahams et. al [15]
estimated ET}, to be the shift 6V in eigenenergy when the boundary conditions are
changed from periodic to anti-periodic, i. e g = 6V/A, and studied the scaling behavior
of g as the system Size L is (exponentially) increased.

In the localized regime, perturbation at the boundaries will make 6 V an exponen-
tially small quantity because eigenstates are exponentially decaying with distance,

and therefore
g(L) 0< eXI)(--L/€)- (2-8)

12

In the metallic regime, on the other hand, one expects Ohm’s law to be valid and

therefore:

g(L) = oLd‘z. (2.9)

One can now attempt to study the rate of change of 9 when two pieces of material of
Size Ld are joined together. The quantity of interest is the logarithmic derivative, fl,

defined as
= dlng
_ dlnL’

 

4 (2.10)

which measures the rate of change of 9 when L changes exponentialy.
On the other hand, 5 can be directly calculated from (2.8) in the localized regime

to be

3 = ln(9/gc), (2.11)

where 9C is a characteristic conductance that separates localized and metallic regimes
9 << go and 9 >> 96 respectively (and therefore fl-function in (2.11) is negative). It

follows from (2.9) that [3 in the metallic regime is
D = d — 2, (2.12)

Abrahams et. at then proceeded to construct the whole fi-function starting from
the asymptotic expressions (2.11) and (2.12). Since only finite-size systems are con-
sidered, one can expect that two limiting cases are connected analytically. Assuming
furthermore that no parameter other than g is relevant, i.e that )6 E C(g), they
obtained the functional form of [3(9) for d = 1, 2, 3 shown in Fig. 2.2.

Since ,8 < 0 in d = 1 and 2, dimensionless conductance g decreases as L increases,
and all states become localized when L -—> oo (arrows in the figure indicate the flow
of In g as L grows exponentially fast). In d = 3 there is a critical point C(gc) = O

(mobility edge) and a continuous phase transition into the metallic regime (where

[3(9) > 0).

13

 

 

 

 

Figure 2.2: Scaling function ,8(g) E d In g/ d In L for d = 1,2, and 3 of Abrahams etal
[15].

The critical exponent V of the Anderson transition can be obtained from the
linear expansion of fl-function around the critical point, C(ln g — In gc) z )6(ln g.) +
D’ (In gc)(ln g — 1n gc), from which one obtains the critical exponent l/ of the transition

to be 1/1/ = E’ (In gc), which is the Slope of the dotted line in Fig. 2.2.

2.3.1 Weak localization corrections

According to the scaling theory of localization, all states in d = 2 disordered system
become localized as L —» 00 at T = O. This result is particularly surprising since
the relation (2.3) gives a simple qualitative understanding of the Anderson transition
in d = 3 as being due to the violation of kpf >> 1, that is due to the enhanced
interference effects of wave-like properties of the electron wavefunction, known as the

Ioffe-Regel criterion [8].

14

To see the effect of quantum interference even in the limit kpf >> 1, let us recon-
sider electron propagation in the potential of impurities as described by Eq. (2.2).
There, the semi-classical Boltzmann approach is expected to be valid, and one, as
a first approximation, can start from the system in which the classical electron of
charge e and mass m performs diffusive motion across the sample. The correspond-
ing probability distribution P(7"', t) to find the electron at a position I" at the time t,

starting from the initial distribution P(F,t = 0) = 6(7") is

_, 1 F2
P(T',t) = Wexp (—fl)—t) . (2.13)

This corresponds to a Gaussian distribution, whose width is slowly increasing (oc
x/th) characteristic of the diffusion process.

Let us now consider the quantum corrections to o. This will be first done using
qualitative arguments from a path-integral representation of the electron propagator
in the semi-classical limit. This approach is not going to give here any quantative re-
sults, but instead only an insight into the origin of quantum interference contributions
to the conductivity.

To that end we can start with the Feyman path integral expression for the electron
propagator [12]

FB(t=tB)=FB

DF(t) exp {-2- fttB dt’ (E2/2m — V(F(t')))} ,

C(FBatBifAit/l) =/
A
(2.14)

r74(t=tA)=r“A
where V(7") is the potential of impurities of Hamiltonian (2.2), and ’DF(t) represents
the integration over all classical trajectories F(t). The number of these trajectories is
greatly reduced by taking into account the semiclassical condition hp! = 27r€//\p >> 1,
which constrains them to only those trajectories that belong to the tubes of width AF
around the lines connecting the impurities. Furthermore, paths undergoing scatter-
ing off two completely different sets of impurities are expected to have uncorrelated

phases, and thus to interfere incoherently.

15

 

Figure 2.3: Two types of pair of paths that interfere coherently in the path-integral
representation of the electron’s Green’s function (2.14) in the semi-classical limit
kpf >> 1.

The paths that will retain their quantum coherence to the largest extent are those
corresponding to electron trajectories undergoing the scattering off the same sequence
of impurities, such as the case (i) presented in Fig. 2.3. The second possibility is (ii)
the interference between a path and a close time—reversed path.

Neither of these two contributions alone, however, is expected to change the prop-
agator significantly, i.e. beyond the order (kp€)‘1. On the other hand, a is related
to a two-particle Green’s function and the localization can be understood to a great
extent at the level of correlation functions. To see this, let us look now at the product
of two Green’s functions (2.14) for a given configuration of impurities. In order to
estimate quantum corrections to the conductivity, we will again restrict ourselves to
the paths of the form (i) and (ii).

The new interference effect, not present at the level of a Single particle Green’s
function, is shown in the Fig. 2.4. Since we study the product of two Green’s functions,

the integrand in the exponent of (2.14) is the sum of contributions from two paths,

16

one from each of the two Green’s functions. We see that whenever two paths are close
to each another, integration along them makes that contributions from two distant
points (marked by the arrows in the figure) add up as the integration is performed for
different t’ . It is these non-local quantum interference effect that affects substantially
transport properties of an electron.

In a much more rigorous approach of the perturbation theory over impurity di-
agrams [12], this is seen in the following way: The Single particle disorder-averaged
Green’s function turns out to decay exponentially fast at a distance of Z. The disorder-
averaged two-particles Green’s function is, in the lowest approximation, the product
of two disorder-averaged Single-particle Green’s functions, which gives the Drude ex-
pression for o from the Kubo—Greenwood formula. The first quantum correction,
however, is obtained from a resummation of an infinite sequence of two—particle dia-
grams involving repeated scattering off the same impurity with electron pr0pagation
described by the disorder-averaged single-electron Green’s function.

We can now estimate quantum contributions from the self-intersecting paths from
Fig. 2.4 to D (and therefore to conductivity due to the relation (2.5)) by considering
the probability that the diffusive electron with probability distribution (2.13) will be

at a distance AF or closer to the initial point after time t:

TO
1’3- = 5’59 ~ — / va%P(t)dt. (2.15)

The lower cutoff is set to 7' since at smaller times there is no diffusion, the upper
cutoff is To, to be discussed below, 1)}? is the Fermi velocity, and P(t) is the return

probability P(f’ = O, t) = (Dt)“’/2. The estimate (2.15) now can be easily evaluated:

 

do 6D 1 2 _ _
a z’D‘N—hmm-dag .142 d)’ (2'16)

where pd is density of states per unit volume and L0 = \/ Dro, the cutoff length scale

corresponding to To, that will also be discussed below. The d = 2 case can be handled

17

 

Figure 2.4: A pair of paths from the product of two single-electron Green’s functions
for a Single configuration of impurities that gives quantum interference contribution
not present at the level of the single-electron Green’s function. Integration of the
exponent of the product includes coherent contributions from two distant points at
integration time t’ (marked by arrows).

by setting 6 = d — 2 to be small and then expanding 93‘ = exp(eln:r) z 1 + elna:
followed by taking 6 —> 0, which, together with relation (2.5), gives an estimate of

the quantum corrections to o as:

—;—h(Lo—£) d=1,
60 ~ -5311} L; d = 2, (2.17)
82
‘275 (Ii; — i) d = 3-

We see that the quantum corrections to conductivity are very large in one- and
two-dimensional systems, diverging in the infinite size system, which thus semi-

quantitatively accounts for the Anderson localization.

18

Upper cutoff L0

While the lower cutoff of the estimate (2.15) was naturally set by the microscopic
length scale 3, the upper cutoff L0 has not been determined yet. This corresponds
to a certain macroscopic length scale, or, correspondingly, to electron dynamics at
some large time scale 7'0. Dependence on this length scale is particularly strong in
the d = 2 case of the weak-localization corrections (2.17). One natural choice for L0
is the system size L. The obtained result from Eq. (2.17) is then easily reformulated
as C(g) = —(7r2g)‘1, which can be interpreted as a rigorously derived result of the
d = 2 case of the scaling theory of localization of Abrahams et. al.

Another possibility is to consider L0 to be the dephasing length L¢ due to var-
ious inelastic processes that destroy quantum coherence. These occur due to the
quasi-elastic scatterings that are consequence of the finite temperature, or due to
the electron-electron interaction, for instance, when electrons can be treated as new

effective scatterer other than impurities [12].

2.3.2 Scaling theory of localization: aftermath

The scaling theory of localization of Abrahams et. al has been formulated based
on a rather general arguments, and since it has many important and unexpected
consequences (as discussed in the previous section), it generated a lot of interest and
was checked by more rigorous approaches.

The hypothesis that change of the g(L) depends only on g has been checked by
several authors using the impurity diagram technique, which is based on the expansion
in (kp€)‘1. Altshuler et al. [8] found that corrections to conductivity in the order
(II:F€)‘2 did not give contribution to g, which thus furthermore strengthened the
idea that kpt’ is irrelevant for the scaling of g. 9 now has the meaning of the average

conductance over all realizations of disorder, g E (g). The asymptotic expression (2.8)

19

now can be more explicitly written as

ln(g(L)) o< —L/§. (2.18)

Kramer and MacKinnon [17, 18] were able to compute numerically [3(9) by study-
ing how localization length of long quasi-1D systems of cross-section L"'1 and length
>> L changes when L grows exponentially fast. This method is also used in this
dissertation, and explained in more details in Sec. 3.3.4.

Although this seemed already enough to establish the validity of one-parameter
scaling, discovery of universal conductance fluctuations [9], showed that g of small
systems is not a self-averaging quantity, and therefore one must address the question
of the whole distribution of g over the ensemble of configurations of disorder. This
in turn seemed to invalidate the basic premise of the one-parameter scaling, which
addresses the behavior only of the first moment of the distribution of 9. To address
this question, Altshuler et. al [19, 20] have developed an extension of the non-linear
o—model (NLSM) allowing a perturbative calculation of moments of the conductance
distribution, which turned out to be a broad function not Simply described in terms
of a few moments, but the method could not address the large L limit.

Anderson et. al [21], on the other hand, pointed out that although distribution
of g can be wide, its scaling still can be described by scaling of an appropriately
chosen single parameter, and addressed the one-dimensional case using the Landauer
formulation of the transport. The important conclusion of this work was that one
possible correct scaling variable in the localized regime is (1n 9) rather than (g), and

therefore the correct asymptotic scaling behavior is
(In g(L)) o< —L/§. (2.19)

Shapiro further extended this idea by developing an approximate renormalization-

group approach to study conductance distribution starting from the quasi-one-

20

dimensional systems [22], and found that for large L in d = 3 there are three limit-
ing distributions corresponding to the localized phase, mobility edge and delocalized
phase whose scaling can be described by a single parameter and gave an interpretation
[23] of results of Ref. [19, 20] in agreement with this finding.

The question of whether indeed the scaling theory of localization is described in
terms of only one-parameter is still under active investigation (see, e. g., work [24] and
references therein). The evidence accumulated so far, however, seems to validate the

one parameter hypothesis.

2.4 Universal conductance fluctuations and meso-
scopic physics

The most important finding in this area after the scaling theory of localization was the
discovery of universal conductance fluctuations [9]. They are direct consequence of
the absence of self-averaging in an ansamble of disordered systems when the quantum
coherence is present. Usually, if the correlation length in a system is lc, one may
average physical quantity X by averaging over NC = (L/lc)d subsystems of Size la.
The relative magnitude of fluctuations is then described by WAX), where
6X = X - (X). Since we assume that small pieces are statistically independent, and
since the variance of the distribution of X over these pieces is finite (due to their finite
size lo), the central limit theorem can be applied which gives that the distribution of

X is Gaussian with the relative variance equal to 1/\/ NC, thus

«6202) I. “ _
(X), “('5) ch 0‘. (2.20)

 

The principal finding, however, is that the variance of g of small disordered con-

ductors is universal and of the order of conductance quantum:

((69)?) ~ (62”.)? (2.21)

21

Since (9) = oLd‘z, we can rewrite Eq. (2.21) as

which shows that the actual fluctuations in quantum case are much larger compared to

 

the classical case (2.20), suggesting non-locality of the correlations inside the system.
This is in a qualitative agreement with our previous discussion of weak localization

corrections, whose source was indeed attributed to long-range interference effects.

2.5 Random matrix theory and universality

Wigner and Dyson proposed an approach for studying complex nuclei that consisted
of approximating the exact Hamiltonian of the system by a matrix drawn from an
ensemble of random Hamiltonians [26]. These ensembles are defined as follows: a
set of random Hamiltonians H can be fully defined by giving a probability density
P(H) of matrix elements hij of H in some basis. Assuming that they are uncorrelated
random variables, P(H) = Hi2]. P(hij).

An ensemble of random H’s in the random matrix theory (RMT) is then defined
by a requirement that P(H) is invariant with respect to a largest group of unitary
transformations that preserves global symmetries of the system (or, more generally, a
given set of constraints). The last requirement means that no preferential basis (that
is, no particular set of microscopic states) is chosen in defining the ensemble. The
necessary condition for this to be fulfilled is that P(H) is also invariant with respect
to the same group of unitary transformations. There is a freedom to choose P(H),
but statistical properties of the eigenspectrum of an ensemble are generally expected
to be independent of a particular choice of P once the energy spectrum is unfolded,
i.e eigenenergies rescaled such that the average density of states after rescaling is

constant [27], and this property is called universality [26, 27, 28].

22

One common choice is P(H) o< exp (Ter/vz), since the trace is invariant with
respect to any unitary transformation, and this class of RMT ensembles is called
Gaussian.

Contrary to the choice of P(H), the choice of a particular group of unitary trans-
formations that leaves P(H) invariant has important consequences on the ensemble
properties. Depending on whether transformations are maximal group of orthogonal
(hij are real), unitary (hij are complex) or symplectic (hij are two component Spin
1 / 2 spinors) transformations, one distinguishes three classes: GOE, GUE and GSE,
respectively.

Since every H in the ensemble can be diagonalized by a suitable change of the basis
set, the probability density P({E,}) to have a given sequences {E,-} of eigenenergies
will be, for Gaussian ensembles, 1'], exp(E, /v2) times the J acobian of the transforma-
tion, which gives [27]

Pa ({E.}) = Ca H lEi - Ejlfi exp(- 2 Elf/v2), (223)

i<j k
where D is the symmetry index, and B = 1, 2, 4 corresponds to GOE, GUE and GSE,
respectively. It is through the Jacobian that the chosen group of unitary transfor-
mations determines properties of the ensemble. It should be noticed that Peon —+ 0
whenever IE,- — E] —> O, which means that eigenenergies in the spectrum “repel” each
other. Probability distribution of level spacings is very accurately described by the

Wigner surmise [27]
7r 3

Dw(s) = a3“}e,xp(—:1—(5)2), (2.24)

where level spacing is s = EH1 — E,, and 6 is the average level spacing.
If the eigenenergies were just independent random variables, their spacing would

have been described by the Poisson distribution

Dp(s) = E’s-exp (—§). (2.25)

23

 

 

 

 

 

 

 

V

Figure 2.5: On the left is a schematical drawing of an experimental setup for studying
the conductance fluctuations and spectral properties of a mesoscopic two-dimensional
disordered conductor. Electron scatters many times before it leaves the system, and
typical trajectory is schematically drawn on the left. On the right is schematically
drawn expected conductance fluctuations as the voltage V is changed.

In this case Dp(0) = 1/6, i.e no level repulsion takes place.

The level spacing distribution can be studied in the context of condensed matter
physics through an experimental setup described in Fig. 2.5, where a small two-
dimensional sample is attached through the tunneling contacts to an electric circuit
[12]. Temperature is assumed to be small enough so that dissipation, together with
the coupling to the leads, introduces level broadening much smaller than the mean
level spacing between energy levels. Furthermore, if the localization length 5 > L,
the sample will conduct electricity whenever the chemical potential coincides with
an eigenenergy near the Fermi level. By changing the voltage V one thus obtains
a set of conductance spikes, and the distribution of spacings between the nearest-
neighboring spikes Should coincide with the distribution of level spacings of spectrum
of the microscopic Hamiltonian.

System considered in such an experiment, however, cannot be too large since

otherwise the individual levels cannot be resolved anymore even in principle [25].

24

On the other hand, in mesoscopic systems, which are of intermediate size between
microscopic and macroscopic, similar experiments can be performed giving novel kinds
of quantum coherence phenomena.

In the Anderson model, In the limit of W/ t -—> oo, eigenstates are localized at var-
ious sites of the lattice with eigenenergies approximately equal to the on-site disorder
strengths 6,. Since the latter are random variables, D(s) is correspondingly approx-
imately equal to Dp. As the disorder strength W decreases, eigenstates overlap in
space and D(s) smoothly interpolates between Dp and Dw [29]. We illustrate this in
Fig. 2.6, where a crossover between random (for large disorder strengths) and chaotic
(for small disorder strengths) spectra of a L = 40 sites d = 2 system occurs as W/ t

is varied.

2.5.1 Quantum Chaos

The main feature of the electron transport through the disordered medium is the
appearance of the long-range quantum interference effects that are responsible for
quantum corrections to the conductivity (2.17), and, for sufficiently large systems
in d = 2, Anderson localization. In one particular case, however, these corrections
disappear. That corresponds formaly to the case when L0 x Z in (2.17), which can
be achieved by setting the system size L to be of the Size of the mean free path 3.

Although this means that there are practically no impurities at all in the system,
one class of disordered systems satisfies this criterium, and those systems are small
metallic grains. Disorder here is present through the roughness of grain boundaries.
In this case, however, r z rm, and therefore one would expect that their spectral
properties would be in full accordance with RMT, not only for small level spacings as
it was the case with small but finite d = 2 systems discussed in the previous section.
Gorkov and Eliashberg [30] proposed this conjecture, which has been since proved by
Efetov [31].

25

 

Figure 2.6: 50 levels of unfolded energy spectrum of the Anderson model on a square
lattice L = 40 sites as a function of the strength of the disorder W/t, keeping the
disorder configuration fixed. For large W/ t, the spectrum exhibits characteristics of
a random spectra, characterized by large fluctuations in the number of levels in a
given energy range and a large number of avoided level crossing, as opposed to the
low W/ t regime in which both the quantities are smaller, and the spectrum exhibits
level repulsion and level rigidity characteristic for chaotic spectra.

26

Since the long—range quantum interference effects are absent, metallic grains are
quantum analogs of a classical particle performing motion in a cavity. This classical
motion is generally expected to be chaotic (even for many regularly-shaped cavities),
and for this reason the properties of an analog quantum systems are called quantum
chaos. Quantum chaos generally refers to the properties of quantum systems whose
classical analogs have chaotic dynamics [26].

In the case of a small metallic grain, spectral prOpertieS are described by the
RMT. The hope is therefore that RMT alone will be enough to describe completely
quantum chaos through classifying all the quantum chaotic systems in different RMT
ensembles. This hope, however, is still far from realization since, even in the case of a
diffusive electron in Fig. 2.5, whose classical dynamics is chaotic, quantum interference
effects give corrections to the RMT results which have to be accounted for by using the
specific features of the particular disordered system [12]. Nevertheless, so far RMT
approach to the quantum chaos has had considerable success in studying systems that
are effectively zero-dimensional [27, 12], such as small magnetic grains [31], nuclei

[26, 27 , 28], and even in the quantum chromo-dynamics [32].

27

CHAPTER 3

Two dimensional chiral disordered

systems on a square lattice

This chapter is devoted to the study of Hamiltonians of the form (2.4) on two dimen-
sional square lattice. The main feature of this class of disordered models is the chiral
symmetry, due to which all eigenenergies come in the opposite-energy pairs for each
realization of disorder.

Anderson transition, as discussed in Sec. 2.3, is a continuous quantum phase
transition, defined by a diverging correlation length 5, which in the localized phase
is identified with the localization length characterizing the exponential decay of the
eigenstates at infinity, as a parameter in Hamiltonian is being changed.

Scaling theory of localization predicts absence of extended states in two dimen-
sional disordered systems. If spin-scattering is present, however, picture changes
and two dimensional systems from symplectic ensemble exhibit localization transi-
tion even in two dimensions [8]. Presence of strong magnetic field in two dimensional
disordered systems, on the other hand, leads to a completely different behavior - the
integer quantum Hall effect (IQHE) — where critical states are present at the middle
of each of disorder-broadened Landau levels [11].

Another class of models exhibiting localization properties different from the sys-

28

tems mentioned above are systems with chiral (particle-hole) symmetry. Such systems
are defined on a bipartite lattice with only hopping (off-diagonal or bond-) disorder.
Wegner was first to realize the importance of this symmetry in disordered systems
[40, 38, 39], and even one-dimensional systems with this symmetry are known to have
peculiar properties, such as diverging DOS at the band center [33], where the eigen-
state decays as exp(—'y\/1—") [34, 35], in contrast to one-dimensional site-disordered
systems which have DOS bounded [36] and all states localized.

There are several models with chiral symmetry that have been extensively studied.
The simplest two, in the sense that only one orbital per site and nearest neighbor
hopping are included, time-reversal symmetry present and the spin not relevant, are
the Anderson bond-disordered model (ABD) [33, 34, 35, 44, 45] and the random Dirac
fermion model (RDF) [41, 42, 43]. The main difference between these two models is
that, in the non-disordered case, ABD model has a line of points as the Fermi surface
at half-filling while RDF has a point Fermi surface and linear dispersion of energies.

Four particular models, defined by four particular distributions of hopping energies
t,,- in Eq. (2.4), are considered: In Sec. 3.2 the models studied are the double-exchange
model with (DE) and without (|DE|) complex phase (to be discussed below); Sec-
tion 3.3 is devoted to a study of Anderson bond-disordered model (ABD), where
tiJ‘,S are uniformly distributed between —1 and 1; In Sec. 3.4 the random flux model
(RFM) is studied, describing an electron on a square lattice with randomly distributed

magnetic flux per plaquette.

3. 1 Chiral symmetry

Suppose that the lattice is composed of two sublattices A and B with, respectively, N A

and N 3 sites. The corresponding bond-disordered Hamiltonian with chiral symmetry

29

can be most generally written as:
H = Z (t,,jc[cj + H.c) . (3.1)
iEAJEB
The diflerence between this Hamiltonian and Hamiltonian (2.4) is that here no par-
ticular restriction on what sites are connected by h0pping is imposed.

It is easy to Show that for every eigenstate lib) of (3.1) with energy E there is an
eigenstate with energy —E with a wavefunction that has the opposite Sign at each
Site of one of the two sublattices.

If the total number of cites N = N A + N B is odd and open boundary condition is
applied (in order to keep the symmetry), then, since all eigenstates come in opposite
energy pairs, there will be exactly one state with eigenenergy 0. This can be further-
more generalized, and if m = N A — N B > 0, then there exist exactly m zero-energy
eigenstates that have vanishing amplitude on the sublattice B [39].

On the other hand, if m = 0, the electron has equal probability of occupying each
of the two sublattices [93]. To Show this, (3.1) is represented in the basis where the
first and second half of the basis vectors are eigenstates of the position operator on

sites of sublattice A and B, respectively. The Hamiltonian is then represented as

0 M
H = , (3.2)

MlO

where M is a square matrix of hopping elements from one sublattice to the other.

Eigenstate lib) = WA) satisfies:
We)
E = <wIHiw> = 2Re<¢AIMIwBi (3.3)
On the other hand,
M E
HM) = [$13) = WA) . (3.4)
M’WA) Ell/)3)

30

From (3.3) and (3.4) now follows that, for the real M, (ibAlibA) = 1/2.

In this work only even L finite size systems on a square lattice with periodic
boundary conditions are studied, because one of the main goals of this work is to
understand the behavior of eigenstates in the vicinity of the critical point of ABD
model on the infinite square lattice, which in turn has m = O. In contrast, the limit

L —-> 00 for odd L and open boundary conditions has m = 1.

3.2 Double-exchange model
The double exchange model (DE) studied in this section is described by the Hamil—
tonian (2.4) on a square lattice, with the hopping elements

7'. m o 0n o 6m ' .
tnm = cos — cos — + sm — sm 7e’(4’"‘¢m) = |tnm|e'”’""‘, (3.5)

2 2 2

where (6", (in) define the orientation of Spins Sn at the Site 7",, of the lattice. The model
corresponds to an electron hopping on a square lattice in the presence of randomly
oriented local classical spins to which it couples with an infinitely strong ferromagnetic
coupling.

Here we focus on the structure of the electronic wave functions as a function of
energy in the presence of localized spin disorder in d = 2. In particular we study how
the the density of states and the nature of the wave function compares with another
studied disordered model, as well as what is the influence of the complex phase (Dam
of tum on the eigenstates.

Thermodynamic prOpertieS of this model have been studied by Yunoki et al. [51]
using a Monte Carlo (MC) procedure where the MC moves of the localized classical
spins are determined by the free energy of the itinerant spinless fermions. Also a
similar procedure has been used by Yunoki and Moreo [52] to look at the dynamic

properties in one dimension.

31

Localization properties of the d = 3 DE model have been studied by Li, Bishop,
and Soukoulis [55] in the same limit, where they have found that the mobility edge lies
too close to the band edge and consequently the paramagnetic phase is not insulating
for realistic conduction electron (hole) concentrations. They argue that additional
disorder mechanisms such as strong on—site disorder have to be included to understand
the experimental observations in La1_xCaIMnOg and similar 3D systems.

AS mentioned above, we are interested in understanding the nature of the elec-
tronic wave functions in the presence of localized spin disorder (which is self-
consistently determined by the electronic free energy) in d = 2 (2D). quasi-2D and 2D
DE systems such as LaermanO-l , La1_xSr1+anO4 have been discovered [53, 54]
and our results should be applicable to these systems. Here we only study the ex-
treme disorder limit when the localized spins are randomly oriented. Physically, this
corresponds to the high-temperature limit of the spin system.

Our study has been carried out by the exact numerical diagonalization of DE,
[DE], and ABD Hamiltonians for systems of N sites on a torus in 2D, for system sizes
up to L = 70. For L even, the lattice is bipartite with nA = 713 , while for L odd,
lattice is not bipartite due to the torus geometry of the system, so that in either case

states found in Ref. [39] do not appear in systems we have studied.

3.2.1 Density of states

The normalized average density of states (DOS) pL(E) E (2,6(E — E,-))/L2 of a
system of size N = L x L was calculated by binning all the eigenenergies for large
number M of different realizations of disorder.

Density of states for all three models in 2D are presented in Fig. 3.1. In the absence
of disorder, the DOS is flat near the band edges (E = :E4) and has a singularity at
the band center. In the presence of hopping disorder, the DOS near the band edge is

dramatically reduced. Going from ABD to [DE], there is a suppression of the strength

32

 

0.6 r

0.5 -

 

   
  

0'4 ' ABD N= 302

[DE | N = 402
0-3 r DE N: 302, 402
0.2 -

0.1 -

 

 

 

0.0 U I I l I I 1

Figure 3.1: Density of states for the two—dimensional DE, [DE[, and ABD models.

of the central peak, and states closer to the band center are now pushed toward the
band edges. The complex phase of the hopping amplitude in the DE model further
reduces the central peak, and also the DOS near the band edges. DOS for the ABD
model calculated here is in agreement with the one presented in Ref. [45].

DOS is continuous even for small size systems due to the averaging over disorder,
and depends on N for smaller system sizes. p L(E), however, converges with increasing
N to a limiting DOS p(E) in 2D which we illustrate below using the example of the
DE model in 2D.

In passing I would like to mention that calculations of the DOS for two electrons
on a four-site plaquette for [DE] and DE were presented in Ref. [56], as an illustration

of some general DOS properties of DE and [DE]. The four-sites DOS of two electrons

33

presented there, however, is decreasing near the band center and vanishing at the
band center. Such DOS does not have a large peak at E = O, which contains 3/4
of all the eigenstates. There is no obvious physical reason why this peak Should be
absent, and, therefore, the interpretation of the DOS presented in Ref. [56] may not

be very meaningful, even qualitatively.

3.2.2 Scaling properties of eigenstates

Localization properties were studied by calculating the participation number PN of a

normalized energy eigenfunctions ‘1! [57]:

PN(‘I’)_1 = Z: I‘I’O‘ll", (3-6)

where the summation goes over all N sites. The participation number gives the
number of sites at which \II significantly differs from zero. For example, PN of a state
localized at one site is 1, while in the case of a Bloch wave, PN equals N.

W is delocalized if PN(\II) ——> N when N —> 00, or, equivalently, PN/N -—> p > 0
for N —> 00, p being the fraction of the total volume significantly occupied by ‘1!
(participation ratio). Although this distinguishes localized and extended states, wave
functions also can have fractal character, for instance at the mobility edge in d = 3
[58, 59]. This critical behavior of the wave functions can be characterized by a scaling
of the participation number (for a given energy averaged over disorder as described

below) with system size, i.e.,
PN oc N3,O < [3 <1. (3.7)

The quantity P151 (commonly referred to as the inverse participation number)
is the second moment of the electron probability distribution pe(f') E [\II(F)|2, and
Eq. (3.7) reflects the scaling properties of the density-density correlation function.

Phrther insight into the nature of the eigenstates can be obtained from the higher

34

moments of p607), whose anomalous scaling properties indicate that the eigenstates
are multifractals. Multifractality of eigenstates is studied in Sec. 3.3.3, where further
discussion with definitions and references is given.

On length scales smaller than the localization length, ,8 characterizes the fluctu-
ations of the wave function, and the closer the value of B is to 1, the more similar is
the wave function to a delocalized one.

In our numerical work we have calculated the average PN for a given energy E, by

averaging over all eigenstates that belong to an energy bin 5 E (E—AE / 2, E + AE / 2),

 

for M realizations of disorder wl, . . . ,wM; as follows:
1 M
P E AE E P \II,~ , 3.8
N( t ) N(€);EZ€€ N( J) ( )

where ‘Ilij is the i-th eigenstate with the eigenenergy Eij for the j-th realization of
disorder wj, and, N (e) is the total number of states in the bin 6 for all M realizations
of disorder. The value of PN(E) is now obtained from the M —> 00, AE —> 0 limits.
We emphasize this because near the band center PN is particularly sensitive to E, as
discussed in Sec. 3.2.3.

The energy dependence of the participation ratio pN = PN/N is given in Fig. 3.2
for system sizes L = 20, 30, 40, 50, 60, for all three models. All calculations presented
in this section are done for E < —0.02 (results are even functions of E, as dis-
cussed before), while the results for the band center states are discussed separately
in Sec. 3.2.3. Averaging is performed over all eigenstates within consecutive energy
intervals of width AB = 0.04. We find that pN decreases with N, indicating that
states get localized (PN = const) in the limit N -—> oo. pN for the ABD model is
smaller than for the other two, since the ABD model has more disorder than DE and
[DE], and it appears that there is a cusp at the band center in all the three models.
An interesting effect in Fig. 3.2 is that there exists an energy range about the band

center inside which pN increases by turning on the complex phase of the hopping

35

 

 

 

 

 

 

 

II‘IIIII'IITUIIIII

-3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0

Figure 3.2: Participation ratio my 5 PN/N for DE, [DE], and ABD models as a
function of energy for E < —0.02. For each model the system size is changing from
L = 20 (uppermost curve) to L = 60 (lowermost curve) in steps of 10. pN is symmetric
about E = 0. The decrease of pN with L suggests localization.

elements tnm, and outside which pN decreases.

Now we discuss in detail the scaling behavior of PN. Figure 3.3 shows the depen-
dence of PN on the system size in the indicated energy range, in log-log plot, for all
three models. There are three different regimes. The first regime occurs for energies
near the band edges, where the participation number is independent of the system
size. The second regime corresponds to intermediate energies for which ln(PN) de-
pends nonlinearly on ln(N), and, finally, the third regime corresponds to eigenstates
for which PN o< N H. We have extracted the exponent B from a least square linear fit
of the data in Fig. 3.3, and the result is presented in Fig. 3.4.

To measure the goodness of the fit, we also did a least—square fitting of the same
data to a quadratic function, lnPN(E) = C(E) lnN + 7(E)(ln N)2 + C(E), and the

strength of the quadratic term (the curvature) 7 is presented in Fig. 3.5 for all three

36

 

103

P 102 /

N

 

 

 

 

101

 

103

P 2
N10

 

 

 

 

 

 

101

 

103

P 2
N10

 

101
102 103 104
N

Figure 3.3: Dependence of the participation number PN on the system size N = L x L
for DE, [DE], and ABD models as a function of energy.

37

 

1.0
. DE .

0.9 - - 1+: ..
0.8 - ./ IDEI
0.7 - .0 -

I O

 

...-......‘.-.-...-....III-.

0-6 ‘ .° ° ABD
[3 0.5 - ° - '
'] o o .-
0.4 - . . ,-
0.3 - ° .° -
0.2 - ,-
0.1 - ° . .-

0.0 -—V°‘-‘h—¢

-O.1 IIjiUlIIUI'IIIIIIIIIIIIIT'ITII'IIII

-3.5 -3.0 -25 -20 -1.5 -1.0 -0.5 0.0
E

I I I I I 1 r l I I I l I l I I

 

I

 

 

 

Figure 3.4: Energy dependence of exponent [3 describing the power-law scaling of the
participation number PN. for the DE, [DE], and ABD models. Full curve corresponds
to the case when PN can be fitted by NB for the data in Fig. 3.3. The dotted lines
correspond to the situation when there was significant nonlinearity of the fit (see text
for the discussion and Fig. 3.5 for a measure of the goodness of the linear fit).

38

 

 

0.05 ‘E I DE ]

 

 

 

 

 

Figure 3.5: Goodness of the linear fit of ln PN vs lnN for all three models. 7 = 0
corresponds to the power-law dependence of PN on N. Fluctuation of 7 near the
band edge are due to low DOS (i.e., small number of states that were averaged over).
In this range, however, '7 = 0 because the participation number is independent on N.

models. 7 a: 0 therefore indicates the validity of the power-law behavior, and the
thick lines in Fig. 3.4 correspond to the energy ranges in which the power-law (3.7)
is obeyed.

The curvature of the fit also clearly distinguishes the three different energy regimes
introduced above in the case of DE and [DE] models. For example, in Fig. 3.5, states
of the DE model with energy E > —1.88 are well described by a power-law dependence
(the third regime), while those between —3.00 and —1.88 are not (the second regime).
States with E < —3.00 can be again described by N 5 dependence with B = 0, i.e.,
PN = const, and fluctuations in both 3 and 'y in this regime are due to the low DOS.

For the ABD model, however, power-law dependence of PN is not seen for the system

39

sizes studied. From Fig. 3.4, we also see that the exponent D has the same value for
DE and [DE] model for the energy EC = —1.49 :l: 0.11.

In order to obtain '7 accurately, averaging over a large number of disorder
configurations is needed. In our calculations we used: for the DE model and
L = 510, 20, 30,40, 50, 60, respectively, M = 104, 104, 4000, 1200,310, 199; for the
[DE] model and same system sizes, M = 5 x 104, 104, 6000, 2000, 500, 295; for
the ABD model and L = 10, 12, 14, 16, 18, 20, 22, 30, 40, 50, 60, respectively, M =
5 x 104,5 x 104, 5 x 104,5 x 104, 2 x 104, 104, 104,6000, 2000, 500,300.

The observed dependence of PN on N (for all three models) can be understood by
assuming the existence of a localization length g (E) for any eigenstate with energy E
for an infinite system. If L > 5 (E), a wave function with energy E is determined only
by local realizations of disorder in the part of the system where the wave function is
localized, and PN is consequently independent of the system size. We refer to this
case as the first regime, in which D(E) = 0 and 7(E) = 0. Fluctuations of ,6 and '7
presented, respectively, in Figs. 3.4 and 3.5 in this regime are due to the small number
of states that were averaged over (low DOS).

The second regime corresponds to eigenstates for which the linear system sizes
studied L a: g (E) This is a crossover regime in which PN is determined from the
interplay of both localization and finite-size effects, and in this regime we find 7(E) <
O. The third regime, however, corresponds to states that have localization length
much larger than the system sizes studied L < é (E) In this case, the participation
number can be remarkably well described by a power-law dependence on the system
size, and, therefore, D(E) > O and 7(E) z 0.

All these discussions are summarized in Fig. 3.6, where we plot schematically
ln(PN) vs ln(N), for three different energies E1 < E2 < E3 corresponding to the
three regimes described above, and for the range of system sizes studied in our work

(indicated by the two vertical dotted lines). We can clearly identify the three regimes,

40

 

 

L<< A

 

 

 

ln(PN)

 

 

 

 

 

 

 

ln(L)

Figure 3.6: Schematic al representation of ln PN vs 1n N for different energies sharing
a common behavior. The curves 1, 2, and 3 correspond to three different energies
E1 < E2 < E3 . Two long vertical dotted lines symbolically represent the range of
system sizes studied in this dissertation. A is the localization length 5.

and curves 1, 2, and 3, correspond to In PN = const, nonlinear dependence of In PN
on In N, and PN = N a , respectively. For the third energy regime we also indicate in
Fig. 3.6 a crossover region (between the two short vertical dashed lines) separating
regions L << 5 and L > 5.

Critical behavior described by Eq. (3.7) is present only for system sizes much
smaller than 5, which in turn never happens in the ABD model away from the vicinity
of the band center and in the DE and [DE] models near the band edges, due to too
small localization lengths. This is in the agreement with the calculation of localization
lengths by the transfer-matrix method [18] in the case of the ABD model presented
in Ref. [45], which shows that the localization length for [E] = 0.1 is slightly above
100 sites, and thus comparable to the system sizes studied here.

As can be seen in Figs. 3.4 and 3.5, the power-law dependence in the [DE] model

begins at lower energies compared to the DE model; in the energy region —-2.08 < E <

41

—1.88 the critical behavior is present only in the [DE] model. This strongly suggests
that in this energy range the localization length is also larger in the [DE] model
compared to that in the DE model. Once the power law is obeyed in both models, the
exponent need not be simply related to the localization length and thus fiDE(E) >
B[DE|(E) for E > Ec does not necessarily imply 5DE(E) > 5IDE[(E). Nevertheless,
eigenstates in this energy range become geometrically more similar to extended states,

within the localization length, if the complex phase of hopping elements is present.

3.2.3 Band center states

As seen in Fig. 3.2, pN(E) appears to develop a cusp at E = 0. This suggests some
sort of singular behavior of the band center states. To analyze this we give in Fig. 3.7
the dependence of participation number on system size in a log-log plot, for E = 0
and various bin sizes for all three models. 5 (E = 0) is calculated from the least-square
linear fit of the data presented for the smallest bin size studied, as a slope of straight
lines in Fig. 3.2. The number of disorder configurations M was chosen such that the
relative error of PN is smaller than 1.4% for the ABD and [DE] models, and 0.55% in
the case of the DE model. For the [DE] model we find ,B(E = 0) = 0.668 :l: 0.002. On
the other hand, if we extrapolate the D(E) given in Fig. 3.4, we estimate D(E —> 0) =
0.85 i002. In the case of the DE model, ,8(E = 0) = 0.757:l:0.002 and extrapolation
from Fig. 3.4 gives ,B(E —+ 0) = 0.90 :l: 0.02. In the case of the ABD model, the
dependence on DE is much stronger than for the above two models, and therefore we
had to perform calculation for AE = 0.0001 for L = 10, 20, 30, 40 and AE = AE2(L)
for L = 50, 60 (where AE2(L) is defined below). The fit gives D(E = 0) = 0.41 :E0.01,
Which is again much lower than the values presented in Fig. 3.4 (in this case power-law
is only approximately satisfied). These results suggest that ,B(E = 0) and [3 (E —> 0)
are not the same, i. e., they suggest a nonanalytical behavior of PN at the band center.

To further support this claim, we give now more detailed analysis of the results

42

 

103

U
Ed

I J l l
->

 

 

.0 o AE=0.01
...-r <> AE=0.001
«2 AE<0.0001

101' T I I IIIIII I I I IIIIII
102 103 104

 

 

Figure 3.7: Dependence of the participation number on system size in log-log plot for
the DE, [DE], and ABD models for band center states (E = 0) and various bin sizes
AE. Solid lines are the least-square linear fits, while dotted lines are for guidance
only. Arrows indicate independent estimates of the system sizes at which a deviation
from the power-law begins to show up.

43

obtained for the band center states. The curves in Fig. 3.7 for larger AES have slight
upward curvature, i. e., there is a deviation from the power-law dependence which is
found for the smallest AE for each model. Furthermore, sensitivity to AE depends
on the model and increases with the increase of the system Size N.

We argue now that the two states closest to the band center for a given disorder
configuration have different scaling properties than the rest of the states. The bending
therefore begins for system sizes when the averaging takes more than two eigenstates
per disorder configuration for a given energy bin. In the limit N —> 00 this corresponds
to a discontinuous change of scaling properties of two eigenstates at the band center
going from E ¢ 0 to E = 0. The two states differ only in the Sign of the phase at all
sites of one of two interpenetrating lattices, due to the bipartite nature of 2D square
lattice.

The average number of states, N0(AE), in the interval AE about the band center

(for which PN is calculated) equals

0
N0(AE) = 2L2] pL(E)dE. (3.9)
—AE/2
The linear size L2 of a system for which every realization of disorder gives two eigen-

states per disorder configuration in the given energy bin AB is obtained from Eq. (3.9)

by equating No to 2:

MI:-

L2(1_\.E)=(/O pL,(E)dE)— . (3.10)

—AE/2
We use this quantity to estimate the system size for which the bending begins to
Show up, and arrows in Fig. 3.7 indicate the estimates. Let us discuss the ABD
model first. This model has the largest DOS in the vicinity of the band center (see
Fig. 3.1), thus we expect this model to have the most pronounced bending of all
the three models studied, which is indeed the case in Fig. 3.7. The ABD model has

L2 (0.001) = 41 :L- 1, which very accurately estimates the system size for which the

44

bending of the corresponding curve begins to Show up in Fig. 3.7. For the [DE] model
we find that L2(0.01) = 22 :l: 1 which again accurately estimates the value (L z 20)
at which the system size dependence on AE occurs in Fig. 3.7. Excellent power-law
dependence of PN on N for AE = 0.001 for the same model is compatible with value
of L2(0.001) = 59 :l: 1. Finally, for the DE model, which has the lowest DOS near
the band center, we find L2(0.01) = 26 :i: 1. This model has the least pronounced
bending, and the value of B (E = 0) for this model is obtained from the fit of data for
system sizes up to L = 26.

Scaling of the participation number with system size has been studied by Eilmes,
Rbmer, and Schreiber [45] for a bond disorder model in which the hopping matrix
elements are uniformly distributed random numbers between -1/2 and 1/2 (ABDH
model) for three different energies. The participation number was found in Ref. [45]
to scale as a power-law with B = 0.50 :E 0.06 for E = 0. The scaling was done for
system sizes from L = 50 up to L = 200, using the bin size AB = 0.0005. The ABDH
model, however, has higher DOS near the band center compared to the ABD model
[45], and, therefore, has higher sensitivity on the bin size (smaller L2(AE)) when one
calculates PN. Our analysis applied to the ABDH model gives L2(0.0005) = 49 :l: 1.5,
and, thus, we expect the upward curvature of PN vs N to start from L z 50. Analysis
of data for system sizes L > 50 therefore will tend to overestimate the value of B, so

that we expect )8 to be considerably smaller than the value quoted above.

3.3 Anderson bond-disordered model

This work is concerned with the ABD model on a square lattice of size L and periodic

boundary conditions, defined by the Hamiltonian:

H = —60 Z (ti’jCICj + H.C) , (3.11)
(M)

45

where the angular brackets denote neighboring sites on the lattice, c,- (c[) is armihi-
lation (creation) operator for the electron at site i, and t’s are uniformly distributed
random variables tm- E (1 — 2w, 1), with 0 < w 3 1. They represent random hOpping
energies between nearest neighbors, expressed in units of energy so, which is set to 1
hereafter.

Interest in this model mainly comes from its unusual scaling properties at the
band center, where Soukoulis et. al [44] have found a critical state at the band center
using Green’s function[17] and transfer matrix method (TMM) [47, 48]. More recent
TMM calculation by Eilmes et. al[45] confirmed this result with a higher accuracy
and showed the validity of one-parameter scaling not too close to the band center.
Nevertheless, Miler and Wang [46] have found in their study of two models with chiral
symmetry an apparent band of extended states near the band center, and it remained
unclear what is the fate of these states in the infinite 2D system. Another study
[50] yet showed that scaling exponent of the average participation number changed
discontinuously near the band center, and the explicit dependence of this energy
on the system size proposed by authors implied the existence of another diverging
length scale in the problem. The last effect is rather subtle to calculate and led to a
different participation number scaling exponent of the ABD model at the band center
in Ref. [45] compared to the one calculated here, as discussed in detail in Sec. 3.3.3
below. Furthermore, Brouwer et al. [49] have calculated the conductance distribution
of quantum wires described by (3.11), and showed its non-universality and pointed
out the necessity to introduce an additional microscopic parameter.

It is thus the goal of this section to present a detailed study of the scaling of
localization length on the approach to the band center for an infinite 2D square
lattice, and test the validity of one-parameter scaling, as well as to calculate the
multifractal properties of the electron probability density on length scales smaller than

the localization length. Also a new analytical expression for the DOS of the infinite

46

two-dimensional system near the band center is proposed. Some general properties
and exact results were already presented in Sec. 3.1, while this section is organized as
follows: Calculation of DOS is presented and analyzed in Sec. 3.3.1; Sec. 3.3.2 contains
analysis of level spacing distributions between the nearest neighbors; Multifractal
properties of eigenstates are studied in Sec. 3.3.3; and, finally, scaling of the Lyapunov
exponents of transfer matrices of long strips and the scaling of localization length are

studied in Sec. 3.3.4.

3.3.1 Density of states near the band center

Density of states (DOS) is calculated by exact numerical diagonalization of finite
size Hamiltonians for various L for many configurations of disorder, and binning of
the eigenenergies. Obtained disorder-averaged DOS for each system size, pL(E), is
normalized to 1. The L—dependent parts of such obtained pL(E) are then removed
leaving the L independent DOS p(E), which is therefore expected to be correct in
the L —> oo limit. The removal of finite size dependence is based on an observation
that DOS converges quickly away from the band center with increasing L. Thus,
only a small number of eigenenergies (up to 20) closest to the band center and the
corresponding DOS histograms has been calculated for each L. The calculated pL(E)
plotted on a single graph revealed that three bins closest to the band center are where
the system size dependence sets in. Their removal thus led to the DOS p(E) of the
infinite system.

Results for pL(]E]) are given in Fig. 3.8, for system sizes L = 10,20, ...,60 and
number of disordered configurations ranging, respectively, from 160000 to 4100. They

are fitted to a power-law pL(E) = CLlEl‘aL, as well as to the Gade’s result [38]:

1
pL(E) = CLmexp (mm—1mm). (3.12)
All the calculations were done for several bin sizes, and the obtained values of fitting

47

 

 

 

 

 

 

 

 

 

 

 

10-3 102 10«3 10-2 104 10-3 102
E E E

Figure 3.8: Density of states pL(E) of the ABD model for w = 1 near the band center
for system sizes L = 10, ..., 60 in log-log plot. Full lines are fits to the Gade’s form
(3.12), while dotted lines are fits to the power-law. Fit to the Gade’s form is done for
all the points except the one closest to the band center, while fit to the power-law is
done for all the bins except the three bins closest to the band center, which are the
L—dependent parts of pL(E).

48

parameters for a given L were the same within error bars.

Figure 3.8 shows that Gade’s result (3.12) describes pL very accurately for
L 2 40, including the sizedependent part. Gade has obtained this result from
the renormalization-group analysis of a non-linear sigma model description of the
problem, and the calculation could not give the value of the EL, which is now calcu-
lated for the first time. Despite this, the expression found to best fit the obtained

L—independent DOS, given in Fig. 3.9, is

1
= Cm exp (—m/—ln ]E[), (3.13)

with n = 1.345 i 0.005 and C = 1.30 :l: 0.03, represented by the full line in the

p(E)

same figure. The observed range in which (3.13) is accurate is for all energies studied

smaller than 6 x 10”.

3.3.2 Distribution of the nearest-neighbor level spacings

In the localized regime, an eigenstate is determined mainly by local configurations
of disorder where the wavefunction is localized, and two eigenstates close in energy
are spatially far apart. Level repulsion is therefore absent and the distribution of the
nearest neighbor level spacings s E E,“ — E,- is Poissonian Dp(s) [29], as discussed
in Sec. 2.5. In the delocalized phase, on the other hand, eigenstates are extended
throughout the system and the level repulsion becomes significant for eigenstates
with close energies and the level spacing distribution is that of GOE, very accurately
described by the Wigner surmise Dw(s) [29]. In finite-size systems, localized states
are on average at a distance L rather than infinitely far apart. This leads to a repulsion
between adjacent energy levels and non-universal distribution DL(3) .

Shklovskii et. al. [29] have shown that 01(3) of the d = 3 Anderson site-disordered
model (ASD) exhibits linear dependence on s, characteristic of Dw(s) for small s

and exponential tail characteristic of Dp(s) for large 3. They were able, from the

49

 

p(E)

 

 

 

 

10'3 10‘2 10'1
E

Figure 3.9: Density of states of the ABD model for w = 1 near the band cen-
ter. The graph is obtained from the data in Fig. 3.8 by removing the three bins
closest to the band center, leaving the L-independent p(E). The fit is p(E) =

Cexp (—m/— 1n |E|) /,/|E|, with n =1.345 :1: 0.005 and C = 1.30 :E 0.03.

50

finite size scaling analysis of the tail, to accurately determine the critical point and
also the exponent. In the infinite size limit, they have recovered not only Dp(s)
in the insulating phase and Dw(s) in the conducting phase, but also a system size
independent non-universal distribution at the critical point which was furthermore
Shown by Braun et. al. [60] to be dependent on boundary conditions.

To see the effect of the chiral symmetry of (3.1) on the distribution of level
spacings, let us for a moment consider the i-th eigenenergy E,- of the ASD model.
Upon averaging over disorder, E,- will be distributed between E?“ and Er”, ac-
cording to some distribution. Some of the eigenenergies, for i close to N / 2, will have
EF’" < 0 < Eff“. This is, however, forbidden for eigenstates of the ABD model
since every E,- of (3.1) is negative for i < N / 2 and positive for i > N / 2. This means
that eigenenergies of (3.1) close to the band center are effectively pushed away from
it due to symmetry.

If L is much smaller than the localization length, states will be repelled due to
their large spatial overlap. But the two states closest to the band center, being simply
related to each another by the symmetry, will not repel at all, i.e. the state closest to
the band center is at the (high energy) end of the spectrum. Thus, these two states
are distributed around zero, where distributions of all other individual levels go to
zero.

To explore consequences of this simple analysis, the level spacing distribution is
calculated between each pair of adjacent levels separately. Let us denote by D,(s)
level spacing distribution between i and i + 1 energy level after unfolding of the
spectrum [27], i.e. expressing level energies in units of mean level spacing, where
0 < E1 < E2 < < EN/g. Figure 3.10 shows 01(3),...,D5(s), for L = 20,40 and,
150000,120000 configurations, respectively. It can be seen that D1(s) is distinctly
different than D2(s), ..., D5(s) in both cases.

The same effect was also present for L = 10 and 30, while for L = 50 and 60 the

51

 

 

 

 

 

 

 

 

 

 

r L=20 .-.- [4:40
D(S) , i=1 L i=1
‘ -------- i=2 """" i=2
0.5_ ............ i=3 0.5’ ............ i=3
i ”w” i=4 ............. i=4
W i=5 ............... i=5
00 1 . 4 0.0 A ' '
0 2 3 0 1 2 3
8 8

Figure 3.10: Distributions of level spacings (after unfolding of the spectrum) D;(s)
between i and i + 1 level, counted from the band center. D1(s) is distinctly different
than other distributions due to the presence of the symmetry.

number of disorder configurations was insufficient for accurate enough determination
of individual D,(s). This illustrates how the presence of chiral symmetry can pro-
foundly influence the spectral characteristics near the band center, despite the fact
that DOS of ABD and ASD models seem to have the same shape for appropriately
chosen pairs of disorder parameters 21) and W away from the band center (and after

rescaling of 60) [45].

3.3.3 Multifractality of eigenstates

Eigenstate of an electron in a random potential fluctuates from site to site and it was
proposed that the eigenstate at the mobility edge in disordered systems in general
should have fractal structure [58, 59]. It was shown that even localized states in
d = 1, 2 exhibited fractal character on length scales smaller than 5 [61, 62, 63].
Inverse participation numbers Zq (IPN) are particularly convenient quantities to

describe scaling properties of the probability distribution of an electron. IPN of an

52

eigenstate \I’ for different q’s are defined as:

LxL
2.010 a Z Ina-M2? (3.14)
i=1
Intuitively, their meaning can be seen by looking at the participation number
Z2(\Il)‘1 : it is equal to 1 for a state localized at one site and to N = L2 for plane-
waves. Participation number thus gives generally the number of Sites at which the
wavefunction is significantly different than zero. Participation numbers Zq(\II)"1 then
generalize this by giving the number of sites where probability distribution of electron
is very high (for large positive q’s), very low (for large negative q’s), and in between
these extrema, continuously parameterized by q.
More convenient, with the advantage of being defined as averages over disorder at

a given energy E, are IPN defined as functions of E and system size L,

Zq(L, E) 2 <2 Zq(\II,-)6(E,- — E)>, (3.15)

where the brackets denote averaging over disorder. Zq(E, L) can be numerically
calculated by averaging (3.14) over all eigenstates from M configurations of disorder
belonging to an energy interval of width AE around E, and studying the limit AE —+
O for large M [50].

Wegner [64] pioneered this kind of investigation, and Castellani and Peliti [65]
proposed that eigenstates near the critical point were multifractal on length scales
smaller than 5. The most important feature of IPN of eigenstates is their scaling with

system size and energy [64, 65]:

2.,(L, E) ~ L‘Tq, (3.16)

Zq(L,E) ~ [E—ECI’W, (3.17)

where Ec is the critical energy. The former scaling is present at any E for L < 5 (E),

while the latter holds in the critical region of the transition.

53

Within the framework of multifractality [66, 67, 68, 69], electron probability den-
sity is characterized by several quantities that can be derived from rq — the gen-
eralized dimension D, and the Singularity strength aq of the q—th Singularity (q-th

fractal) with the fractal dimension fq :

dr

—3 , fq(aq) E aqq — rq. (3.18)

(q-1)DQET<II aqE dq

Dq represents a generalization of the fractal dimension, and it is constant and equal
to the fractal dimension for ordinary fractals, while fq(aq) is the singularity strength
Spectrum describing multifractal as an interlaced set of fractals with fractal dimen-
sions fq, where the measure on the q-th fractal scales as a power-law with exponent
(19. These quantities have several general properties: Do is the fractal dimension of
the support containing probability distribution (two in this case); D1 is called in-
formation dimension since it describes scaling of the entropy of the measure (Dq=1
is defined from the limit q -+ 1 of the definition (3.18 [11]), and there exist finite
Dmin = anoo and Dmax = Dq_._oo.

This work is concerned mainly with the spectrum of generalized dimensions Dq
characterizing the spatial structure of eigenstates on the length scales smaller than
the localization length. Before a detailed discussion of the results, an overview of
the main results of this part of the section is given. Scaling of IPN at the band
center is calculated first, and I Show that scaling properties (that is, whole spectrum
of generalized dimensions 0.) changes discontinuously near the band center, at an
energy E’(L) for a range of L studied. It is then shown that E’ can be quite accurately
identified with half of the width of the energy range around band center within which
two states occur on average in the ensemble of disordered systems. Existence of this
energy reveals the existence of a new length scale 5’(E), diverging when E —-> 0, that
is the system size at which IPN’s change their scaling properties from one power-

law dependence on L to another. This change is then explained as a finite size

54

Within the framework of multifractality [66, 67, 68, 69], electron probability den-
sity is characterized by several quantities that can be derived from rq — the gen-
eralized dimension 0,, and the singularity strength 01., of the q-th singularity (q-th

fractal) with the fractal dimension fq :

fl

dq , fq(aq) E aqq — rq. (3.18)

(q—1)Dq_=.rq, 0.15

D, represents a generalization of the fractal dimension, and it is constant and equal
to the fractal dimension for ordinary fractals, while fq(aq) is the singularity strength
spectrum describing multifractal as an interlaced set of fractals with fractal dimen-
sions fq, where the measure on the q—th fractal scales as a power-law with exponent
aq. These quantities have several general properties: Do is the fractal dimension of
the support containing probability distribution (two in this case); D1 is called in-
formation dimension since it describes scaling of the entropy of the measure (Dq=1
is defined from the limit q —> 1 of the definition (3.18 [11]), and there exist finite
Dmin = Dq_.oo and DW = Dq—b—oo-

This work is concerned mainly with the spectrum of generalized dimensions Dq
characterizing the spatial structure of eigenstates on the length scales smaller than
the localization length. Before a detailed discussion of the results, an overview of
the main results of this part of the section is given. Scaling of IPN at the band
center is calculated first, and I Show that scaling properties (that is, whole spectrum
of generalized dimensions Dq) changes discontinuously near the band center, at an
energy E’(L) for a range of L studied. It is then shown that E’ can be quite accurately
identified with half of the width of the energy range around band center within which
two states occur on average in the ensemble of disordered systems. Existence of this
energy reveals the existence of a new length scale 5’ (E), diverging when E —> 0, that
is the system size at which IPN’s change their scaling properties from one power-

law dependence on L to another. This change is then explained as a finite size

54

manifestation of the critical point, and the critical exponent calculated by a finite-
Size analysis.

Calculation of Zq(E, L) starts with calculation of Zq(E, L, AE), which is just
IPN averaged over all eigenstates from an energy interval (E — AE / 2, E + AE / 2),
taken from N9 realizations of disorder, followed by studying the limit AE —+ 0 [50].
Results at the band center for system L = 80 and several different q’s, are presented
in Fig. 3.11. The error bars in the figure are taken to be the standard deviation of
the average value.

The figure suggests the existence of an energy E’ independent of q (and therefore
defined by the whole multifractal measure) such that decreasing AE below E’ does
not change Zq(E, L, AE) significantly. Decrease of Zq by a smaller extent, however,
is still present for AE < E’, and the main source of this is the mismatch between
average and typical value of IPN at a given energy. Thus, the effect should become
smaller as the number of disorder configurations that are averaged over is increased,
and Zq(E, L, AE < E’) z Zq(E, L) up to the corresponding statistical error.

This can be seen in Fig. 3.11 and 3.12, where values of Z,, for AE < E’ are
approximately constant within the error bars (as indicated by the horizontal dashed
lines), while for AE > E’ there is approximately a linear dependence of Z, on the bin
size AE. In this sense Figure 3.12 suggests that E’ exists for all the systems studied
(and can be shown to be independent of q for each of them analogously as shown in
Fig. 3.12).

Similar analysis is carried out for energies away from the band center and Zq(E, L)
determined accordingly, where it turns out that the convergence for these energies is
slightly easier to establish and occurs at larger AE than at the band center. Such
obtained energy intervals used in calculations of IPN for different E were small enough
so that there was practically no overlap among them.

The vertical dashed lines in Fig. 3.11 and Fig. 3.12 represent half of the energy

55

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

25500-""' I 35- ''''' - .— 1 f
' q = 0.5
L - 80
. 34 .
25000 '
Z
q I 33 -
24500 » T1
. l_ ,7
I 32
t
24000 - »
0.105 ~ q = 1.5 0.022 - q - 2.0
b T L-80 un- .... ........... L—8o
0.100 1111111111 0.020 "
Z 0.095 0,013 . ‘ 'I’
‘1
0.090 - 0.016 -
0.035 - 0.014 -
10'5 10‘4 10‘5 10'4
AE/2 AE/2

Figure 3.11: Dependence of the average IPN for several different q’s and L = 80 on
the bin size AE. The vertical dashed lines represent energies E2(L) / 2 (see text for
definition), while horizontal dotted lines Show the error bars of Zq(L, E, AE = E2(L)).

56

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.90- . _
L: 100 . L=90 ’
q=09. ' q=03 . l
E 1.85” -i
Z
q I
1'85” T ernr ' ’ ,. J.
L .
' l". TT'PI E
l l 1
. 1.80-
1301...... - - Mu.-. - -- ....Lu . . l-
1w.nnq .. .nm .. ..nm - fir r r,
’ 14:80 ‘ 1.30.. L=70
q=09 q=09
Z .
q ' d
l
1.80 - d . ........................................
9T W _ 1"
i " .
r . 1.75~
1.77....... . . 4 . . 1.73”“... . - ...Tm
1.76- L = 60 . 1.72 _ L = 50
q=09 q=O9
1.75- l
171- v".I’T ........
Z 1.74 wt; . ........................................... . .............
q ' "'1 1.70- l ..
1.73- .
172- 1.69.
1.71“““' ‘ ‘44---1 . . - 1.68“““' . . ....1‘1
10'5 104 10-6 1041
AE/2 AE/2

Figure 3.12: Dependence of the average IPN for several different L and same q = 0.9
on the bin size AE. The vertical dashed lines represent energies E2(L) / 2 (see text
for discussion), while horizontal dotted lines Show the error bars of Zq(L, E, AB =
E2(L))

57

E2(L), defined as the width of the energy interval around zero in which every system
from the ensemble of disordered systems has two states on average [50],
0
1 = L2 / pL(e) dc. (3.19)
E2/2

Ptom the figures one can see that E’ z E2 / 2 for all system sizes studied except L = 40
(the smallest system studied, not shown in Fig. 3.12), where the convergence seems
to be somewhat slower. Equation (3.19) defines a new length scale 5’ (E) that can be
described as the system size L for a given energy E such that the number of states
within the energy interval between —E and E is two on average. This can be defined

as

€’(E) = E;1(E), (3-20)

where E2"1(E) is the inverse of E2(L), which is defined by (3.19). Since E2(L) goes
to zero when L —+ 00, the new length scale diverges when E —> 0.

It is tempting to integrate results for m, from See. 3.3.1 to obtain 5’ (E) explicitly.
This, however, does not give the correct result since the whole analysis of DOS from
Sec. 3.3.1 is done for energies larger than the width of the distribution of two states
closest to the band center, which in turn defines 5’ in (3.20). In other words, there
is an energy cutoff, vanishing when L —> 00, below which the fits are not accurate,
most obviously seen by noticing that all of the assumed analytical forms of pL are
diverging at the band center, while the actual pL is not.

Nontrivial feature connected with the existence of the new length scale 5’ is that
scaling exponents rq(E) are different for L S, 5’ (E) and L 2, 5’ (E) This is a general-
ization of findings from Ref. [50], where a deviation from the power-law scaling of the
average participation number for different AE at the band center appeared whenever
AE exceeded E’ (L) in several models with chiral symmetry. It can be straightfor-
wardly shown that the same thing happens with scaling of Zq(E = 0, L, AE). This

suggests that a new scaling characteristic should be attributed to the two states clos-

58

 

 

* +1031: = .300
I —o—iogE = -3.25 ,
0.65 - +logE= -3.50 % .,
' + logE = -3.76

-3.6 -

   
    

U : —0—- logE = 4.00
N 4.2 . —+— logE = 4.25 .
.5 0°60 ' —X— logE = 4.50 4
C + E = 0.00 1
-4.8 -

 

0.55 ~ //
t A

40 50 60 70 80 90100 0'50 40 50 60 70 80 90100

L L

q=2.0

L

 

 

 

 

 

 

Figure 3.13: Dependence of the average IPN on system Size for several energies near
the band center for q = 0.9, 2.

est to the band center, and that E’ scales as E2, with a coefficient of proportionality
close to one. This also implies that the corrections to the constant Zq(E, L, AE) for
AE < E’ are small. Therefore, the length at which change of scaling occurs should be
close to and depend on the energy proportionally to 5’(E), while the change of scaling
should be a narrow crossover, as opposed to much broader crossover from power-law
to constant IPN that occurs at the termination of multifractal scaling (3.16) for L z 5.

This is compatible with the results presented in Fig. 3.13, which gives calculated
Zq(L, E), for q = 0.9, 2 for different values of L and E studied. Smaller q’s allow for
more accurate determination of the scaling properties, and it can be seen from the
data that for both q’s there are two characteristic scaling behavior — one occurs at
the band center and nearby energies for smaller L, while the other scaling holds at
energies further away from the band center, and for larger L at energies close to the
band center.

Scaling exponents rq are determined form the linear regression of the data for

system sizes L = 50,60, ..., 100. Goodness of the power-law fit is quantitatively

59

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-0.11 r """"""" Ti . 1.2 g """ i VVVVV .1
. q = 0.9 . 5 I I f a
-0.12 - - 1.0,- 1
-0.13_- f . O.8:- 2 0 I!
Tq -0.14: f 2 0.6:- q — °
-0.15 ~ i i 2 ‘ 0.4} f ‘3
-O.16~ ’ ' 0.2% i
. : § : .L : c : : c c : ‘ 1.5 f 4 : t 4. : : : : r : c ‘
F 1 t I
0.05} { ”E i i
: : 0.5:- ] I i -1
' _ ‘ E i i
Yq 0007—} i I : 0.0E 1 1 3
: ; -0.5Z- i
-005 . . : :
; ; 4.0;- } -‘

-1x10'3 -5x10'4 0 -1le3 -5x10‘4 0

E E

Figure 3.14: Scaling exponent of IPN, rq, and the goodness of fit coefficient 7.],
calculated from the data in Fig. 3.13 (except L = 40 results).

characterized by a coefficient 7q(E) next to the quadratic term from an additional
quadratic fit of the data. Such obtained r, and 7,, from the data in Fig. 3.13 are
presented in Fig. 3.14.

Results Show the existence of three different cases: (i) away from the band center,
7,, z 0 and 7,, is independent of E, indicating power-law dependence of IPN on L; (ii)
approaching the band center, 7., becomes different than zero, indicating that power-
law is not obeyed. This is due to the emergence of the new scaling for system sizes
L ,2 5’ (E) discussed above and present in Fig. 3.13; and (iii) for E = 0 power-law is
obeyed again (7,] z 0), but with a different 7,, than for energies away from the band
center.

Discontinuity of rq(E) for all the other q’s studied naturally leads to two different

spectra of generalized dimensions, and Fig. 3.15 Shows calculated D, for all energies

60

 

 
  
  

 

 

 

 

2.5 IIII IIITIIIIIIIIIIIIII—IIIIIIIIIIIIWTTI‘IIITIIIIII-
----- ~---- log|E| =-3 .
I ------« log|E| =-3.252
2.0 - A
: ---'v---- log|E| =-3.5 :
: ......... log|E| =-3.75:
1-5 r "r" loglEl =-4 1
D - .
q ; I
1.0 - :-
05:- -
i
0.0 1111 11111111111111111111111111111]11111“In 1111

0 1 2 3 4

‘1

Figure 3.15: Multifractal spectrums Dq at the band center and nearby energies.

except the three nonzero energies closest to the band center (which cannot give 7},
from the fitting procedure used here due to the change of scaling properties discussed
above).

In particular, the participation number grows with the number of sites L2 as a
power-law with exponents ,B(E = 0) = 0.25 :L- 0.02, and D(E aé 0) = 0.55 :l: 0.05. This
should be compared with the result of Ref. [45] for the same ABD model, D(E = 0) =
0.50:l:0.06. It is easy to explain this discrepancy since the bin size AE = 4 x 10‘4 used
in Ref. [45] was, depending on the system size, roughly over an order of magnitude
too large to detect the correct scaling behavior, and therefore 3 (E as 0) was obtained
instead.

Dq is calculated for only three negative values of q, mainly because IPN for nega-
tive q’s are determined mostly by the parts of eigenstates with the smallest probability

to find the electron, which in turn acquire the highest relative error during numerical

61

diagonalization of the Hamiltonian. Difficulties in calculating Dq in this regime even
arose suspicion that multifractality might break down for negative q [73]. It is thus
important to Show that D, is defined for negative q’S as well as for positive ones.
Accuracy of all the calculations of IPN done in this section can be straightforwardly
improved by increasing the number of configurations of disorder that was averaged
over. This would lead to smaller error bars of all the quantities calculated, as well as
to the wider range of q’s for which Dq can be calculated.

Results of this section give the following picture of the scaling of IPN with system
size: for any energy E close enough to the band center, there exist two power-law
scalings: one for L0 < L S, 5’ (E) described by the set of exponents rq(E = 0),
and another one for 5’(E) S, L < 5 (E), described by a difi'erent set of exponents
rq(E 75 0), which leads to the two different spectra of generalized dimensions Dq.

Dependence of the new length scale on energy, 5’ (E), can be easily determined
from (3.20) by integrating the actual numerical data for pL(E), and the result, ob-
tained from Fig. 3.16, gives 5’(E) oc [E]‘"’, with

u’ = 0.35 i 0.01 . (3.21)

The meaning of this new length scale and corresponding exponent u’ can be under-
stood by assuming that the additional scaling of the two states is due to the finite size
effect of “smearing” of the E = 0 critical point of the infinite system, because exactly
the two states closest to the band center become critical when L —» 00. The critical
energy then changes by A5 which is, in a finite system of size L (in units of the lattice
spacing), equal to AE2(L) o< L“‘5 = L‘l/V' . If furthermore 5 oc ]E|"" in the infinite
system, the shift Ace of the critical energy in the system of size L is Aec oc L‘l/V,
from the general theory of the finite size scaling [3]. Therefore, V’ equals V, the critical
exponent of the correlation length of the infinite two—dimensional system. It should

be noticed that the finite size scaling applied here is somewhat different than usual,

62

 

'U T I r 'I‘ITI I

3 0C [El—035(1)

 

100,

 

 

 

10'3 10’4 10‘5

Figure 3.16: Dependence of the new length scale on energy near the band center, as
determined from the Eq. (3.20).

where Aec E [E — Eel/EC [3]. Here, E6 = 0 and A56 E [E — E], where both energies
are expressed in units of 60, as discussed in the Introduction. Thus, energy 60 appears

in the denominator of Ace rather than the critical energy itself.

3.3.4 Localization length near the band center

In order to calculate the localization length of the infinite two dimensional system,
finite-size scaling analysis of MacKinnon and Kramer [17] (FSS) was applied to the
transfer matrix method (TMM) of Pichard and Sarma [47].

The analysis is based on rewritting the Schrb dinger equation H tb = Ed) of the

system on a strip of length L and width M in the following form:

’1)“ =7". ’1)" ,n=1,...,L; (3.22)

wn 7)bn-l

63

where #2,, is M dimensional vector representing the part of 112 at the n-th “slice” and

t;,1.+.(E — H.) —t;1.+.t‘

n—1,n

Tn E ' , (3.23)
1 0

where H, is an M x M dimensional matrix representing the projection of H at n—th
Slice, tum“ is a diagonal matrix of the same dimensionality with h0pping elements
tij from the n-th to the n + 1-st slice (due to Hermicity of H, tn," = tfnm).

In the limit L —> 00, the product
1/2L
1“,, 2 (TL . - .TIT,’ - - erg) , (3.24)

is a self-averaging matrix of dimension 2M x 2M, with eigenvalues {exp(i)\M,m), m =
1, . . . , M}, and AMm’s are called Lyapunov exponents [10]. The smallest of these
exponents, AM, corresponds to the inverse largest characteristic length scale of the
system.

The scaling analysis consists of assuming that the change of the largest ILE,
A(E, M) _=_ 1//\M(E)M, due to rescaling M —> bM can be compensated by an ap—
propriate change of energy, after which A will remain the same. This implies that
[17]

A(E, M) = A(M/5(E)). (3.25)

Figure 3.17 shows the calculated A(E, M) for several energies E close to the band
center and for various strips up to 128 sites wide. The figure also gives the second
largest (dashed line) renormalized ILE, A2, for the two lowest non-zero energies stud-
ied (E = 10‘5 and 10‘s) and for M 5 20. All values are obtained with relative error
of 1% or better.

At E = 0 and for M even, all ILE become doubly degenerate due to the presence
of chiral symmetry [46, 49], and they scale linearly with M for M 2 16, reflecting the
scale invariance of A characteristic of a critical state. To see what effect this degener-

acy has on scaling prOperties of ILE, we should recall that ILE of transfer matrices of

64

 

 

 

 

 

2 I I r I I I I I
1 — ..
m“ .-
2“ t '
- -0- - A2( 5:10-6, M) .\0
- -A—- — A2(E=10-5,M) ~ '
E = 10-1
0.2 . . m1 . l 1 I l m 1 . I I 1 l 1 l
l 10 100

Figure 3.17: The largest renormalized inverse Lyapunov exponents A(E, M), calcu-
lated for long stripes of width M and various energies near the band center (full line).
The energies range from loglE] = -1,—1.5, ..., -5, and additional log [E] = —6, as
well as E = 0. Dashed lines connect points of the second largest ILE, A2(E, M), for
energies E = 10‘5,10’6, and M _<_ 20.

65

disordered systems repel each other in general [13] and become self-averaging quan-
tities for sufficiently long stripes [47]. The symmetry now enforces “dimerization”
of pairs of ILE, acting as an effective attractive force between each pair of ILE that
become degenerate at zero energy, as in Fig. 3.17, where, for a fixed M, A and A2 are
closer together for the smaller energy. The largest ILE thus decreases in a strip of
width M on approaching the band center, as in models with chiral symmetry studied
in Ref. [46]. On the other hand, at any given energy, every pair of ILE becomes
more repelling with increasing of M. Such increased number of ILE thus diminishes
the attractive effect of the symmetry, and, depending on the relative strengths of
the two effects, there are two regimes: (i) for smaller M, the largest ILE increases
approximately linearly with M, and, since the slope of the rise changes with energy,
finite-size scaling (FSS) cannot be done; (ii) for sufl'lciently large M, on the other
hand, the repulsion among ILE due to disorder dominates and FSS is possible.

The one-parameter universal function A(M/5(E)), obtained for system sizes M 2
50 for energies [E] > 10'4, and for system sizes M 2 64 for energies 10"5 S [E] 3
10'4, is presented in Fig. 3.18. The obtained localization length 5 (E), in units of the
localization length of the smallest energy studied (E = 0.1), is shown in the inset
of the same figure. Calculated 5 (E) /5 (0.1) is fitted to the power-law for energies
10‘5 S E S 10'3 , and the result suggests power-law diverging localization length at

the band center, with the exponent
u = 0.335 :1: 0.034, (3.26)

in agreement with the result (3.21) of Sec. 3.3.3.

Some additional analysis of the results can be done by introducing Am(M),
defined as the maximal A(E, M) for a given strip width M. Importance of this
quantity comes from the fact that points where A reaches its maximum for various

energies cannot be described by FSS, and, at the same time, Am” limits by its

66

 

 

 

  

 

 

 

 

 

 

1.5 1 T r T 1 I I 1' I I I T I l I ll
AA
A
“.8
1 ... _
- °¢ .
... ...... ..."... ...... ....“ . 3
A _ ' ‘30 ‘q.
’ gm [ EI-0.335(15) ‘ 9+
. 10 a l
£15 on
- — ’5 i
g ’a
0.1
fl
1 :1
10—5‘ ”‘10-4 “10'3- “-10'2‘ A lO‘l‘
O 2 1 1 1 El 1 1 1 1 I 1 1 1 1 I I I 1 I
1 10 100

M/E

Figure 3.18: One parameter universal function A(M/5(E)) for the ABD model. The
four separate points are A(E = 0, M) for M = 16, 32, 64, 128, corresponding to the
band center critical state. The inset shows the calculated localization length 5 (E) in

units of the localization length 5 (E = 0.1), with error bars smaller than the symbol
size.

67

 

definition possible values that Mr) (obtained from the F SS) can have. The main
observation is that Am“ seems to grow slower than linearly in Fig. 3.17. Linear (or
faster) growth would imply existence of an additional energy scale below which FSS
would break down for all large M and A(E, M) would grow linearly and indefinitely
with M in the same figure, implying the existence of a whole band of extended states
where one-parameter scaling would not hold. Slower than linear growth of Am seen,
therefore, furthermore suggests the existence of a single critical point at the band
center of the system on an infinite square lattice.

The obtained value of the critical exponent V < 1 seems to contradict a rigorous
theorem of Chayes and Chayes et.al [74], which states that the critical exponent of
the Anderson transition must satisfy condition V Z 1. The contradiction, however, is
most likely only apparent, because the theorem requires the existence of both localized
and delocalized phase separated by a mobility edge, while in the case of ABD model

in d = 2 the delocalized phase does not exist, only localized phase and a critical point.

3.4 Random flux model

Problem of a quantum particle moving in a random magnetic flux has attracted con-
siderable interest due to its relation to the problem of high Tc superconductivity
[75, 76, 77], theory of the half-filled lowest Landau level [78, 79, 80], as well as its un-
resolved localization properties within the context of the scaling theory of localization
(Sec. 2.3).

The model of a quantum particle hopping on a square lattice studied here is de-
scribed by the Anderson model (2.1) in which the hopping amplitudes t,,,- E exp(i¢,-,j),
where the phases 05,-0- are chosen such that the total flux per plaquette, (I),- E (-p7r,p1r),
is a uniformly distributed random variable, parameterized by 0 < p S 1. The uni-

formly distributed 6, E (-W/ 2, W/ 2) representing the on—site disorder are set to 0

68

unless explicitly stated differently. This is the random flux model (RFM).

Localization prOperties of the RFM remained only partially understood despite a
significant effort of the community: and three incompatible conclusions (with varia-
tions) have emerged: (i) all states are localized and the model is in the unitary class
[81, 86, 88]; (iia) there is a critical point at the band center of the RFM [46, 90] due
to the chiral symmetry [36, 37, 35, 38, 39] of the model (discussed below), while all
the other states are localized, (iib) with an additional Kosterlitz-Thouless-Berezinskii
(KTB) transition into a phase of critical states around the band center for 0 < p < 1
[89]; or (iiia) a full metal-insulator transition [82, 83, 84, 85, 87], (iiib) even in the
presence of (sufficiently weak but finite) on-site disorder [87].

This dissertation presents some further evidence that the the scenario (iia) is the
correct one and that the localization length 5 diverges as ]E["’ at the band center
with a critical exponent V, which has been calculated for the first time for p = 0.5
and 1.

The main source of difficulties in numerical studies of the RFM is very large 5,
growing approximately exponentially from the band edges toward the band center
[81, 86, 90]. The standard numerical method to determine the scaling of 5 is from
the finite-Size scaling analysis of either the smallest Lyapunov exponent AM of the
transfer-matrix of a long strip of width M (TMM), or the logarithm of the norm
of retarded Green’s function connecting the two far ends of the strip (GF) [18]. In
numerical studies of systems with small M, the exponential rise of 5 may give a false
impression of a genuine continuous phase transition, as illustrated by Xie et. al [89].
Particularly difficult is to distinguish between the scenarios (i) and (iia), because the
eigenstate at a given energy with very large 5 in a finite size system looks like a critical
state even if it is localized in the infinite system.

Another approach was to study the Chern number of eigenstates, by analyzing the

boundary-phase averaged Hall conductance in the system with generalized periodic

69

boundary conditions [91, 92]. The conductance is quantized by an integer, the Chern
number, and states with zero Chern number are localized. The necessity of averag-
ing the conductance over boundary conditions limits numerical studies to systems
of smaller sizes than in other methods, although one may expect that the averag-
ing somewhat reduces the boundary effects due to the finite Size. Sheng and Weng
[85] have given numerical evidence of scenario (iiia) using this method. Additional
calculations on larger systems were done by Yang and Bhatt [87], who studied the
total number of states with a nonzero Chern number, NC, as a function of L and W,
and then used the scaling analysis to determine WC. They found a critical disorder
strength, E; > 0, such that for 0 S W < Wc a band of states with nonzero Chern
number opens up around the band center. What is, however, puzzling about this
result is that the calculated NC, even deep inside the ”metallic” region, appear to
grow as N a, with a < 1, approaching 1 when W —) 0. Hence, numerical data from
Ref. [87] seems to be compatible with the conclusion that the ratio NC / N goes to zero
as N -—> 00, suggesting that there is no band of non-zero Chern number states for any
W Z 0 in the infinite system. This is indirectly supported by the TMM calculation
of Xie et al. [89] who found that the states are localized at the band center for W
significantly smaller than We.

The most extensive numerical study so far of the RFM was done by Furusaki
[90]. He found the one-parameter scaling to be obeyed for [E] > 2.55, all states
localized for [E] 2 0.1, and gave some further evidence for the scenario (iia) and
against the scenario (iib). Let us therefore first investigate the scaling properties of
AM 5 l/AMM in the energy region [E] g 0.1 and [E] < 0.1. The TMM is applied
and scaling of the A M is studied for the RFM in a gauge where all the complex phases
are zero along the strip of width M and length up to 2 x 107 Sites, the largest systems
studied so far.

Furusaki [90] found that AM decreases with M for [E] > 0.1 [90], as typical for

70

 

TjIIIII I I IIIIIII I I IIIITII

3- . a ., .. .
. 5 . . —A—Iong=-1.0

—v-— log10 E = -1.5

AM —<>— loglo E = -2.0
-<l— log10 E = -2.5
2 —l>— log10 E = -3.0

—0— log10 E = -3.5
—*— loglo E = -4.0

 

—o— 1/M AM(E = 0)
-----o 1/MAM.1(E = 0)

 

 

1 1111111 1 1 111111l 1 1 lllLLJl

10 100 1000
M

 

Figure 3.19: Largest renormalized inverse Lyapunov exponent, A M, of the transfer
matrix of RFM model near the band center. AM initialy increases with strip width
M but eventually start to decrease for sufficiently wide strips, indicating localized
behavior.

71

localized states. Figure 3.19 gives the result of my extensive calculations AM for
E = 10”, 10'”, ..., 10‘4,0, and M = 4,...,128. For E z 0.1, there is a sudden
increase of A M for small system sizes before the decrease. This, previously thought
of as being a finite-size artifact [90], marks the beginning of a different scaling regime
that is characteristic for the chiral disordered systems near the band center. In this
regime, AM increases with increasing of M, implying the existence of some sort of
extended states in the infinite system. Figure 3.19 suggests, however, that this rise
is not indefinite and that, for sufficiently wide strips at any given E, AM reaches
maximum and then decreases with further increase of M, and numerical evidence
for this behavior is seen for [E ] 3 10’3. This is very reminiscent of the TMM
calculation from Sec. 3.3.4 for uniformly distributed real t,”- E (—1/2,1/2), W = 0
model (ABD), where a similar asymptotic behavior was found. Due to much smaller
localization lengths for small [E] in the ABD model, it was possible to analyze in
more detail the behavior of AM when M is increased. In this case a restoration of
one-parameter scaling, A M(E) = A(M / 5 (E)), was found for sufficiently wide strips,
leading to 5 ~ ]E]“’, with V = 0.335 :l: 0.015. Performing analogous calculation in
the case of RFM, however, would require about two orders of magnitude wider strips,
which seems to be beyond the present-day computational means.

Despite this severe limitation, an alternative approach from Sec. 3.3.3 will be
used for calculation of V of the RFM model. The link between the two models is
the already discussed chiral symmetry. As already mentioned, V has been obtained
numerically for the ABD model in the Sec. 3.3, using TMM (Sec. 3.3.4) as well as
from a finite-size scaling analysis related to the two states closest to the band center
of a system on a square lattice of linear size L. Since they become zero energy states
when L —> co and therefore critical, it was conjectured in Sec. 3.3.3 that one should
be able to apply a finite size scaling analysis to the new length scale 5’ obtained from

the width E2(L) of the energy interval that they occupy in a system of size L, and

72

 

I I I I T I I T V II

100-" L
r °° , 0 p = 1.0, V = 0.481(4) ’
.. . El ,9 = 0.5, v = 0.494(5) ‘

L‘ O T

 

 

20’...1 1 .. ...1 3"”
103 I? 10'

Figure 3.20: Scaling of the 5’ length as determined from the Eq. (3.19), for p = 0.5
and 1. Power-law is obeyed well in both cases, and the V’s are determined from two
linear fits, represented in the figure by lines.

from there obtain V.

Figure 3.20 presents the calculated 5’ for p = 0.5 and p = 1, W = 0, on a square
lattice with periodic boundary conditions for various even L between 20 and 90 with
the number of configurations of disorder ranging from 2 x 103 to 100, respectively.
The power-law is well obeyed for 1.5 decades in energy, and the extracted values
of the critical exponent V are given in the figure. Although the difference between
exponents is small and both are close to 1 / 2, the values nevertheless support the
conclusion that Vp=1 < Vp=0,5 < 1 / 2, implying also a weak dependence of V on the
disorder distribution. The second inequality comes from the fact that pL depends
on L near the band center [90]. In the ABD model this dependence is much more

pronounced [93], so that V is even smaller.

73

We now discuss the possibility of the KTB transition (scenario (iib)). Xie et al.
[89] have found from TMM calculations that AM exhibits a scale—invariant behavior
around the band center for 0 < W < Wc(p), leading to the conclusion (iib). On
the other hand, Furusaki [90] has given numerical evidence that the variance of the
conductance distribution of RFM quantum wires follows one of the unitary class for
W > 0 and a different one for W = E = 0, emphasizing the importance of the
chiral symmetry and implying that all the states are localized for W > 0, and for
[E] 7E 0,W = 0, for p 2 0.2.

Figure 3.21 gives further evidence that there is no KTB transition (iib) in the
RFM model by redoing one of the calculations from Ref. [89] with higher precision.
The figure shows AM at the band center of the RFM, p = 0.5 model for different
strengths W of the on-site disorder, and M = 16,32, 64, 128. The length of the strips
used in the calculation is, however, 2 x 106 sites, which is about an order of magnitude
larger than that used in Ref. [89]. This gives more accurate values of AM and it can
be noticed that, for W 2 2.5, AM is systematically smaller for M = 128 than for
M S 64, as shown in the inset, indicating clearly localization at the band center for
W around 4, and thus the absence of KTB transition seen in Ref. [89].

Let us now turn to the analysis of the E = 0 states. In the ABD model there
exists an additional symmetry among Lyapunov exponents at E = 0 due to which all
of them come in degenerate pairs [93]. This property of chiral disordered models, first
noticed by Miler and Wang [46], is naturally expected to hold for the RFM model as
well, as it was expected by these authors. Figure 3.19, however, clearly shows that
AM < AM_1, implying that the additional symmetry at the band center is broken in
the RFM model, in contrast to the ABD model where the two A’s are equal to a high
accuracy.

This has its analog in the language of eigenenergies, where the two states closest

to the band center are repelled only in the case of RFM. This important difference is

74

 

 

 

 

 

 

 

 

 

0.011L1.1.1.1.1.4.1
7

Figure 3.21: Scaling of the largest inverse renormalized Lyapunov exponent at the
band center, for various strengths of the on-site disorder W. The inset shows mag-
nified region where the KTB transition (vertical dashed line) should be expected
according to Ref. [14]. Due to higher accuracy of the present calculation, it is clearly
seen that the transition is absent, since A M decreases with M, signaling localization.

75

 

Density of individual eigenenergies

 

 

 

 

 

I
-0.01 0.00 0.01
Energy

Figure 3.22: Distributions of 8 individual eigenenergies, enumerated —4, —3, ..., 4,
closest to the band center of the ABD and RFM systems of size L = 40 for 2 x 10"
configurations of disorder. Both systems have symmetric spectrum around E = 0
due to the chiral symmetry, but the central two states are repelled only in the RFM
system.

76

shown in Fig. 3.22, where the distribution of several eigenergies closest to the band
center is presented. The power-law behavior shown in Fig. 3.20 and discussed above
then suggests that the distribution of central two states may be scale invariant and
thus p1,...» equal to 0 in the case of the RFM, as opposed to a divergence in the case
of the ABD model.

In the language of eigenstates, the additional symmetry at the band center of the
ABD model is 0(2), since any linear combination of the two E = 0 states is also
an eigenstate at E = 0. In the RF M then the broken additional symmetry means
that two states closest to band center in the infinite system are separated by an
“infinitesimal gap”, which means that they have, as all other states E 74 0, equal
probability density at each of the two sublattices [93]. Nevertheless, it is still possible
to have two critical states at E = :l:0 in the infinite system, because the FSS from

Ref. [93] that was used to interpret the scaling seen in Fig. 2 implies
5(32) ~ E'lEzl N [E2l—V "V L- (3-27)

Thus, eigenstates that are not localized in a system of size L for [E] S, E2(L) will
not get localized when L is increased. Since these states exhibit multifractal behavior
[86, 93], they will remain multifractal as L —> 00, and this novel kind of the critical
behavior, present in both the ABD and RMF, we refer to as anomalous critical.

In summary, the RPM model has been investigated by a high precision TMM and
exact numerical diagonalization. The scenario (iia) was found to be valid in the limit
of an infinite 2D square lattice, i.e there are two anomalous critical states at the crit-
ical point E = :l:0, with all the other states localized. Also some differences between
the two chiral systems with and without time-reversal symmetry are identified and
discussed, most striking of which is the broken chiral symmetry at the critical point

in the absence of time-reversal symmetry.

77

CHAPTER 4

Monte-Carlo studies of cluster

magnetization

Studies of atomic clusters over the last decade have shown that there are significant
changes in the behavior of matter as one reduces the size from bulk to small clusters
[106]. The changes arise due to the finite size, preponderance of low coordinated
surface sites as well as geometrical arrangements which are different from those in
bulk. These can profoundly modify the electronic structure which affects various
properties. One of the most intriguing developments occurs in the magnetic behavior
[107]. For example, clusters of conventional itinerant ferromagnetic metals like Fe, Co,
and Ni have been found to exhibit superparamagnetic relaxations [108, 109, 110]. The
magnetic moments are generally enhanced over the bulk, but exhibit large variations
with size for smaller clusters [111]. It is also found that clusters of Rh [112, 113]
which is non-magnetic in bulk, and those [114] of Mn which has an antiferromagnetic
bulk phase, exhibit ferromagnetic order with large moments (ferromagnetic order here
refers to the ferromagnetic alignment of the atomic moments).

An interesting case of the novel magnetic order in clusters is that of Gd. Bulk Gd
is ferromagnetic with a moment of 7.55 [13 per atom [115]. A Gdg molecule has been

found to have a spin moment of 8.82 pg per atom with ferromagnetic coupling [116].

78

The Gdfl clusters containing 10 to 30 atoms were however found to have magnetic
moments between (2.94:E0.35—3.88i0.47)u3 per atom [117, 118], considerably below
the bulk value. A Gd atom has 7 unpaired f-electrons and one unpaired d—electron. To
obtain a moment of less than 6MB per atom in a cluster, either some of the f-electron
spins (in a given atom) have to be paired or the magnetic coupling between the atoms
modified. Note that the f-spins are unpaired in the molecule and in the bulk where
the exchange splitting is around 12 eV. It is therefore highly unlikely that they could
be paired in clusters. This puzzle was resolved by Pappas et al. [119], who argued that
like bulk Gd, small clusters are marked by indirect Ruderman-Kittel—Kasuya—Yosida
RKKY type exchange interactions [120] which oscillate as a function of the distance
between two magnetic atoms. Assuming a nearest neighbor ferromagnetic interaction
and a non-nearest neighbor antiferromagnetic coupling within a Heisenberg model,
they studied the spin configuration in Gd13 clusters. These studies show that for
a range of interaction strengths, the atomic spins assume configurations where the
surface Spins are canted (See Fig. 4.1). This leads to lower net magnetization and
accounts for the observed anomalously low moments of the clusters.

Experiments on Gd clusters also show another interesting new feature. Studies of
clusters as a function of temperature Show that the measured magnetization increases
with temperature. This feature has also been recently seen in Gd nanoparticles [121]
where the measured magnetization Shows an increase with temperature for certain
range of temperatures. These latter experiments on a range of particle sizes further
show that while the magnetization as a function of field saturates at high temper-
atures, it does not saturate at low temperatures. It was suggested [121] that these
effects could be understood by assuming that the spins near the surfaces were not in
alignment with the interior spins. While such a picture would be consistent with the
ground state spin configurations for Gd13 found by Pappas et al. [119], no study of

the temperature dependence of magnetization is available.

79

In this Chapter, we discuss our theoretical investigations of the effect of tempera-
ture on the magnetic properties of small Gdn clusters. The key issues we wish to focus
are: (1) how does the magnetization change with temperature and how is this related
to the details of the coupling? (2) at a given temperature, what is the distribution of
magnetic moments amongst different clusters? We Show that the presence of compet-
ing ferro-magnetic and anti-ferromagnetic couplings can lead to situations where the
magnetization initially increases with temperature and then decreases. This behavior
is linked to the excitations from the canted configuration to the ferromagnetic ex-
cited states. We also show that the clusters can have a wide distribution of magnetic
moments. This could lead to a broadening of the cluster beams in Stern-Gerlach
deflection experiments. We would like to point out that while our studies are carried
out on Gd13 clusters, we believe that the results have a more wider applicability.

The magnetic interaction on the cluster is described by a Hamiltonian H,

H=J _zs.§,+7§y§..§, , (4.1)

(52') (M)
where Sis are the spins of the Gd atoms, the prime Sign in the second sum denotes
summation over all pairs of Spins that are not nearest neighbors, and, 7 _>_ 0. J
determines the energy scale, and, thus, we assume J = 1 from now on. The quantity

that is calculated is the magnetization per atom M, defined as

If all the moments are oriented parallel to each other M = 1. As in previous studies
[119], the cluster is assumed to gave an hexagonal closed packed structure as shown
in Fig. 4.1.

For 0 S 7 < 0.29, the ground state is found to be ferromagnetic, i.e. all spins are
aligned in the same direction. When the value of 7 is further increased beyond 0.29,

competition between ferro— and antiferro-magnetic interactions leads to spin canting.

80

 

 

Figure 4.1: Gd13 cluster. Spins are oriented to have D3,, symmetry. Atoms belonging
to the lower, middle and upper planes are connected (except the middle one).

The symmetry of the ground state is consequently reduced from the symmetric S13 to
the canted C3,, symmetry [119]. The new symmetry, together with invariance under
simultaneous rotation of all spins, is still sufliciently high to allow very accurate
determination of the ground state. By rotating all spins such that the central Spin
is oriented along the z-axis of the cluster, and rotating all spins around the z-axis
until spins outside the middle plane are in vertical symmetry planes, the energy of
the cluster (4.1) can be expressed as a function of only three angles. First is the angle
6,, between spins outside the middle plane and the central spin, and other two angles
(0, d) that determine the orientation of one of the spins in the middle plane. Values
of other 21 angles are then fully determined from the symmetry. The energy of the

cluster in units of J becomes

H(6u,0, gt) = 3 [1 — 7 — (1+ 7) cos(2d>) + c0326 (—2 + 47 + (1 + 7) cos(2¢)) +
—2 c050,, + 3 (—1 + 7) cos2 6,, + 2 cosfl (—1+(—2 + 47) 0036,) +
—2 (1 + 7) sind (1/3 0034) — sin (t) sin 0”]. (4.3)

Calculation of the exact minimum of (4.3) is a simple computational task for any

given value of 7. However, in order to confirm the symmetry of the ground state,

81

we have minimized (4.1) first without assuming any symmetry, by using simulated
annealing [122]. N ear-global minima were obtained and compared with the results of
minimization of (4.3).The latter always had slightly lower energy, while the former
always had almost 03,, symmetry. Thus, the symmetry of the ground state can
be understood as a fine modification of the near-global minima into the exact one.
This holds when the ground state is ferromagnetic as well, and, it can be used for
an accurate determination of the transition point between ferromagnetic and canted
ground states.

We notice that if the ground state had only 03,, symmetry, additional angle 6.1
would have been needed to define the angle between one of the lower three spins and
the middle spin, while 0,, would have been used only for the upper three Spins. We
have assumed that ground state has 0,, = 6,1. We have ensured that this assumption
is correct by minimizing H(0u,0, (15,21). This means, however, that the full symmetry
of the ground state is D3,, instead of C3”.

When 7 > 0.37, the ground state remains canted but the symmetry changes
to a mirror symmetry defined by one plane (Oz group), and, the magnetization of
the ground state changes discontinuously at the transition point [119]. Number of
relevant degrees of freedom (angles) now becomes 12, and, therefore, finding exact
ground states by using symmetry considerations is a considerably more difficult task

than in the previous case.

4. 1 Monte-Carlo method

Let classical system be described by a set of N parameters 8 E {51, ..., SN} and a
Hamiltonian H (S ) We are interested in calculating the thermodynamic average of

an observable 0(8),
(0) E N [as 0(5) exp(—5H(s)), (4.4)

82

 

where B = (103T)‘1 is the inverse temperature and the normalization constant N E
Z ‘1, Z being the partition function.

In the problem of magnetic clusters studied in this dissertation, S,’S are classical
spins parameterized in polar coordinates by the two angles (15,- and 0,: representing the

orientation of the Spin at the i-th lattice Site,
Sf = cos 6,, Sf = sin 6,- cos d,, Sf = sin 6,- sin at, (4.5)

which defines the measure d8 E dSl...dSN.

To calculate (4.4) we consider a random sequence of states (8M, ...,81), defined
by a joint probability density P(Sk[5k_1, ..., 81) to have a state 8,, after the sequence
(Sk_1, ...,81) of k — 1 previous states. If P is the function of only two last states in

the sequence,

P(SkISk-1,-~-,51) = P(SkISk—I): (45)

the sequence is called the Markov chain.

The essence of the Monte Carlo method is to choose the transition probability
density P(Sk]8k_1) such that the probability density P(Sk) of all individual system
configurations is equal to the wanted density f (8),), independently of the choice of
the initial configuration 81. In the spin problem, f is the Boltzmann distribution of
Eq. (4.4).

Once we have such a Markov chain M = (SM, ..., 81), the average (4.4) is easily

computed as
M
1
(0) = M k§=1:0(5k) + AOM, (47)

where the last term is the error of the first term, which, as will be discussed below,
goes to zero as M —> 00.

In order for M to exist, the first condition that M must satisfy is the ergodicity,
that is, for any finite volume 52 of the configuration space, the probability p(S;c E

9) > 0 for some k. If this is not satisfied then (I is never visited by the chain, and

83

it is therefore impossible in general to estimate contribution of system configurations
from Q to the average (0).

P(Slek-1), being a probability density, satisfies the following identity:
/ dsP(S|s,.) = 1, (4.8)
which can be rewritten by inserting (1 — 6(8 — 81,) + 6(8 — 81.)) under the integral as
f 480 — 6(8 — 8.))P(8Is..) = 1— P(sklsk). (4.9)
Also, from the definition of Markov chain (4.6), the probability distribution of 8;, is
P(s.) = / dSP(S)P(Sk|S), (4.10)
which after inserting (1 — 6(8 — 5k) + 6(8 — 8;,» under the integral can be written as
P(Sk)(1— P(Sklskll =fd5P(5)(1— 5(3 — Sk))P(Sk[S)- (4-11)
Inserting the expression for (1 — P(SkISk» frOm Eq. (4.10) into Eq. (4.11), one gets
/ d8 {P(sanIS.) — P(S)P(Sk|8)} = 0. (4.12)

which is the condition of detailed balance.

Now we can obtain an equation for the transition probability by requiring that
P(Sk) = fish),
/ 48 05010515.) — mama} = o. (4.13)

4. 1 . 1 Metropolis algorithm

Sufficient condition for Eq. (4.13) to be satisfied is to set its integrand to zero:
f(Sk)P(S[Sk) = f(S)P(Sk[8)° (4-14)

The choice for the transfer probability density

 

puns.-.» = min (1, $32)) , (4.15)

84

is known as the Metropolis method for sampling of an arbitrary distribution f.
Proof that this choice of the transfer probability indeed satisfies the condition of
details balance can be done by a direct insertion of Eq. (4.15) into the Eq. (4.14) (or
Eq. (4.13)). This gives
f(8k) min (1,%) = f(S) min (1,;(3?) . (4.16)
If f(8) Z f(Sk), the Eq. (4.16) becomes f(Sk) = f(S)f(Sk)/f(8), while when
f(S) < f(8k) the Eq. (4.16) becomes f (8),) f(8)/ f(8;.) = f(8), which proves that

 

the condition of detailed balance is satisfied.

It is now possible to formulate the Metropolis algorithm that generates M. Start-
ing from an arbitrary initial configuration 81, at a k-th step (k 2 1), (i) a new
configuration 8’ is chosen, (ii) if the ratio f (S) / f (8),) is greater than 1 or smaller
than a randomly chosen number r between 0 and 1 then the new configuration 8H1

is chosen to be 8’, (iii) otherwise the new configuration 8H1 is chosen to be 8),.

4. 1.2 Thermalization

The Metropolis algorithm will therefore generate a sequence of configurations with
wanted probability distribution f for any ergodic transfer probability P(Sllsg), in
the limit of infinitely many computational steps performed. The problem now is how
to ensure that the convergence of the Markov chain toward the desired probability
distribution is achieved after finitely many steps of a particular computer simulation.

This issue is commonly referred to as the problem of thermalisation. In general,
there is little we can tell about the speed of convergence of the distribution of sam-
pled configurations toward the desired ones. Near the second order phase transitions,
for instance in the case of the Ising model, the necessary number of sampled con-
figurations can be exponentially large in the system Size if the transfer probability

P(81 [82) is chosen poorly, for instance if the new configuration is chosen by flipping

85

 

of individual spins. In such cases the cluster algorithms, where whole sets of spins are
simultaneously updated, can achieve exponential gain in the speed of convergence.

In simulation performed in this work, the spin updates are performed by small
random changes of angles 6,, d,- of a randomly chosen i-th spin. The amount by which
angles are changed is chosen such that in approximately 50% of changes the new spin
configuration is accepted. This ensures that the newly chosen spin configurations
diffusively explore the configuration space, but since the updates are local, the spin
configurations chosen in the process remain correlated.

Thermalization was monitored by recording the distributions of all 26 angles for
each of the configurations one system goes through at the given temperature during
the simulation. Due to the finite size of the system, global SO(3) symmetry of spin
configurations cannot be broken, and, thus, every global orientation of the spins has
to be equally probable once the thermal equilibrium has been reached. Furthermore,
since the spin updates are local, the uniform distribution of all of the Spin angles is a

good indication that the thermalisation has been achieved.

4. 1.3 Parallel tempering

Apart from carefully choosing transition probability, further improvements in the
convergence properties of Monte-Carlo simulations can be made by analyzing the
numerical data collected in the course of simulation [124], or even constructing several
Markov chains and combining the information obtained from each of them in choosing
new members of chains.

One such technique is used in this dissertation to calculate finite temperature
magnetization, that of the parallel tempering Monte-Carlo (PTMC) method [123,
124]. The algorithm consists of applying usual Metropolis MonteCarlo (MC) steps
(spin flips in our case) on a set of systems at different temperatures fir, followed by

a sequence of additional swaps of the i-th and the j-th system with a probability

86

 

 

 

 

 

0.55- -
- 7:040
0.547 —O— 106 PTMC
- —0— 5x106 PTMC
053‘ —+— 107 PTMC
‘ —o— 108 PTMC
0'52" —x— 106 x 10 MC
M ‘ —A— 107 x 10 MC
051- .
0.50-
0.49-
0.48 . . . . . . . , .
0.0 0.1 0.2 0.3 0.4 0.5
T [kB/J]

Figure 4.2: Results of PTMC and MC simulations for various number of time steps. In
the case of MC simulations, the result is an average over ten independent Simulations.

min(1,exp(—(,8j — ,8,)(E, — E,))). The swaps are performed sequentially between
two systems at the neighboring temperatures, beginning from the pair at the high
temperature end.

Systems at higher temperatures explore more system configurations in the course
of simulation, and, form the point of view of a system at some lower temperature, a
swap allows the latter system to make a jump from one configuration to the another
one with similar energy but which is much harder to access through usual MC steps.

As a result, the system samples configuration space more efliciently.

87

4.2 Magnetization of Gd13 clusters

4.2.1 Results of the numerical simulations

Figure 4.2 Shows M (T) for several MC and PTMC simulations. In the case of MC,
results are averages over ten independent simulations. We have used these results to
estimate error, which is about 0.005 for low T, and is smaller for higher T.

The temperature dependence of magnetization is given in Fig. 4.3, and distribu-
tions of magnetization for three different values of 7 are presented in Fig. 4.4. Results
are calculated performing 107 PTMC steps per point. When the ground state is fer-
romagnetic (7 < 0.29), magnetization decreases with the temperature, as expected.

Once the ground state has only D3,, symmetry, magnetization increases with tem-
perature when the cluster is heated from T = 0. This can be understood by noting
that the thermal fluctuations lead to excitations of the system to other configurations
with higher magnetization than the ground state. For 7 = 0.36, highest relative in-
crease in magnetization (compared to the ground-state value) is about 13%. Similar
behavior occurs when the ground state has mirror symmetry, for the lower values of
7, but the magnitude of this increase is smaller, due to reduction in the level of com-
petition between ferro and canted (low moment) configurations which get thermally
excited. Finally, for large values of 7 this effect disappears, and the magnetization
again slowly decreases with the temperature.

As pointed out above, the observed low magnetic moments of Gdn clusters were
earlier explained as due to canting of surface spins. Using the experimentally mea-
sured magnetic moments on Gd.“ clusters and the calculated magnetization as a func-
tion of 7 in [119] , one obtains a value of 7 in the range 0.35—0.4 . It is interesting
to note that for these values of 7, Fig. 4.3 shows that the magnetization increases
with temperature. This is in agreement with the observed increase of magnetization

seen in Gdn clusters [117, 118]. Note, however, that the quantitative change depends

88

 

 

 

 

 

 

 

 

 

 

 

M fi‘ I "I

I
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
T [kB/J]

Figure 4.3: Magnetization vs. temperature for different values of 7. Difference in 7
for two neighboring curves is 0.01. The error in the low T regime is of the size of the
circles (0.005) and is smaller for higher T.

89

on 7 which depends on the electronic structure. This shows that the temperature
dependence would depend on the cluster size.

What is most surprising is that the clusters have a distribution of moments (see
Fig. 4.4) at a given temperature. The experimental measurements of moments in
clusters are carried out by passing size selected clusters through the gradient field
in Stern Gerlach magnets [111]. Variation of moments amongst clusters would lead
to broadening of the beams. The effect would be pronounced if the passage time
through the magnet is shorter than the relaxation times for the moments [111]. This
shows that the earlier interpretation of the observed broadening of beams in Stern
Gerlach experiments [111] as being due to variations in cluster temperatures may not
be necessarily correct.

To summarize, we have shown that the observed increase in magnetization in Gd
clusters can be understood using a Simple classical Heisenberg model with competing
interactions. The variation of magnetization with temperature are linked to the na-
ture of the excited states and the presence of ferromagnetic excited states can lead to
the observed increase. We would like to point out that although the present studies
are carried out on Gd13 clusters, the results are applicable to other systems as well.
For example, Hehen et al. [125] have recently obtained an array of cobalt dots on
A1203 substrate and have observed the presence of canted configurations including
vortices. This is similar to the case of Gd clusters and it will be interesting to study

the temperature variations of magnetization in these magnetic arrays.

4.2.2 Results of experimental studies

Gerion et. al [126] have studied temperature dependence of magnetization of Gdn
(n = 13, 21, 22) clusters, and this section is devoted to discussion of their experimental
findings.

The clusters are created using the laser evaporation, thermalized in the nozzle

90

 

 

 

7:025

'I'I‘lrl'l'l

 

S3

o

cc
'l‘l'l'l'l'l
III'IIITI'

 

 

 

'ITI'I'I'I'T

 

 

bI'III'I'I'T'Il

 

Figure 4.4: Distribution of magnetization for three different values of 7. Lowest
temperature is 0.025, and, the temperature step is 0.025. Broader curves correspond
to higher temperatures.

91

at temperatures T = 70—600K and their colimated beam studied in the Stern-
Gerlach experiement in magnetic fields B up to 9.8kG. All the information about
the cluster are then extracted from the study of deflection patterns of the molecular
beam. In the case of Gd13, the authors find that the measured averged magnetization
M is proportional to B. This suggests the superparamagnetic behavior (SP), which
means that the measurement duration is long enough so that the cluster explores the
whole phase space during the interaction with the external magnetic field B. When
[18 << kBT, where u is the magnetic moment of the cluster, the SP is characterized

by

a:
kBT’
where M is the projection of the cluster magnetization in the field, M = [M - B] / B.

M = (4.17)

The linear dependence on B found by the authors therefore allows one to determine
experimentaly the cluster magentization per atom u/N.

The obtained p. for the Gd13 is 5.4:E0.3 in the units of Bohr magneton [13, which is
considerably smaller than the bulk value value 7.55, which suggests that some sort of
spin canting takes place in the cluster. Furthermore, observed a slowly decreases with
T for higher temperatures, which is in agreement with the results of MC simulations
preseneted in the previous Section of the Heisenberg model (4.1) of Pappas at. al
[119]. Furthermore, Gerion et. al [126] were able to extract 7 = 0.365 from the ratio
of u at low-T and the one caluclated in Ref. [119], and fit their experimental findings
with the results of the simulation from the previous section by adjusting J (Fig. 3(b)
of Ref. [126]), which was found to be 4.3 :l: 0.4meV. This value is larger that the bulk
one, abulk z 3meV, in agreement with somewhat larger Curie temperature observed
in the cluster compared to the bulk [126].

We notice that the MC numerical simulations provide the distribution of cluster
magnetizations as well as p(T), and it Should be possible to relate the former more

directly to the measured distribution of intensities. To this end, and in the light of

92

 

the experiment from Ref. [126] discussed above, it seems that MC Simulations should

include also the external magnetice field of the Stern-Gerlach experimental setup in

the model (4.1).

93

 

CHAPTER 5

Conclusions

In conclusions, main contributions of this dissertation are summarized. These findings
were published in three papers [50, 93, 127], and in another one that is about to be
published [128]:

Properties of several models from the class of quantum disordered systems with
chiral symmetry on a two-dimensional square lattice are studied using numerical sim-
ulations. The first one, the Anderson bond disordered model (ABD), is known to have
a logarithmically divergent density of states (DOS) and the parameter determining
this divergence is calculated for the first time in this dissertation. It is also shown that
the known asymptotic expression of DOS is inadequate for describing of divergence
of the infinite size system and a new asymptotic expression is proposed. Multifractal
analysis associated with the scaling of inverse participation numbers is carried out
for eigenstates near the band center and it is found that scaling at the band center
is described by a significantly different set of generalized dimensions than of those
at nearby energies. This discontinuity at the band center is shown to be a property
of the two central states of the spectrum. It is known from previous studies that
ABD model has a critical point at the band center, and this dissertation studies the
corresponding divergence of the correlation length, which is shown to be accurately

described by a power—law, and the critical exponent is calculated for the first time.

94

 

This calculation is done using a standard transfer-matrix method (TMM) as well as
a new finite-size scaling analysis (FSS) related to the two central states of the spec-
trum, and both methods give the same exponent within error bars. Another model in
the same class is the random flux model (RFM), describing an electron moving on a
square lattice in the presence of random magnetic flux through each plaquette. RFM
has attracted a lot of interest and previous studies have proposed several mutually
incompatible solutions to the problem of its scaling properties. The model is stud-
ied using TMM, and the high-precision of calculations allowed us to conclude that
scaling of the correlation length near the band center has similar behavior to that of
the ABD model. Therefore the FSS related to the two central states of the Spectrum
is applied to obtain the critical exponent. The differences between ABD and RFM
are also discussed, the main one being a broken symmetry at the critical point in the
former model.

Another problem studied in this dissertation is the temperature T dependence of
the magnetic properties of a small classical magnetic cluster of Gd13 at finite T. It
was known from previous T = 0 studies that the cluster has a hexagonal close-packed
structure and its magnetic properties could be described by a classical Heisenberg
model with competing ferromagnetic (among the nearest) and antiferromagnetic in-
teractions (among the next-nearest neighbors). This model has been used here to
probe finite T properties of the cluster. Magnetization is calculated from the parallel
tempering Monte-Carlo numerical simulations, and it exhibits an anomalous increase
with T in a certain range of T and the ratio of fero— to antifero—couplings. Our

simulation results are in good agreement with experimental findings.

95

 

BIBLIOGRAPHY

96

BIBLIOGRAPHY

[1] W. Kohn, Rev. Mod. Phys. 71, 1253 (1999).

[2] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum
Field Theory in Statistical Physics (Dover Publications, Inc., 1963).

[3] J. Cardy, Scaling and Renormalization in Statistical Physics (Cambridge Uni-
versity Press, 1996).

[4] R. B. Laughlin, Rev. Mod. Phys. 71, 863 (1998).

[5] P. W. Anderson, The Theory of Superconductivity in the High- Tc Cuprate Su-
perconductors (Princeton Series in Physics, 1997).

[6] P. W. Anderson, Phys. Rev. 109, 1492 (1958).
[7] D. J. Thouless, Phys. Lett. C 13, 93 (1974).
[8] P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. 57, 287 (1985).

[9] B. L. Altshuler, P. A. Lee, R. A. Webb, Mesoscopic Phenomena in Solids, (North-
Holland, 1991).

[10] B. Kramer and A. MacKinnon, Rep. Prog. Phys. 56, 1469 (1993).
[11] B. Huckestein, Rev. Mod. Phys. 67 , 357 (1995).

[12] B. L. Altshuler and B. D. Simons, Universalities: from Anderson Localization to
Quantum Chaos (Proceedings of the Les Houches School, 1996).

[13] C. W. J. Beenakker, Rev. Mod. Phys. 69, 731 (1997).

[14] S. L. Sondhi, S. M. Girvin, J. P. Carini, and D. Shahar, Rev. Mod. Phys. 69,
315 (1997).

97

[15] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan,
Phys. Rev. B 42, 673 (1979).

[16] D. C. Licciardello, D. J. Thouless, Phys. Rev. Lett. 35, 1475 (1975).
[17] A. MacKinnon and B. Kramer, Phys. Rev. Lett. 47, 1546 (1981).
[18] A. MacKinnon and B. Kramer, Z. Phys. B 53, 1 (1983).

[19] B. L. Altshuler, V. E. Kravstov, and I. V. Lerner, Sov. Phys. JETP 64, 1352
(1986).

[20] B. L. Altshuler, V. E. Kravstov, and I. V. Lerner, Phys. Lett. A 134, 488 (1989).

[21] P. W. Anderson, D. J. Thouless, E. Abrahams, and D. S. Fisher, Phys. Rev. B
22, 3519 (1980).

[22] B. Shapiro, Phys. Rev. B 34, 4394 (1986).

[23] B. Shapiro, Phys. Rev. Lett. 65, 1510 (1990).

[24] K. Slevin, P. MarkoS, T. Ohtsuki, Phys. Rev. Lett. 86, 3594 (2001).

[25] L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon Press Ltd, 1958).

[26] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wrong,
Rev. Mod. Phys. 53, 385 (1981).

[27] F. Haake, Quantum Signatures of Chaos (Springer-Verlag, 1991).
[28] V. Zelevinsky, Annu. Rev. Nucl. Part. Sci. 46, 237 (1996).

[29] B. I. Shklovskii, B. Shapiro, B. R. Sears, P. Lambrianides, and H. B. Shore,
Phys. Rev. B 47, 11487 (1993).

[30] L. P. Gorkov and G. M. Eliashberg, Sov. Phys. JETP 21 940 (1965).
[31] K. B. Efetov, Adv. in Phys. 32, 53 (1983).

[32] J. J. M. Verbaarschot and I. Zahed, Phys. Rev. Lett. 70, 3852 (1993).
[33] G. Theodeorou and M. Cohen, Phys. Rev. B 13, 4597 (1976).

[34] L. Fleishman and D. C. Licciardello, J. Phys. C 10, L125, (1977).

[35] C. M. Soukoulis and E. N. Economou, Phys. Rev. B 24, 5698 (1981).

98

[36] F. Wegner, Z. Phys. B 44, 9 (1981).

[37] R. Gade and F. Wegner, Nucl. Phys. B360, 213 (1991).

[38] R. Gade, Nucl. Phys. B 398, 499 (1993).

[39] M. Inui, S. A. Trugman, and E. Abrahams, Phys. Rev. B 49, 3190 (1994).
[40] F. Wegner, Phys. Rev. B 19, 783 (1979).

[41] Y. Hatsugai, X.—G. Wen, M. Kohmoto, Phys. Rev. B 56, 1061 (1997).

[42] Y. Morita and Y. Hatsugai, Phys. Rev. Lett. 79, 3728 (1997).

[43] Y. Morita and Y. Hatsugai, Phys. Rev. B 58, 6680 (1998).

[44] C. M. Soukoulis, I. Webman, G. S. Crest, and E. N. Economou, Phys. Rev. B
26, 1838 (1982).

[45] A. Eilmes, R. A. Rdmer, and M. Schreiber, Eur. Phys. J. B 1 29, (1998).
[46] J. Miler and J. Wang, Phys. Rev. Lett. 76, 1461 (1996).

[47] J. L. Pichard and G. Sarma, J. Phys. C 14, L127 (1981).

[48] J. L. Pichard and G. Sarma, J. Phys. C 14, L617 (1981).

[49] P. W. Brouwer, C. Mudry, and A. Furusaki, Nucl. Phys. B 565, 653 (2000).

[50] V. Z. Cerovski, S. D. Mahanti, T. A. Kaplan, and A. Taraphder, Phys. Rev. B
59, 13977 (1999).

[51] S. Yunoki, J. Hu, A. L. Malvezi, A. Moreo, N. Furukawa, and E. Dagotto,
Phys. Rev. Lett. 80, 845 (1998).

[52] S. Yunoki and A. Moreo, Phys. Rev. B 58, 6403 (1998).

[53] Y. Morimoto, A. Asamitsu, H. Kuwahara, and Y. Tokura, Nature (London) 380,
141 (1996).

[54] T. G. Perring, G. Aeppli, and Y. Tokura, Phys. Rev. Lett. 78, 3197 (1997).
[55] Q. Li, J. Zang, A. R. Bishop, and C. M. Soukoulis, Phys. Rev. B 56, 4541 (1997).
[56] E. Miiller-Hartmann and E. Dagotto, Phys. Rev. B 54, R6819 (1996).

[57] R. J. Bell, P. Dean, Discuss. Faraday Soc. 50, 55 (1970).

99

 

[58] H. Aoki, J. Phys. C 16, L205 (1983).

[59] H. Aoki, Phys. Rev. B 33, 7310 (1986).

[60] D. Braun, G. Montambaux, and M. Pascaud, Phys. Rev. Lett. 81, 1062 (1998).
[61] C. M. Soukoulis and E. N. Economou, Phys. Rev. Lett. 52, 565 (1984).

[62] M. Schreiber, J. Phys. C 18, 2493 (1985).

[63] M. Schreiber, Phys. Rev. B 31, 6146 (1985).

[64] F. Wegner, Z. Phys. B 36, 209 (1980).

[65] C. Castellani and L. Peliti, J. Phys. A 19, L429 (1986).

[66] T. Hasley, M. Jensen, L. Kadanoff, I. Procaccia, and B. Shraiman, Phys. Rev. A
33, 1141 (1986).

[67] M. H. Jensen, L. P. Kadanoff, and I. Procaccia, Phys. Rev. A 36, 1409 (1987).
[68] G. Paladin and A. Vulpiani, Phys. Rep. 156, 147 (1987).

[69] M. Janssen, Int. J. of Mod. Phys. B 8, 943 (1994).

[70] M. Schreiber and H. Grussbach, Phys. Rev. Lett. 67 , 607 (1991).

[71] V. I. Fal’ko, K. B. Efetov, Eur. Phys. Lett. 32, 627 (1995).

[72] V. I. Fal’ko, K. B. Efetov, Phys. Rev. B 52, 17413 (1995).

[73] S. N. Evangelou, Physica A 167, 199 (1990).

[74] J. T. Chayes and L. Chayes, D. S. Fisher, and T. Spencer, Phys. Rev. Lett. 24,
2999 (1986).

[75] L. B. Ioffe and A. I. Larkin, Phys. Rev. B 39, 8988 (1988).

[76] N. Nagaosa and P. A. Lee, Phys. Rev. Lett. 64, 2450 (1990).

[77] P. A. Lee and N. Nagaosa, Phys. Rev. B 45, 966 (1992).

[78] B. I. Halperin, P. A. Lee, and N. Read, Phys. Rev. B 47, 7312 (1993).
[79] V. Kalmeyer and S. C. Zhang, Phys. Rev. B 46, 9889 (1992).

[80] J. K. Jain, Phys. Rev. Lett. 63, 199 (1989).

100

[81] T. Sugiyama and N. Nagaosa, Phys. Rev. Lett. 70, 1980 (1993).
[82] Y. Avishai, Y. Hatsugai, and M. Kohmoto, Phys. Rev. 47 , 9561 (1993).
[83] V. Kahneyer, D. Wei, D. Arovas, and S. Zhang, Phys. Rev. 48, 11095 (1993).

[84] D. Z. Liu, X. C. Xie, S. Das Sarma, and S. C. Zhang, Phys. Rev. B 52, 5858
(1995).

[85] D. N. Sheng and Z. Y. Weng, Phys. Rev. Lett. 75, 2388 (1995).

[86] K. Yakubo and Y. Goto, Phys. Rev. B 54, 13432 (1996).

[87] K. Yang and R. N. Bhatt, Phys. Rev. B 55, R1922 (1997).

[88] J. A. Vergés, Phys. Rev. B 57 , 870 (1998).

[89] X. C. Xie, X. R. Wang, and D. Z. Liu, Phys. Rev. Lett. 80, 3563 (1998).
[90] A. Furusaki, Phys. Rev. Lett. 82, 604 (1998).

[91] Q. Niu, D. J. Thouless, and Y.-S. Wu, Phys. Rev. B 31, 3372 (1985).

[92] D. P. Arovas, R. N. Bhatt, F. D. M. Haldane, P. B. Littlewood, and R. Rammal,
Phys. Rev. Lett. 60, 619 (1988).

[93] V. Z. Cerovski, Phys. Rev. B 62, 12775 (2000).

[94] S. Jim, T. H. Tiefel, M. McCormack, R. A. Fastnacht, R. Ramesh, and L. H. Chen,
Science 264, 413 (1994).

[95] G. H. Jonker and J. H. van Santen, Physica (Amsterdam) 16, 337 (1950).
[96] J. H. van Santen and G. H. Jonker, Physica (Amsterdam) 16, 599 (1950).
[97] P. W. Anderson and H. Hasegawa, Phys. Rev. 100, 675 (1955).

[98] P. G. de Gennes, Phys. Rev. 118, 141 (1960).

[99] K. Kubo and N. Ohata, J. Phys. Soc. Jpn. 33, 21 (1972).

[100] N. Furukawa, J. Phys. Soc. Jpn. 65, 1174 (1996).

[101] S. Yunoki, J. Hu, A. L. Malvezi, A. Moreo, N. Furukawa, and E. Dagotto,
Phys. Rev. Lett. 80, 845 (1998).

101

[102] J. Zang, H. Rbder, A. R. Bishop, and S. A. ngman, J. Phys: Condens. Matter
19, L157 (1997).

[103] T. A. Kaplan and S. D. Mahanti, J. Phys: Condens. Matter 19, L291 (1997).
[104] S. Yunoki and A. Moreo, Phys. Rev. B 58, 6403 (1998).
[105] P. Horsch, J. Jaklic, and F. Mack, Phys. Rev. B 59, 6217 (1999).

[106] Small Particles and Inorganic Clusters, ed. H. H. Anderson, (Springer Verlag,
New York, 1997).

[107] J. Shi, S. Gider, K. Babcock, and D. D. Awschalom, Science 271, 937 (1996).
[108] W. A. deHeer, P. Milani, and A. Chatelain, Phys. Rev. Lett. 65, 488 (1990).

[109] J. P. Bucher, D. C. Douglass, and L. A. Bloomfield, Phys. Rev. Lett. 66, 3052
(1991).

[110] S. N. Khanna and S. Linderoth, Phys. Rev. Lett. 67, 742 (1991).

[111] S. E. Apsel, J. W. Emmert, J. Deng, and L. A. Bloomfield, Phys. Rev. Lett.
76, 1441 (1996).

[112] B. V. Reddy, S. N. Khanna, and B. I. Dunlap, Phys. Rev. Lett. 70, 3323 (1993).

[113] A. J. Cox, J. G. Louderback, and L. A. Bloomfield, Phys. Rev. Lett. 71, 923
(1993).

[114] M. R. Pederson, F. Reuse, and S. N. Khanna, (Submitted for publication).

[115] L. M. Falicov, D. T. Pierce, S. D. Bader, R. Gronski, K. B. Hathaway, H.
Hopster, D. N. Lambeth, S. S. Parkin, G. Prinz, M. Salamon, I. K. Schuller, and
R. H. Victora, J. Mater. Res. 5, 1299 (1990).

[116] R. J. Van Zee, S. Li, and W. Weltner, Jr., J. Chem. Phys. Rev. 100, 4010
(1994).

[117] D. C. Douglass, J. P. Bucher, and L. A. Bloomfield, Phys. Rev. Lett. 68, 1774
(1992).

[118] J. P. Bucher, Int. J. of Mod. Phys. 7, 1079 (1993).

[119] D. P. Pappas, A. P. Popov, A. N. Anisimov, B. V. Reddy, and S. N. Khanna,
Phys. Rev. Lett. 76, 4332 (1996).

102

[120] D. Mattis, The Theory of Magnetism (Harper and Row, New York, 1965).
[121] N. V. Shevchenko (Private Communication).

[122] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical
Recipes in C : The Art of Scientific Computing (Cambridge University Press,
1998).

[123] E. Marinari, and G. Parisi, Report No. hep-lat/9205018, (1992).
[124] E. Marinari, Report No. cond-mat/9612010, (1996).

[125] M. Hehn, K. Ounadjela, J. P. Bucher, F. Rousseaux, D. Decanini, and B. Chap-
pert, Science 272, 1782 (1996).

[126] D. Gerion, A. Hirt, and A. Chatelain, Phys. Rev. Lett. 83, 532 (1999).

[127] V. Z. Cerovski, S. D. Mahanti, and S. N. Khanna, Eur. Phys. J. D 10, 119
(2000).

[128] V. Z. Cerovski, (to appear in Phys. Rev. B Rapid Comm.)

103

lllllilltll