Includes bibliographical references (leaves 224-234).

Additional Physical Form:

Also available online.

General Note:

Typescript.

General Note:

Vita.

Statement of Responsibility:

by Brian Gregory Wilson.

Record Information

Source Institution:

University of Florida

Holding Location:

University of Florida

Rights Management:

Copyright Brian Gregory Wilson. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.

A DISSERATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA 1987

TABLE OF CONTENTS

PAGE

ABSTRACT................................................. ....... iv

CHAPTER

I INTRODUCTION.. ............................. ................. 1

II BACKGROUND THEORY................................... 7

The System......................................... 8
Basics of Linear Response........................... 9
The Dynamic Density Response ....................... 11
Relation to Density Correlation Functions............ 16
Further Relations .................................. 17
Distribution Functions.............................. 19

III QUANTAL FUNCTIONAL EXPANSIONS ........................... 27

The QOZ Relation................................... 27
Classical Generating Functionals.................... 32
Quantal Generating Functionals...................... 38
Density Functional Approach......................... 45

IV THE BASIC QHNC EQUATION.................................. 51

The QHNC Effective Potential ....................... 54
Numerical Calculations.. ........................... 57

V LOCAL FIELD CORRECTIONS.... ............................. 69

Mean Field Theories ................................ 71
The STLS Method.................................... 73
The Relaxation Function............................. 78
Scalar Products of Operators ....................... 81
Microscopic Theory ................................. 86

VI NUMERICALLY EVALUATION S(Q). ............................ 93

The Lindhard Function............................... 97
Convergence Acceleration............................ 103

VII THERMODYNAMIC PROPERTIES................................. 111
Thermodynamics with Slater Type
Effective Potentials ............................. 116
Interaction Energy Results.......................... 118
Kinetic Energy Results.............................. 120

VIII RESULTS OF THE EXTENDED QHNC EQUATIONS .................. 131

Method of Solution................................. 133
Quantal Hartree Results............................. 141
The Zwanzig Equation ............................... 144
Blending........................................... 147

IX DYNAMICS OF THE LFCF.................................... 189

Calculating the Memory Function..................... 197
Results............................................ 200

X CONCLUSIONS............................................. 219

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

THE ELECTRON LIQUID AT ANY DEGENERACY

Brian Gregory Wilson

May 1987

Chairman: Charles F. Hooper, Jr. Major Department: Physics

The dielectric formulation of the many-body problem is applied to the study of the static correlations in electron liquids at non-zero temperatures. The intermediate coupling effects arising from the exchange and Coulomb correlations are taken into account through a local-field correction factor obtained from quantal extensions of classical fluid integral equations. This results in a unified theory which can describe the thermodynamic and dielectric properties of the

-iv-

electron liquid covering the limits of the zero temperature jellium to the classical one-component plasma.

Detailed numerical calculations are provided for the

self-consistent set of equations describing the static structure factor. The variation of interaction energy per particle versus temperature at constant density is shown to exhibit a peaked structure near absolute zero. This feature is not observed using effective potential approaches in classical fluid equations. A comparison of zero temperature results with those obtained from the Greens function Monte-Carlo method is good; however, agreement worsens when dynamical effects of the local-field correction factor are included using the Mori continued fraction method.

-v-

CHAPTER I

INTRODUCTION

The one component plasma (OCP) is a system consisting of a single species of charged point particles together with a uniform oppositely charged rigid background to ensure charge neutrality. It was first introduced by Sommerfeld as an approximation to metals when band structure effects are negligible. By coupling the standard adiabatic (Born-Oppenheimer) approximation with an assumption of a weak electronion pseudopotential interaction, it can be shownl-3 that the thermodynamics of a simple metal is the sum of that due to a one component electron plasma and that due to a system of classical ions interacting via a state-dependent pair potential. Viewed as a model fluid, the OCP exhibits the characteristics which distinguish real coulomb systems, such as plasmas and electrolyte solutions, from neutral fluids. These include the phenomena of screening and plasma oscillations, which arise from the long-range nature of the coulomb interaction. Thus the OCP is a standard approximation for both astrophysical and laboratory fusion plasmas,4 where the electrons are treated as a polarizable background for the gas of ions. Further applications lie within the local density approximation of the density functional theory of electronic structure, where the excess chemical potential of a one-component electron plasma plays the role of an

-1-

-2

exchange-correlation potential in the one-electron Schroedinger equation.

It is thus not surprising that there has been a large amount of work expended in obtaining accurate theoretical descriptions of both the zero temperature fermion OCP (degenerate electron liquid) and the high temperature classical OCP (electron gas).5-9 Active research is ongoing in three main areas. First, in the course of study of linear dispersion of plasmons in the simple organic polymers polyacetylene and polydiacetylene, it has proved important to contrast the properties of the exchange hole, and its screening, in quasi-one-dimensional solids with those in isotropic three-dimensional jellium (classical OCP).10 Secondly, one has been able to determine experimentally, using inelastic X-ray scattering, the dynamical structure factor S(q, w) for the electrons in a few metals.11-13 It is of particular interest that the experimental S(q, w) show a tendency of two-peak structure for q within the particle-hole continuum, a feature believed to be a property of the uniform electron liquid not presently accounted for accurately.14-19 Lastly renewed attempts are being made to quantify the intermediate degeneracy regime of the fermion plasma.20-22

There are, in fact, many physical systems whose electrons exist in partially degenerate states and are inadequately described by either zero temperature or classical formalisms. For example, the finite temperature thermodynamic and dielectric properties of the electron liquid are needed for a proper treatment of the equation of state of high temperature liquid metals forund in shock wave experiments23 and experiments on liquid metals expanded to near their critical

-3

points,24-25 of high density inertially confined fusion plasmas, and of the finite temperature exchange-correlation potentials used in density functional calculations of atomic properties at high pressure and temperature.26 Semidegenerate electron liquids are also found in semiconductors where the low densities involved lead to correspondingly low Fermi temperatures. The finite temperature equation of state of the electron liquid is also needed for a quantitive explanation of the miscibility gap in solutions of alkali metals in their alkali-halide melts.27

The problem of the OCP at intermediate degeneracy has only
recently received attention because of complications not present in the extremes of a classical or a ground state quantum description. For example at zero temperature diagrammatic,28 self-consistent distribution function methods29-31 and pseudo integral equation methods,32-37 together with recent Quantum Monte Carlo38-40 calculations, have combined to present an accurate description of the paramagnetic electron liquid.41-42 But at finite temperature we no longer deal with a unique ground state (invalidating the Monte Carlo approach) nor may we utilize the ground state energy variational principle (invalidating Fermi Hypernetted Chain approaches).* Diagrammatic approaches may be generalized to finite temperatures43

Fermi Hypernetted Chain calculations of quantum fluids9 are
actually variational calculations of the ground state using an adjustable Slater-Jastrow type wavefunction and should not be
confused with the classical integral equation method per se.Rather
the name stems from a diagramatic expansion of the ground state
distribution function which is regrouped into a structure of the
same form as the classical HNC equation.

-4

but are perturbation expansions reliable only in the regime of weak coupling.44
On the other hand the classical limit may be handled by a variety of techniques, including the Monte Carlo method,45-50 the molecular dynamics method,51 and various approximate methods involving integral equations, most notable amongst these the Hypernetted-Chain equation52-58 and modifications thereof.59-62 Here quantal effects introduce two complications. First a quantal calculation of the pair distribution function is hampered fundamentally by the fact that the canonical partition function can no longer be factored into a solely configuration-integral part and a part involving solely momentum coordinates. Secondly the thermodynamic description will also depend on the momentum distribution (or off diagonal single body density matrix)63-64 and this information is not trivially contained in the pair distribution function.

Various attempts have been made to overcome these complications within the framework of the successful methods available for classical fluids. Wigner-Kirkwood expansions65-66 may be used for those cases where quantum corrections to the classical partition function are perturbative, but their convergence as a series in temperature or Plancks constant is slow (if existent) and is restricted to short ranged potentials;67 it is, in fact, an ill defined expansion for the OCP.68 Published work often centers around diffraction effects and assumes Boltzmann statistics.

Previous calculations of the OCP at intermediate degeneracy and coupling have been of several types. The most obvious course has been

-5

to approximate the true quantum mechanical Slater sum which appears in the canonical partition function with the corresponding classical expression involving a temperature and density dependent effective potential. Pokrant69 and others70-71 have developed a finite temperature variational principle to approximate the Slater sum as a product of an ideal fermion Slater sum and that of pairwise additive potentials. The energy is calculated by approximating three-body distribution functions in terms of pair distribution functions, which are obtained from the HNC equation using the effective potentials. Zero temperature correlation energies from this method are, however, in substantial disagreement with the parameterized fittings of the quantum Monte-Carlo results.41-42

Several other groups have studied the thermodynamic21-22'72-75 and dielectric76,77 properties of the OCP by using finite temperature perturbation theory within the random phase approximation. In particular Perot and Dharma-Wardana21 and Kanhere et al.22 have presented closed form parameterizations of the exchange-correlation energy and chemical potential within the random phase approximation (RPA) for the paramagnetic and spin polarized cases. However the RPA is a weak coupling theory that results in considerable error at metallic densities.

Dandrea et al.20 attempt to overcome the limitations of the RPA by including a static local field correction factor (LFCF) within the framework of the dielectric formulation. They use an approximate form for the LFCF which has two adjustable parameters. The first is fixed by assuming a value for the pair distribution function at the radial

-6

origin. The second is obtained from an appropriate interpolation between known Monte-Carlo results at zero temperature41-42 and in the classical regime47-50 for the bulk modulus. Tanaka et al.78 approximate the LFCF at all temperatures by the static zero temperature LFCF of Singwi-Tosi-Land-Sjolander (STLS) self-consistency scheme.29-31

This thesis concerns itself with a unified theory capable of

describing the dynamic density response function of the OCP throughout a broad range of densities and temperatures: from a fully degenerate quantum plasma at absolute zero temperature to a nondegenerate classical gas. It differs from the approach of Dandrea and Ashcroft or Tanaka et al. in two important respects. First, that the theory be particularly applicable to electron liquids at metallic densities, which are not weakly coupled, perturbative local field corrections to the random phase approximation are replaced altogether by a selfconsistent field determined by a quantal extension of classical integral equations. Formally we can define a LFCF that, unlike Dandrea et al., requires no known results, and unlike Tanaka et al., is temperature dependent; indeed it may be viewed as a generalization of the STLS LFCF of zero temperature and reduces to the proper classical result. Secondly, dynamical dependencies of the LFCF can be rigorously included. This has been shown to be a mechanism for modeling the double peaked shape of the dynamic structure factor observed in metals.11-19 It is also an essential theoretical featurestatic LFCFs cannot simultaneously satisfy the compressibility sum rule and the third frequency moment sum rule for the density response.79

CHAPTER II

BACKGROUND THEORY

This thesis will concern itself solely with the OCP as a Fermi liquid. The most basic characteristic of a liquid is that it possesses short range order as opposed to the long range periodicity of a crystalline solid. Since the structure of a crystal is determined experimentally by observing the Bragg reflection of X-rays, it is natural to seek a quantitative description of the liquid structure via the intensity of X-ray, thermal neutron, or light scattering.3,80-81 Of particular interest then are the fluctuations of space and time dependent densities which describe, at long wavelength, the cooperative motion of many particle systems. For this purpose define the particle density correlation function as

S (r ', t) : <6n(, t) Sn(6, o)>
nn eq

where 6n(r,t) is the excess particle density operator

Sn(r, t) = n(r, t) <6n(r, t)>eq

and the equilibrium average is indicated by the brackets <...> eq Specifically, even though eq = 0, there will be spontaneous, usually small, fluctuations on a local scale. The density correlation function is a measure of these fluctuations.

-7-

-8

In this chapter we will begin our investigation into the

structure of quantum fluids by reviewing the general properties of the time correlation functions related to the particle density correlation function. These include 1) various symmetries (for example time reversal), 2) positivity properties which are related to the dynamical stability of the system, 3) fluctuation dissipation theorems which connect spontaneous fluctuations and energy dissipation in a thermally equilibriated system, 4) Kramers-Kronig relations which express causality, like the well known one between the index of refraction and the absorption coefficient, and 5) sum rules which provide restrictions that any approximate microscopic theory must fulfill.

The System

Because we will be dealing with a many particle system from a quantum statistical point of view, we will employ the second quantization representation of quantum mechanics even though we will mostly be working in a canonical ensemble where the number of particles is fixed.

For our equilibrium "system" we shall define in the Schroedinger (time dependent state) representation the unperturbed Hamiltonian

H = T + V

in terms of field creation and annihilation operators.43 The kinetic energy operator is

-9

h2 82
T = J dI ~I+(I) {- 2nn -2 (I) 8r1

while the interaction energy operator is given by

2
S f dl dlI k'(I) (I) { e '(I) t(I) Irl r21

Roman numerals were used as shorthand notation for labels; that is "I" stands for the position coordinate as well as the discreet spin label of particle 1.

If there is an external (possibly time dependent) field present, there is in addition a perturbation to the hamiltonian of the form

SH(t) = f drn(r) 6V(r, t)

where we introduced the electron density operator

n(X) = i '+(I) T(I)
spins

Basics of Linear Response82-85

In the external field the hamiltonian is explicitly time dependent and given by

H(t) = H + 6H(t)

in the Schroedinger representation where the operator n(r) is time independent. The time dependance is carried by the density matrix, or ensemble operator p(t) which describes the state of the system such that the average of n(r) is

= tr p(t) n(r) tr p(t) = 1

-10

What follows is entirely parallel to traditional derivations of time dependent perturbation theory in elementary quantum mechanics texts. We have to solve the time evolution equation for the density matrix

i at p(t) = [H(t), p(t)] = [H, p(t)] + [6H(t), p(t)]

(note the commutator is reversed from the Heisenberg equation of motion for operators) subject to the initial condition

p(t = m) = Peq [H, p eq] = 0

The initial condition expresses the fact that the system is stationary before the external field is turned on--we also require that the external field decay sufficiently rapidly as time reaches infinity.

For manipulations it does not matter what p eq is but since the system starts from thermal equilibrium we take

Peq = P-BH/tr p-BH

If we assume a time dependance to the density matrix of the form

p(t) = Peq + 6p(t)

then the solution to the time evolution equation to first order in 6p and linear in the external field is

Sp(t) = dr e-iH(t r)/f [6H(), eq e+iH(t )/

(This is easily verified by differentiation remembering that time appears both in the limits and in the integrand.) From this equation the induced change in the average density

-11

6 tr p(t) n(r) tr Peq n(r) is given by

t
6 = f dt' f dr' < [n(? t), n(F' t')]>eq 6V(6' t')

t
zi f dt' J d'' x"(rt; r' t') 6V(6' t') (2.A)
-w

where [A, B] = AB BA is the commutator bracket. In the above equation are the Heisenberg operator for the unperturbed system

n(r' t) = et/ n(r) e-iHt/

and we have introduced the dynamic density response function (please note the double prime superscript) defined by

x"(rt; 't') = < En(t), n('t')]>eq

This equation is the fundamental result of linear response theory. It shows that the density response to a small external potential is the averaged commutator rather than the correlation function S(r,t) as one might expect. However we shall see that the two are intimately related.

The Dynamic Density Response Function

In equilibrium the system is space translationally invariant. We adopt the convention that forward spatial Fourier transforms have a minus one signature and write

-12

x"(E) = d(r -r) e x"(r r')

= < 1 [n (t), nj(t')]> k k

Here we have used the definition of the spatial transform of the density operator

-ik-r
n (t) = S dr e n(rt)

which can be expressed in terms of plane wave creation and annihilation operators as

n (o) = f dP a+ 'a+
K a (21) P+K,a from which it follows that

n n -K K

Time translation invariance of the dynamic density response

function is easily established from the definition of a Heisenberg operator and the cyclic nature of the trace. We adopt the notation that forward temporal Fourier transforms have a plus one signature:

x"(Kw) = f d(t t') ei(t t') x"(K, t t')

It is easy to show that this space and time Fourier transform is real (since the function is a commutator of hermitian operators), an odd function of frequency (since the equilibrium state is invariant under time reversal and parity), and depends only on the magnitude of k (spatially isotropic). It can also be shown that

ex"(Kw) > 0

-13

in a stable system.85 This is just a statement that a dissipative many body system takes more energy out of the external field than it gives back.

In addition to the density-response function it is convenient to introduce certain interrelated functions. The complex response function (please note the tilda) is defined in transform space as

x(Kz) = f S x"(kw)
n w-Z

This is an analytic function of the complex variable z as long as imaginary z does not equal zero (on the real axis it has a branch cut). For z above the real axis it reduces to the temporal Laplace transform of the density response:

x(Kz) = zi S dt eizt x"(Kt)
0

while below the real axis

o izt
= (- zi) f dt e x"(Kt) This follows from the identity

iwt izt
co L e e t>o nn
f 21i w z 0 t < o

o t>o
t >o Inn z < o
-izt
e t
which can be proven using the Cauchy integral formula in the complex omega plane. (For t>O we can close the contour on top-exp(iwt) being bounded-and the contour is in the positive sense. For t

-14

negative sense. We then just sum the residues, the pole being included depending on the sign of imaginary z.)

The physical response may be defined in transform space as the limit of the complex response as we approach the real frequency axis from above

x(Kw) = lim x(K, z = w + ic) = lim f do x"(Kw')
C0 CO _. ir w' -(w + ic) Using the identity (P denoting Cauchys principal value)

= P + ir 6(x)
X+ic X

the physical response can be expressed in terms of its real and imaginary parts as

x(kw) = x'(Kw) + ix"(Kw)

where the imaginary part is just the previously defined density response and the real part (note the single prime) is

x' (Ko)) = P f d x"(Kw)
-wnw'an even function of frequency, is a manifestation of the KramersKronig relations.86

The utility of the physical response function is that it links

the induced change in the time-space Fourier component of the density with the time-space Fourier component of an arbitrary external potential

6 = x(Kw) 6(Kw)

-15

One may also note that this same relation may be easily obtained via eq. 2.A using the definition of the physical response in real (untransformed) space, namely

x(rt) = zi e(t) x"(r,t)

The response to a static external potential may be derived from linear response theory by assuming that the potential is adiabatically switched on starting from the remote past

-E(t)
SV(Ft) = V() e(t) t < o o t > o At t = 0 the external potential is at full strength and we find

6 = zi J d fJ dI' x"(r r', -) e-EZ 6n(')
0

or

6 = x(K) 6V(K)

where the static response is given by zero frequency limit of the physical response

x(K) E x(K, = o) = lim x(K, z = ic) C-o
A note on this last equality. Formally we have lim x(K, z = ic) = lim (Pf x"(Kw) + ix"(K, w= o) s0o E-o -C W -but because x" is a real odd function of frequency the imaginary term vanishes and we can drop the "principal part" condition.

-16

Relation to Density Correlation Functions

At this point we should explicitly confirm the statement that the density response is intimately related to the density correlation function, the latter being the experimentally pertinent quantity. Now because the Fourier transforms of correlation spectra are related by the equation83

S eq e -it dt ePO fc e-it dt it follows directly that

< e-i"t 1 e-Of -wt
1 < [A, B(t)]> e dt (1 -e ) eq e dt and so in particular (using the fact that x"(Kr) is odd in frequency) x"(Kw) ~ (1 e-f ) S(Kw)

where S(k,w) is the time-space Fourier transform of the time dependent density correlation function. We note that S(r,t) is Fourier transformable because when r and/or t are very large, then n(r,t) and n(O,O) are statistically independent

eq eq eq

and so S(r,t) as defined as a correlation of excess density vanishes.

The above relation between the density response and density correlation functions is an explicit example of what is called the fluctuation dissipation theorem84 which relates two physically distinct quantities of fundamental experimental significance: 1)

-17

Spontaneous fluctuations, which arise even in the absence of external forces from the internal motion of the constituent particles. Described by S(k,w) these fluctuations give rise to scattering of neutrons or light. 2) Dissipative behavior--all or part of the work done by external stirring forces--is irreversibly disseminated into the infinitely many degrees of freedom of thermal systems.

The fluctuation-dissipation theorem shows that S(k,w) is not

quite symmetric in w, being a little stronger at positive frequencies than at negative. Indeed, since wX"(Kw) is always even in w, it is generally true that

S(K, -w) = e-fK S(K, w)

This result makes sense from a neutron scattering point of view. Positive frequency means the neutron has lost energy to the system (by creating an excitation of energy hw) while negative frequency describes a process in which the neutron has picked up energy from the system (by destroying an excitation). Of course, to destroy an excitation you must first have one, and their relative abundance is given by exp[-bhw]. This dissymmetry of the scattering intensity, proportional to S(k,w) is only pronounced at low temperatures, kt
-- it is absent classically.

Further Relations79

Besides its relation to the time dependent density correlation function the density response satisfies certain sum rules, usually presented in the form of moments of the density response:

-18

<<)>> E f J { x"(Kc)} dw

The j = 0 or zeroth moment follows from analyticity of x"(k,w).

D x"(Kw) = x(K, o) = x(K)

All odd order moments vanish because x" is odd in frequency. The remaining higher order moments can be derived from the expression

ia8 1
(t) x"(rt; 't') = < ) n(, t) n('t')]> taken at equal times t = t'; this means that

I dw J x"(Kw) I d( e') ei(r r')

<(i)J+1 [[...En(rt), H] ..., HI, n(r't)]> The right hand side contains a sequence of equal time commutators which can in principle, and sometimes in fact, be exactly calculated. For example the second order moment yields the Thomas-Reiche-Kuhn or F-sum rule equivalent

CK2
f d wX"(Kw) = eK
Ira m

For the electron gas the fourth moment has also been calculated.87-89

The sum rules are often used to provide coefficients for an expansion of x(k,z) for large z from the definition

x(Kz) = J d xK( <> <>
(Kz)f Z z 2 4+

-19

From its derivation, which expands (1 w/z)-1 = 1 + (w/lz) + (w/z)2 + ... it is clear this expansion can only be asymptotic, i.e. valid only when Izl is large compared to all frequencies in the system, which means all frequencies for which x"(k,w) is not effectively zero.

An interesting feature of the sum rules is their existence. There is no reason why the thermodynamic average of the multiple commutator should exist to all orders, and indeed although the sixth moment has been calculated in the classical fluid limit90 it may diverge for the degenerate electron gas in certain regions.91 As the sixth moment determines the w-6 term in an inverse frequency expansion, its divergence would herald the presence of a term of the form k4w-11/2 in S(k,w) at long wavelength for the high frequency behavior of the spectrum. Great caution is thus indicated in the use of moment sum rules beyond the first three mentioned.

Distribution Functions

It is natural to describe the structure and thermodynamic

properties of liquids in equilibrium by employing static distribution functions. Most treatises on fluids introduce such quantities first and later generalize to include time dependance as a prelude to exploring non-equilibrium statistical mechanics. Here we have taken a reverse path; we have already related time dependent density correlation and density response functions and it is our intent now to relate these with the static distribution functions common to liquid theory.

-20

The most basic static distribution function is the single

particle distribution function. If we define the number density operator in N particle configuration space as

n
n(r) = X 6C- ii)
i=l

(the ri are quantum operators) then the single particle distribution is given by

p(r) = N f
p() = = N J W(r 1, r2 ". rN) dr2 .. dN

where the configurational part of the partition function is

Q =f W(r .. rN) drl .. drN

and the configuration probability W(rl, ... r n) is given by the diagonal slater sum

W(rl, ... rnlrl, ... rn)

where, in general, the off diagonal Slater sum is defined as

3N -Hop W(r rN I ...rN) = N 3N rr" N) e -So ;i r N )

In the above equation the index "i" denotes any complete set of properly symmetrized N-particle basis states and we have introduced the DeBroglie thermal wavelength

2 1/2
T m

-21

In the classical limit the probability function "W" is just the Boltzmann factor of the potential energy of the configuration; at zero temperature it consists solely of the ground state wavefunction.

The information contained in the single particle density distribution is minimal, for in the case of a translationally invariant Hamiltonian it is a constant equal to the average number density of the system. However we introduce it here because relevant thermodynamic properties we will consider require the information of the off-diagonal single particle distribution function63-64

p( I ') = N (r r2... N .. ) d ... d Q 1 2"' Ni r2 rNdN The physical significance of the off-diagonal distribution function is that its spactial Fourier transform provides us with the momentum distribution of the interacting system.

The pair distribution function is defined by
2 I
p g(? ) = N(N 1) f H d-r3 .. rN = < 6(r r) 6(r' r>
p g(r Q 3 N i
lsi

and for a translationally invariant system is of the form g(r-r'). The pair distribution function in turn is related to the static structure factor through the relation

S(K) = I + p f ei1' (g(r) 1) dr

Within a factor of density the static structure factor is simply the time independent density correlation function.

S( <(n() ) (n(') )>
S(F ') =P

-22

Note that if we had defined the static structure factor in terms of the density operator instead of the excess density we would replace g(r) 1 by g(r) in the formula relating the two. This amounts to including the experimentally unmeasurable "forward scattering."

The definition of the static structure factor is immediately generalized to form the dynamic structure factor, essentially the previously introduced time dependent density correlation function (again within a constant factor of density) by replacing the quantum operators ri with their time dependent equivalents in the Heisenberg representation. Bearing in mind that ri(t) and rj(t) are operators which in general do not commute (except at t = 0), we easily obtain several important properties of the dynamic structure factor for translationally invariant systems:

1) S(-r, t) = S*(r,t) where S* is the complex conjugate of S.
2) In the classical limit the imaginary part of S(r, t) vanishes.
3) S(Kw) iwtik
3 2n dt e f d e-ikr S(rt)
is real in both the classical and quantal cases.
4) The elastic sum rule

S(K) = J S(Kw) dw

is an immediate consequence of the definitions.

All the relationships amongst the various functions discussed are best conveyed diagrammatically (See Fig. 2.1). Note that in the upper

-23

left corner we have introduced a new function C(k,w) called the Kubo or Relaxation function, about which we will have much to say later. For now it simply expresses the difference between the (dynamic) physical response and the static response function.

The point of our review of the response, correlation, and
distribution functions, so concisely summarized in Fig. 2.1, is that if we wish to generalize the classical many-body methods of calculating g(r) into the quantal regime, there are two possible relations open to us. First there is Scl(k) as the classical limit of Sqm(k) (obviously) or, as

F 6n(r I U) 1
6 u(') x(K)

where the 10 denotes evaluation at constant density/zero external potential, and the F denotes the Fourier transform with respect to the spactial variable r-r'.

Let us contrast the two alternatives. The quantum mechanical static structure factor Sqm(k) is essentially statistical in nature; it is strongly dependent on the dynamics of the particles, depending on S(k,w) Imag x(Kw) over the full range of frequency. On the other hand x(K) depends on the single zero frequency limit of x(Kw), and is essentially mechanical, arising from the Schroedinger equation by our linear response derivation, which makes but implicit reference to statistical questions. This is clearly indicated by the F-sum rule-which holds in any stationary

-24

x(kw)

x(kw) = X(K) + iwM3C(k) y(Kw) Im x(Kw) lim W o S(kw) =2 (1 + cot h ) y(kw)

x(K) = J y(Kw) dw

x(K) S(K) = J S(Kw) dw Figure 2.1

Relations between the linear density-density response function and the structure factor.

-25

ensemble--expressing nothing but the conservation of particles and the elementary commutator of position and momentum.

This dissertation will concern itself with the integral equation methods of solving the many body problem. This route seems particularly promising for treating quantal fluids in that for the classical limit there exists a wide variety of integral equations for g(r) which have been investigated and give results in good agreement with experimental results.93'3 The purpose here is to extend this procedure into the quantal regime. A fundamental ansatz of this treatise is that one would expect the more accurate integral equation approximation to correspond with the less sensitive function. The extension of classical integral equations into the quantal regime via the density response function is presented in Chapter II.

With this fundamental ansatz, however, there is a fundamental drawback: From x(K) one obtains an integral equation not for gqm(r) (which is used in the computation of thermodynamic quantities) but rather n(r/v), that is the distribution function in the presence of an external potential which has the form of the interparticle interaction. While it is true that in the classical limit n(r/v) = nog(r), this is not true in the quantal case, since we are neglecting exchange and recoil effects. The external potential can be thought of as arising from a distinguishable particle of infinite mass and therefore fixed in position. An extreme example occurs for the case of very weak interactions between fermions. In such a case n(r/v)/n0 is approximately unity everywhere, whereas g(r) has approximately the ideal fermion value of 1/2 at the origin.

-26

Thus our task is twofold: first to generate an accurate quantal integral equation for n(r/v)/no and second to generate an approximate g(r) from knowledge of n(r/v)/n In essence our second task is to go from the bottom left of Fig. 2.1 to the bottom right along the paths drawn in this figure. We cannot go along the diagonal, as this path is one directional (we lose information when we integrate y(k,w) over w in order to obtain x(k)). We must first traverse over the upper left path of Fig. 2.1; we must obtain an accurate approximation to the Relaxation function. This is taken up in subsequent chapters.

CHAPTER III

QUANTAL FUNCTIONAL EXPANSIONS

In this chapter we will establish integral equations93 for the quantity n(rlv)/n (the density distribution in the presence of an external potential of the form of the interparticle potential) through the functional expansion method.94-95 Because in the classical limit these equations reduce to the familiar Hypernetted-Chain (HNC) or Percus-Yevick (PY) equations for the pair distribution function, we will review the classical derivation with the aim of establishing corresponding quantal generalizations.

The OOZ Relation

In classical fluids the Ornstein-Zernicke (OZ) relation plays an important role in treating integral equations for g(r). Thus in this section we derive a quantal extension of the OZ relation (QOZ).

We start from the known functional dependance of the classical canonical partition function on pairwise additive potentials. This allows us to write the functional derivative of the density distribution as96

n( n 6(r r') + n2 {g(Ir r') 1} = F {n S(K)}
6 Ur') o

-27-

-28

The I denotes evaluation when the arbitrary static external potential U(r) is turned off, that is at n(r) = no (n ; the average density). The (negative power) script F denotes the (inverse) Fourier transform on the spatial variable r r', and S(k) is the previously defined static structure factor

S(K) = 1 + no f dr e- ik {g(r) }

Now consider the functional

6 BU(r')
6 N(r")

This is a vacuous quantity unless such a functional of density can actually be constructed.97 If such a creature does not exist, we define it through the chain rule of functional calculus:

f dr" 6 BU(F) 6n(F") =
6 n(") 0 6 U(r')I

It follows that the inverse functional derivative is

S- BU(r') F-1 1 ( )
Sn(r") o o

where we have split up the expression for the derivative into two pieces (this should be considered an arbitrary process) with the first piece chosen to be the non-interacting (g(r) = 1) result and the second piece (c(r)), called the direct correlation function, representing the deficit from from the non-interacting result. The functional chain rule as expressed in terms of c(r) takes the form

-29

h(r) = c(r) + no f c(Ir r'l) h(r') dr'

This is the classical Ornstein-Zernicke (OZ) relation, where h(r) = g(r) 1 is called the pair correlation function.

Furthermore, in the classical limit where there are no exchange or recoil effects (the momentum integrations factor from the configurational piece of the partition function) one can rigorously derive the bootstrap relation96

n(f I U = v) = no g(F)

where the external potential U(r) has the form of the interparticle interaction V(r). This bootstrap relation allows us to reformulate the OZ relation in terms of the inhomogeneous density distribution as follows:

f d?' 6 BU(f) n(F' I U = v) no = c(F) (3.1)
6 n( ')

A quantal extension of the OZ relation is straightforward;98 we will follow many of the above steps.

First, from linear response theory, we saw that in the quantal case

6n(r I U) = F1 {n x(K)}
6 BU(r) 0 K

where x(k) differs by constant factors from the static response function previously defined in Chapter II. Defined here it is conveniently dimensionless; note that only the product Bx(k) remains finite as temperature goes to zero.

-30

Again the chain rule of functional differentiation requires that

6 KU(r') F {1 L 1 { 1 "
6n(" I U) = K n x(K) n- ) (3.2)

but here we have decomposed the inverse derivative into two new (and as yet unspecified) functions. These are determined as follows.

As we wish c(r) to reduce to the classical direct correlation function in the proper limit, we need only require J(k) = 1 in the classical limit. We furthermore make the ansatz that

f d' 6 U(F) n(1' I U = v) no = c(r) (3.3)
6n(r') o

also holds for the quantal case.

We have seen that this equation is strictly true in the classical case; by combining the two above equations we constrain the function c(r), and so simultaneously J(k), by the relation

FK{n(? I U = v) n } = 1 (3.4)

In hindsight we now see the motivation of requiring eq. 3.3 be fulfilled: as J(k) tends to unity in the classical limit we have merely imposed the condition that we recover the bootstrap relation in the classical limit.

The above equation also determines J(k) (and so c(r)) in the general quantum limit, for if the system described by n(r/v) and x(k) were that of non-interacting fermions, then in the limit that U(r) vanishes we see that 3(k) must equal xo(k). This is an appealing result, for it amounts to defining the quantal direct

-31

correlation function by separating out the non-interacting result, just as we did in the classical case. Henceforth we will implicitly use this result, which will be referred to as the quantal bootstrap relation (QBS).

As a consequence of this now specified decomposition of the

inverse functional derivative eq. 3.2, the functional derivative chain rule expressed in terms of the quantal direct correlation function becomes

( 1) = M{c(F)} + n f d?'(M{c(I ?'I)}) (n((V) 1)
o o where M{...} is an operator defined by

FK{Mf(F)} = xo(R) FK{f(F)}

This is the quantal generalization of the OZ relation (denoted QOZ for short). Note that this equation can be cast into the classical form by defining analogue functions:

gan(F) n( Iv) han) an

can(?) E M{c()}

The classical equivalents of these functions all have "physical" interpretations that are only loosely identified with the analogue. (For example gan should not be confused with the true or quantum mechanical pair distribution function.) To distinguish classical analogues we shall henceforth employ the superscript "an."

-32

Classical Generating Functionals

The heart of the functional expansion method of deriving integral equations lies in the ability to construct a functional W[n] of the density distribution n(rlu) in a nonuniform system which is reasonably insensitive to the function n(rlu).99 (For other points of view, see Bakshi100 and Choquard.101) If such is the case we can approximate the functional W[n] by making a functional "Taylor series" expansion, for example about the average density no, and truncate after the first (linear) term

W[n] = W[n ] + df'( 6w (n(F' I U) no) Sn(&')

The OZ relation then forms a closed set of two equations in two unknown functions. As an illustration we will consider a set of functionals and the integrals equations they generate in classical statistical mechanics.

First consider the functional

n(F I U)
W[n] = ln
no eu()

The denominator forms what is called the reference distribution. In order to make the functional variationally insensitive, it is constructed to approximate the numerator as close as possible. Here we have taken the classical non-interacting Boltzmann factor. The first functional derivative is easily found to be

-33

6w 6(?-?') 6 U(?)
6n(?') n(r I U) 6 n(?')

and from the definition of the direct correlation function:

6 BU(r) 6r- r') c
6n(F') no we find

6w i c( ')
6n(F')

The functional expansion of W[n] as the potential U(r) is changed from the intitial value 0 to its final value

In n(r I V) = J d?'(n(' I v) n ) 6(- I )
n e- )V(?) 6n(')

along with the bootstrap relation n(r/v)/no = g(r) and the OZ relation then combines to give the HNC equation

g(F) = e- RV(F) ey(7)

where y(r) = h(r) c(r) is called the nodal function.

The PY equation

g(f) = e- RV(F) [1 + y(F)]

follows from the same steps by starting with the functional

W[n] = n( I U)
no e- BU(r)

Because the HNC functional is the logarithm of the PY functional, it is the more insensitive functional of the two from variations. One might be tempted to conclude that it would therefore produce the

-34

superior approximate integral equation. However the validity of the HNC equation (or any other derived by a functional expansion method) depends on whether the cumulative effects of the higher order derivatives in the functional expansion are negligible. In actuality their effect is to push the form of the equation towards that of the PY (note how the form of the PY equation is the small y(r) expansion of the exponential y(r) term of the HNC). This is evidenced by the fact that for purely repulsive potentials the PY and HNC equations of state bracket the "exact" Monte-Carlo or molecular dynamics computer simulation results.102 What we do know from diagrammatic analysis52-58 is that the exact closure equation to the OZ relation for g(r) is of the form

g(?) = e- RV(r) + y(F) + B(F)

where the unknown "bridge function" B(r) is the sum of bridge diagrams of two point functions. The evaluation of the bridge diagrams in terms of higher functional derivatives103'61 or other attempts to account for bridge effects will be discussed in Chapter V. Here our intent is to generate zeroth order integral equations from variationally insensitive functionals, thereby minimizing the effects of the higher order terms that we neglect in the Taylor functional expansion.

It is well known that the PY equation has no solution for a one component plasma system.104 Numerically this arises from the fact that in the PY equation for g(r) the long range tail of the potential is not compensated by that of the nodal function y(r), as it is in the HNC equation.105 This is symptomatic of a physically meaningful

-35

inadequacy, namely; in real systems the particles do not feel the "bare" long ranged potential but a collectively screened or short ranged potential.

An improvement on the PY generating functional would be to

replace the bare potential appearing in the Boltzmann factor of the reference distribution by an appropriately screened one. For example

n(r I U)
W[n] =
RUH(r)
n e

where the Hartree potential is defined as

BUH(F) = BU() + no f dr' (1U(? ?') (g(&') 1)

As the Boltzmann factor now more closely describes the physical distribution n(r/v) this functional is variationally less sensitive than the usual PY

6w ( ') n(r I U) 6 UH(F)
6n(r') BUH(r) BUH(r) Sn(r')
no e no e

Using the bootstrap relation in the definition of the Hartree potential we see that

(n(S) no) J dS BU(P 5) 6(?' s) which when evaluated in the limit v = O(n(rlv) = n ) and using the definition of the direct correlation function (eq. 3.4) yields

-36

6 PUH(r) s(F r')
S (r = c( ') BU(r r')
6n(r') no
0

and so

w = c(r ') + BU(r r')
6n(F')

The linear expansion of the functional W[n] thus generates the equation

n(P I U)
n(? I U) = 1 + f d?' (n(F' I U) n ) (c(F F') + BU(F F'))
BUH(r)
n e

which can be re-expressed with the help of the bootstrap equation and the OZ relation as

IVH(r)
g(r) = e {1 + IBVH(F) + y(F) V(F)}

Note that in this equation the long range part of the potential is cancelled by the long range part of the nodal function, and all other functions appearing in this equation are screened. This equation, which we shall term the PYH equation (for Percus-YevickHartree) is equally applicable to long and short range potentials, unlike the PY equation.

If we try to improve the HNC equation by introducing the screened Hartree potential into its generating functional

W[n] = ln n(T I U)
BUH(F)
no e

-37

we find that its first variation is the same as in the above PYH derivation and (following many of the same steps) the functional expansion leads to the equation

BVH () y(r) + RVH(r) BV() e- RV(F) ey()
g(?) =e e = e

namely the original HNC equation. In other words the HNC equation has the striking property of being invariant to (Hartree) screening of the reference distribution in its generating functional.

From a related point of view it has also been shown106 that the HNC equation is the limit in a series of integral equations that can be systematically generated starting from the PY equation, among others.

Finally we should note that the HNC equation can be derived from an even more interesting point of view. We ignore higher order functional derivatives and assert that the best functionals for our closure purpose are insensitive to variations, and this property shows up in the magnitude of its derivative. This can be minimized by using an even more realistic screened potential for the reference distribution. In fact it is easy to show that if we modify the Hartree potential by replacing the convoluted bare potential by the direct correlation function

Uc() 7 BU(r) f dr' c(P ?') (g(F') 1)

then the derivative of both the HNC and PY like functionals

S= ln n(r -rU) n(r I U)

n e- BU Tr) n e- U (T)
no e n

-38

actually vanishes. They in turn both yield the same equation

g(f) = e (F)

which one will now notice, via the OZ relation, is simply the HNC equation.

Quantal Generating Functionals

Corresponding quantal equations can be obtained from the above classical generating functionals if we replace the reference distributions--the Boltzmann factor n exp{-bU} describing classical non-interacting particles in the presence of an external potential U(r)--by its quantal extension n*(rlu)--the density distribution of non-interacting fermions in an external potential. (Henceforth non-interacting fermion system quantities will be distinguished by a star.) That is we replace the reference distribution by a single particle Hamiltonian approximation. As in the classical presentation of the above section the functionals are then linearly expanded and evaluated in the limit where the external potential is set equal to the interparticle interaction potential.

We have previously derived a quantal extension of the pair correlation function

FQ{no c(F)} = 1

from the ansatz that the exact classical relation

f6 BU() ) (n(F' IU = V) n ) di' = c(F)
6n(r')

-39

remains valid in the quantal case, alternatively expressed by the ansatz that the classical bootstrap relation is extended as

FQ{n(F(U = n) n} = x( 1

By making use of these relations and defining the analogue potential

SUan(r) n*(P I U)
e = a

a closed set of equations is formed with the QOZ relation. Expressed in terms of previously introduced analogue functions, these integral equations bear a close correspondence with their classical counterparts, and indeed reduce to those counterparts in the appropriate limit.

As a first example we consider the functional

n(r I U)
W[n] = ln n*(r I U)

Its first variation is

Sw n*(r I U) 6(? F') n(r I U) Sn*(r I U)
n(r I U) {n*(r I U) [n*(r I U)]2 6n(r')
Sn(r') In*(r U)] This can be simplified in the following manner:

F6n*(r I U) (F{Sn*(r I U))(F{ 6 BU(r")) (n x*(Q))( Sa n(F') 6 RU(F") 6n(F') n x(Q

= 1 n Fa{c an

where it should be remembered that the analogue pair correlation function differs from the quantally extended "physical" pair correlation through the operator

-40

M{f(F)} = F~Q{x*(Q) F{f(i)}}

From the above considerations we see that

6n*(r I U) nan
6n(F') 0

and so

6Sw an
= c (r-r')
Sn(?')

By Taylor expanding the functional W[n] about U(r) O to v(r), we find

n n*(r I V) = f d?' (n(r I v) n ) ( 6 Sn(r')

or

gan e- Van () ean()

This equation we will call the Quantal Hypernetted Chain or QHNC equation.

Following similar steps with the functional

n(r I U)
W = n*(r I U)

we can arrive at

gan (F) = e- an) + yan()}

or what shall be referred to as the QPY equation.

It is of no trivial importance that the QOZ, QPY and QHNC

equations are of the exact same form as their classical counterparts (this in addition to the property of reducing to the classical equation in the limit). This is because phenomenological ways of incorporating the effects of truncating classical functional

-41

expansions at the linear term may still be of use. [In the classical HNC case such a truncation is equivalent to neglecting bridge diagrams in the Mayer expansion of non-ideal gases.] This will be considered further in Chapter V.

It should also be stressed again that these integral equations, even if they could be made exact, are for gan(r) = n(r/v)/n0 and not the pair distribution function. This is evidenced by the fact that if we turn off interparticle interactions, Van(r) vanishes and the solutions of either the QHNC or the QPY equation yield gan (r) =

1. This is correct for the density distribution but the fermion pair distribution has an exchange hold. This difference is taken up in the next chapter.

The exact correspondence between quantal analogue and classical integral equations is lost when one considers screening effects on long range potentials. Rigorously incorporating screening effects is complicated in the quantal case by symmetry effects. It can be shownl07 that the best single particle Hamiltonian (best in the sense that the statistical average over the square of the difference of the exact and single particle Hamiltonian is minimized) is provided by a Hartree-Fock potential, not the Hartree. On the other hand we will see in the next section that as far as density distributions are concerned there exists an exact effective potential single particle Hamiltonian. In the quantal integral equation approach we will consider here screening of the reference potential will ignore exchange effects; we will see in the next section that they are actually incorporated in the closure of the linear functional expansion, namely the QOZ relation.

-42

When we screen the reference distribution of the QPY generating functional

W = n( U)
n*(r I UH)

with the Hartree potential

OUH(?) = BU(?) + J d?'(MU(F' F) (n(?' I U) no) its linear expansion requires the variation

S r 6n*(r I UH)
Sn(n') 1 n0 no 6n(r' I U)
o o

This can be evaluated by invoking the chain rule over the Hartree potential, which in Fourier space becomes

Sn*(? I UH) Sn*(F I UH) 6 BUH(F")
FQ{- F } Fa
6n(r' I U) 6 RUH(") a 6n(F' I U)

6 OUH(F")
nx*(Q) F{ } 6n(F' I U)

Now the variation of the Hartree potential yields

6BUH() I 6 U)
6n(F') 6n(F') or

6 3BUH() 1
Fat Sn(F') n x(Q) Fa{BU(? F')} from which a few simple manipulations yield

Plugging into the functional expansion of W[n] we obtain

n(F I V) = 1 + J d?'(can(r F') (n(F' I v) n0)
n*(r I VH)

+ f dr'(M{U(r r')}) (n(r' I v) no) In terms of analogue functions and employing the QOZ relation this becomes

V(an
gan = e [1 + an () + F Q 1 {x*(Q) F{n( I v) no} F{BV()}}

SVan
= e H + yan(T) + M{f dT'(RV(T T')) (n(F' I V) no)}]

BRan
= e [1 + yan(F) + M{BVH(F) BV(F)}]

This equation will be referred to as the QPYH equation.

Unlike the classical case, if we screen the QHNC generating functional, we do not recover the QHNC equation. A now familiar procedure yields instead the result

g (r) y(?) + M{IVH(?) BV(?)}
gan() = e e

to be known as the QHNCH equation. We see that for both the QPYH and QHNCH equations long range tails cancel but the potential enters in a nontrivial manner.

Lastly, in the classical case we saw how the HNC equation could also be derived by constructing a screened potential such that the first variation of the generating functional vanished. The same

-44

procedure can be applied to the quantal case. Define an effective potential

BUc(F) = BU(F) f d?'(c(F F')) (n(F' I U) no)

(where c(r) is the quantal pair correlation function and not the analogue quantity); then both the QHNC- and QPY-like generating functionals

W = In n(r I U) W n(T I U)
n*(F I Uc) n*( I Uc)

have vanishing first variations and both yield the same equation:

g- Van(r)
an(F) = e c

This equation is distinct from the QHNC equation but like the QHNC equation reduces to the classical HNC equation in the limit. For this reason it shall be termed the QHNC2 equation.

Interestingly enough the QHNC2 equation can be obtained from an entirely different approach. Instead of assuming the quantal bootstrap relation and utilizing Percus' method we can arrive at the QHNC2 result directly from Kohn-Sham-Mermin density functional theory.107-110 From there we can recover the quantal bootstrap relation (our starting ansatz) and so justify our variant equations by making what we will see as a reasonable assumption concerning the effective potentials of many body systems.

-45

The Density Functional Approach

Following the zero temperature approach of Hohenberg and
Kohnl08 and Kohn and Shaml09 for the ground state properties of an electron gas in an external potential U(r), Mermin showed that110 i) the Helmholtz free energy A of an inhomogeneous electron gas in thermal equilibrium is a unique functional A[n] of the one electron density n(r), and ii) the grand potential

0[n] = A[n] p J n(r) d

is minimum for the equilibrium density n(r). We employ the grand canonical rather than a canonical ensemble because macroscopic quantities derived from a non-interacting Fermi system are more easily calculated and the fundamental variable, namely the density distribution, is a physical observable that should be the same from either ensemble. The free energy is customarily decomposed as A[n] = I J v( F') n(F) n(F') d? d?' + J U(F) n(r) dF + A*[n] + Axc[n]

to separate out respectively 1) solely classical electron-electron interaction contributions, denoted by v(r), 2) the external potential contribution, 3) the contribution A.En] that a non-interacting Fermi system with density n(r) would still exhibit, and 4) all other effects are grouped to define the unkown exchange and correlation functional Axc[n]. The minimum property of the grand potential leads to the functional equation72-73

6A,[n]
--- + Ueff() = (3.5)
6n(r)

-46

where u is the chemical potential and Ueff(r) is defined by

Equation 3.5 is equivalent to the equation for non-interacting electrons in this effective potential; that is we have n(rIU) = n*(rlIUeff). The non-interacting fermion distribution n*(rlU) is given by the following system:

(- 2 A + Ueff) i = ci = M = 1 units

n*(F I Ueff) = f(Ci) I ti(r) 12

where the sum over the index "i" refers to all bound and continuum eigenstates and f(e) is the Fermi-statistics occupation number

f() = {l + exp[B(ci p)]}-1

dependent on the temperature and the chemical potential. The latter is determined as usual from the total number of electrons by spatially integrating over n*(r).

For the jellium model of the electron gas the external potential must not only include the source of the external potential Uo but the effect of the neutralizing background as well

U(F) = U o() I dF' no V(F F') so that the effective potential is

6A [n]
Ueff(F) = Uo(F) + f d?' V(r F') (n( I Uo) n + xc 6n(r)

-47

At large radii, the displaced electron charge vanishes, and any long range tail of the source potential is cancelled on the assumption of perfect screening, so that the assymptotic value of the effective potential is SAxc
Ueff(r n sn

Applying eq. 5 at large radius gives

6A, 6A
S6n0 =- 6n or, equivalently

no = / 3/2 1/

where the standard definition of the Fermi-Dirac functions has been used:

I () i y* dy
S o 1 + ey-n

It is convenient to shift the zero point of our energy scale to define a new effective potential which goes to zero at infinity

Ueff(F) = U o() + f d?' V(i F) (n(? I Uo) n )

A+ ( xc 6Axc n no

Up to this point our equations are exact; the exchange

correlation part of the free energy functional is, however, quite unknown. But if we make a functional expansion to first order

-48

6Axc SAxcI
Sn n S ) + f K(S ?') (n(?'l U ) n ) dr'
n n

K(F T') ( xc

Sn(r) Sn(F')

we obtain an approximate integral equation in terms of the unknown kernal K(Ir-r'I):

n(? I V ) no = n*( I Ueff) n

Ueff = Uo + r dT'{V(F r') + K(r r')} {n(r' I Uo) n} Since the unknown kernal is independent of the external potential U(r), we can determine this function by considering the case of a very weak external potential where the linear response formula can be employed on both sides of the first of the above equations with impunity. The left hand side yields

6n(Q) = FQ{n(r I Uo) n} n0 x(Q) U (0) while the right hand side gives

FQ{n*(r I Ueff) no} noBx*(Q) Ueff(Q)

= n Bx*(Q) {U (0) + [V(Q) + K(Q)] Sn(Q)} Using the last equation we can show that the Fourier transform of the unknown function K(r) is

K(Q) = (- -( -1-_) V(Q)
no0 x(Q) x*(Q)

-49

Due to this assignment the effective potential reduces in form to the previously defined screened potential

U (F) = Uo () f d?' ( c(? -?') (n(r I Uo) nO) and we obtain the relation

n(r I U ) = n*(r I Uc)

Note that the source external potential is completely arbitrary, and so this result is a generalization of the QHNC2 equation, which follows immediately when the external potential is taken to be the interparticle interaction. From this generalized result we can work backwards and obtain the quantal bootstrap relation, and so justify the other (QPY, QHNC, etc) equations we have obtained.

This is done by now making the assumption that the average one body potential felt by a constituent electron due to a fixed test charge Se, namely

Ve-t =Vtest d'( ')) (n(' I Vtest) n) where

V e*6e
test r

is equal to the average one body potential felt by a free test charge 6e, arising from a fixed electron in an electron gas, namely

Vt-e = (6e) W(i)

V2 W(F) = 4n{ eS6() + e[n(? I U = e2/r n ]}

-50

In the limit of small Se we can again invoke linear response with impunity and solve for

FQ{n(? I U = V) n } = -1 x*(Q)

thus recovering the quantal bootstrap relation.

We see that our density functional approach, aside from

justifying the quantal bootstrap ansatz of the Percus method, gives physical insight to its meaning. Even though it is not an exact approach, it clearly points out its limitation, as we approximated the variation of the exchange-correlation free energy functional by a first order functional Taylor expansion. This type of truncation is a property shared by the Percus method and ultimately may be responsible for the self-consistency of the two approaches.

CHAPTER IV

THE BASIC QHNC EQUATION

An alternative ansatz for the extension of classical integral equations into the quantal regime lies in the construction of effective pair potentials. These incorporate the bulk of quantal effects; for example in multicomponent plasmas quantal effects prevent the collapse of particles in an attractive coulomb field. More generally the effective potential yields finite values for the pair distribution factor g(r) in the limit of zero interparticle separation

r.

Such classical pseudopotentials have been constructed by

requesting an approximate account of the N-body Slater sum11-112 or, in an approach which can in principle be implemented exactly but is more limited in its physical range of applicability, by reproducing the two-body Slater sum113-118

BVeff(r) Hel
e (r) = (42 /m)3/2 (4.1)

Hrel f 2 + V(r)
2 m

It can be shown that this Binary Slater Sum approximation correctly represents the derivative of the classical pseudopotential for the N-body problem in the limit r tends to zero.119

-51-

-52

As shown in the previous chapter, the QHNC equation for n(rlv) is also equivalent to the classical HNC equation with an effective potential defined by

BVeff(r) n*(r I V)
e
no

where n*(rlv) is the density distribution of non-interacting fermions in the presence of an external potential of the form of the interparticle interaction--the bare coulomb potential. Normally this entails the solution of the Kohn-Sham-Mermin equations.108-110 That is calculating n*(rlv) is the solution of the following system:120-121

(- 2V + V) i = ci i

n*(r) = f(ci) i 1'i (4.2)

where f(ei) is the Fermi-statistics occupation number

f(ci) = {1 + exp[B(ci p)]}depending on the temperature and the chemical potential fixed by the total number of particles. The index "i" runs over all continuum and bound (if existent) states.

The treatment of the delocalized states proposed by Friedel for impurities at zero temperaturel2 can be applied to the nonzero temperature case too. The eigenfunctions F KLM of eq. 4.2 for energy

sK = 2m

-53

are

TKlm(r) = AKim *Kl(r) Ylm(r)

Kl(r) is normalized by requiring that it have the asymptotic form

K1(r) sin[Kr e i + n(K)]

n1(K) being the phase shifts in the potential V(r) and 1n(K = = 0. Overall normalization is accomplished by requiring that the wavefunctions vanish at a spherical wall of very large radius R. We assume that R is so large that most of the integrand is approximated by its asymptotic form:

1=1 AK1 2 j dr sin2[Kr e i + 1]
0

with the boundary conditions

KR + 1l(K) = e = m m = 0, 1, 2, 3, ....

incrementing m by unity for a given L yields the density of states (including spin multiplicity) as

2R 1 a
91(K) = (1 +R aK nl(K))

while the normalization factor for very large R is given by

AKI = v2/R (1 -0( ))

(These same results can be derived for a long range coulombic potential.123 cf Eq. 36.23) A straightforward derivation gives for the density

n*(r) n = Xf(cb) b b ~ f f(c) dK I (2e + 1) [ 2K1- K 2 (Kr)]
b Ir 0 1

-54

where subtracting off the uniform density result accelerates convergence of the integral and summation. We note that in general a potential may support a number of bound states and we have included them also.

The OHNC Effective Potential

We thus see that the QHNC pseudopotential is essentially a

one-particle Schrodinger equation calculation. This is a feature in common with the Binary Slater Sum. It is important to contrast the two however. The Binary Slater sum (BSS) is used to calculate the pair distribution function g(r) while the QHNC equation yields n*(rlv); this difference is important in consideration of quantum effects. It should be noted that in the BSS approach the Schrodinger equation is solved with the use of the reduced mass. The matrix element in eq. 4.1 can be evaluated by inserting a physical basis set of properly symmetrized or antisymmetrized wavefunctions, depending on the statistics of the two particles. Parallel spin versus antiparallel spin contributions can be calculated separately. In the BSS approach the effective potentials are density independent. In the QHNC approach the Kohn-Sham system of equations can also be written in the form of a matrix element

n*(r I V) =

where the operator "f" is defined as

f = (n + exp B[H ull-1

-55

H is the coulomb hamiltonian operator, mu the chemical potential and we have included a variable "eta" to be set equal to unity to describe Fermi statistics. In the QHNC approach the external potential has no spin attribute, and the matrix element is evaluated using a complete basis set with no regard to statistics; the statistics are built into the functional form of the matrix element. The QHNC effective potentials are implicitly dependent on the density through the chemical potential.

All these subtle differences are manifestations of the fact that in the QHNC approach g(r) must be obtained from n*(rlv) by a further calculation. One obvious advantage into this splitting of tasks is that the BSS pseudopotential breaks down at zero temperature, whereas the QHNC remains valid.

A surprising feature of the QHNC effective potential is that it is analytically calculable. This allows one to bypass the computationally time consuming Kohn-Sham-Mermin system of equations. The starting point for the solution is the fact that the effective potential is related to the I 2 limit of the off diagonal density matrix for fermions in an external coulomb potential:107

P(r I 1 = < l I f I 2

Spherical symmetry reduces the functional dependance of the off diagonal density matrix does from the full six variables of r1 and 2 to only the magnitudes r1, r2 and the angle between rl and 2. This is easily seen if we expand the matrix element using a complete set of spherical coulomb wavefunctions:

-56

SKI(r) 2K2
mKlm = Y (r) =2m

2 2
d S + [K2 2m (1 + 1) S = o
dr 6 r

Using the addition theorem of spherical harmonics then yields (2 1 + 1) PI(rl I r2)
P(11 I 4n r1 r2 P1(cos()

where the 1-wave radial density matrix components are defined by

l(rl I r2) = (spin dK) 22 SK1(rl) SK1(r2) e 2m +

Furthermore, because the potential is coulombic, there is an additional symmetry: the matrix element "f" also commutes with the Runge-Lenz vector.124 This reduces the functional dependance of the off diagonal density matrix to the two variables "x" and "y":

x = (r1 + r2)/2 y = (r2 + r2 2rl r2 cose) /2/2

The solution then follows the steps of Hostler125-126 and Storer127 except that the Bloch equation satisfied by the off-diagonal canonical density matrix is replaced by a Fermi statistics generalization

[(H p) (1 +n a + a] P(r1 I r2) = o

(Here H acts only on the r1 coordinates and we set eta equal to unity at the end of the calculation.) Because of the aforementioned symmetries, the solution to this equation reduces to a functional form identical with the s-wave radial density matrix component (however in the variables x and y):

-57

P(rI 1 r2) 1 ay Po(x + y I x -y) The diagonal limit rl = r2 = r is given by

n*(r) = p( ) 8 ay2 Po(x + ) x y) )

The second derivatives of the radial coulomb wavefunctions can be eliminated in favor of the radial wavefunctions themselves by employing the radial coulomb differential equation. He thus obtain

n*(r) = sUn f dK (K, r)
212 eB I2K2 r
22m ] + 1

2 2m
J(K, r) { dp SKo(r) (2 K2) I SKo(r)2 2rfir

Numerical Calculations

The effective potential, obtained from the ideal fermion density distribution in a coulomb potential, will be calculated from the analytic formula derived above. For comparison with the classical HNC equation, we will scale lengths in units of the interparticle separation x r/a where

no = (ra3 )Temperature/density points will be expressed by the classical coupling parameter

r =
a

-58

along with a quanticity temperature (two divided by the temperature in rydbergs)

B* E me4/fi2

In the classical limit all functions are solely dependent on the parameter r6; we expect any dependance on the quanticity temperature to indicate the pervasiveness of quantum effects.

The calculations were performed corresponding to the temperature density points of tabulated results provided by Pokrant69 for later comparison. These consisted of values at constant density, as measured by the dimensionless ground state coupling parameter

rs = ala

where ao is the bohr radius, as well as a dimensionless degeneracy temperature

t = Kt/EF

where Ef is the ground state Fermi energy. The results are presented in the following graphs.

In Fig. 4.1 we present plots of the ideal fermion density in the presence of an external coulomb potential. For comparison these are plotted together with two approximations. The first is merely the classical limit

n*(r I V)
no

The second is the Thomas-Fermi approximation.107 It is obtained by expanding the diagonal density matrix element "f" using a complete set

-59

of plane waves. If we then neglect the non-commutativity of the potential and kinetic energy operators in the exponential (the neglected terms can be systematically calculated by Trotter formulas128 and shown to vanish for large r) we then obtain

n*(r I V) F 1/2 (3p r/x)
n F1/2(3P)

where the Fermi-Dirac integral of order p is defined by the integral

Fp r r(P + 1) f dy -x o 1 + e-x

(In this formula the symbol r denotes the gamma function.) It is seen that the exact result quickly merges with the Thomas-Fermi approximation.

By using the relation

dF (x)
dx Fp-1(X)

we find that asymptotically the Thomas-Fermi approximation has the form F 1 (Op)
n*(r I V) r 2
n x F (1P)

This is indeed the correct asymptote of the analytic solution as can be verified by looking at the appropriate expansion of the spherical coulomb wavefunctions.129 This tells us that the effective potential in the QHNC approximation has a long range l/r tail--but with a strength modified from the classical theory. The QHNC tail contains the additional factor

-60

F 1 (p)/F P)
2 2

not found in the Binary Slater Sum approach. In Fig. 4.2 we present plots for the QHNC effective potential. For comparison the corresponding effective potentials of the Binary Slater Sum approximation are presented in Fig. 4.3.

The approach to solving the HNC system of equations with long tail potentials is well known--we follow the convention of Ng105 and renormalize the constituent functions by subtracting out an analytically transformable long range function of the form

z erFc(ar)
r

We employed the Ng method of solution (essentially a picard iteration with an accelerated convergence procedure) although it should be mentioned that two new algorithms have recently been published130-131 that employ a hybrid Picard/Newton-Raphson iteration procedure which promises faster convergence.

The results of the integral equation solver, that is the analogue pair distribution function n(r/v)/no, are presented in Figs. 4.4-4.8. This distribution is of theoretical interest in its own right, for example in the study of impurities in metals or positron annihilation rates.132 However we will limit our concern to the approximate calculation of the pair distribution function. A method of calculating g(r) based in part on the information contained in n(rlv)/no, and a comparison of results with the Binary Slater Sum approximation and that of Pokrant will be taken up in the next two chapters.

-61

0.

storer

.... boltzman .1 -. rs- 2.

'0. tau 2g.8 OD .6-,. J"
4.)- tau- i 1.68

tau- 6.5e

tau- 6. t

-1.0 i
0. .5 1.0 1.5 2.0 2.5 3.0 3.5 x-r/a

Figure 4.1 Plots of n*(rlv)/no-1 versus x = r/a for a density of rs = 2 at temperatures of kt/Ef of .1, .5, 1.0, 2.0, 20.0. Storer denotes the analytic solution, following Russian literature the Thomas-Fermi approximation has been denoted Hartree, and Boltzmann is the classical approximation.

-62

1.0

storer

.8 har tree

S r-. 2.88 tau-20.88

4 tau- 2.88
x
tau- 1.88

2 -- tau- 8.58 tau- 8. 18

0 .I I I
0. .5 1.0 1.5 2.0 2.5 3.0 3.5

x-r/a

Figure 4.2

Plots of the QHNC effective potential versus x = r/a for a density of rs 2 at temperatures of kt/Ef of .1, .5, 1., 2., 20. We have plotted the potential in the form

xr Veff(x)

A purely classical coulombic interaction would correspond to a constant values of unity. Storer represents the analytic solution; Hartree is the Russian literature denotation of the Thomas-Fermi Approximation, which is truncated near the origin to prevent divergences.

-63

effective potential

1.2

slater sum

1.00
tau-5.00

tau-2.00 tau- 1.00

tau-0. 50 .2
tau-*. 10

0. .5 1.0 1.5 2.0 2.5 3.0 x-r/a

Figure 4.3 Plots of the Binary Slater Sum effective potential versus x r/a for a density of rs = 2 at temperatures of kt/Ef of .1, .5, 1.,
2., 3., 5. We have plotted the potential in the form

x V ff(K)

A purely classical coulombic interaction would correspond to a constant values of unity. Note that as we approach classical temperatures the curve reaches this constant value quickly, and in fact first overshoots this value.

-64

static structure factor

hnc

-. free chihera
teu-.18
-. 4 rs-2. e

-.8

-1.0I
0. 2. L. 6.

k-q*a

Figure 4.4

Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rs = 2 at temperatures of kt/Ef of .1. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

-65

static structure factor

hnc

free
-.2

Chhare
tau-a.5 S
4 rs-2. 8
-,

-.6

-.8

-1.0
0. 2. q. 6.

k-q-a

Figure 4.5

Plots of density times the Fourier transform of h(r) g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rs = 2 at temperatures of kt/Ef of 0.5. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

-66

static structure factor

0.

hnc

-.2 ---- f ree

ch her a
tau l.89
S4, rs 2. 88

-.6

-.8

-1.0 -0
0. 2. 4. 6.

Figure 4.6

Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rs = 2 at temperatures of kt/Ef of 1.0. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(rlv)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

-67

static structure factor

-.2- tr ee cahira
tau-2.88
.9 rs-2. e

-.

-.6

-.8

-1.0
0. 2. L. 6.

k-q*a

Figure 4.7

Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rs = 2 at temperatures of kt/Ef of 2.0. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

-68

static structure factor

---------------------hnc

-.2- rree Chi hara
tu-8. 88
rs-2.00

0

-.6

-.8

-1.0 i I
0. 2. 4. 6.

k-q-a

Figure 4.8

Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rs = 2 at temperatures of kt/Ef of 8.0. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

CHAPTER V

LOCAL FIELD CORRECTIONS

We have seen in the previous chapter how classical integral equations generalized into the quantal region can yield accurate results for the quantity n(rlv), namely the density distribution of interacting particles in the presence of an external potential of the form of the interparticle interaction. This in effect is equivalent to knowledge of the static response function x(k). If from this information the dynamic response x(k,w) could be extracted, then we could easily calculate the static structure factor or pair distribution function via the fluctuation dissipation theorem

S(K) = f S(Kw) dw = f Im Bx(Kw) dw
e 1

These operations would constitute the quantal generalization of the classical bootstrap assignment g(r) = n(rilv)/nO.

To demonstrate that such a program is plausible, we note that the analytical properties of the dynamic susceptibility make it possible to write an exact "dispersion" relation for x-1(k,z) of the form133

x (KZ) = x- (K) + 2 [Z2 + iZK2 G(KZ)] pK

-69-

-70

This particular relation, which is one of many possible representations, introduces an unknown function e(k,z) which is analytic in the upper half plane. The physically meaningful dynamic response is given in the limit z = w + ie, and in this limit the real and imaginary parts of

O(K, Z = w + ic) = O'(Kw) + iQ"(Kw)

are, like the dynamic susceptibility itself, related by the Kramers-Kronig relations. The function e'(k,w) is known to be a positive, even function of frequency, as follows from the known parity of x'(k,w) and x"(k,w) and the fact that S(k,w)> 0. The utility of the function e(k,z) is that it interpolates smoothly between the known results for x'(k,z) in the limits of large and small z, and being less sensitive than x(k,z) itself, is more amenable to approximation.

By subtracting out a corresponding expression for a noninteracting Fermi system we can express the dynamic susceptibility as xo(KZ)
x(KZ) = (5.1)
1 Veff (KZ) x (KZ)

Here we have grouped all of the actual physics into a dynamic effective potential Veff(k,z), actually a functional of the static response x(k) and the unknown function D(k,z). We will show through the use of the Mori projection operator technique134 that under reasonable assumptions the dynamic effective potential can be approximated by the static (z = 0) limit of the above equation

1 1
Veff(KZ) (x(K) xK)

-71

It should be noted that although x (k,z) in eq. 5.1 is

usually taken to be the free response function, this choice is not a unique one; xo(k,z) plays a role of a reference system in terms of which the expression for x(k,z) is constructed.135-136 A possible improvement in the procedure may be achieved by using an adjusted choice of x (k,z) that includes self-energy effects, for example.14-19'137 This method of improvement will not be taken up in this thesis in lieu of other avenues of extension. We will retain the free response as our reference system; this will allow us comparison with certain mean field theories.

Mean Field Theories

Mean Field Theories bypass the dependence of veff(k,z) on x(k) and directly try to relate Veff to the bare interparticle potential by solving approximately the many body problem. (Something we have already done via our integral equations.) The efective potential is usually taken to be "static," that is independent of z, and is commonly expressed in the form

Veff(K) V (K) [1 G(K)]

where V (k) is the Fourier transform of the bare interparticle potential and G(k) is called the local field correction factor. A review of such local field function theories can be found in Kugler.138 The many body part of the problem is generally solved by the Vlasov kinetic equation139 or, what amounts to the same thing, appropriately decoupling the equations of motion for the density

-72

fluctuations.140 Whatever route, at some point collisions between particles are neglected and the particles are assumed to move in a mean or average field based on the collective motions of all the other particles. The physical significance of the correction factor is thus to incorporate some of the effects of short range correlations between individual particles not taken care of by the average field.

Mean field theories for the most part have been limited to the ground state quantal electron gas. These have disagreed on the form of the local field correction factor on three main points.

First some workers141-145 have been led to the conclusion that G(k) is a universal function of k/kf (kf the Fermi wave vector), while others146-149'29-31 contend that G(k) is density dependent.

Second, there is wide diversity of opinion regarding the value of G(k) in the limit of large wavevector k.150-151 At this limit, a constant 1/3, irrespective of density, has been given by Geldart and Taylor,142 Rajagopall144 and others,143,145-147 while that of Togio and Woodruffl48 and Kugler38 use 2/3. On the other hand Singwi et al.29,30 derive this value to be l-g(O) in terms of the pair distribution function g(r) at the origin, and that of Vashista and Singwi31 is

1 g(o) ap(a )

(where a is an adjustable parameter determined self-consistently) in contrast with the result of Niklasson149 who gives 2/3(1 g(O)) to this limit.

Thirdly, there are several authors138,146-147,142 who have shown the local field factor G(k) has a maximum around k = 2kf.

-73

Kugler138 among them further indicates the slope of G(k) has a logarithmic singularity at k = 2kf in addition to a maximum. On the contrary, G(k) of Singwi and collaborators29-31 has neither maximum nor singularity in its slope around k = 2kf. It should be noted that the behavior of G(k) in the large wave vector region is quite extremely sensitive to small differences in x(k) so that it offers a stringent test for the distinction of various theories.

In the case of the degenerate electron gas, there is neither

computer simulation nor experiment to be used as a criterion to check these diverse conclusions concerning G(k). However, at present the prescription of Singwi et al.29-30 is regarded as fairly successful in treating a degenerate electron gas and is relatively easy to program numerically. It encompasses in the classical limit, as do the integral equation methods of Chapter II, the Hyper-Netted Chain equation, which is known to give the best fit to the pair distribution function of computer simulations. For this reason the derivation of the Local field correction factor in the STLS theory will be presented here as a basis for later comparisons.

The STLS Method

Following Singwi8 we will derive the density response function by following the time evolution of the off diagonal single particle density function

p (x I x" t) -a < ( t) t(',t)>

-74

If we employ a Heisenberg equation of motion approach with a Hamiltonian containing an external potential Vext(x,t), we obtain

{i mA A + VH(it) VH(i't)} pa( I '; t)

S d [V(R ") V(' ")] <,F; (xt) n(R"t) T (R' t)> where the Hartree potential

VH(it) = Vext(it) + f di' V ( i") p(i", t)

depends on the interparticle interaction potential Vo and the average single particle density distribution

p(xt) = E < I +(xt) T (Xt)>

We were able to regroup the terms of the equation of motion into its present form by defining the cumulant bracket

=
Sa c

a a a a Our equation of motion is equivalent to an infinite set of equations of motion for observable physical quantities with classical equivalents. This follows from the fact that the Wigner distribution function96,152

f (Rt) = J d? e- ipr/i <+ ( 1 t) (R + t> pa at) (R+ r, t)> 2

-75

has properties analogous to the classical phase space distribution f(r,p;t). In particular it can be used to construct physical observables such as the single particle density distribution

p(Rt) = I f (Rt)
pa

the particle current density

(Rt) = m- P f (Rt)
po
po
and the kinetic stress tensor
'(Rt)=m- P P f (Rt)
pa
pa

By expanding the equation of motion for the off diagonal density distribution about its diagonal, that is in powers of r = x' x, the coefficients of the first few powers of r yield the beginning of a hierarchy of one particle equations. First we have the usual continuity equation

at p(Rt) + *= (Rt) = o

second the equation of motion for the current density

m t J(Rt) = ~.k(t) p(Rt) VH(t) I d IV(R R)

The last term in the above equation describes through the

correlation function c all the complicated effects of the Paull and Coulomb hole surrounding each electron (in the presence of the external field). In lieu of an exact evaluation of this term or main aim is to extract a local field correction. Noting that

c = n(Rt) n(it) {g(R, i; t) 1}

-76

where g(r,x;t) is the non-equilibrium pair correlation function, the simplest approximation is to replace this function by its equilibrium value g(r x), thus accounting for the local depletion of charge density but neglecting the dynamics of the hole.

This approximation replaces the Hartree potential by a local effective potential given by

Veff(O) = VH(Kw) G(K) p(Ew) with

G() = (23 q[S -) 11 (5.2)
(2r)3 q

Relating ii (k,w) to thas esectave potential An the same manner as a free particle calculation, namely as ( P f (p K/2) fo(p + K/2) ir (Kw) = 2 V (Kw)
OC pa p*R/m + i e
then the continuity and current density equations together yield

x(Kw) = x (Kw)/l 4re2 (1 -G(K)) x (Kw)

where xo(k,w) is the Lindhart polarizability153 of an ideal Fermi gas

f (p K/2) f (p + K/2)
xo(K) = X o
pa W p R/m + ic

Although G(k) contains the unknown structure factor S(k) of the electron fluid, this can be determined by requiring consistency with the structure factor obtained through the use of the fluctuation dissipation theorem. Of use numerically is the little known direct relation between G(q) and the pair distribution function g(r)154

-77

G(Q)= 1 Q f dr g(r) jl(Qr)
o

and by using orthogonality of the spherical bessel functions

2r2
g(r) = 1 f o dQ Q G(Q) jl(Qr)
o

These results are extremely useful because they enable one to generate a G(q) directly from the pair correlation function and to test G(q) by constructing the corresponding g(r).

In the STLS scheme both exchange and correlation corrections to the Hartree field are automatically taken into account through its self-consistency. Moreover it is straightforward to recover the Hubbard141 result

1 K2
G(K) 2 2
K2 + K2

by substituting for the structure factor in eq. 5.2 its Hartree-Fock value. The Random Phase Approximation, which is the dynamic extension of the Debye-Huckel theory, can be viewed in this scheme as the approximation that g(r,x;t) = 1, and leads to the result G(k) = 0.

At zero temperature a further constraint can be used to check the validity of any approximate microscopic theory. The compressibility "sum rule" states that the long wavelength limit of the static response should agree with the result obtained for the compressibility from differentiation of the ground state energy. (This and other ground state properties of the quantal electron gas are covered in Ichimaru.5 Both the Hubbard and STLS equations violate this

-78

constraint. This has led to modifications of the Hubbard correction
factorl155-157

1 K2
2
G(K) =
K + aKF

and of the STLS scheme by Vashista and Singwi31

G(K) = (1 + a ) {- 1f [S(K Q; p) 11]} (5.2)

where the additional free parameter "a" is to be varied such that compressibility sum rule is nearly satisfied. However it has been shown79 that if the correction factor G(k,w) is static, it is impossible to satisfy the compressibility sum rule and the third frequency sum rule simultaneously (the f sum rule being denoted as the first moment).

The Relaxation Function

As an alternative approach to calculating the correction factor G(k,w), which unlike the methods covered so far easily generalizes to finite temperatures, we shall consider approximations to the Kubo function provided by the Mori continued fraction method.134

In Chapter I the Kubo function was identified from the difference of the dynamic and static density response functions:

x(KZ) = x(K) + iRZ C(KZ)

In this section we shall demonstrate that this follows from the definition of the Kubo function as the equilibrium averaged density auto-correlation function

-79

Co XH -XH
C,(r r', t t') j dX

where <.... > stands for the thermal expectation average.

It is easily shown that the Kubo function vanishes at large

times; spatial Fourier transformability is insured by the subtracted term . The intimate relationship of the Kubo function with the density correlation function can by exposed by writing the Kubo function in the form

C(r F', t t') = where we defined the Kubo transform of an operator as

1H
J fJ dX e J eH

By taking matrix elements diagonal in the Hamiltonian of the Kubo transformed operator, it is easily seen that the Kubo function reduces in the high temperature limit to the density correlation function. (The temperature = 0 limit however is ill-defined). Other notations abound in the literature; in lieu of the Kubo transform we shall employ the semicolon bracketl4

THE ELECTRON LIQUID AT ANY DEGENERACY By BRIAN GREGORY NILSON A DISSERATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1987

PAGE 2

TABLE OF CONTENTS PAGE ABSTRACT iv CHAPTER I INTRODUCTION 1 II BACKGROUND THEORY /. 7 The System 8 Basics of Linear Response 9 The Dynamic Density Response 11 Relation to Density Correlation Functions 16 Further Relations 17 Distribution Functions [ 19 III OUANTAL FUNCTIONAL EXPANSIONS 27 The QOZ Relation 27 Classical Generating Functionals 32 Ouantal Generating Functionals 33 Density Functional Approach 45 IV THE BASIC QHNC EQUATION 51 The QHNC Effective Potential 54 Numerical Calculations 57 V LOCAL FIELD CORRECTIONS 69 Mean Field Theories 71 The STLS Method '.'.'.'.'.'.'.'.['.'.'.'.]'. 73 The Relaxation Function 78 Scalar Products of Operators 81 Microscopic Theory 86 VI NUMERICALLY EVALUATION S(Q) 93 The Lindhard Function 97 Convergence Acceleration 103

PAGE 3

VII THERMODYNAMIC PROPERTIES Ill Thermodynamics with Slater Type Effective Potentials 116 Interaction Energy Results 118 Kinetic Energy Results 120 VIII RESULTS OF THE EXTENDED QHNC EQUATIONS 131 Method of Solution 133 Quantal Hartree Results 141 The Zwanzig Equation 144 Blending 147 IX DYNAMICS OF THE LFCF 189 Calculating the Memory Function 197 Results 200 X CONCLUSIONS 219 REFERENCES 224 BIOGRAPHICAL SKETCH 235

PAGE 4

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy THE ELECTRON LIQUID AT ANY DEGENERACY Brian Gregory Wilson May 1987 Chairman: Charles F. Hooper, Jr. Major Department: Physics The dielectric formulation of the many-body problem is applied to the study of the static correlations in electron liquids at non-zero temperatures. The intermediate coupling effects arising from the exchange and Coulomb correlations are taken into account through a local-field correction factor obtained from quantal extensions of classical fluid integral equations. This results in a unified theory which can describe the thermodynamic and dielectric properties of the -i V-

PAGE 5

electron liquid covering the limits of the zero temperature jellium to the classical one-component plasma. Detailed numerical calculations are provided for the self-consistent set of equations describing the static structure factor. The variation of interaction energy per particle versus temperature at constant density is shown to exhibit a peaked structure near absolute zero. This feature is not observed using effective potential approaches in classical fluid equations. A comparison of zero temperature results with those obtained from the Greens function Monte-Carlo method is good; however, agreement worsens when dynamical effects of the local-field correction factor are included using the Mori continued fraction method. -V-

PAGE 6

CHAPTER I INTRODUCTION The one component plasma (OCP) is a system consisting of a single species of charged point particles together with a uniform oppositely charged rigid background to ensure charge neutrality. It was first introduced by Sommerfeld as an approximation to metals when band structure effects are negligible. By coupling the standard adiabatic (Born-Oppenheimer) approximation with an assumption of a weak electronion pseudopotential interaction, it can be shown^"^ that the thermodynamics of a simple metal is the sum of that due to a one component electron plasma and that due to a system of classical ions interacting via a state-dependent pair potential. Viewed as a model fluid, the OCP exhibits the characteristics which distinguish real coulomb systems, such as plasmas and electrolyte solutions, from neutral fluids. These include the phenomena of screening and plasma oscillations, which arise from the long-range nature of the coulomb interaction. Thus the OCP is a standard approximation for both astrophysical and laboratory fusion plasmas,^ where the electrons are treated as a polarizable background for the gas of ions. Further applications lie within the local density approximation of the density functional theory of electronic structure, where the excess chemical potential of a one-component electron plasma plays the role of an

PAGE 7

-2exchange-correlation potential in the one-electron Schroedinger equation. It is thus not surprising that there has been a large amount of work expended in obtaining accurate theoretical descriptions of both the zero temperature fermion OCP (degenerate electron liquid) and the high temperature classical OCP (electron gas).^"^ Active research is ongoing in three main areas. First, in the course of study of linear dispersion of plasmons in the simple organic polymers polyacetyl ene and polydiacetylene, it has proved important to contrast the properties of the exchange hole, and its screening, in quasi-one-dimensional solids with those in isotropic three-dimensional jellium (classical OCP).^'^ Secondly, one has been able to determine experimentally, using inelastic X-ray scattering, the dynamical structure factor S(q, w) for the electrons in a few metals. ^^"^^ It is of particular interest that the experimental S(q, w) show a tendency of two-peak structure for q within the particle-hole continuum, a feature believed to be a property of the uniform electron liquid not presently accounted for accurately. ^^"^^ Lastly renewed attempts are being made to quantify the intermediate degeneracy regime of the fermion plasma. There are, in fact, many physical systems whose electrons exist in partially degenerate states and are inadequately described by either zero temperature or classical formalisms. For example, the finite temperature thermodynamic and dielectric properties of the electron liquid are needed for a proper treatment of the equation of state of high temperature liquid metals forund in shock wave experiments^^ and experiments on liquid metals expanded to near their critical

PAGE 8

-324-25 points, of high density inertially confined fusion plasmas, and of the finite temperature exchange-correlation potentials used in density functional calculations of atomic properties at high pressure and temperature. Semidegenerate electron liquids are also found in semiconductors where the low densities involved lead to correspondingly low Fermi temperatures. The finite temperature equation of state of the electron liquid is also needed for a quantitive explanation of the miscibility gap in solutions of alkali metals in their alkal i-hal ide melts. The problem of the OCP at intermediate degeneracy has only recently received attention because of complications not present in the extremes of a classical or a ground state quantum description. For example at zero temperature diagrammatic, self-consistent distribution function methods^^"^^ and pseudo integral equation methods, ^^"^'^ together with recent Quantum Monte Carlo^^"^ calculations, have combined to present an accurate description of the paramagnetic electron 1 iquid.'*^"'*^ But at finite temperature we no longer deal with a unique ground state (invalidating the Monte Carlo approach) nor may we utilize the ground state energy variational principle (invalidating Fermi Hypernetted Chain approaches).* Diagrammatic approaches may be generalized to finite temperatures^^ Fermi Hypernetted Chain calculations of quantum fluids^ are actually variational calculations of the ground state using an adjustable Slater-Jastrow type wavefunction and should not be confused with the classical integral equation method per se. Rather the name stems from a diagramatic expansion of the ground state distribution function which is regrouped into a structure of the same form as the classical HNC equation.

PAGE 9

-4but are perturbation expansions reliable only in the regime of weak 44 coupling. On the other hand the classical limit may be handled by a variety of techniques, including the Monte Carlo method/^"^ the molecular dynamics method, and various approximate methods involving integral equations, most notable amongst these the Hypernetted-Chain equation^^"^^ and modifications thereof .^^"^^ Here quantal effects introduce two complications. First a quantal calculation of the pair distribution function is hampered fundamentally by the fact that the canonical partition function can no longer be factored into a solely configuration-integral part and a part involving solely momentum coordinates. Secondly the thermodynamic description will also depend on the momentum distribution (or off diagonal single body density matrix)^^"^^ and this information is not trivially contained in the pair distribution function. Various attempts have been made to overcome these complications within the framework of the successful methods available for classical fluids. Wigner-Kirkwood expansions^^"^^ may be used for those cases where quantum corrections to the classical partition function are perturbative, but their convergence as a series in temperature or Plancks constant is slow (if existent) and is restricted to short ranged potentials;^^ it is, in fact, an ill defined expansion for 68 the OCP. Published work often centers around diffraction effects and assumes Boltzmann statistics. Previous calculations of the OCP at intermediate degeneracy and coupling have been of several types. The most obvious course has been

PAGE 10

-5to approximate the true quantum mechanical Slater sum which appears in the canonical partition function with the corresponding classical expression involving a temperature and density dependent effective potential. Pokrant^^ and others^^"^^ have developed a finite temperature variational principle to approximate the Slater sum as a product of an ideal fermion Slater sum and that of pairwise additive potentials. The energy is calculated by approximating three-body distribution functions in terms of pair distribution functions, which are obtained from the HNC equation using the effective potentials. Zero temperature correlation energies from this method are, however, in substantial disagreement with the parameterized fittings of the quantum Monte-Carlo results. '^^"^^ Several other groups have studied the thermodynamic^^"^^*^^"^^ and dielectric^^'''^ properties of the OCP by using finite temperature perturbation theory within the random phase approximation. In particular Perot and Dharma-Wardana^^ and Kanhere 22 et al. have presented closed form parameterizations of the exchange-correlation energy and chemical potential within the random phase approximation (RPA) for the paramagnetic and spin polarized cases. However the RPA is a weak coupling theory that results in considerable error at metallic densities. 20 Dandrea et al attempt to overcome the limitations of the RPA by including a static local field correction factor (LFCF) within the framework of the dielectric formulation. They use an approximate form for the LFCF which has two adjustable parameters. The first is fixed by assuming a value for the pair distribution function at the radial

PAGE 11

-6origin. The second is obtained from an appropriate interpolation between known Monte-Carlo results at zero temperature^^ ""^^ and in the classical regime^'^"^^ for the bulk modulus. Tanaka et al.'^^ approximate the LFCF at all temperatures by the static zero temperature LFCF of Singwi-Tosi-Land-Sjolander (STLS) self-consistency scheme. ^^"^^ This thesis concerns itself with a unified theory capable of describing the dynamic density response function of the OCP throughout a broad range of densities and temperatures: from a fully degenerate quantum plasma at absolute zero temperature to a nondegenerate classical gas. It differs from the approach of Dandrea and Ashcroft or Tanaka et al in two important respects. First, that the theory be particularly applicable to electron liquids at metallic densities, which are not weakly coupled, perturbative local field corrections to the random phase approximation are replaced altogether by a selfconsistent field determined by a quantal extension of classical integral equations. Formally we can define a LFCF that, unlike Dandrea et al requires no known results, and unlike Tanaka et al., is temperature dependent; indeed it may be viewed as a generalization of the STLS LFCF of zero temperature and reduces to the proper classical result. Secondly, dynamical dependencies of the LFCF can be rigorously included. This has been shown to be a mechanism for modeling the double peaked shape of the dynamic structure factor observed in metals. ^^"^^ It is also an essential theoretical featurestatic LFCFs cannot simultaneously satisfy the compressibility sum rule and the third frequency moment sum rule for the density response.

PAGE 12

CHAPTER II BACKGROUND THEORY This thesis will concern itself solely with the OCP as a Fermi liquid. The most basic characteristic of a liquid is that it possesses short range order as opposed to the long range periodicity of a crystalline solid. Since the structure of a crystal is determined experimentally by observing the Bragg reflection of X-rays, it is natural to seek a quantitative description of the liquid structure via the intensity of X-ray, thermal neutron, or light 3 80Â—81 scattering. Of particular interest then are the fluctuations of space and time dependent densities which describe, at long wavelength, the cooperative motion of many particle systems. For this purpose define the particle density correlation function as S^pCr ?'. t) = <6n(r, t) 5n(r. o)>^^ where 6n(r,t) is the excess particle density operator 6n(r, t) = n(r, t) <6n(r, t)> eq and the equilibrium average is indicated by the brackets <. > eqSpecifically, even though gq = 0, there will be spontaneous, usually small, fluctuations on a local scale. The density correlation function is a measure of these fluctuations. -7-

PAGE 13

-8In this chapter we will begin our investigation into the structure of quantum fluids by reviewing the general properties of the time correlation functions related to the particle density correlation function. These include 1) various symmetries (for example time reversal), 2) positivity properties which are related to the dynamical stability of the system, 3) fluctuation dissipation theorems which connect spontaneous fluctuations and energy dissipation in a thermally equilibriated system, 4) Kramers-Kronig relations which express causality, like the well known one between the index of refraction and the absorption coefficient, and 5) sum rules which provide restrictions that any approximate microscopic theory must fulfill. The System Because we will be dealing with a many particle system from a quantum statistical point of view, we will employ the second quantization representation of quantum mechanics even though we will mostly be working in a canonical ensemble where the number of particles is fixed. For our equilibrium "system" we shall define in the Schroedinger (time dependent state) representation the unperturbed Hamiltonian H = T + V in terms of field creation and annihilation operators. The kinetic energy operator is

PAGE 14

-9T = J dl ^""(I) {} ^(I) while the interaction energy operator is given by V = ^ J ; dl dll f^il) = tr p(t) n(f) tr p(t) = 1 spins Basics of Linear Response 82-85

PAGE 15

-10What follows is entirely parallel to traditional derivations of time dependent perturbation theory in elementary quantum mechanics texts. We have to solve the time evolution equation for the density matrix ifi ^ p(t) = [H(t), p(t)] = [H, p(t)] + [6H(t), p(t)] (note the commutator is reversed from the Heisenberg equation of motion for operators) subject to the initial condition P(t = Â— ) = [H. p^q] = 0 The initial condition expresses the fact that the system is stationary before the external field is turned onÂ— we also require that the external field decay sufficiently rapidly as time reaches infinity. For manipulations it does not matter what p is but since eq the system starts from thermal equilibrium we take Peq = P ^^'^ P If we assume a time dependance to the density matrix of the form p(t) = pgq + 6p(t) then the solution to the time evolution equation to first order in 6p and linear in the external field is 5p(t) = ^ / dx e-^"^t -^^/^ [6H(x). p ] e-^^'H^t _oo eq (This is easily verified by differentiation remembering that time appears both in the limits and in the integrand.) From this equation the induced change in the average density

PAGE 16

-116 = tr p(t) n(r) tr p n(r) is given by 5 = ; df J dr' <1 [n(r t), n(r' t')]> 6V(r' f ) _00 Q t = zi ; dt' ; df' x"(ft; f' t') 6V(r' t') (2. A) where [A, B] = AB BA is the commutator bracket. In the above equation are the Heisenberg operator for the unperturbed system n(?' t) = e^"*/^ nCr) e"^"*/^ and we have introduced the dynamic densitv resp onse function (please note the double prime superscript) defined by x"(ft; r'V) [n(ft), n(f't')]>gq This equation is the fundamental result of linear response theory. It shows that the density response to a small external potential is the averaged commutator rather than the correlation function S(r.t) as one might expect. However we shall see that the two are intimately related. The Dvnamic Densitv Response Fiinrtinn In equilibrium the system is space translational ly invariant. We adopt the convention that forward spatial Fourier transforms have a minus one signature and write

PAGE 17

-12x"(k) = J d(? -?) e"^"^*^"^ "'''^ x"(r V < k f^^t)' n!(t')]> k k Here we have used the definition of the spatial transform of the density operator n^(t) = ; dr e"^"^*^ n(?t) k which can be expressed in terms of plane wave creation and annihilation operators as njo) = i; aj ^a^ ^ from which it follows that n ^ = n^ -K K Time translation invariance of the dynamic density response function is easily established from the definition of a Heisenberg operator and the cyclic nature of the trace. We adopt the notation that forward temporal Fourier transforms have a plus one signature: 00 x"(Kto) = ; d(t t') e^"^* ^'^ x"(K. t t') 00 It is easy to show that this space and time Fourier transform is real (since the function is a commutator of hermitian operators), an odd function of frequency (since the equilibrium state is invariant under time reversal and parity), and depends only on the magnitude of k (spatially isotropic). It can also be shown that x"(Kco) > 0

PAGE 18

-13QC in a stable system. This is just a statement that a dissipative many body system takes more energy out of the external field than it gives back. In addition to the density-response function it is convenient to introduce certain interrelated functions. The complex response function (please note the tilda) is defined in transform space as ^(Kz) = f ^ ^^^^ ^ ir CO z -.00 This is an analytic function of the complex variable z as long as imaginary z does not equal zero (on the real axis it has a branch cut). For z above the real axis it reduces to the temporal Laplace transform of the density response: 00 x(Kz) = zi ; dt e^^* x"(Kt) 0 while below the real axis 0 = (zi) J dt e^^^ x"(Kt) This follows from the identity M i"t dco e Ziri co-z~ 0 t 0 J 7 = 4. ^ Inn z > 0 Â—00 0 t > 0 T ,'-,4. Inn z < 0 e"^^^ t < 0 which can be proven using the Cauchy integral formula in the complex omega plane. (For t>0 we can close the contour on top-exp(iwt) being bounded-and the contour is in the positive sense. For t<0 we must close the contour on the bottom, the contour now being in the

PAGE 19

-14negative sense. We then just sum the residues, the pole being included depending on the sign of imaginary z.) The physical response may be defined in transform space as the limit of the complex response as we approach the real frequency axis from above 00 x(Kco) = lim x(K. 2 = CO + ie) = limf Â— /"i'^"'^ e-o e-o -co ^ + Using the identity (P denoting Cauchys principal value) = P 1 iir 5(x) the physical response can be expressed in terms of its real and imaginary parts as x(ku) = x'(Kco) + ix"(Ku) where the imaginary part is just the previously defined density response and the real part (note the single prime) is x'(K(o) = PJ dcoLx;^ _^ ir u 0) an even function of frequency, is a manifestation of the KramersKronig relations. The utility of the physical response function is that it links the induced change in the time-space Fourier component of the density with the time-space Fourier component of an arbitrary external potential 5 = x(Ka)) 6(Kco) #

PAGE 20

-15One may also note that this same relation may be easily obtained via eq. 2. A using the definition of the physical response in real (untransformed) space, namely x(rt) = zi e(t) x"(r,t) The response to a static external potential may be derived from linear response theory by assuming that the potential is adiabatical ly switched on starting from the remote past 6V(rt) = ^<'^> I I 0 t > 0 At t = 0 the external potential is at full strength and we find 00 6 = zi J dx ; dr' x"(r r', x) e"" 6n(r') 0 or 5 = x(K) 6V(K) where the static response is given by zero frequency limit of the physical response X(K) = x(K, CO = 0) = lim x(K, z = ie) e-O A note on this last equality. Formally we have lim x(K. z = ie) = lim f ^ = (P f ^co iHlMi^ ^ ^ ^ ^ e^O e"*o -00 but because x" is a real odd function of frequency the imaginary term vanishes and we can drop the "principal part" condition.

PAGE 21

-16Relation to Density Correlation Functions At this point we should explicitly confirm the statement that the density response is intimately related to the density correlation function, the latter being the experimentally pertinent quantity. Now because the Fourier transforms of correlation spectra are related by 83 the equation r _ e"^"* dt = e"^^ f e"^"* dt -00 _oo it follows directly that J <|:^ [A. B(t)]> e-^"* dt = j:^ (1 e"!^^) J e"^"^ dt and so in particular (using the fact that x"(K(o) is odd in frequency) x"(Kco) = ^ (1 e""^^) SCKco) where S(k,w) is the time-space Fourier transform of the time dependent density correlation function. We note that S(r,t) is Fourier transformable because when r and/or t are very large, then n(r,t) and n(0,0) are statistically independent Â„ = Â„ *-M "q eq and so S(r,t) as defined as a correlation of excess density vanishes. The above relation between the density response and density correlation functions is an explicit example of what is called the fluctuation dissipation theorem^^ which relates two physically distinct quantities of fundamental experimental significance: 1)

PAGE 22

-17Spontaneous fluctuations, which arise even in the absence of external forces from the internal motion of the constituent particles. Described by S(k,w) these fluctuations give rise to scattering of neutrons or light. 2) Dissipative behaviorÂ— all or part of the work done by external stirring forcesÂ—is irreversibly disseminated into the infinitely many degrees of freedom of thermal systems. The fluctuation-dissipation theorem shows that S(k,w) is not quite symmetric in w, being a little stronger at positive frequencies than at negative. Indeed, since tox"(Ku) is always even in w, it i s general ly true that S(K, -(o) = e"''^'^ S(K, co) This result makes sense from a neutron scattering point of view. Positive frequency means the neutron has lost energy to the system (by creating an excitation of energy hw) while negative frequency describes a process in which the neutron has picked up energy from the system (by destroying an excitation). Of course, to destroy an excitation you must first have one, and their relative abundance is given by exp[-bhw]. This dissymmetry of the scattering intensity, proportional to S(k,w) is only pronounced at low temperatures, kt
PAGE 23

-1800 (o^ E ; X"(K{o)} do) _00 The j = 0 or zeroth moment follows from analyticity of x"(k,w). f ^ = X(K. 0) = X(K) _00 All odd order moments vanish because x" is odd in frequency. The remaining higher order moments can be derived from the expression (^)'' x"(rt; ?'t') = <^ [(i 1^)"^ n(r. t) n(r't')]> taken at equal times t = t'; this means that r x"(Kco) = J d(? f') e^^*^'' '''^ _00 <(j;)^''^ [[...[n(rt), H] .... H] n(r't)]> The right hand side contains a sequence of equal time commutators which can in principle, and sometimes in fact, be exactly calculated. For example the second order moment yields the Thomas-Rei che-Kuhn or F-sum rule equivalent 00 2 r do) Â„ii/L^ \ e K J ^ wx (Kco) = Â— Â— -CD ^ "> For the electron gas the fourth moment has also been calculated. ^^"^^ The sum rules are often used to provide coefficients for an expansion of x(k,z) for large z from the definition vri^y'* r" x"(K(o) (/ -Â•'^Trco-z -^ 2 4~ + Â— ^ _oo 2 z

PAGE 24

-19From its derivation, which expands (1 w/z)~ = 1 + (w/z) + (w/z) + ... it is clear this expansion can only be asymptotic, i.e. valid only when Izl is large compared to all frequencies in the system, which means all frequencies for which x"(k.,w) is not effectively zero. An interesting feature of the sum rules is their existence. There is no reason why the thermodynamic average of the multiple commutator should exist to all orders, and indeed although the sixth 90 moment has been calculated in the classical fluid limit it may 91 diverge for the degenerate electron gas in certain regions. As the sixth moment determines the w~^ term in an inverse frequency expansion, its divergence would herald the presence of a term of the 4 -11/2 form k w in S(k,w) at long wavelength for the high frequency behavior of the spectrum. Great caution is thus indicated in the use of moment sum rules beyond the first three mentioned. Distribution Functions It is natural to describe the structure and thermodynamic properties of liquids in equilibrium by employing static distribution functions. Most treatises on fluids introduce such quantities first and later generalize to include time dependance as a prelude to exploring non-equilibrium statistical mechanics. Here we have taken a reverse path; we have already related time dependent density correlation and density response functions and it is our intent now to relate these with the static distribution functions common to liquid theory.

PAGE 25

-20The most basic static distribution function is the single particle distribution function. If we define the number density operator in N particle configuration space as n n(r) = I 6(f r.) i = l ^ (the r. are quantum operators) then the single particle distribution is given by p(r) = = N ^ J" W(rp .. r^^) dr2 .. df^^ where the configurational part of the partition function is 0 = J W(r^ rj^) df ^ drj^ and the configuration probability WCr^ ... r^^) is given by the diagonal slater sum W(r^, ... r^Ir^, ... r^^) where, in general, the off diagonal Slater sum is defined as N(r^ r^ I r| ...r') = N X^'^ Z ?t(r^ .. f^) e" P
PAGE 26

-21In the classical limit the probability function "W" is just the Boltzmann factor of the potential energy of the configuration; at zero temperature it consists solely of the ground state wavefunction. The information contained in the single particle density distribution is minimal, for in the case of a translational ly invariant Hamiltonian it is a constant equal to the average number density of the system. However we introduce it here because relevant thermodynamic properties we will consider require the information of the off-diagonal single particle distribution function^"^"^^ p(r I r') = N 1 J W(r^ r^... r,^ I r j r^ r,^) df^ ... dr^ The physical significance of the off-diagonal distribution function is that its spactial Fourier transform provides us with the momentum distribution of the interacting system. The pair distribution function is defined by g(r ?') = N(N 1) 1 J N d?. .. r. = and for a translational ly invariant system is of the form g(r-r'). The pair distribution function in turn is related to the static structure factor through the relation S(K) = 1 + p J e^^*'' (g(r) 1) dr Within a factor of density the static structure factor is simply the time independent density correlation function. S(r -?')=! <(n(r)

PAGE 27

-22Note that if we had defined the static structure factor in terms of the density operator instead of the excess density we would replace g(r) 1 by g(r) in the formula relating the two. This amounts to including the experimentally unmeasurable "forward scattering." The definition of the static structure factor is immediately generalized to form the dynamic structure factor, essentially the previously introduced time dependent density correlation function (again within a constant factor of density) by replacing the quantum operators r. with their time dependent equivalents in the Heisenberg representation. Bearing in mind that r.(t) and rj(t) are operators which in general do not commute (except at t = 0), we easily obtain several important properties of the dynamic structure factor for translational ly invariant systems: 1) S(-r, t) = S*(r,t) where S* is the complex conjugate of S. 2) In the classical limit the imaginary part of S(r, t) vanishes. 3) S(K(o) E^f dt e-^"^ ; dr e"^''*^ S(rt) _00 is real in both the classical and quantal cases. 4) The elastic sum rule 00 S(K) = J S(Ka)) dto Â—CO is an immediate consequence of the definitions. All the relationships amongst the various functions discussed are best conveyed diagrammatical ly (See Fig. 2.1). Note that in the upper

PAGE 28

-23left corner we have introduced a new function C(k,w) called the Kubo or Relaxation function, about which we will have much to say later. For now it simply expresses the difference between the (dynamic) physical response and the static response function. The point of our review of the response, correlation, and distribution functions, so concisely summarized in Fig. 2.1, is that if we wish to generalize the classical many-body methods of calculating g(r) into the quantal regime, there are two possible relations open to us. First there is s''^(k) as the classical limit of $^""(1^) (obviously) or. as P 6n(r I U) 6 pU(f') = ^ X(K) where the Iq denotes evaluation at constant density/zero external potential, and the F denotes the Fourier transform with respect to the spactial variable r-r' Let us contrast the two alternatives. The quantum mechanical static structure factor S^'"(k) is essentially statistical in nature; it is strongly dependent on the dynamics of the particles, depending on S(k,w) ~ Imag x(Kw) over the full range of frequency. On the other hand x(K) depends on the single zero frequency limit of x(Ku), and is essentially mechanical, arising from the Schroedinger equation by our linear response derivation, which makes but implicit reference to statistical questions. This is clearly indicated by the F-sum rule-which holds in any stationary

PAGE 29

-24xCkco) X(kco) = x(K) + itoPCCku) yiKta) E Im x(Kco) ir (0 lim (0 0 SCkco) = ^ (1 + cot h ^) Y(kco) x(K) = ; y(Kw) d(o X(K) S(K) = J S(Kto) dco Figure 2.1 Relations between the linear density-density response function and the structure factor.

PAGE 30

-25ensembleÂ— expressing nothing but the conservation of particles and the elementary commutator of position and momentum. This dissertation will concern itself with the integral equation methods of solving the many body problem. This route seems particularly promising for treating quantal fluids in that for the classical limit there exists a wide variety of integral equations for g(r) which have been investigated and give results in good agreement With experimental results. The purpose here is to extend this procedure into the quantal regime. A fundamental ansatz of this treatise is that one would expect the more accurate integral equation approximation to correspond with the less sensitive function. The extension of classical integral equations into the quantal regime via the density response function is presented in Chapter II. With this fundamental ansatz, however, there is a fundamental drawback: From x(K) one obtains an integral equation not for g^"^(r) (which is used in the computation of thermodynamic quantities) but rather n(r/v), that is the distribution function in the presence of an external potential which has the form of the interparticle interaction. While it is true that in the classical limit n(r/v) = n^qir) this is not true in the quantal case, since we are neglecting exchange and recoil effects. The external potential can be thought of as arising from a distinguishable particle of infinite mass and therefore fixed in position. An extreme example occurs for the case of very weak interactions between fermions. In such a case n(r/v)/nQ is approximately unity everywhere, whereas g(r) has approximately the ideal fermion value of 1/2 at the origin.

PAGE 31

-26Thus our task is twofold: first to generate an accurate quantal integral equation for n(r/v)/nQ and second to generate an approximate g(r) from knowledge of n(r/v)/n^. In essence our second task is to go from the bottom left of Fig. 2.1 to the bottom right along the paths drawn in this figure. We cannot go along the diagonal, as this path is one directional (we lose information when we integrate y(k,vi) over w in order to obtain x(k)). We must first traverse over the upper left path of Fig. 2.1; we must obtain an accurate approximation to the Relaxation function. This is taken up in subsequent chapters.

PAGE 32

CHAPTER III QUANTAL FUNCTIONAL EXPANSIONS In this chapter we will establish integral equations for the quantity n(r\y)/n^ (the density distribution in the presence of an external potential of the form of the interparticle potential) through the functional expansion method. ^'^"^^ Because in the classical limit these equations reduce to the familiar Hypernetted-Chai n (HNC) or Percus-Yevick (PY) equations for the pair distribution function, we will review the classical derivation with the aim of establishing corresponding quantal generalizations. The OOZ Relation In classical fluids the Ornstein-Zernicke (OZ) relation plays an important role in treating integral equations for g(r). Thus in this section we derive a quantal extension of the OZ relation (QOZ). We start from the known functional dependance of the classical canonical partition function on pairwise additive potentials. This allows us to write the functional derivative of the density distribution as^^ = n^ 6(? ?) + n^^ ^g(,^ r'l) 1} = F"^ {n^ S(K)} 0 6n(r I U) 6 3U(r') -27-

PAGE 33

-28The denotes evaluation when the arbitrary static external potential U(r) is turned off, that is at n(r) = n (n ; the 0 0 average density). The (negative power) script F denotes the (inverse) Fourier transform on the spatial variable r r', and S(k) is the previously defined static structure factor S(K) = 1 + n^ ; dr e"^'^*'^ {g(r) 1} Now consider the functional S BU(r') 6 N(f") This is a vacuous quantity unless such a functional of density can actually be constructed.^^ If such a creature does not exist, we define it through the chain rule of functional calculus: 6 n(r") Sn(r") 6 3U(r') = 6(r r') It follows that the inverse functional derivative is 6 BU(r') 6 n(r") 1 K ^n S(kT> = r Cdr' ?"l) 0 0 where we have split up the expression for the derivative into two pieces (this should be considered an arbitrary process) with the first piece chosen to be the non-interacting (g(r) = 1) result and the second piece (c(r)), called the direct correlation function, representing the deficit from from the non-interacting result. The functional chain rule as expressed in terms of c(r) takes the form

PAGE 34

-29h(r) = c(r) + J c(lr r' I) h(r') dr' This is the classical Ornstein-Zernicke (OZ) relation, where h(r) = g(r) 1 is called the pair correlation function. Furthermore, in the classical limit where there are no exchange or recoil effects (the momentum integrations factor from the configurational piece of the partition function) one can rigorously derive the bootstrap relation^^ n(r I U = v) = n^ g(r) where the external potential U(r) has the form of the interparticle interaction V(r). This bootstrap relation allows us to reformulate the OZ relation in terms of the inhomogeneous density distribution as fol lows: 6 n(r') n(r' I U = v) n^ = c(r) (3.1) 0 A quantal extension of the OZ relation is straightforward;^^ we will follow many of the above steps. First, from linear response theory, we saw that in the quantal case Sn(r I U) 6 3U(r) ^ = I^K^ {% X(K)} where x(k) differs by constant factors from the static response function previously defined in Chapter II. Defined here it is conveniently dimensionl ess ; note that only the product (5x(k) remains finite as temperature goes to zero.

PAGE 35

-30Again the chain rule of functional differentiation requires that 6 BU(r') 6n(r" I U) = ^K^ ^=k^K ^J(K)> ^"'> (3.2) but here we have decomposed the inverse derivative into two new (and as yet unspecified) functions. These are determined as follows. As we wish c(r) to reduce to the classical direct correlation function in the proper limit, we need only require J(k) = 1 in the classical limit. We furthermore make the ansatz that J ^ 6n(r') n(r' I U = v) n^ = c(r) (3.3) also holds for the quantal case. We have seen that this equation is strictly true in the classical case; by combining the two above equations we constrain the function c(r), and so simultaneously J(k), by the relation F,{n(r I U V) n^) = xili ^^ ^^ In hindsight we now see the motivation of requiring eq. 3.3 be fulfilled: as J(k) tends to unity in the classical limit we have merely imposed the condition that we recover the bootstrap relation in the classical limit. The above equation also determines J(k) (and so c(r)) in the general quantum limit, for if the system described by n(r/v) and x(k) were that of non-interacting fermions, then in the limit that U(r) vanishes we see that J(k) must equal This is an appealing result, for it amounts to defining the quantal direct

PAGE 36

-31correlation function by separating out the non-interacting result, just as we did in the classical case. Henceforth we will implicitly use this result, which will be referred to as the quantal bootstrap relation (QBS). As a consequence of this now specified decomposition of the inverse functional derivative eq. 3.2, the functional derivative chain rule expressed in terms of the quantal direct correlation function becomes ^ 1) = M{c(?)} + n ; dr'(M{c(lr r'l)}) (^i^^^ i) 0 where M{...} is an operator defined by F|^{Mf(r)} = Xq(K) F^{f(?)} This is the quantal generalization of the OZ relation (denoted QOZ for short). Note that this equation can be cast into the classical form by defining analogue functions: gan(-) = nirj_vi ^an^-^ ^ ^an^-^ ^ 0 c^"(r) E M{c(r)} The classical equivalents of these functions all have "physical" interpretations that are only loosely identified with the analogue. (For example g^" should not be confused with the true or quantum mechanical pair distribution function.) To distinguish classical analogues we shall henceforth employ the superscript "an."

PAGE 37

-32Classical Generating Functionals The heart of the functional expansion method of deriving integral equations lies in the ability to construct a functional W[n] of the density distribution n(rlu) in a nonuniform system which is reasonably insensitive to the function n(r|u).^^ (For other points of view, see Bakshi^ and ChoquardJ^^) If such is the case we can approximate the functional W[n] by making a functional "Taylor series" expansion, for example about the average density n^, and truncate after the first (linear) term W[n] = WCn^] + ; dr'( 6w 6n(r') ) (n(r' I U) n^) The OZ relation then forms a closed set of two equations in two unknown functions. As an illustration we will consider a set of functionals and the integrals equations they generate in classical statistical mechanics. First consider the functional W[n] = In U> e-PU(r) The denominator forms what is called the reference distribution. In order to make the functional variational ly insensitive, it is constructed to approximate the numerator as close as possible. Here we have taken the classical non-interacting Boltzmann factor. The first functional derivative is easily found to be

PAGE 38

-336w ^ 6(r r') 6 SU(r) 6n(r') n(r I U) 6 n(r') and from the definition of the direct correlation function: 6 m?) 6n(r') 6(r r') = ^ c(r r' ) 0 we find Sw 6n(r') = c(r r' ) 0 The functional expansion of W[n] as the potential U(r) is changed from the intitial value 0 to its final value In "^^ = ; d?'(n(?' I V) nj (6w ) along with the bootstrap relation n(r/v)/nQ = g(r) and the OZ relation then combines to give the HNC equation g(r) = ee^^"^^ where yd) = h(r) c(r) is called the nodal function. The PY equation g(?) = ePV^^^ [1 Y(-r)] follows from the same steps by starting with the functional n^ eBecause the HNC functional is the logarithm of the PY functional, it is the more insensitive functional of the two from variations. One might be tempted to conclude that it would therefore produce the

PAGE 39

-34superior approximate integral equation. However the validity of the HNC equation (or any other derived by a functional expansion method) depends on whether the cumulative effects of the higher order derivatives in the functional expansion are negligible. In actuality their effect is to push the form of the equation towards that of the PY (note how the form of the PY equation is the small y(r) expansion of the exponential y(r) term of the HNC). This is evidenced by the fact that for purely repulsive potentials the PY and HNC equations of state bracket the "exact" Monte-Carlo or molecular dynamics computer 1 02 simulation results. What we do know from diagrammatic 52Â—58 analysis is that the exact closure equation to the OZ relation for g(r) is of the form g(r) = e" '^^^''^ Y^r) + B(r) where the unknown "bridge function" B(r) is the sum of bridge diagrams of two point functions. The evaluation of the bridge diagrams in terms of higher functional deri vati ves^^*^^ or other attempts to account for bridge effects will be discussed in Chapter V. Here our intent is to generate zeroth order integral equations from variationally insensitive functional, thereby minimizing the effects of the higher order terms that we neglect in the Taylor functional expansion. It is well known that the PY equation has no solution for a one component plasma system. Numerically this arises from the fact that in the PY equation for g(r) the long range tail of the potential is not compensated by that of the nodal function yd), as it is in the HNC equation. This is symptomatic of a physically meaningful

PAGE 40

-35inadequacy, namely; in real systems the particles do not feel the "bare" long ranged potential but a collectively screened or short ranged potential An improvement on the PY generating functional would be to replace the bare potential appearing in the Boltzmann factor of the reference distribution by an appropriately screened one. For example W[n] = "^'^ l5UÂ„(r) "o ^ where the Hartree potential is defined as 3U^(r) = 3U(r) + n^ / dr' (3U(r ?') (g(r') 1) As the Boltzmann factor now more closely describes the physical distribution n(r/v) this functional is variational ly less sensitive than the usual PY _6w_ ^ Sir ?) n(r I U) ^ ^^h^"^^ 6n(?') 3U.(r) 3Um(?) 6n(r') Using the bootstrap relation in the definition of the Hartree potential we see that LlJMIl ^ s m?) ^ ^ ^5 ( S m? Sn(r') 6n(r') 6n(r') (n(S) n^) ; dS fJU(? 5) 6(?' s) which when evaluated in the limit v = 0(n(r|v) = n^) and using the definition of the direct correlation function (eq. 3.4) yields

PAGE 41

-366 3U^(r) = ^^'^n "^'^ c(r r') 3U(r ?') 0 6n(r') and so = c(r ?') + 3U(r ?') 6n(r') Q The linear expansion of the functional W[n] thus generates the equation "^^ = 1 + ; d?' (n(?' I U) n^) (c(? ?') + mr ?')) which can be re-expressed with the help of the bootstrap equation and the OZ relation as BV.(?) g(r) = e {1 + pV^C?) + y(?) pvCr)} Note that in this equation the long range part of the potential is cancelled by the long range part of the nodal function, and all other functions appearing in this equation are screened. This equation, which we shall term the PYH equation (for Percus-Yevi ckHartree) is equally applicable to long and short range potentials, unlike the PY equation. If we try to improve the HNC equation by introducing the screened Hartree potential into its generating functional W[n] = In U) 3U^(r) "o

PAGE 42

-37we find that its first variation is the same as in the above PYH derivation and (following many of the same steps) the functional expansion leads to the equation namely the original HNC equation. In other words the HNC equation has the striking property of being invariant to (Hartree) screening of the reference distribution in its generating functional. From a related point of view it has also been shown^^ that the HNC equation is the limit in a series of integral equations that can be systematically generated starting from the PY equation, among others. Finally we should note that the HNC equation can be derived from an even more interesting point of view. We ignore higher order functional derivatives and assert that the best functional for our closure purpose are insensitive to variations, and this property shows up in the magnitude of its derivative. This can be minimized by using an even more realistic screened potential for the reference distribution. In fact it is easy to show that if we modify the Hartree potential by replacing the convoluted bare potential by the direct correlation function W^ir) = (3U(r) / dr' c(r ?') (g(?') i) then the derivative of both the HNC and PY like functional H = In -^U) ,, n(? I IJ) g(r) = e BV^(r) Y(r) + 3V^(r) 3V(r) n 3Ujr) n 0 e 3U (?)

PAGE 43

-38actually vanishes. They in turn both yield the same equation 3V ^r) g(r) = e ^ which one will now notice, via the OZ relation, is simply the HNC equation. Quantal Generating Functionals Corresponding quantal equations can be obtained from the above classical generating functionals if we replace the reference distributionsÂ— the Boltzmann factor n^expf-bU} describing classical non-interacting particles in the presence of an external potential U(r)Â— by its quantal extension n*(rlu)Â— the density distribution of non-interacting fermions in an external potential. (Henceforth non-interacting fermion system quantities will be distinguished by a star.) That is we replace the reference distribution by a single particle Hamiltonian approximation. As in the classical presentation of the above section the functionals are then linearly expanded and evaluated in the limit where the external potential is set equal to the interparticle interaction potential. We have previously derived a quantal extension of the pair correlation function V"o ^("^^> -rk)-^ from the ansatz that the exact classical relation J ( 6 SU(r) 6n(r' ) ) (n(r' I U = V) n^) dr' = c(r)

PAGE 44

-39remains valid in the quantal case, alternatively expressed by the ansatz that the classical bootstrap relation is extended as By making use of these relations and defining the analogue potential g3U^"(?) ^ n*(? I U) "o a closed set of equations is formed with the QOZ relation. Expressed in terms of previously introduced analogue functions, these integral equations bear a close correspondence with their classical counterparts, and indeed reduce to those counterparts in the appropriate limit. As a first example we consider the functional W[n] = In "^"^ I n*(r I U) Its first variation is _-5w_ n*(r I U) r 6(? r') n(r I U) Sn*(r I U) n(r I U) ^n*(r I U) ^^.^^ ^^^2 6n(r") > This can be simplified in the following manner: 6n(r') 5 3U(?") ^ 6n(r') ^'''n^ xCQ) ) = 1 FJc^"} where it should be remembered that the analogue pair correlation function differs from the quantal ly extended "physical" pair correlation through the operator

PAGE 45

-40-1 M{f(r)} = F-'{x*(0) F{f(r)}} From the above considerations we see that 6n*(r I U) 6n(r') = 6(r r') n^c^^Cr r') and so 6w 5n(r') = c'"(? r') By Taylor expanding the functional W[n] about U(r) = 0 to v(r), we find In "(r I V) p n*(r I V) = J ^"^^ I "o^ ^ 6w 6n(r') or This equation we will call the Quantal Hypernetted Chain or QHNC equation. Following similar steps with the functional u n(r I U) ^ n*(r I U) we can arrive at or what shall be referred to as the QPY equation. It is of no trivial importance that the QOZ, QPY and QHNC equations are of the exact same form as their classical counterparts (this in addition to the property of reducing to the classical equation in the limit). This is because phenomenologi cal ways of incorporating the effects of truncating classical functional

PAGE 46

-41expansions at the linear term may still be of use. [In the classical HNC case such a truncation is equivalent to neglecting bridge diagrams in the Mayer expansion of non-ideal gases.] This will be considered further in Chapter V. It should also be stressed again that these integral equations, even if they could be made exact, are for g^"(r) = n(r/v)/nQ and not the pair distribution function. This is evidenced by the fact that if we turn off interparticle interactions, V^"(r) vanishes and the solutions of either the QHNC or the QPY equation yield g^"(r) = 1. This is correct for the density distribution but the fermion pair distribution has an exchange hold. This difference is taken up in the next chapter. The exact correspondence between quantal analogue and classical integral equations is lost when one considers screening effects on long range potentials. Rigorously incorporating screening effects is complicated in the quantal case by symmetry effects. It can be shown^^'^ that the best single particle Hamiltonian (best in the sense that the statistical average over the square of the difference of the exact and single particle Hamiltonian is minimized) is provided by a Hartree-Fock potential, not the Hartree. On the other hand we will see in the next section that as far as density distributions are concerned there exists an exact effective potential single particle Hamiltonian. In the quantal integral equation approach we will consider here screening of the reference potential will ignore exchange effects; we will see in the next section that they are actually incorporated in the closure of the linear functional expansion, namely the QOZ relation.

PAGE 47

-42When we screen the reference distribution of the QPY generating functional ^ ^ n(? I U) n*(r I U^) with the Hartree potential 3U^(r) = (5U(r) + J dr'(pU(r' r) (n(r' I U) n^) its linear expansion requires the variation 6w 6n(r') 6(? ?) 1 V "o "o 6n(r' I U) This can be evaluated by invoking the chain rule over the Hartree potential, which in Fourier space becomes 6n*(r I U^) 6n(r' I U) 6n*(r I U^) 6 pU^(r") S 3U^(r") 6n(r' I U) = "oX*(0) Fq{ 6 3U^(r") 6n(r' I U) Now the variation of the Hartree potential yields 6 3U^(r) 6n(r') 6 BU(r) 6n(r') |iU(r r') or 6 3U^(r) Sn(r') } = n^x(Q) 'a K{mr ?')} from which a few simple manipulations yield

PAGE 48

-436w = c^"(r-r') + Fq^ x*(0) F{(iU(r-r')} = c^"(r-r') + M{|3U(r-r')} 0 ^ 5n(r') Plugging into the functional expansion of W[n] we obtain "^"^ = 1 + ; d?'(c^"(? r') (n(?' I v) n ) n*(? I V^) 0 + ; dr'(M{flU(r ?')}) (n(r' I v) n^) In terms of analogue functions and employing the QOZ relation this becomes g (r) = e [1 + Y (r) + Fq {x*(0) F{n(? I v) nj F{(JV(?)}} = e [1 + Y (r) + M{; dr'(pV(r ?')) (n(r' I V) n^)}] = e [1 + Y (r) + M{3V^(r) 3V(?)}] This equation will be referred to as the QPYH equation. Unlike the classical case, if we screen the QHNC generating functional, we do not recover the QHNC equation. A now familiar procedure yields instead the result g^"(?) e' '^""^''^ ^^""^ ^ M{3V^(r) pvc?)} to be known as the QHNCH equation. We see that for both the QPYH and QHNCH equations long range tails cancel but the potential enters in a nontrivial manner. Lastly, in the classical case we saw how the HNC equation could also be derived by constructing a screened potential such that the first variation of the generating functional vanished. The same

PAGE 49

-44procedure can be applied to the quantal case. Define an effective potential (5U.(r) = (5U(r) J dr'(c(r r')) (n(r' I U) n^) I. 0 (where c(r) is the quantal pair correlation function and not the analogue quantity); then both the QHNCand QPY-like generating functional s W = In U> H^_ni^_LUi_ n*(r I U^) n*(r I U^) have vanishing first variations and both yield the same equation: an This equation is distinct from the QHNC equation but like the QHNC equation reduces to the classical HNC equation in the limit. For this reason it shall be termed the 0HNC2 equation. Interestingly enough the 0HNC2 equation can be obtained from an entirely different approach. Instead of assuming the quantal bootstrap relation and utilizing Percus' method we can arrive at the QHNC2 result directly from Kohn-Sham-Mermi n density functional theory. From there we can recover the quantal bootstrap relation (our starting ansatz) and so justify our variant equations by making what we will see as a reasonable assumption concerning the effective potentials of many body systems.

PAGE 50

Â•1 -45The Density Functional Approach Following the zero temperature approach of Hohenberg and Kohn^*^^ and Kohn and Sham^^^ for the ground state properties of an electron gas in an external potential U(r), Mermin showed that^^ i) the Helmholtz free energy A of an i nhomogeneous electron gas in thermal equilibrium is a unique functional A[n] of the one electron density n(r), and ii) the grand potential Q[n] = A[n] vi J" n(r) dr is minimum for the equilibrium density n(r). We employ the grand canonical rather than a canonical ensemble because macroscopic quantities derived from a non-interacting Fermi system are more easily calculated and the fundamental variable, namely the density distribution, is a physical observable that should be the same from either ensemble. The free energy is customarily decomposed as A[n] = 1 J v(r ?') n(r) n(r') dr dr' + J U(r) n(r) dr + A*[n] + A [n] X c to separate out respectively 1) solely classical electron-electron interaction contributions, denoted by v(r), 2) the external potential contribution, 3) the contribution A*[n] that a non-interacting Fermi system with density n(r) would still exhibit, and 4) all other effects are grouped to define the unkown exchange and correlation functional A^^[n]. The minimum property of the grand potential leads to the functional equation 6A*[n] Â— + U f^(r) }i = 0 (3.5) 6n(r)

PAGE 51

-46where u is the chemical potential and U^^^d) is defined by SA rn] Ug^^(r) = U(r) + J dr V(r ?') n(r) + 6n(r) Equation 3.5 is equivalent to the equation for non-interacting electrons in this effective potential; that is we have n(rlU) = n*(rlUg^^). The non-interacting fermion distribution n*(rlU) is given by the following system: (^ A + Ug^^) = e. fi = M = 1 units n*(r I Ug^^) = I f(e.) I ^.(r) 1^ where the sum over the index "i" refers to all bound and continuum eigenstates and f(e) is the Fermi-statistics occupation number 1-1 f( e) = {1 + exp[p(e. li)]}' dependent on the temperature and the chemical potential. The latter is determined as usual from the total number of electrons by spatially integrating over n*(r). For the jellium model of the electron gas the external potential must not only include the source of the external potential but the effect of the neutralizing background as well U(r) = U^Cr) J dr' n^ V(r ?') so that the effective potential is Ueff(r) = U^Cr) + ; dr' V(r r') (n(? I U^) n^ +

PAGE 52

-47At large radii, the displaced electron charge vanishes, and any long range tail of the source potential is cancelled on the assumption of perfect screening, so that the assymptotic value of the effective potential is 6A W^^">= 6n xc Applying eq. 5 at large radius gives 6A* 6A. 6n. = ^ xc 6n or, equivalently v2 Â„3/2 nÂ„ = Â— 3 ^ I^/^Pii) 0 2 ir where the standard definition of the Fermi-Dirac functions has been used: I (n) E J y* 0 1. e^-^ It is convenient to shift the zero point of our energy scale to define a new effective potential which goes to zero at infinity Ugff(r) = Uq(?) + J d?' V(? r) (n(r I U ) n ) 6A + (xc 6n 6A. ) (xc 6n ) Up to this point our equations are exact; the exchange correlation part of the free energy functional is, however, quite unknown. But if we make a functional expansion to first order

PAGE 53

-486A. xc 6n 6A ) = (xc 5n ) + ; K(r r') (n(r' I U^) n^) dr' 5 A K(r r') E (xc 6n(r) 5n(r') we obtain an approximate integral equation in terms of the unknown kernal K(lr-r'l): n(? I V^) n^ = n*(-r I U^^^) n^ ^eff = ^0 -f" dr'{V(r r') + K(r ?)} {n(?' I U^) n}^ Since the unknown kernal is independent of the external potential U(r), we can determine this function by considering the case of a very weak external potential where the linear response formula can be employed on both sides of the first of the above equations with impunity. The left hand side yields 6n(0) E FQ{n(? I U^) n^} ^ n^p x(Q) U^CQ) while the right hand side gives FQ{n*(r I Ug^^) n^} : n^3x*(Q) U^^^CQ) = n^^x*(Q) {Uq(0) + [V(0) + K(0)] 8n(Q)} Using the last equation we can show that the Fourier transform of the unknown function K(r) is K(0) = (V x(0) x*(0) ) V(Q)

PAGE 54

-49Due to this assignment the effective potential reduces in form to the previously defined screened potential Uj.(?) = Uq(?) ; d?' (1 c(? -?') (n(r I U^) n^) and we obtain the relation n(r I U^) = n*(r I U^) Note that the source external potential is completely arbitrary, and so this result is a generalization of the QHNC2 equation, which follows immediately when the external potential is taken to be the interparticle interaction. From this generalized result we can work backwards and obtain the quantal bootstrap relation, and so justify the other (QPY, QHNC, etc) equations we have obtained. This is done by now making the assumption that the average one body potential felt by a constituent electron due to a fixed test charge 6e, namely Vt =Vtest ^ d-r'(^ c(-r r')) (nCr' I V^^^^) n^) where is equal to the average one body potential felt by a free test charge Se, arising from a fixed electron in an electron gas, namely V^-^g = (5e) l^(r) W(r) = 4Tr{ e8(r) + e[n(r I U = e^/r n ]}

PAGE 55

-50In the limit of small Se we can again invoke linear response with impunity and solve for F^{n(r I U = V) n^} = 1 thus recovering the quanta! bootstrap relation. He see that our density functional approach, aside from justifying the quanta! bootstrap ansatz of the Percus method, gives physical insight to its meaning. Even though it is not an exact approach, it clearly points out its limitation, as we approximated the variation of the exchange-correlation free energy functional by a first order functional Taylor expansion. This type of truncation is a property shared by the Percus method and ultimately may be responsible for the self-consistency of the two approaches.

PAGE 56

CHAPTER IV THE BASIC QHNC EQUATION An alternative ansatz for the extension of classical integral equations into the quantal regime lies in the construction of effective pair potentials. These incorporate the bulk of quantal effects; for example in mul ti component plasmas quantal effects prevent the collapse of particles in an attractive coulomb field. More generally the effective potential yields finite values for the pair distribution factor g(r) in the limit of zero interparticle separation r. Such classical pseudopotential s have been constructed by requesting an approximate account of the N-body Slater sum^^^~^^^ or, in an approach which can in principle be implemented exactly but is more limited in its physical range of applicability, by reproducing the two-body Slater sum^^^"^^^ (r) (JH^^^ e = (A-^'^^/my'^ (4.1) "2'^ = ^ v2 + V(r) It can be shown that this Binary Slater Sum approximation correctly represents the derivative of the classical pseudopotential for the N-body problem in the limit r tends to zero.^^^ -51-

PAGE 57

-52As shown in the previous chapter, the QHNC equation for n(rlv) is also equivalent to the classical HNC equation with an effective potential defined by P^eff^^^ n*(r I V) e = Â— Â— Â— "o where n*(r|v) is the density distribution of non-interacting fermions in the presence of an external potential of the form of the interparticle interactionÂ—the bare coulomb potential. Normally this entails the solution of the Kohn-Sham-Mermin equations That is calculating n*(r|v) is the solution of the following system:120-121 (^ v2 V) ^. = c. f. n*(r) = I f(e.) (4.2) i 111 where f(e.) is the Fermi-statistics occupation number f(e.) = {1 + exp[3(c. \i)]y^ depending on the temperature and the chemical potential fixed by the total number of particles. The index "i" runs over all continuum and bound (if existent) states. The treatment of the delocalized states proposed by Friedel for impurities at zero temperature can be applied to the nonzero temperature case too. The eigenfunctions i^^^^ of eq. 4.2 for energy

PAGE 58

-53are (|>j^^(r) is normalized by requiring that it have the asymptotic form ^^^^(.r) ^ sin[Kr e J + n^CK)] n-jCK) being the phase shifts in the potential V(r) and n^(K = >) = 0. Overall normalization is accomplished by requiring that the wavefunctions vanish at a spherical wall of very large radius R. We assume that R is so large that most of the integrand is approximated by its asymptotic form: 1 = I A^^ 1^ ; dr sin^CKr e J + n^] with the boundary conditions KR + T]^iK) = e J = mir m = 0, 1, 2, 3. .... incrementing m by unity for a given L yields the density of states (including spin multiplicity) as 9i(K) =f ^1 ^ rIkHiCK)) while the normalization factor for very large R is given by A^^ = v27r (1 -0(^)) (These same results can be derived for a long range coulombic potential .^^^ 36.23^ ^ straightforward derivation gives for the density n*(r) = I f
PAGE 59

-54where subtracting off the uniform density result accelerates convergence of the integral and summation. We note that in general a potential may support a number of bound states and we have included them al so. The OHNC Effective Potential We thus see that the QHNC pseudopotential is essentially a one-particle Schrodinger equation calculation. This is a feature in common with the Binary Slater Sum. It is important to contrast the two however. The Binary Slater sum (BSS) is used to calculate the pair distribution function g(r) while the QHNC equation yields n*(r|v); this difference is important in consideration of quantum effects. It should be noted that in the BSS approach the Schrodinger equation is solved with the use of the reduced mass. The matrix element in eq. 4.1 can be evaluated by inserting a phvsi cal basis set of properly symmetrized or anti symmetrized wavefunctions depending on i the statistics of the two particles. Parallel spin versus antiparallel spin contributions can be calculated separately. In the BSS approach the effective potentials are density independent. In the QHNC approach the Kohn-Sham system of equations can also be written in the form of a matrix element n*(r I V) = where the operator "f" is defined as f = {n + exp m p]}"^

PAGE 60

-55H is the coulomb hamiltonian operator, mu the chemical potential and we have included a variable "eta" to be set equal to unity to describe Fermi statistics. In the QHNC approach the external potential has no spin attribute, and the matrix element is evaluated using a complete basis set with no regard to statistics; the statistics are built into the functional form of the matrix element. The QHNC effective potentials are implicitly dependent on the density through the chemical potential. All these subtle differences are manifestations of the fact that in the QHNC approach g(r) must be obtained from n*(r|v) by a further calculation. One obvious advantage into this splitting of tasks is that the BSS pseudopotential breaks down at zero temperature, whereas the QHNC remains valid. A surprising feature of the QHNC effective potential is that it is analytically calculable. This allows one to bypass the computationally time consuming Kohn-Sham-Mermin system of equations. The starting point for the solution is the fact that the effective potential is related to the r^ = limit of the off diagonal density matrix for fermions in an external coulomb potential : p(f^ I r^) = Spherical symmetry reduces the functional dependance of the off diagonal density matrix does from the full six variables of r^ and 12 to only the magnitudes r^ and the angle between r^ and T^. This is easily seen if we expand the matrix element using a complete set of spherical coulomb wavefunctions :

PAGE 61

-56dL c ri^2 2m (1 + 1) ^ X S + [K ~^ ~ I 0^ S = 0 Using the addition theorem of spherical harmonics then yields where the 1-wave radial density matrix components are defined by 00 Pl(r^ I r2) = ; (spin J dK) Â— ^-^ S*^(r^) S^^^r^) ^ ^ 2m ^ ^ Furthermore, because the potential is coulombic, there is an additional symmetry: the matrix element "f" also commutes with the 1 24 Runge-Lenz vector. This reduces the functional dependance of the off diagonal density matrix to the two variables "x" and "y": 7 7 ^ 17 X = (r^ + r2)/2 y = (r^ + r2 2r^ cose)"^/2 The solution then follows the steps of Hostler^^^"^^^ and Storer^^^ except that the Bloch equation satisfied by the off-diagonal canonical density matrix is replaced by a Fermi statistics generalization [(H 1^) n + n + fj^] p(r^ I r2) = 0 (Here H acts only on the r^ coordinates and we set eta equal to unity at the end of the calculation.) Because of the aforementioned symmetries, the solution to this equation reduces to a functional form identical with the s-wave radial density matrix component (however in the variables x and y):

PAGE 62

-57The diagonal limit r] = r2 = r is given by n*(r) p(? 1 7) = ^ p (X + y) x ^ ay y) y=o The second derivatives of the radial coulomb wavefunctions can be eliminated in favor of the radial wavefunctions themselves by employing the radial coulomb differential equation. We thus obtain n*(r) = ^ J dK ( 1 2ir' J(K, r) = { 3 rfi] + 1 ) J(K. r) d? ^Ko^^^ 2m 2 K^) I S^^(r)|2} Numerical Calculations The effective potential, obtained from the ideal fermion density distribution in a coulomb potential, will be calculated from the analytic formula derived above. For comparison with the classical HNC equation, we will scale lengths in units of the interparticle separation x = r/a where n^ = u a3)-l Temperature/density points will be expressed by the classical coupling parameter

PAGE 63

-58along with a quanticity temperature (two divided by the temperature in rydbergs) 13* E fime^/fi^ In the classical limit all functions are solely dependent on the parameter r^; we expect any dependance on the quanticity temperature to indicate the pervasiveness of quantum effects. The calculations were performed corresponding to the temperature density points of tabulated results provided by Pokrant for later comparison. These consisted of values at constant density, as measured by the dimensionl ess ground state coupling parameter r^ = a/a^ s 0 where a^ is the bohr radius, as well as a dimensionl ess degeneracy temperature t = Kt/Ep where is the ground state Fermi energy. The results are presented in the following graphs. In Fig. 4.1 we present plots of the ideal fermion density in the presence of an external coulomb potential. For comparison these are plotted together with two approximations. The first is merely the classical limit n*(r I V) The second is the Thomas-Fermi approximation. It is obtained by expanding the diagonal density matrix element "f" using a complete set

PAGE 64

-59of plane waves. If we then neglect the non-commutativity of the potential and kinetic energy operators in the exponential (the neglected terms can be systematically calculated by Trotter 1 23 formulas and shown to vanish for large r) we then obtain 'j n*(r I V) ^1/2 (Pv^r/x) I "o ^1/2^^^^ where the Fermi-Dirac integral of order p is defined by the integral 0 1 + e-' (In this formula the symbol r denotes the gamma function.) It is seen that the exact result quickly merges with the Thomas-Fermi approximation. By using the relation dF (X) we find that asymptotically the Thomas-Fermi approximation has the form F ^ (3n) n*(r I V) ~ 1 I 2 ^ 1 ^^V"^ ^ 2 This is indeed the correct asymptote of the analytic solution as can be verified by looking at the appropriate expansion of the spherical 1 29 coulomb wavefunctions. This tells us that the effective potential in the QHNC approximation has a long range 1/r tailÂ—but with a strength modified from the classical theory. The QHNC tail contains the additional factor

PAGE 65

-60F ^ ((5p)/F ^ "2 "^2 not found in the Binary Slater Sum approach. In Fig. 4.2 we present plots for the QHNC effective potential. For comparison the corresponding effective potentials of the Binary Slater Sum approximation are presented in Fig. 4.3. The approach to solving the HNC system of equations with long tail potentials is well known Â— we follow the convention of Ng^^^ and renormalize the constituent functions by subtracting out an analytically transformable long range function of the form ^ erFc(ar) We employed the Ng method of solution (essentially a picard iteration with an accelerated convergence procedure) although it should be mentioned that two new algorithms have recently been published^^^"^^^ that employ a hybrid Picard/Newton-Raphson iteration procedure which promises faster convergence. The results of the integral equation solver, that is the analogue pair distribution function n(.r/yj)/r\^, are presented in Figs. 4.4-4.8. This distribution is of theoretical interest in its own right, for example in the study of impurities in metals or positron annihilation 1 32 rates. However we will limit our concern to the approximate calculation of the pair distribution function. A method of calculating g(r) based in part on the information contained in n(r|v)/n and a 0 comparison of results with the Binary Slater Sum approximation and that of Pokrant will be taken up in the next two chapters.

PAGE 66

-610. XT/a Figure 4.1 Plots of n*(r|v)/no-l versus x = r/a for a density of rg = 2 at temperatures of kt/Ef of .1, .5, 1.0, 2.0, 20.0. Storer denotes the analytic solution, following Russian literature the Thomas-Fermi approximation has been denoted Hartree, and Boltzmann is the classical approximation.

PAGE 67

tau8.10 0. -fT 1.0 1 Â— ,5 2.0 x-r/a 2.5 3.1 3.5 Figure 4.2 Plots of the QHNC effective potential versus x = r/a for a density of rg 2 at temperatures of kt/Ef of .1, .5, 1., 2., 20. We have plotted the potential in the form I Veff('^> A purely classical coulombic interaction would correspond to a constant values of unity. Storer represents the analytic solution; Hartree is the Russian literature denotation of the Thomas-Fermi Approximation, which is truncated near the origin to prevent divergences.

PAGE 68

-63Figure 4,3 Plots of the Binary Slater Sum effective potential versus x = for a density of rj = 2 at temperatures of kt/Ef of .1, .5, 1., 2., 3., 5. We have plotted the potential in the form I Veff
PAGE 69

-64static structure factor m o Figure 4.4 Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rc 2 at temperatures of kt/Ef of .1. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

PAGE 70

-65static structure factor hnc free ctitnara tau-a.se r3-2.ee Figure 4.5 Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rj = 2 at temperatures of kt/Ef of 0.5. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

PAGE 71

-66static structure factor k-q-a Figure 4.6 Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rg = 2 at temperatures of kt/Ef of 1.0. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

PAGE 72

-67static structure factor k-q-a Figure 4.7 Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rg 2 at temperatures of kt/Ef of 2.0. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

PAGE 73

-68static structure factor k-q-a Figure 4.8 Plots of density times the Fourier transform of h(r) = g(r) -1 versus k = q*a (q the reciprocal space to r) for a density of rj = 2 at temperatures of kt/Ef of 8.0. Free denotes the ideal fermion pair distribution function, HNC the results of the QHNC calculation where g(r) represents the analogue pair distribution n(r/v)/no. Chihara denotes the results of a yet to be discussed construction of g(r) from the analogue g(r).

PAGE 74

CHAPTER V LOCAL FIELD CORRECTIONS We have seen in the previous chapter how classical integral equations generalized into the quantal region can yield accurate results for the quantity n(rlv), namely the density distribution of interacting particles in the presence of an external potential of the form of the interparticle interaction. This in effect is equivalent to knowledge of the static response function x(k). If from this information the dynamic response x(k.,w) could be extracted, then we could easily calculate the static structure factor or pair distribution function via the fluctuation dissipation theorem S(K) = ; S(Kco) d
PAGE 75

-70This particular relation, which is one of many possible representations, introduces an unknown function e(k,z) which is analytic in the upper half plane. The physically meaningful dynamic response is given in the limit z = w + ie, and in this limit the real and imaginary parts of G(K, Z = to + ie) = 0'(Kco) + ie"(K(o) are, like the dynamic susceptibility itself, related by the Kramers-Kronig relations. The function e'(k,w) is known to be a positive, even function of frequency, as follows from the known parity of x'(k,w) and x"(k,w) and the fact that S(k,w)> 0. The utility of the function e(k,z) is that it interpolates smoothly between the known results for x'(k,z) in the limits of large and small z, and being less sensitive than x(k,z) itself, is more amenable to approximation. By subtracting out a corresponding expression for a noninteracting Fermi system we can express the dynamic susceptibility as Here we have grouped all of the actual physics into a dynamic effective potential Vg^^(k,z), actually a functional of the static response x(k) and the unknown function D(k,z). We will show through the use of the Mori projection operator technique^^^ that under reasonable assumptions the dynamic effective potential can be approximated by the static (z = 0) limit of the above equation x(KZ) = (5.1) 1 V eff (KZ) x (KZ) V eff (KZ) = ( X(K) Xq(K) )

PAGE 76

-71It should be noted that although x^{k,z) in eq. 5.1 is usually taken to be the free response function, this choice is not a unique one; x^Ck.z) plays a role of a reference system in terms of which the expression for x(k,z) is constructed. A possible improvement in the procedure may be achieved by using an adjusted choice of y.^(k,z) that includes self-energy effects, for 14-19 137 example. This method of improvement will not be taken up in this thesis in lieu of other avenues of extension. We will retain the free response as our reference system; this will allow us comparison with certain mean field theories. Mean Field Theories Mean Field Theories bypass the dependence of v^^^Ck.z) on x(k) and directly try to relate V^^^ to the bare interparticle potential by solving approximately the many body problem. (Something we have already done via our integral equations.) The efective potential is usually taken to be "static," that is independent of z, and is commonly expressed in the form Vg^^(K) = Vq(K) [1 G(K)] where V^(k) is the Fourier transform of the bare interparticle potential and G(k) is called the local field correction factor. A review of such local field function theories can be found in 1 og Kugler. The many body part of the problem is generally solved by the Vlasov kinetic equation or, what amounts to the same thing, appropriately decoupling the equations of motion for the density

PAGE 77

-72fluctuations J'*^ Whatever route, at some point collisions between particles are neglected and the particles are assumed to move in a mean or average field based on the collective motions of all the other particles. The physical significance of the correction factor is thus to incorporate some of the effects of short range correlations between individual particles not taken care of by the average field. Mean field theories for the most part have been limited to the ground state quanta! electron gas. These have disagreed on the form of the local field correction factor on three main points. 141_145 First some workers have been led to the conclusion that G(k) is a universal function of k/k^ (k^ the Fermi wave vector), 146_i4g 29-31 while others contend that G(k) is density dependent. Second, there is wide diversity of opinion regarding the value of 1 50-1 51 G(k) in the limit of large wavevector k. At this limit, a constant 1/3, irrespective of density, has been given by Geldart and Taylor. Rajagopal^^^ and others, ^^^'^^^"^^^ while that of Togio and Woodruff^^^ and Kugler^^^ use 2/3. On the other hand 29 30 Singwi et al derive this value to be l-g(O) in terms of the pair distribution function g(r) at the origin, and that of Vashista 31 and Singwi is 1 g(o) ap(^) (where a is an adjustable parameter determined self-consistently) in contrast with the result of Niklasson^'*^ who gives 2/3(1 g(0)) to this limit. Thirdly, there are several authors^^^'^^^"^^^'^^^ who have shown the local field factor G(k) has a maximum around k = 2k^.

PAGE 78

-73Kugler among them further indicates the slope of G(k) has a logarithmic singularity at k = 2k^ in addition to a maximum. On the 29-31 contrary, G(k) of Singwi and collaborators has neither maximum nor singularity in its slope around k = 2k^. It should be noted that the behavior of G(k) in the large wave vector region is quite extremely sensitive to small differences in x(k) so that it offers a stringent test for the distinction of various theories. In the case of the degenerate electron gas, there is neither computer simulation nor experiment to be used as a criterion to check these diverse conclusions concerning G(k). However, at present the 29-30 prescription of Singwi et al is regarded as fairly successful in treating a degenerate electron gas and is relatively easy to program numerically. It encompasses in the classical limit, as do the integral equation methods of Chapter II, the Hyper-Netted Chain equation, which is known to give the best fit to the pair distribution function of computer simulations. For this reason the derivation of the Local field correction factor in the STLS theory will be presented here as a basis for later comparisons. The STLS Method o Following Singwi we will derive the density response function by following the time evolution of the off diagonal single particle density function P (X I X' ; t) E < ^^(xt) t (x't)> a a CT i 1

PAGE 79

-74If we employ a Heisenberg equation of motion approach with a HamiUonian containing an external potential Vg^^Cx.t), we obtain 3t 2m ^" ^H^^^^ V^Cx't)} p^(x I x'; t) = r d [V(x X") V(x' X")] <^''^(xt) n(x"t) t (x' t)> where the Hartree potential V^(xt) E Vg^^(xt) + ; dx' V^(x x") p(x". t) depends on the i nterparti cl e interaction potential and the average single particle density distribution p(xt) = E < T 'f'^Cxt) t (xt)> S We were able to regroup the terms of the equation of motion into its present form by defining the cumulant bracket <'f'^(xt) n(x"t) ^ (x't)> E a a C <^'"^(xt) (x't)> a a a a Our equation of motion is equivalent to an infinite set of equations of motion for observable physical quantities with classical equivalents. This follows from the fact that the Wigner distribution function^^'^^^ fp^(Rt) = ; d? e" '^^'^'^ <

PAGE 80

-75has properties analogous to the classical phase space distribution f(r,Â£;t). In particular it can be used to construct physical observables such as the single particle density distribution p(Rt) = I f (Rt) po the particle current density J(Rt) = m"^ ^ P f (Rt) ^ pa pa and the kinetic stress tensor r(Rt) = m"M P P f (Rt) pa pa ^ By expanding the equation of motion for the off diagonal density distribution about its diagonal, that is in powers of r = x' x, the coefficients of the first few powers of r yield the beginning of a hierarchy of one particle equations. First we have the usual continuity equation 1^ p(Rt) + VJ (Rt) = 0 second the equation of motion for the current density m 1^ J(Rt) = V^(Rt) p(Rt) VV^(Rt) J dx VV(R x) ^ all the complicated effects of the Pauli and Coulomb hole surrounding each electron (in the presence of the external field). In lieu of an exact evaluation of this term or main aim is to extract a local field correction. Noting that ^ = n(Rt) n(xt) {g(R, x; t) 1}

PAGE 81

-76where g(r,x;t) is the non-equilibrium pair correlation function, the simplest approximation is to replace this function by its equilibrium value g(r x), thus accounting for the local depletion of charge density but neglecting the dynamics of the hole. This approximation replaces the Hartree potential by a local effective potential given by 2 Vg^^(R) = V^(Ko)) ^ G(R) p(Ku) wi th G(R) = i ; [S(K q) 1] (5.2) Relating Tr^p(k,w) to thas eaaectave potentaal an the same manner as a free particle calculation, namely as P P. f (p K/2) f (p + K/2) pa to p.^/m + ie then the continuity and current density equations together yield xCKco) = X-(Kco)/l ^ (1 -G(K)) x^CKco) 0 where Xo(k,w) is the Lindhart polarizabi 1 i tyl53 of an ideal Fermi gas f^(p K/2) f (p + K/2) XqCKco) = I : 5 pa CO p^/m + ie Although G(k) contains the unknown structure factor S(k) of the electron fluid, this can be determined by requiring consistency with the structure factor obtained through the use of the fluctuation dissipation theorem. Of use numerically is the little known direct relation between G(q) and the pair distribution function g(r)^^'*

PAGE 82

-7700 G(0)= 1 0 J dr g(r) j^CQD 0 and by using orthogonality of the spherical bessel functions 2r 2 0 g(r) = 1 ; dO 0 G(0) j^(Qr) 0 These results are extremely useful because they enable one to generate a G(q) directly from the pair correlation function and to test G(q) by constructing the corresponding g(r). In the STLS scheme both exchange and correlation corrections to the Hartree field are automatically taken into account through its self-consistency. Moreover it is straightforward to recover the Hubbard^^^ result by substituting for the structure factor in eq. 5.2 its Hartree-Fock value. The Random Phase Approximation, which is the dynamic extension of the Debye-Huckel theory, can be viewed in this scheme as the approximation that g(r,x;t) = 1 and leads to the result G(k) = 0. At zero temperature a further constraint can be used to check the validity of any approximate microscopic theory. The compressibility "sum rule" states that the long wavelength limit of the static response should agree with the result obtained for the compressibility from differentiation of the ground state energy. (This and other ground state properties of the quantal electron gas are covered in Ichimaru.^ Both the Hubbard and STLS equations violate this G(K) =

PAGE 83

-78constraint. This has led to modifications of the Hubbard correction factor'"-'" G(K) = -ir^ X r + aKp and of the STLS scheme by Vashista and Singwi^l G(K) = (1 + a |-) {^ ; ^ ^ [S(K Q; p) 1]} (5.2) where the additional free parameter "a" is to be varied such that compressibility sum rule is nearly satisfied. However it has been 79 shown that if the correction factor G(k,w) is static, it is impossible to satisfy the compressibility sum rule and the third frequency sum rule simultaneously (the f sum rule being denoted as the first moment) The Relaxation Function As an alternative approach to calculating the correction factor G(k,w), which unlike the methods covered so far easily generalizes to finite temperatures, we shall consider approximations to the Kubo function provided by the Mori continued fraction method. In Chapter I the Kubo function was identified from the difference of the dynamic and static density response functions: X(KZ) = x(K) + ipz C(KZ) In this section we shall demonstrate that this follows from the definition of the Kubo function as the equilibrium averaged density auto-correlation function

PAGE 84

-79, > >^H -XH C^j^Cr r' t t') E ^ ; dX 0 where < > stands for the thermal expectation average. It is easily shown that the Kubo function vanishes at large times; spatial Fourier transformability is insured by the subtracted term . The intimate relationship of the Kubo function with the density correlation function can by exposed by writing the Kubo function in the form C(r ?'. t t') = where we defined the Kubo transform of an operator as J E 1 ; dX e^" J e"^" 0 By taking matrix elements diagonal in the Hamiltonian of the Kubo transformed operator, it is easily seen that the Kubo function reduces in the high temperature limit to the density correlation function. (The temperature = 0 limit however is ill-defined). Other notations abound in the literature; in lieu of the Kubo transform we shall employ the semicolon bracket^"^^ P XH -XH E ^ ;^ dX Tr {p^^ e A e B} Using time translation invariance and the cyclic property of the trace one finds that the Kubo function has the property ^ ft "^^^ = p ^"^"^