v.3. .ur 5. ”7‘: $1.] . : L . a. $9.". 1:374. OJ .t..!o.1. .- 1 . ‘u NIV 9 Wm «WMWW l l 060 WMianaeafig7 4016 LIBRAA . “lg." Sta 3‘ University This is to certify that the dissertation entitled Structure, Dynamics and Superconducting Correlations of Layered Solids presented by Wei J in has been accepted towards fulfillment of the requirements for Ph.D. degree in Physics / Major professor Date August 8, 1989 MCIIi-nn Am... .1"... l ‘ 7‘ In, - v - . 0.12771 PLACE IN RETURN BOX to remove We checkout from your rooord. TO AVOID FINES return on or before one we. DATE DUE DATE DUE DATE DUE MSU Is An Affirmative Action/Equal Opportunity Institution Structure, Dynamics and Superconducting' Correlations in Layered Solids By Wei J in A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Physics and Astronomy 1989 (00400:; ABSTRACT Structure, Dynamics and Superconducting Correlations in Layered Solids By Wei J in In this thesis I study three aspects of layered solids. These are ( 1) the struc- ture of intercalated layered compounds, (2) structure and dynamics of molecular monolayer adsorbed on a substrate, and (3) superconducting correlations in re- cently discovered high-Tc layered oxide superconductors as probed by neutron. To study the structural properties of randomly intercalated layered solids, I have set up a harmonic spring model that describes both the layer rigidity and the size and stiffness of the intercalant species. In certain limiting cases, the model can be solved exactly. I also give an effective medium solution that reproduces all the known exact results, and agrees well with numerical simulations in other cases. These simulations are performed for both one and two-dimensional systems. If the two intercalant species have the same spring constant, the usual linear behavior, the well known Vegard’s law, is recovered. I compute the probability distribution of various interlayer distances and apply the results to two-dimensional alloys of ternary graphite intercalation compounds. Using a constant-pressure molecular-dynamics simulation, I have investigated the structure and dynamics of a two-dimensional molecular monolayer undergoing ferroelastic and order-disorder phase transitions. Phonons and librons have been studied below and above the ferroelastic phase transition by both self-consistent lattice dynamics and molecular-dynamics. Comparison between phonon frequen- cies for certain symmetry directions obtained by using these two methods clearly shows the presence of large anharmonicity effects in the paraelastic phase. Time dependent correlation functions and dynamic structure factors obtained in the molecular-dynamics simulations are given. Softening of phonon frequencies near the ferroelastic phase transition is discussed. Dynamic response and neutron scattering characteristics for paired fermion and paired spinless charged boson superconductors are presented. For the paired fermion system, I show that both orbital and spin currents associated with the Cooper pairs make equally important contributions to the inelastic scattering, the former dominating the later in the regime of small wave vector transfer. I suggest that a direct measurement of the the current-current response function by neutron scattering may be a diagnostic tool for distinguishing the quasiparticle statistics as well as determining the superconducting parameters of the new high-Tc layered oxide superconductors. To My Parents iv ACKNOWLEDGMENTS First of all I would like to thank my advisor Professor S.D. Mahanti for his careful and patient guidance during the last four years and for his support, un- derstanding, and encouragement to me. I am very inspired by his approach to research. His influence on my graduate study is greately appreciated and acknowl- edged. I would like to express my special thanks to Professor M.F. Thorpe for his ad- vice and collaboration on the problem of random systems, to Professor S.A. Solin for his advice and collaboration on the problem of intercalated layered solids. I would like to express my sincere thanks to Dr. A.K. Rajagopal of Naval Re- search Laboratory for his advice and enlightening collaboration on the problem of superconductivity and response function characteristics of the new high-Tc su- perconductors. His knowledge of theoretical physics and mathematical techniques has been a real inspiration to me. I would also like to thank Dr. R.K. Kalia of Argonne National Laboratory for his encouragement. I am grateful to Professor B.A. Brown, Professor W. Repko, and Professor S.A. Solin for their service on my Ph. D guidance committee. Also I must offer special thanks to Mrs. Janet King for her help, to Professor J. Kovacs for his assitance. .Many fellow graduate students, in particular those worked in Professor Mahanti’s group over the years deserve my special thanks for interesting general discussions. I would also like to thank two good friends, Mr. Yuming Chen and Mr. Hongming Xu for their friendship and humor. I would like to thank the Department of Physics and Astronomy at Michigan State University for providing me a great study and research environment and for honouring me with the Sherwood K. Haynes Graduate Physics Award. Finally, financial supports from the National Science Fundtion, from Michigan State University, and from the Center for Fundamental Materials Research at Michigan State University are greatfully acknowledged. Table of Contents List of Tables List of Figures Chapter 1 Introduction Chapter 2 Rigidity of Randomly Intercalated Layered Solids 2.1 Introduction 2.2 The Model 2.3 Limiting Cases and Exact Solutions 2.4 Effective Medium Theory 2.5 Numerical Simulations 2.6 Summary Appendix List of References for Chapter 2 Chapter 3 Dynamical Properties of Two-Dimensional Molecular Solids 3.1 Introduction 3.2 Internal Stress Tensor of Anisotropic Molecular Solids 3.3 Static Structure Factor 3.4 Self-Consistent Phonon and Translational-Rotational vii ix 13 27 37 52 55 61 65 65 73 85 Coupling in the Ferroelastic Phase 3.5 Density Correlation Function and Dynamic Structure Factor 3.6 Orientational Disorder and Phonon Anomaly 3.7 Summary List of References for Chapter 3 Chapter 4 Theory of Neuton Scattering from Paired Fermion Superconductor 4. 1 Introduction 4.2 Inelastic Scattering Cross Section from Superconductors 4.3 Dynamic Response in Paired Fermion Superconductor 4.4 Summary Chapter 5 Response Function Characteristics of Paired Spinless Charged Boson Superconductor 5.1 Introduction 5.2 Paired Spinless Charged Boson Superconductivity 5.3 Dynamic Response Functions 5.4 Quasi elastic Scattering 5.5 Summary List of References for Chapters 4 and 5 viii 95 111 117 130 132 137 137 146 153 195 197 197 199 204 212 218 219 Table 3.1 Table 3.2 Table 3.3 Table 3.4 Table 4.1 Table 4.2 Table 4.3 Table 4.4 Table 5.1 Table 5.2 List of Tables Parameters used in the MD simulation Coordinates of neighboring molecules Equilibrium atomic distances Lattice parameters from MD Limits of the k-integral Response functions for a 3D normal fermi liquid Response functions for a 2D normal fermi liquid Quasi elastic scattering cross section Comparison of fermi and hose systems Comparison of neutron scattering characteristics ix 86 104 105 106 168 169 170 194 202 217 Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17 List of Figures Top view of a layer and intercalants Side view of the intercalants Watson integral for a linear chain Average and fluctuation in height and average energy Correlation length for a linear chain Watson integral and its derivate for various lattices Equivalent spring in the effective medium theory Local rigidity in the effective medium theory Effective coupling constant for a triangular lattice A configuration of the relaxed layer A configuration of the relaxed layer Average and fluctuation in height and average energy Average and fluctuation in height and average energy Average and fluctuation in height and average energy Average and fluctuation in height and average energy Average and fluctuation in height and average energy Probability distributions obtained from simulation 12 18 23 24 26 29 30 36 39 40 43 44 45 46 47 49 Fig. 2.18 Probability distributions obtained from simulation 50 Fig. 2.19 Probability distributions obtained from simulation 51 Fig. 2.20 Effective medium for a complex lattice 54 Fig. 3.1 Ground state structure of a 2D molecular solid 72 Fig. 3.2 A diatomic molecule 79 Fig. 3.3 Rotational contribution to the internal stress tensor 84 Fig. 3.4 Radial distribution function at T‘ = 0.36 89 Fig. 3.5 A snapshot of molecular configuration at T" = 0.36 90 Fig. 3.6 Radial distribution function at T" = 0.40 91 Fig. 3.7 A snapshot of molecular configuration at T“ = 0.40 92 Fig. 3.8 Radial distribution functions 93 Fig. 3.9 Radial distribution function for a quenched configuration 94 Fig. 3.10 The first Brilliouin zone for an isosceles triangular lattice 107 Fig. 3.11 Phonon and libron dispersion curves obtained from SOLD 108 Fig. 3.12 Phonon and libron dispersion curves obtained from SOLD 109 Fig. 3.13 LA phonon dispersion curves in the ferroelastic phase 79 Fig. 3.14 Density-density correlation functions 114 Fig. 3.15 Dynamic structure factor 115 Fig. 3.16 Dynamic structure factor 116 3.17 3.18 3.19 3.20 3.21 3.22 3.23a 3.23b 3.24 3.25 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 Density-density correlation functions LA phonon dispersion curves LA phonon dispersion curves Density-density correlation functions Comparison of F(q,t) and S(q,w) at two temperatures Comparison of F(q,t) and S(q,w) at three temperatures Density-density correlation function Density-density correlation function LA phonon dispersions in the disordered phase LA phonon velocity Green’s functions for paired fermion superconductor Diagrams that contribute to the response functions Orbital current response function of normal fermi liquid Spin current response function of normal fermi liquid Orbital current response function of normal fermi liquid Orbital current response function of normal fermi liquid Sipn current response function of normal fermi liquid Orbital current response function of normal fermi liquid Density response function for a 2D normal fermi liquid 120 121 122 123 124 125 126 127 128 129 163 164 171 172 173 174 175 176 177 Fig. Fig. 4.10 4.11 4.12 4.13a 4.13b 4.14 4.15 4.16 4.17 4.18a 4.18b 5.1 5.2 5.3 5.4 5.5 5.6 Density response functions for a 2D normal fermi liquid Inelastic orbital scattering of fermion superconductor Inelastic spin scattering of fermion superconductor Inelastic orbital scattering of fermion superconductor Inelastic spin scattering of fermion superconductor Inelastic orbital scattering of fermion superconductor Inelastic spin scattering of fermion superconductor Quasi elastic orbital scattering of fermion superconductor Quasi elastic spin scattering of fermion superconductor Quasi elastic orbital scattering Quasi elastic spin scattering Chemical potential and order parameter Coherence factor of boson superconductor Current response function of paired boson superconductor Density response function of paired boson superconductor Quasi elastic scattering of paied boson superconductor Inelastic orbital scattering in the paired fermion model xiii 178 184 185 186 187 188 189 190 191 192 193 203 206 210 211 215 216 Chapter 1 Introduction A basic aim of theoretical condensed matter physics is to understand and predict structural, dynamic, transport, electrical and magnetic properties of the solids starting from a fundamental microscopic basis. There are a large number of solids such as graphite, boron nitrate, layered dichalogenides, artificial super- lattices, and the new high-Tc oxide superconductors which have layered structure. These systems show a large anisotropy in their physical properties due to this layered structure. A common feature of these solids is that the interactions of atoms or molecules in the same layer is usually much stronger than the interlayer interaction. Thus a quasi two-dimensional model is a good starting point to probe different physical properties of these solids even though the interlayer interaction can not be completely ignored. The scope of this thesis spans the study of the structure, dynamics and superconducting properties of different types of layerd solids. This thesis is divided into three nearly independent parts, the common link being that each part deals with one particular property of layered solids. The first part consists of chapter 2 where I study the structure of randomly intercalated layered solids like ternary intercalation compounds. These systems have recently become extremely interesting both from experimental and theoretical point of view 2 since the transverse rigidity of the layers can span a broad range, starting from a rather small value (floppy) for graphite to an extremely large value (rigid) for complex layered oxides. The second part which consists of chapter 3 deals with the dynamics of layered solids with particular focus on physisorbed monolayer molecular systems. My major interest here is to probe the nature of the dynamics of two-dimensional molecular solids which undergo order-disorder and ferroelastic phase transitions. In this chapter, I have investigated various dynamic properties of a model two-dimensional molecular solid which closely resembles oxygen molecular on graphite surface. In particular, I study the dynamic density correlation function, dynamic structure factor, phonons and librons in this two-dimensional diatomic molecular solid both by constant-pressure molecular dynamics simulation and self- consistent phonon theory. Finally the third part consists of chapter 4 and chapter 5 where I study the dynamic response and neutron scattering characteristics of paired fermion and paired boson superconductors. Iapply this theory to the recently discovered high-Tc oydde superconductors, almost all of which have layered structure as far as the electronic transport is (concerned, the only exception being Ba1_,K, 13120; which has a three-dimensional structure. Rather than reviewing the background materials associated with these three different topics here in chapter 1, each of the following chapters contains a rather detailed introduction to the problem under investigation. Sirnilary, the references of different parts are given at the end of each chapter excepting chapters 4 and 5 which have a common reference section, given at the end of chapter 5. Chapter 2 Rigidity of Randomly Intercalated Layered Solids 2.1 Introduction All crystalline solid solutions show a composition dependence of the aver- age unit cell volume (V) which increase with the concentration of the largest constituent.l Simple examples of such systems are: binary alloys of the type 111-33, and ternary systems of the type .41.,B,.,L.l The linear variation of (V) with :c is the well known Vegard’s law,2 although most solid solutions do exhibit a complex nonlinear (superlinear, sublinear, sigmoidal) behavior. The x-dependence of (V) depends, at the microscopic level, on the competition between local and global energies associated with forming a solid solution. These energies depend upon the relative size and compressibility of the different atomic species and over- all rigidity of the system. To address the question of the x-dependence of both the average volume (V) and fluctuations in the unit cell volume from site to site, it is somewhat simpler theoretically to consider either systems of reduced dimensionality or those with large anisotropy in physical properties. For example, there is a considerable variety 4 of ternary layered intercalated compounds of the form A1_,B,L, with 0 < a: < 1, [B(A) is the larger (smaller) intercalant, L denotes the layer] where the major expansion takes place along the direction perpendicular to the layers, denoted as the c-axis. 3"“ In the past, severalattempts have been made to study the nonlinear x- dependence of the c-axis expansion within two quite distinct types of models. One, 15'7 in which the layers are assumed to be flat referred to as the rigid layer mode i.e. perfectly rigid against transverse distortions. In this model the nonlinearity in the c-axis expansion arises from the finite but different compressibilities of the two intercalants A and B. The second, referred to as the layer rigidity model,”—16 in which the layers are assumed to be deformable but the intercalants are taken to be perfectly incompressible but having different sizes. This model has been solved in the case where there are only two allowed heights and gives a simple functional form that seems to fit the observed layer thickness well in some compounds like 0.9sz1_,Vermiculite.12'”'15 The second model has also been solved using an elas- tic continuum approximationm'17 in the low-defect concentration limit (i.e. a: close to 0 or 1). Simultaneous inclusion of the effects of both the local ionic compress- ibility and the finite layer rigidity on the c-axis expansion is the main purpose of this chapter.6 In this chapter, I set up a spring model that describes both the layer rigidity and the size and stiffness of the intercalant species. In certain limiting cases, the model can be solved exactly. I also give an effective medium solution that reproduces all the known exact results, and agrees well with numerical simulations in other cases. These simulations are performed for both one and two dimensional systems. If the two intercalant species have the same spring constant, Vegard’s law is recovered. I compute the probability distribution of the various interlayer distances and apply the results to two dimensional alloys of Li and vacancies in graphite, and to K and Rb in graphite. The outline of this chapter is as follows. Secction 2.2 introduces the spring model and in Sec. 2.3, I present results of exact calculations in certain limiting cases. In Sec. 2.4, I develop an effective medium theory and test its accuracy by comparing the results with the exact calculations. In Sec. 2.5, numerical simula- tions are compared with the effective medium theory results and finally in Sec. 2.6, I make some general remarks and compare the simulation and effective medium theory results with available experiments in graphite intercalation compounds. 2.2 The Model Consider a layered ternary system with composition A1_szL, where L rep- resents the host layer such as graphite, dichalcogenide or vermiculite. A and B are two different types of intercalants which are assumed to occupy (randomly) a set of well defined latticeosites (see Fig. 2.1). The total energy of the system can be approximated by a sum of two major contributions; one associated with the in- teraction between the intercalants and the host atoms and the other between host atoms themselves. Since I assume the intercalants to be frozen (i.e. quenched), the direct interaction energy between the intercalants does not play any role in the layer distortion. In a real ternary which consists of many layers, the experimentally measured average interlayer spacing will, to some degree, depend on the interlayer correlation. The latter results from large intercalant ions in the adjacent galleries repelling each other,5 and other similar effects. The qualitative effect of Fig. 2.1 Top View of a hexagonal lattice host layer where the interaclants form a triangular lattice 7 interlayer correlation has been addressed in Ref. 5 and will not be discussed further in this thesis. In the present analysis, I neglect this correlation and consider a single gallery bounded by two host layers with intercalants randomly occupying lattice sites inside the gallery. The total energy of the host layer-intercalant system can be written as E=%ZK(h;-— h°)2+— Z KT( (hi—hi+6)2 +- :ZKF [EUR-final] 1 i 2, 6 (2.1) where hr is the gallery height at the site i where an intercalant (either A or B) sits. The angular brackets in the second summation indicate that each bond 6 is only counted once. The terms involving the spring constants K7 and K p describe respectively the transverse and bending rigidity of the layers. In a continuum model, these terms would give an energy. density “5 —K—T(v m2. (g—Zrép—(VW . (22) i2—'( where a is the area per atom, z is the number of nearest neighbors, 6 is the distance between the nearest neighbor intercalants, and d is the dimension. This leads to an equation of motion, 82 2'62 {pV ,I “a? ”(261)? (26d) a Vh’ (2'3) p where h(r,t) is the height, V is the d—dimensional gradient operato, 6 is the lattice constant, and p is the mass density. From Eq. (2.3), we see that the isotropic layer dispersion relation at long wavelengths is 1 2KT K__p 2 2+ 2 q4 = _ 2.4 w. p (ii,— 2.) q +(—— )7, a . < ) where q is the wavevector. This phonon has propagation direction in the basal plane but displacements out of plane. When KT = 0, the positive q4 term is the source of anomalous dispersion wq ~ q2 which is seen in pristine graphite.9 The host-intercalant interaction is approximated by a harmonic spring of strength Kr(= KA or K3) and equilibrium height h;(= h0 or ho B.) In princi- 18,19 Here ple, K A and K 3 can be calculated by using density- functional theory. I regard them as two parameters characterizing the compressibilities of the inter- calant atoms (ions). The energy associated with host layer deformation has two types of contributions common to layered solids. The first one, proportional to K7, is the transverse layer rigidity,5 and the second one, proportional to K p is the flexural rigidity.16 For binary systems, i.e. AL or BL, all the h's are equal and the gallery height (or equivalently the interlayer spacing) is he, or ’12,. I define a variable 0'; at each site i by assigning it values +1(—1) if the inter- calant at that site is B(A). The corresponding probability distribution is P(cr;) = 2:6(a; — 1) + (1 — r)6(a'; + 1). (2.5) {0'} is a set of random variables whose joint probability distribution ”P( {0}) is given by I], P(a’;). K; and h? can be expressed as Kr = %(K3 + KA)+-21-(KB — KA)O’i, and h? = %(h23+ haw —(h% —- hay!“ 9 To obtain various averaged quantities of interest, such as (h) , (hA), (hg) where =1—tf';hiy (ha) : Niaghipf‘v Na : nglzv (2'6) and where pfia = A, B) is the projection operator defined by PiB’A = :0 i O'l)- (2.7) We have to minimize the energy E given in Eq. (2.1) with respect to the heights h; for a given realization of the random variables Ki(= K A,K 3) and h?(= he“ hog). I can then calculate (h), (hA), (hp), the average energy per site 8 = (E)/N, the fluctuations in height ((h — (h))2) etc., as functions of ax, KA, KB, h?“ hoB, KT, KF and the structure of the host lattice by averaging over different configurations, i.e., integrating over ’P({a})d{0'}. This last configuration average will be denoted simply as ( )av. Obviously, (0).", = 2:1: — 1, (62).", = 1, and (03),", = (0)“. One simplifying aspect of the harmonic model with uniaxial displacements, which I will refer to as the scalar case, is that I can scale the heights and the energy by (hog -— ’12,) and (hi)? — ’13,)2 respectively. If I define5 a local dimensionless height di such that h, =‘hf’4 + (1,022, -— ha), 0 g d; g 1, (2.8) then the energy E can be rewritten as 10 E = $0233 -— hi)”{ 2 KT(di — di+6)2 + Z KF [263% - dual] i +ZKAd§pf+ZKBu —d;)2p,3}. (2.9) Here the superscripts A and B in the last two summations mean that the sums go over only the A sites or over the B sites. In the harmonic model the actual value of the difference (h?1 -— ’13,) is not relevant for the physical properties of the system as it scales out. This is because all the displacements are in the same direction. The diagram accompanying this scaling transformation is shown in Fig. 2.2. The particles are attached to the lower or upper surfaces depending on whether they are A or B. These surfaces are rigidly separated by unit distance and the particles are constrained to move only in a direction perpendicular to the planes. In more complex geometries, such as those in three-dimensional mixed semiconductors,6 the arguments leading to Eqs. (2.8) and (2.9) are only correct for small values of the natural length differences h‘}, — hit The average height given in Eq. (2.6) becomes 1 1 a (d) - TV— ;dt, (dd) = N—agdi Pi- (2-10) Note that by defination (d) = (1 — $di) + mg). (2.11) There is no net macroscopic force on the system in equilibrium which leads to the relation (1 - 2)KA(dA) = 23KB(1 —- ((13)). (2.12) 11 Thus if any one of the three quantities (d), ((1,4) and (d3) is known, the other two can be found using Eqs. (2.9) and (2.10). The fluctuations in d and do, are related by ((d-(d))2) = (1*$)((dA—>)) . ' (2.14) with 77(0) = «<1 —- (d))’>- 12 (a) 7:7 7”, r'v‘V/V r Vf’I/CI’I/ DC 77"- V (b) d=1 d=0 Fig. 2.2 Side View of the intercalants fi'om Fig. 2.1, showing a larger B in- tercalant and two smaller A intercalants. In (a), the separation between the layeres are given by the variables hi, whereas the diagram associated with the corresponding dimensionless variables d.- is shown in (b) 13 2.3 Limiting Cases and Exact Solutions CompleteILFloppy and Infinitely Rigid Limits If the layers are completely floppy towards transverse distortions ( K T 2 K1? 2 0), then the energy E is minimized trivially. In this limit, I have (d) = :c; the usual Vegard’s law2 with ((1,4) = 0 and (d3) = 1. The corresponding fluctuations are given by ((Ad)2) = (d2) — (cl)2 = :r(l — 2:), and the site-specific fluctuations are zero because all the A sites have dimensionless height zero and all the B sites have dimensionless height 1. In the limit of infinite layer rigidityl3 (i.e. when KT and/or KF —+ 00), the energy E is minimized by having all the d; equal (d A 2 d3 = d). Here the energy per site 6 = $0 — and? +%2K3(1 — d)2 (2.15) is minimized by 13 = («M = = x + (1 gnu/KB (2.16) giving an energy 1 KAKB .. 6:2x(1—x)xKA+(1—at)K3 (2.11) and fluctuations in all the three heights vanish in this limit. In earlier theoretical studies dealing with the gallery expansion using discrete lattice models, all the di’s 14 were assumed to be the same (rigid layer modell) even for systems with finite layer rigidity. In this approximation, the non—Vegard’s law behavior was determined completely by the compressibility ratio K A/ K 3 of the two intercalants. At small 2:, Eq. (2.16) becomes (d) z 712—22. If the larger ion B is highly incompressible compared to A i.e., K B >> K A, one finds a large deviation (superlinear) from Vegard’s law behavior which continues for all ac. On the other hand, if the larger ion is more compressible than the smaller ion, then ((1) shows a deviation from Vegard’s law but with a sublinear behavior. If K A = K B, then one observes a Vegard’s law behavior even for an infinitely rigid layer. Later, I will show that within this model, Vegard’s law is obtained for arbitrary layer rigidity as long as KA 2K3. One Dimensional Chain To calculate the average quantities exactly for the 1D chain, it is slightly more convenient to work with the h variables, rather than the dimensionless (1 variables. One reason for studying the 1D chain problem is that the effective 'medium theory to be discussed in Section 2.4 can be worked out analytically for the 1D case and compared with the exact calculations to test its accuracy. The energy for the 1D chain is given by 1 1 1 E1. = 5 Z: K.(h.- — h?)2 + 5 2: Km. - hi+1)2 + 5 Z Kp(2h.- — hi+1— hi—1)2- (2.18) Minimizing EM with respect to the hi, I find Mh = «r, (2.19) 15 where M is a tridiagonal matrix with random diagonal matrix elements, i.e. Mii = Ki + ZKT + 6KFa Mini] 2 —-KT — 4KF, (2.20) Mi,ii2 = KF 01' I{1 +2KT+6KF —KT-4Kp 0 —KT-4KF Ifz+2KT+6KF 0 M = KF -KT—4KF 0 , 0 O KN+2KT+6KF . (2.21) and h and Q are two coluom vectors h, K1 h‘,’ h2 thg = , = , (2.22) hN Km), In Eqs. (2.18)-(2.22), all quantities are real. For a given random configuration of {0’}, let the eigenvalues and eigenvectors of the matrix M be Ag and (fig. Define a N x N (random) matrix Q by N Q(’\: {0}) Z Z bl" — AQ)¢:¢Q: (2'23) «1:1 then I have h = M’% = d—AA-QQ, {a}). (2.24) 16 The linear chain case can be regard as a special realization of a general class of one—dimensional disorder problem,20 such as, electron localization and electron “'22 random harmonic transport properties of one-dimensional disordered system, chain system.23 This kind of problem can be solved by Green’s function method24 which needs the solution of certain integral equations,2°'21'23'24 When both K ,- and h,- are random variables, which is the case when K3 75 K A, the problem is difficult to solve analytically. Our numerical solutions in this case will be discussed in Section 2.4. However when K 3 = K A = K, one can obtain the hi’s by diagonalizing the M matrix, which is no longer random. The results can be worked out analytically. In this limit, the eigenvalues and eigenvectors of the matrix M are given as Aq = K + 2KT[1 — cos(q6)] + 4Kp[1 — cos(q6)]2, (2.25) and ¢q(n) = 1— eiq" (2.26) with q = 27rr/N with r = 0,1,2,. . .N — 1. Then I have — h>-§(Z(M-‘).—n ) = 7},— f some»: 2 @4026»? on K : AoN Zav : 33,108 + (1 - 33)}131 (2.27) Thus I find ((1) = a: for all values of KT and K )7. However the site-specific heights (d3) and MA) do depend on the layer rigidity parameters KT, K p. Consider Using I obtain and similary (ha) 27v— 1 N = __ h 3 N$ av 2 1, / «was» 2 Z {imwnx 1 + ‘" >h2. q 1,12 q = 3,1; 2 Z §¢.(z>¢;(n)(<1:”')ha) , q 1,71 q an (0103,)“, = (2:1: — 1)2 + 4:1:(1 — 2:)6zn, (d3) = 1 - (1 - av)[1 - W(K)!. (dA) = a:[1 — VV(K)]. 17 (2.28) (2.29) (2.30a) (2.13%) The Watson integral W(K) which is a measure of the layer rigidity is discussed and calculated in the appendix. Figure 2.3 shows W'(K) for linear chain for various values of K p/ [(7. Note that the expressions in Eq. (2.30) obey the general relation in Eq. (2.12). As the general expression for VV(K) given in the appendix is rather complex, I give the results here for K p = 0, i.e., 18 L0 'vv'1*"'I*'*'I""I"v > (1). K'/Kf = O 0 i : (2). Kf/Kf = 1.0 l 0.5 L. (3), Kr/Kf 3 10.0 _ ~ (4). mm, = 100.0 . . I I L 1 0.6 - .. l . " . z I 0.4 1 L C 3 ‘ 0.2 .. . 4 : 0.0 *“l““““““*"”*l 0 1 3 3 4 5 K/xr Fig. 2.3 The Watson integral W(K) for the linear chain, for various values of Etc/KT. 19 1 (d3) = 1 — (1 — 2:) [1 — W] , (2.310.) and 1 (2,.) = 2: [1 — W] . (2.315) Clearly when KT 2 KF = 0, i.e. in the floppy layer limit, the results in Eqs. (2.31) become ((1,4) = 0 and (d3) = 1 and when KT or K}: 2 oo, i.e. in the infinitely rigid layer limit, (d3) = ((1,4) = 2:. Thus, in the limit KA = K3, when ((1) shows a Vegard’s law behavior the transverse layer rigidity effects are seen in the z-dependence of (dA) and (d3). This point will be discussed in more detail later. The fluctuations in d, dA and d3, can also be calculated exactly in this limit. I have 1 N (hangs), (2.32.) 1 K2 ,. , 1'+ 0' = 7v? / d{0}7’({a})§ W[a¢,(z)¢,.(z)¢q(n)¢q,(m)( 2 Utah?" 1 K2 ur- a 1 + 0 = E Z M,1;".¢.(t)¢¢(t)¢.(n)¢.r(m)<(———2—‘flwnh?.)w- (29%) Using (alanam)av : 6ln6lm<03)av + [6111(1 - 6(m)(1 " 6rim) +6lm(1 '" 6ln)(1 - 612m) +6...,.(1 — 5..)(1 — 61m)](a)a,, +(1 — 5,.)(1 — 5,...)(1 — 6m)3. : (223 —1)3 + 43(1 — $)(2$ — 1X6"; ‘i‘ 6lm + 6am _ 261n6lm)9 I obtain ((da - («1302) = ((dA — ((12492) = 3(1 - a=) [W1(K) - W(K)2] 23(1— 1 [1+2KT/K _1] 33 )1+4KT/K mm In addition, I have ((61 - ((9)2) = 3(1 - $)WI(K), I. + ZKT/K : ”(I - ”)(1 + 4KT/K)3/2’ and the average energy 3 = %K$(l — 33)“ " Ill/(KN: éKxfl—x) [1— 1 ], V1 +4KT/K 20 (2.33) (2.340.) (2.34b) ° (2.35a) (2.355) (2.360.) (2.366) where the explicit results in Eqs. (2.34)-(2.36) have used Kp = O for simplicity. From the above equations we see that the fluctuations are maximum when 2: = 0.5 and K p = 0 as would be expected. Note that the expressions for the fluctuations in Eqs. (2.35) and (2.36) obey the general result Eq. (2.13). The results for a 21 linear chain with KA 2 K3 2 K = KT and Kp = 0 are shown in Fig. 2.4 together with simulation results obtained by the method described in section 2.5. To calculate the intrachain correlation, I first consider the average <-11\7; hihl+R> = — N] d{a}r<{a})ZK , A X 209.0 + R>¢;(n)¢;.(m)h3.h:. q l,,n m 1 K2 e-iqR Ir- 0 0 I N 2A: — [Z ¢q(n)¢9(m)60 ' (2'37) Using Eqs. (2.29) and (2.14), I obtain 77(3) = 2:0 — 01220912), (2.38) where (2.39) .>:.| a 1 W2(K,R) = N296 For the linear chain with K p = 0, I obtain 1 1I' If2 . W R = — "IR 2”" ) 21r [quf + 2KT(1 - cosq)]2 e ’ 2 (1 + 4KT/K)3/2 1 + ZKT/K ZKT ’ (2. 40a.) —R 1 2K K R,/1 4K K , + 7” [1+R + T/ J [1+——-—+———\/1+4RT/K 22 = W K 0 1 R V /5. . 2( ,)[+ 1+2KT/K]e (24Gb) Then the correlation length 6 can be identified as K K -1 : ____ ~ __ ,_ - { log [1 k 2KT + ZKT‘A ( 4KT/I1]. (2.41) Figure 2.5 plots the correlation length for the linear chain. One sees that when KT —+ 0,6 ——2 0 and when KT ——> 00,6 —+ 00. 6 can also be thought of as a healing lengthlz'1° whose value is controlled by the transverse layer rigidity. 23 1.0 T I I (t) 0.5 - . I" 0.8 P .1 L (d) 0.‘ r- .. 0.8 0.0 0.15 " 0. 10 [(4.5 - ']"' Energy Fig. 2.4 Showing (a) average height, (b) fluctuations in height and (c) the energy in units of %K(h"a - ’12,)2 {or a one dimensional chain with KA = K3. The solid lines are exact results and the solid dots represent simulation results. 30 _ 20 -— - ‘u‘o a .2 a .2 3 O E O 10 - - o . w 10’310'310'1100 1o1 102 103 KT/K Fig. 2.5 Correlation length for a one dimensional chain model. 24 Two-Dimensional Lattice The average energy, average heights and the corresponding fluctuations can also be calculated exactly for the two-dimensional systems in the limit K A = K 3 = K. One has to replace Aq in the Watson integrals W(K) used above by Aq = K + KTz(1 — 7g) + Kp[2(1 — 7.1)]2, (2.42) where 1 . 7. = - Z 6"", (2-43) and z is the number of nearest neighbors (z = 4 for the square lattice, z = 6 for the triangular net). Clearly the geometry of the lattice will determine Ac, and hence the layer distortion characteristics. The results are formally the same as for the linear chain, except that the appropriate Watson integrals from the appendix for two dimensions must be used. We plot the integral VV(K) and its derivative W1(K) for various lattices as a function of K/(zKT) in Fig. 2.6 for Kp = 0. For the 2D system, I obtain (d3) = 1 _ (1 — .)[1 — W(K)], (2.30.) (.A) = .(1 — W(K)], (2.30.) ((dB — (‘13))2) = ((dA — WAD?) = “3(1 — 2:) [W1(K) — W(K)2], (2.34.1) ((d — (an?) 2 .(1 — 3:)W1(K), (2.35.) e = 312—1941 —- .)[1 — W(K)]. (2.36.) The average height obeys Vegard’s law in this limit for all the lattices. I will compare these exact results with numerical simulations in section 2.5. 1'0 rfrT""I"rfjvvvavr. I w\ , o.e - /' I ~- ; L I, x" 1 f ("I ' f /' . .3 » I 3" ( I EOJT II I. \ ‘ 3 , I I w, ‘ .5 1‘ e ' l f ‘ g 'l f . £04 I t .1 d . . t l l Solid Lines: 2 =- 2 l l Dashed lines: 2 = 4 , i Dotted Lines: 2 =- 6 i or i d . l ‘ t l l . ' l L-LLLA-AL1-AAALLA‘A1LL‘L 0.0o . ‘ e e m K/ (21%) Fig. 2.6 The W(K) integral and its derivate W1(K) are shown for various lattices as a function of K / (2K7) and for K p = 0, where z is the number of nearest neighbors. Here 2 = 2 is the linear chain, 2 = 4 is the square net and z = 6 is the triangular net. 26 27 2.4 Effective Medium Theory In this section I develop an effective medium theory. The work in this section follows the general ideas of Garboczi.25 The effect of the layer on the intercalant atom is contained within an effective local spring constant Kc. This can be found by applying a force F to a single site 0 in the non-random system where K = K A = K B as shown in Fig. 2.7(a). The equation of motion is F = —g,% where the energy E is given in Eq. (2.1). The right hand side of this equation is linear in the {hi} and can be inverted to give the diagonal term ho = W(K)F, (2.44) where the effective spring constant for this kind of displacement Kc is given by K “=va (2.45) and W'(K) is the same Watson integral that I have used in the previous sec- tion. The whole system is now replaced by a single spring Kc as shown in Fig. 2.7(b). The problem is now reduced to just two springs in parallel as shown in Fig. 2.8. One of these springs is K0,, where a can be either A or B with probability 1 —- a: or a: respectively. The other spring is K; = K. - K, (2.46) formed by removing the spring K. Later I will identify K as the effective medium spring constant which will be determined self-consistently. The total energy per site e for a single impurity is given by 1 1 e = 5K;(h — h...)2 + §Ka(h — ha)”, (2.47) 28 where he is the effective height which like the effective spring constant Ke is to be determined self-consistently. We are considering only a single impurity spring K in Eq. (2.47). We minimize the energy he with respect to h and obtain the local height h Kéh. + Kah‘; h = K' + K . (2.48) Substituting this back in Eq. (2.47) gives the energy for a single site a (”impurity”) KgKa m. (2.49) 1 6a 2 Elhe - ’12.)2 Because the springs have different lengths, the postfactor in Eq. (2.49) comes from adding the two springs in series. A general energy expression can be written down by summing over all the different types of ”impurities” which in the spirit of effective medium theory are assumed to be non-interacting. _ _ 1 o 2 KLKO, ‘3 “2P0“! - 2 2 “(he —-h.) 271?”; (2'50) 0 a=A,B e where 1 "L a = — 9'. 2.51 The above expression of energy Eq. (2.50) contains two unknown parameters, the height he and the spring constant K; (or equivalently Ke or K). These are determined by two conditions. For the first I minimize the energy 6 with respect to he to give 29 (3) F Mex.) 000000000 0000000000 0000000000 000000000 " ' " ‘ >6 : : : : :3 :K :K :K 3K (SK 2 I I E E 5000000000 ‘000000000'Cg‘0000000000.‘000000000 (b) Fig. 2.7 Illustrating how the equivalent spring Ke is determined within effec- tive medium theory. Fig. 2.8 Showing how the local rip'dity is determined within efl'ective medium theory. Here X; = K. - K and a = A,B. 30 31 ZpahaKa/(K; + K.) h. (2.52) 2 ZpaKa/(KL + Kc) This expression can be regarded as variational because a minimization is involved. For the second condition, I use the following argument. Apply a force F at the sin- gle impurity site. This produces a displacement F / (K; + K) which when averaged over all sites is set equal to F / Ke to give the self consistency condition 1 1 a _ = —— 2.r 3 :13" Kg+Ke K.’ (”l which for the case of interest here with two intercalants A, B can be written as 33K_KB KL+KB + (1 — ”XXL—1:121— : 0. (2.54) Various solutions for K from Eqs. (2.45) and (2.54) are shown in Fig. 2.9 for the triangular net and compared to the virtual crystal result K” : 3K3 + (1 — 2:)KA. It can be seen that the results used here, which agree very well with simulation results as discussed in the next section, are very different from the virtual crystal result. Together Eqs. (2.52) and (2.53) provide the effective medium solution to the problem. Combining the two, I can rewrite Eq. (2.52) as 0 1(1(° he _ Ea paha e( L a), (2.55) which for the case of the two intercalants A, B of interest here, can be written in terms of the dimensionless variable d as K B x— (d) _ K g + K — :c—KB +(l—z) KA K; + K 3 K g + K A (2.56) 32 After some manipulation Eq. (2.56) can be put into the form (d) = a: + 23(1 — :6)Fd, (2.570.) where _ KeK;(K3 - KA) F“ _ K(K; + K.)(K; + K.) (2‘57”) Either directly, or using Eqs. (2.11) and (2.12) I can find (dc), i.e., KeK'KA d = 1 — 1 — e , 2.58 < B) ( ”MK; + K.)(K; + K.) ( .) and I ((1,4) = K°K°KB (2.58b) ‘” K(K; + KA)(K; + K.) The energy can also be found by substituting Eqs. (2.54) and (2.55) back into Eq. (2.50), and using the probabilities appropriate to two kinds of sites. I find 6 = %./KAKB .(1 — .)(h‘;, — 1.2,)2141, (2.59.) and I Fe K°K‘ V KAKB (2.596) : K(K.’. + KAXKL + KB). The effective spring constants K, Ke : K; + K are determined by Eq. (2.54). Note that the relation Fd KB—KA —— = —, 2.60 F. m l ) gives a useful relation as it is independent of any of the effective medium parame- ters. Although the derivation of Eq. (2.60) is within effective medium theory, it is 33 in fact an exact relation as can be proved using the Feynman-Hellman theorem.26 In the present context, this theorem states that 85 53(5) _ (a), (2.61) where the energy 6(A) = E / N depends on some parameter A. From Eq. (2.1), one sees that A can be set equal to h); or hi2, which leads to Eqs. (2.11) and (2.12) and also allows us to demonstrate that Eq. (2.60) is an exact result. Of course effective medium results are exact in the limit of single impurities where the linear terms in a: or 1 — a: are given correctly for all quantities that I have calculated. Explicit results in this limit can be found from Eqs. (2.56)-(2.60) and for small a: l‘” ___ ” 1 + (KB/Iiiilfiwmn’ (2'62“) («1.) = a: 1l1(f{:v/irfquyIffvi/fza) (2.6%) (d6) = 1 + E;:;§:):Vl($()KA), (2.62.) e = $10. .- (h‘), — 52,)2 1 + ( KiiKTEK1’31W K.) (2.62d) Similar expressions for small 1 -— a: can be found by replacing a: by 1 — 2:, d by 1 — d, and interchanging A and B in Eqs. (2.62). 34 The effective medium results derived in this section contain the previous, exact results as special cases. For infinitely rigid layers, when KT or K p are infinite as discussed in section 2.3, Eq. (2.54) becomes K = :L'KB + (1 — 2:)KA (see Fig. 2.9) and the effective medium results of this section reduce to Eqs. (2.16) and (2.17). When KA 2 KB 2 K, using Eq. (2.45) that relates Ke to the Watson integral W(K), we recover the results of sections 2.3 for ((1), (da) and (e). The fluctuations cannot be easily calculated within effective medium theory. However they can be obtained using the Feynman - Hellman theorem Eq. (2.61), with A set equal to K A and K 3 respectively. Using the effective medium Eq. (2.54) and differentiating Eq. (2.50) with respect to Ka(a = A, B), I find KAKBKg 2 K(Ké + KA)(Ké + KB)(Ké + K.) W(K)2 (K — KA)(K — KB) ' W1(K) - W(KV (K; + KA)(KL + KB) ((da - (da>)2) == 3(1 - 3) (2.63) Using Eq. (2.63), the fluctuations in d can be found using Eq. (2.13). W1(K) KAKBKe 2 WHO - W(K)2 K(K; + KA)(K; + KB) W(K)2 (K — KA)(K — KB) W1(K) — W(K)2 (K; + KA)(K; + K3) ((4 - (d))2> = a3(1 - 33) (2-64) One notes that in the limit KA 2 K3 2 K, the exact results Eqs. (2.34a) and (2.35a) are recovered. Thus the use of the Feynman-Hellman theorem allows me to extract fluctuations in d. Without the use of this theorem it would not be possible to obtain these fluctuations in an effective medium theory. In the next section, I will test how good the effective medium approximation is in general, by comparing 35 results with numerical simulations. I stress that in the appropriate limits, the effective medium results derived in this section, reproduce all the exact results of section 2.2. In particular, the completely floppy (KT 2 K F = 0), the infinitely rigid (KT and/or Kp —0 00), and the KA 2: K3 cases are all given correctly. 36 1.5 TV‘V lrvvtlvrrrlviTTIVVVV ~ ig/x. - 0.05 4 , (1)0 Icf/KI - 0'0 . (2). Ke/K- - 0-1 (a). WK. - 1.0 ‘ . (4). a/x. - 10.0 "71.0 3‘. 4 ’1‘ S I 3 + d a? N P " "\" :40!” 2 "‘ . l ‘ 0.0 ....L...-L-L..l.m.-L.LL 0 0.8 0.4 0.0 0.0 1 X Fig. 2.9 The effective coupling constant K of the effective medium for a tri- angular lattice is compares to the virtual crystal result 3K3 + (1 - r)KA. 37 2.5 Numerical Simulation In this section I describe the numerical simulation procedures and analyze the results obtained. I will also compare these numerical results with the same quantities obtained from exact solutions when K A = K B and from the effective medium theory when K A 75 K B- The problem of finding the stable structure of the random alloy can be regarded as an optimization problem for the energy function Eq. (2.1). To find the equilibrium structure, the total energy is minimized with respect to the variables {h}. At zero temperature, the minimization procedure I have used is based on the conjugate gradient method?"28 I have found that this method is more powerful in terms of computer time and accuracy compared with other methods such as direct matrix diagonalization or the relaxation method.” However, to study the finite temperature structural properties, the molecular dy- namics simulated annealing”'30 total energy minimization method:“'32 would be more appropriate. I have performed extensive studies for the triangular lattice where typical lat- tice sizes used were N = 50 x 50 = 2500 nodes with periodic boundary conditions. All the simulations were done on a VAX 8650 computer. The values of our pa« rameters are chosen to be close to those of Li and vacancies in graphite,”’9 and K and Rb in graphite.1°'n The details of the numerical procedure consist of the following steps for each concentration x: (a) First an initial configuration is generated. This was done by gener- ating a two dimensional triangular network of N nodes with each node labeled sequentially and assigning vertical springs K A and local heights h‘}, to each node. Then I randomly select an A-type node and replace 38 it by a vertical spring K B and local height h% until the total number of B-type nodes is N B = 3N. All nearest neighbor horizontal springs are uniform and given by KT (and K p). (b) The system is relaxed using the conjugate gradient total energy min- imization method and a final configuration characterized by heights {h;} is obtained. (c) The average height, fluctuations in height and average energy, etc., are obtained for the above equilibrium (minimum energy) configuration. (cl) The above steps are repeated 1000 times to obtain an ensemble aver- age. These numerical simulation procedures were carried out for several values of K A,K 3, and KT (while keeping K F = 0). Figure 2.10 shows a typical configu- ration of the relaxed layer at x = 0.2 and K A = K B- One sees the relaxed layer is not flat, but rather has the appearance of ”waves on the ocean”. Figure 2.11 shows a typical configuration of the relaxed layer for K A 79 K 3 at :0: = 0.2. It also has a lot of fluctuations from site to site. As KT —0 00, the layers do become flat as discussed in section 2.3. The simulation results of the average heights, av- erage fluctuations for both the A and B intercalants, and the energy are shown as solid dots in Figs. 2.12-2.16. These figures are all for the triangular lattice. In Figs. 2.12 and 2.13, the simulation results are compared to the exact solution. The agreement is excellent for all quantities, thereby giving us confidence in our simulations in that we are using large enough lattices etc. In Figs. 2.14, 2.15 and 2.16, the solid lines are the effective medium results for the above quantities. 39 £95325? 50¢ nonmfihc m.— m:=mo.~ 23.. .md H a we ES cA H mv~\.~..< 64 H mvwfwi 5.5. .852 SEMEWWS @3233 was cozeuzwwzg 4. Omé .wmh I’ .‘O‘.’ 00....‘OOplI'; 014.; 505,5000000 z. ..5....~...¢....h 5...»..5 95...... . O s I 0 u 0 . IN;’ Ifi0flowCOOOOOOOOOOOOOOOO/ 010.00” I”.9.l)flb. JO}. MMM’O‘HOOOIOJO: O... 4 s 00 :0 . . 00 00 . 6 O I o I O l O o O 000:; I OIVrb) 5 .0 0. 00090002’00...0000060)0; . . .... 0 I 0 I di0l.. O O :. .0 o... ».o.n.,......3 - . --v...ou....,hoo..n.»...o 'I I, .00 le§0l;nb.O.I. .636»... ...V:...... .A ow»... 40 £63335? 50¢ 32—536 m.— mzsmou oars .Né H H as was fio H mv~\.~.u~ J6 H mM\ aw . 5......» 2...... CI ~O...‘ ’l‘ O o o s . ‘00“.Nsoo I 1.560;, 0010,00 O s 3......“ ........,,..,...... 9.3.... ......ua.w.....,....£.6....... to.) I I. o P u I. ~N.;‘..:o..‘>. \ O ‘ ‘00 0’. ’1‘4‘ §>:“ . .5. 0‘0: u 0 ‘00... 0.0.P0Wu0vu” .‘~ 41 The effective medium results provide excellent fits to the simulation results. Indeed the quality of the fits is almost as good as for the exact results in Figs. 2.12 and 2.13, except perhaps for the fluctuations in da (with a = A, B) where there appear to be small systematic differences. The effective medium results are exact for small a: and small 1 — 2:. However effective medium theory is only an approximation in general.8 Nevertheless the quality of the agreement in Figs. 2.14-2.16 is excellent and means that effective medium theory can be used with considerable confidence in interpreting experimental data. When KA = K3, the average heights ((1), ((1,4), (d3) show straight line be- havior; with (d) obeying Vegard’s law. The fluctuations and average energy are symmetric about a: = 0.5. The difference (.3) — ((1,1) = W(K), (2.65) is independent of :6, as can be seen from Eqs. (2.30) and is illustrated in Figs. 2.12 and 2.13. This difference depends only on the transverse layer rigidity parameter KT and decreases as the rigidity of the layer increases. In the limit of a single B impurity in an A host, we have (d3) = W(K) from Eq. (2.30). Thus for rather stiff layers, as in Fig. 2.12. this intercept is quite small ~ 0.2, whereas for rather floppy layers, as in Fig. 2.13, the intercept increases to ~ 0.65. For infinitely rigid layers W(K) = 0 and for perfectly floppy layers W(K) = 1. Thus we can say that W(K) is, through Eq. (2.60), a direct measure of the ”floppiness” of the lattice with 0 < W < 1. Conversely 1 — W may be regarded as a measure of the transverse rigidity of the lattice. Note that the strain energy in the lattice e, given in Eq. (2.36a), is also directly proportional to the rigidity of the lattice. For a 42 completely floppy lattice, the strain energy is zero as the lattice can accommodate any impurity at no cost in energy. When K 3 > K ,4 (Figs. 2.14, 2.15, 2.16), the average height shows superlinear behavior and the partial heights associated with the A and B intercalants are no longer linear. The average energies and fluctuations are peaked at 1: < 0.5, and no longer symmetric about a: = 0.5. The effective medium theory reproduces the simulation results quite accurately when K A 7E K 3. i; 2 8 L' 2 0 Fig. 2.12 Showing (a) average height, (b) fluctuations in height and (c) the average energy per site in units of K(h°3 - ha)“, for a triangular net with K A = K 3 = KT. The solid lines are exact results and the solid dots repre- sent simulations results. 43 1.0 r .. (a) r r "' ' 0.0 " q .9 0.. .. m _ 1.0 ‘ HI. I 0.1 V as - q as b q an us are m m m 0 0.0 0.0 0.0 0.0 1 2 Fig. 2.13 Showing (a) average height, (b) fluctuations in height and (c) the average energy per site in units of K (hos - ’02,)“, for a triangular net with K A = K 3 3!: K7. The solid lines are exact results and the solid dots repre- sent simulations results. 44 45 82’; (d) p . 1 Fig. 2.14 Showing (a) average height, (h) fluctuations in height and (c) the average energy per site in units of K 301% — hfl)’, for a triangular net with K 4 ¢ K 3. The solid lines are exact results and the solid dots represent simu- lations results. 1.0 T 0.. "' (.) 1 0.. " «I A “U V . 00‘ " q Fig. 2.10 Showing (a) average height, (b) fluctuations in height and (c) the average energy per site in units of K 301% - hfi)’, for a triangular net with K ,4 # K a. The solid lines are exact results and the solid dots represent simu- lations results. («13> — «par/2 p.615 - - — a . t. 0.010 - a O :3 I33 0.005 - .. 0000 L L l l o 0.2 0.4 X 0.6 0.8 1 Fig. 2.10 Showing (a) average height, (h) fluctuations in height and (c) the average energy per site in units of K 3(h"3 — ha)“, for a triangular net with K A 50’: K 3. The solid lines are exact results and the solid dots represent simu- lations results. 48 The nature of the disorder effect in the structure of random alloy can be better determined by the partial probability distribution functions of the heights for the A and B-type intercalants. The height probability distributions P(d), P(d3), P(dA) are defined by N 1 P(d) _ ‘N’ I; 6(d — d;), (2.66) and 1 N d.) = EZMM- (11);); , (a = A,B). (2.67) 1:] I plot these distribution functions obtained from the simulation in Figs. 2.17- 2.19. Figure 2.17 shows the probability distribution at a: = 0.4 for a system with KA 2 K3 2 K7. One sees that P(dA) and P(d3) are symmetric around the average values (dA) and ((13). The line shapes of the two peaks are close to Gaussian, since from the simulation results, the moment ratio ((Ad...)“)/((A.l.,)2)2 = 2.8 :l: 0.1 which is close to 3.0 which would be obtained for a perfect Gaussian. The re- sult cannot be exactly a Gaussian, as the distributions are bounded, whereas the Gaussian distribution is unbounded. When K A = K 3, it is possible to calculate the low moments exactly as shown in section 2.3. Figures 2.18 and 2.19 show the probability distributions at a: = 0.4 for systems with KT 2 K A = 0.1K 3 (Fig. 2.18) and KA = 0.1K3,KT = 0.05K3 (Fig. 2.19). It is clear that P(d3) and especially P(dA) are not symmetric when K A 79 K B. 49 $33 .23 Away .334... 3:: 3.64.2.1 13:59» 948 .mcomgauawwzon. =9: .6 D m a O 06 0.6 #6 «6 .3 s t - — - - - — — —. .6 a... .6 v6 «.9 3 0.9 06 v6 .545 2335.3 3.3 map—31.5.5 53.: 0.3 .3153 punch. .mM H VVN 5...: 09.53 535.3 8 x S e 8.. as see .3. e6 86.33398 5.4338; S." 3F.— - «.9 o d]... .. a .. r. .. H -N in — Izmq'qozd 50 .33 ~23 AVE 32? 8:: 30:33 13.3.5; .58 homagu 33 33.3: 35.3 59¢ an 333.— 0345 2:. .maoueamwag :2: .mk x Q .32 8.33 32.9.3: an x cm a Eu 3. was 1% i .3 maczzamhfii bur-cecal m—é .wfm 6 06 v6 O6 an §I1WI1I14Ifi .6 .6 0.: v6 «.9 3 #7 ”h u a... a... v... a... _ _ .n . ”r .«r .TL. . J... A; m .“n .n a." :m". .m an m L. .. 11.“ 1....” "r Vela u. a gang}. M 1... ~.ou-a\ 04H. .m:o_.a._:mw..ou 2:: ~95 33:23 «:3 macaw—35.3 89¢ 0.3 3:80.- 3048 .QN um > 1. We have not attempted to fit the experimental data precisely but a choice of KA/KB = 0.1, and KT/KB = 0.1 (see Fig. 2.14a) can semi-quantitatively fit the experimental data excepting for :2: == 1 where anharmonicity effects may be important. The potassium-rubidium ternary data, where the gallery expands from 5.47.31 for a: = O to 5.68.71 for a: = 1, shows nearly a Vegard’s law behaviour.10 This can be understood if we assume that K A / K 3 z 1 which seems physically reasonable.5 In this case the strength of the layer rigidity can only be determined from measurements of local quantities such as (CIA) and ((13). We would very much like to see experiments performed that 53 couple to these local structural parameters. This could be done using extend x- ray absorption fine structure spectroscopy (EXAFS)33 or NMR. The experimental situation is more advanced for mixed semiconducting alloys34 like GazIn1_,As, where many compounds have been studied using both X-rays and EXAFS. The results for (d), (dA), ((13) are remarkably similar to those I have obtained in Fig. 2.13. The results obtained in this chapter are more general than the simple geome- tries and force constant models that I have used. These have served our purpose in illustrating the kinds of behavior to be expected. However, effective force con— stants can be defined for more complex geometries as illustrated in Fig. 2.20. Here the intercalant atoms have many springs attached to common points in the layers above and below. By applying a suitable force field, as illustrated, an effective force constant K, can be defined that leads to the same local displacements as in the actual system. This effective spring constant Ke can be written in terms of a suitably generalized Watson integral for the lattice. This integral will contain all the information that is necessary to describe the rigidity of the layers and can be computed if a suitable force constant model for the layers is available. In practice the Watson integral W(K) is not very sensitive to such details as can be seen from Fig. 2.6 which compares W(K) for various lattices. Having obtained the effective spring constant Ke via the Watson integral W(K), the effective medium equations developed in section 2.4 are quite general. (a) F F (KrKF) OH'HOHQHH ONCOOO'OOOI QOOOQO'OHOOH\ titlttttlttl Ottttttlttttttt Ottttttltttttt (b) F E K F Fig. 2.20 A illustration of the side view of a more complex lattice than shown in Fig. 2.4. The equivalent spring K. can be determined within effec- tive medium theory if appropriate forces F are applied as shown. 55 Appendix In this appendix, I calculate the Watson integrals W(K) used in the text. In general IV(K) is defined by6 . _ i f: z _A__ K W (K) _ N 2“: Aq (Z‘Ir)2 quK + KTz(1 — 7g) + KF[z(1 — 7.,”2 ’ (A1) where the q-integral is over the first Brillouin and the factor 7., is given in Eq. (2.43) The dispersion relation for the vibrational modes is may: 2 K + KTz(1 —- syq) + Kp[z(1 — yq)]2, (A2) where m is the mass of an atom in the graphite layer. This dispersion relation corresponds to modes of vibration of the whole sample, containing many layers, in which the wavevector is in the plane and the displacements are all perpendicular to the plane. The springs K A and K 3 are irrelevant for this motion as they are not stretched. These modes could be seen for example in inelastic neutron scattering experiments.35 In the limit q —+ 0, 1 2_ 22 7.41—5;(q.6)—1—6q/2d which gives Eq. (2.4). 56 It is convenient to introduce a Green function 9(2) defined by 9(2) = if; m _ Zfl _ 70‘), (A3) from which one can define a spectral weight p(:z:) by pm = firmgu) = 732 Mm — 2(1 — 7.)], (A4) where a: is real and p(:z:) is non zero and positive definite within the band 0 < a: < 22. The Watson integral Eq. (A1) can now be written as an integral over this spectral weight Kp(:c)d:c WK = . ( ) K+qu~:t:-+-vaa:2 0 (A5) The integral Eq. (A5) is well behaved as 'the denominator contains no poles. Note that for large K >> KT, va W(K) —+ 1, (A6) as we would expect. For small K << KT,K F the Watson integral is dominated by the small a: behavior, depends only on the ratio K / KT, and goes to zero as K goes to zero in the 1D and 2D cases of interest here. Another related quantity W1(K) that I require, can be found directly from the Watson integral __1_ E 2- 22 K’p(x)dx __ 2 3 W(K) W1(K)_NZ?(AQ> —0 (K+KT:12+KF:02)2 _ K 6K[ K ] (A7) 57 However, W2(K) in the correlation function Eq. (2.38) is not simply related to I'V1(K) and has to be calculated directly. One Dimensional Chain In the case of one dimensional chain, 2 = 2 and 7., = cosq6, so that 50(3) = 1 (A8) and the Watson integral is given by W(K) — 1 \/1+ 4KT/K + 16Kp/K 1/2 x 1 + SKF/K . (A9) 1 + ZKT/K + \/1 + 4KT/K+ 16KF/K When K is small this result becomes W(K) = -:-(/ 1g: (A10) independent of K F as discussed above. For the case when K F = O, the result Eq. (A9) simplifies to l I’V K 2: , A11 ( ) \/l + 4K/KT ( ) and when KT 2 0, the result Eq. (A9) becomes 1 1 \/1 16K K ”2 W(K) = + + F/ , (A12) \/1 + IGKT/K 2 58 The quantity W1(K) is found by differentiation of Eq. (A9) and for K F = 0 is given by W,(K) = (1 + 2KT/K)(1 + 4KT/K)‘3/2, (A13) For general values of Kp/KT, I plot W(K) in Fig. 2.3. Square Lattice For two-dimensional square lattice, z = 4, and 7., = —:-[cos(q,6) + cos(q,,6)] , (A14) The spectral density can be written as an elliptic integral and leads to a density of states p(:r) given by p(m) = «£sz 2(4 — 3)], (.415) 2 for 0 < a: < 4. Here [C(z) is the complete elliptic integral of the first kind. In general the one dimensional integral Eq. (A5) must be done numerically. However if I put K p = 0, the integral can be done analytically and I obtain 2 l 4KT/K W K = —- A16 ( ) 1r 1+4KT/K [1+4KT/Ki’ ( ) If KT >> K, the result Eq. (A16) becomes K W(K) = ln(32KT/K), (A17) 47rK T 59 For Kp- = O I also have the result from Eqs. (A7) and (A16) 2 1 4KT/K W' K = -— - A18 1( ) 7r 1+8KT/K [1+4KT/KJ’ ( ) where 8(a)) is the complete elliptic integral of the second kind. Triangular Lattice For the two—dimensional triangular lattice, z = 6 and 1 yq : 3 [cos(2aq,6) + 2 cos(a.q,6) cos(bq,,6)] , (A19) where a. = % and b = 3?. Carrying out the integral for the spectral weight, I obtain x w (2:) - 1 f/dmd 1 (A20) 9 — 7r2 ya: - z[1 — g-(cosZz + 2 case: cosy)] ’ o 0 leading to36 p(:1:) : 5.1—(m _ 3)-1/4IC(./1‘:'s"), (A21) 11' where —— — 3 — I s : [x/a: 3 1] [V13 3+3]. (A22) 16\/a: — 3 In general, the Watson integral Eq. (A5) must be done numerically over the spectral weight. If K p = 0, it can be done analytically to give K G x K -1—G(9 + K/KT)”4 , (A23) WK 2 ( ) 47rKT 2 60 where G = 8[\/§+—K77{—T— 1] W2 [W+ 3T1”. (A24) If KT >> K, the result Eq. (A24) becomes K IV K = —————- l ( ) 47rKT\/§ 71(7sz /K). (A25) The quantity ’1(K) can be found by differentiation. List of References for Chapter 2 10 11 12 13 61 List of References for Chapter 2 F.S. Galasso, Structure and Properties of Inorganic Solids (Pergammon, New York, 1970); J .C. Phillips, Bonds and Bands in Semiconductors (Academic, New York, 1973). L. Vegard, Z. Phys. 5, 17 (1921); C.Y. Fong, W. Weber, and J.C. Phillips, Phys. Rev. B 14, 5387 (1976). M.S. Dresselhaus and G. Dresselhaus, Adv. Phys. 30, 139 (1981) ; Inter- calation in Layered Materials, edited by M.S. Dresselhaus (Plenum, New York, 1987); M.S. Dresselhaus, Phys. Today (March 1984), page 60. S.A. Safran, in Solid State Physics, edited by F. Seitz, D. Turnbull, and H. Ehren- reich (Academic, New York, 1987). W. Jin and S.D. Mahanti, Phys. Rev. B 37, 8647 (1988). M.F. Thorpe, W. Jin, and S.D. Mahanti, Phys. Rev. B (in press, 1989). J.R. Dahn, D.C. Dahn, and R.R. Haering, Solid State Commun. 42, 179 (1982). J.E. Fischer and H.J. Kim, Phys. Rev. B 35, 3295 (1987). S.A. Solin, Adv. Chem. Phys. 49, 455 (1982); S.A. Solin and H. Zabel, Adv. Phys. 37, 87 (1988). RC. Chow and H. Zabel, Synth. Met. 7, 243 (1983); Comments Condens. Matter Phys. 12, 225 (1986); Phys. Rev. B 38, 12837 (1988). D. Medjahed, R. Merlin, and Roy Clarke, Phys. Rev. B, 36, 9345 (1987). H. Kim, W. Jin, S. Lee, P. Zhou, T.M. Pinnavaia, S.D. Mahanti, and S.A. Solin, Phys. Rev. Lett. 60, 877 (1988). W. Jin, S.D. Mahanti, S.A. Solin, and H.C. Gupta, in Microstructure and Prop- erties of Catalysts, edited by M.M.J. Treacy, J .M. Thomas, and J .M. White (Materials Research Society, Pittsburgh, 1987), p. 231. 14 15 16 17 18 19 20 21 22 23 62 S. Lee, S.A. Solin, W. Jin, and S.D. Mahanti, in Graphite Intercalation Com- pounds: Science and Applications, edited by M. Endo, M.S. Dresselhaus, and G. Dresselhaus (Materials Research Society, Pittsburgh, 1988), p. 41. M.F. Thorpe, Phys. Rev. B 39, 10370 (1989). S. Lee, H. Miyazaki, S.D. Mahanti, and S.A. Solin, Phys. Rev. Lett. 62, 3066 (1989) K. Komatsu, J. Phys. Soc. Japan, 6, 438 (1951). See e.g., Theory of the Inhomogeneous Electron Gas, edited by N .H. March and S. Lundqvist (Plenum, New Yprk, 1983); W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965); L. Hedin and B.J. Lundqvist, J. Phys. C4, 2064 (1971); C.T. Chan, D. Vanderbilt, and S.G. Louie, Phys. Rev. B 33, 2455 (1986); N.A.W. Holzwarth, S.G. Louie, and S. Rabii, Phys. Rev. B 26, 5382 (1982) D.P. DiVincenzo, E.J. Mele, and N.A.W. Holzwarth, Phys. Rev. B 27, 2458 (1983); D.P. DiVincenzo and E.J. Mele, Phys. Rev. B 29, 1685 (1984). See: E.H. Lieb and D.C. Mattis, Mathematical Physics in One Dimension: Ez- actly Soluble Models of Interacting Particles (Academic, New York, 1966); see also Physics in One Dimension, edited by J. Bernasconi and T. Schneider (Springer-Verlag, 1981). EN. Economou and M.H. Cohen, Phys. Rev. B 5, 2931(1972); D.C. Lic- ciardello and EN. Economoum, Phys. Rev. B 11, 3697 (1975); C. Papatri- antafillou, Phys. Rev. B 7, 5386 (1973); K. Ishii, Prog. Theor. Phys., Suppl., 53, 77 (1973); T.M. Nieuwenhuizen, Physica A 120, 468 (1983); G. Czycholl and B. Kramer, Solid State Commun. 32, 945 (1979); Z. Phys. B 39, 193 (1980), ibid. 43, 5 (1981). S. Alexander, J. Bernasconi, W.R. Schneider, and R. Orbach, Rev. Mod. Phys. 53, 175 (1981); S. Alexander and R. Orbach, Physica B 107, 675 (1981); RA. Lee and T.V. Ramakrishnan, Rev. Mod. Phys. 57, 287 (1985). F.J. Dyson, Phys. Rev. 92, 1331 (1953); H. Schmidt, Phys. Rev. 105, 425(1957); J.E. Hirsch and T.P. Eggarter, Phys. Rev. B 14, 2433 (1976); ibid. 15, 799 (1977); F.J. Wegner, Z. Phys. B 22, 273(1975). 24 25 26 27 28 29 30 31 32 33 63 B.I. Halperin, Phys. Rev. 139, A104 (1965); Adv. Chem. Phys. 13, 123 (1967); B.I. Halperin and M. Lax, Phys. Rev. 148, 722 (1966), ibid., 153, 802 (1967); I.M. Lifshits, S.A. Gredeskul, and L.A. Pastur, Introduction to the Theory of Disordered Systems (Wiley, New York, 1988), chapters 1 and 2. G. Garboczi, Ph.D thesis, Michigan State University (1985); G. Garboczi and M.F. Thorpe, Phys. Rev. B 31, 7276 (1985); ibid. 33, 3289 (1986). R.P. Feynman, Phys. Rev. 56, 340 (1939). R. Fletcher, Practical methods of Optimization (Wiley, New York, 1980); J. Stoer and R. Bulirsch, Introduction to Numerical Analysis (Springer-Verlag, New York, 1980), p.572. W.M. Press, B.P. Flannery, S.A. Teukolsky, and WT. Vetterling, Numerical Recipes (Cambridge University Press, 1986), p.301. 3. Kirkpatrick, G.D. Gelatt, and MP. Vecchi, Science 220, 671 (1983); S. Kirk- patrick, J. Stat. Phys. 34, 975 (1984); _ S. Kirkpatrick and G. Touleuse, J. Phys. (Pairs) 46, 1277 (1985). . B.P. Feuston, M. Ma, S.D. Mahanti, and R.K. Kalia, Phys. Rev. A 37, 902 (1988). R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985), ibid. 60, 204 (1988) ; Solid State Commun. 62, 403 (1987); D. Hohl, R.O. Jones, R. Car, and M. Parrinello, Chem. Phys. Lett. 139, 540 (1987); ibid. J. Chem. Phys. 89, 6823 (1988); P. Ballone, W. Andreoni, R. Car, and M. Parrinello, Phys. Rev. Lett. 60, 271 (1988). M.C. Payne, J.D. Joannopoulos, D.C. Allen, M.P. Teter, and D.M. Vander- bilt, Phys. Rev. Lett. 56, 2656 (1986); M. Needels, M.C. Payne, and J.D. Joannopoulos, Phys. Rev. Lett. 58, 1765 (1987). RA. Lee, P.H. Citrin, P. Eisenberger, and B.M. Kincaid, Rev. Mod. Phys. 53, 769 (1981); T.M. Hayes and J.B. Boyce, in Solid State Physics, edited by H. Ehrenreich, F. Seitz, and D. Turnball (Academic, New York, 1982), Vol. 37, p. 173. 34 J.C. Mikkelsen and J.B. Boyce, Phys. Rev. Lett. 49, 1412 (1982); Phys. Rev. B 24, 5999 (1981), ibid. 28, 7130 (1983). 64 35 H. Zabel, W.A. Karnitakahara, and RM. Nicklow, Phys. Rev. B 26, 5919 (1982); H. Zabel, Physica B 136, 1 (1986); H. Zabel, in Graphite Intercala- tion Compounds, Topics of Current Physics, edited by H. Zabel and S.A. Solin (Springer-Verlag, New York, 1989), in press. 36 T. Horiguchi, J. Math. Phys. 13, 1411 (1972). 65 Chapter 3 Dynamical Properties of Two-Dimensional Molecular-Solids 3. 1 Introduction Molecular solids are composed of molecular units which are clusters of atoms tightly bound to each other by strong covalent forces. These molecules have both center-of-mass (c.m.) translational degrees of freedom and orientational degrees of freedom. One of the important characteristics of molecular solids is the com— petition between direct and indirect (lattice mediated) intermolecular interactions which leads to different types of orientational ordering.“2 In this chapter, I inves- tigate the dynamic properties of a model two-dimensional molecular solid using self—consistent lattice dynamics theory and constant-pressure molecular dynamics simulation.3 The dynamic correlation functions and the dynamic structure factors have been investigated using the latter.4 Furthermore phonon and libron dispersion relations along different symmetry directions have been calculated and softening of phonons near structural phase transition (SPT) has been carefully studied!“5 66 When anharmonicity in solids is not weak, structural phase transitions can take place in which the solids change their crystallographic structure with changes in temperature.6 A common characteristic of the structural phase transition is sym- metry beraking,7 i.e., with decreasing temperature, the system goes from a high symmetry phase to a low symmetry phase at the transition temperature Tc. Struc- tural phase transitions can in general be classified into three categories;6 these are displacive, order-disorder, and ferroelastic. In the case of displacive SPT, the atoms or molecules in the ordered phase are displaced away from the high- temperature equilibrium positions. A typical example is the cubic-tetragonal tran- sition in strontium titanate (SrTi03, antiferrodistortive).8 The order-disorder type SPT has an ordering of some degrees of fredom which are disordered in the high-temperature phase. Typical examples are the phase transitions in N a0; (Ref. 9) and N (LN 02 (Ref. 10) where the appropriate degree of freedom is the orien— tation of the molecules. The unit cell does not change shape in these two types' of SPTs. Thus there is no strain associated with the pure displacive or the pure order-disorder SPTs. However, in the ferroelastic SPT, elastic deformation is as- sociated with the transition.‘*11 Most SPTs in solids are ususlly combinations of the above three.6 In displacive SPTs, there are anomalies in phonons of specific symmetry whose frequency decreases (or approaches zero) when the temperature approaches some characteristic temperature (not necessarly the transition temperature) from above.12 In ferroelastic SPTs, elastic constants of specific symmetry soften.1'13 Because of strong anharmonicity near the structural phase transition, theoretical understanding of the dynamic properties of molecular solids near the SPT is a rather difficult problem. In the past, mean—field type theories,6’7 renormalization 67 group method,14 and molecular dynamics computer simulations15 have been used to study the physical properties of solids near the structural phase transitions. The SPTs and dynamics of three-dimensional (3D) molecular solids have been widely studied.16 Typical systems are solid nitrogen (N2)," solid oxygen (02),” alkali cyanides (CN‘),19 alkali superoxides (02- )9 and alkali nitrates (NOE).10 These systems generally exhibit SPTs in conjunction with orientational order- disorder transitions where both translational and rotational degrees of freedom are involved.20 For example, in three dimensional K C N , at low temperatures, the solid has an orthorhombic structure in which all the CN‘ molecules are oriented in one direction. When temperature is raised, it transforms into a cubic structure in which the CN“ molecules are orientationally disordered.21 However, similar study of the dynamic properties in two-dimensional (2D) molecular solids have not been carried out systematically until recently?"5 In the remaining part of this section, I will discuss some of the basic features of the two-dimensional molecular solid,20 introduce our modelz'22 and sumarize the ferroelastic phase transition in this model system obtained from previous MD simulations by Tang, Mahanti and Kalia (TMK).22 Molecular overlayers adsorbed on a crystalline surface can form the two- dimensional (2D) molecular solid. Molecular nitrgen (N2) or oxygen (02) ph- ysisorbed on graphite surface serve as typical examples of such a 2D system.” In particular, physisorbed oxygen molecule on graphite exhibits a number of different phases depending on the coverage and temperature.24 There are extensive experi- 25 mental investigations of this system. X-ray diffraction, magnetic susceptibility,26 heat capacity,” neutron diffraction28 and low-energy electron diffraction (LEED)"9 68 measurements have discovered a diverse phase structure in this system. In particu- lar, there is a low—coverage low-temperature phase (6-phase) in which the molecular axes are collinear and parallel to the substrate surface. The oxygen lattice in this 6-phase is a centered rectangular (or, isosceles triangular, IT) structure which is incommensurate with the substrate graphite hexagon. This phase was predicted by Etters et al.30 by pattern search energy minimization calculation. The zero tem- perature energy minimization calculation3o wich ignored the substrate corrugation predicted an IT lattice structure for this low coverage 6-phase. Expermentally, this phase was first found by Heiney et al.25 and later careful LEED work29 has given a great deal of information about this phase. Finite temperature properties of the 6-phase have been studied using Monte Carlo simulations.31 The molecular dynamics simulation is another direct method to study the structure and phase transitions and in addition it gives information about the real time dynamics.15 In particular, the constant-pressure molecular dynamics is an excellent method to study physical properties of systems exhibiting structural change with temperature or pressure.” Structural phase transitions in a diatomic molecular monolayer system has been studied recently by TMK.22 The model in this study consists of homonuclear diatomic moleculars interacting through an atom-atom potential of Lennard-Jones (LJ) type, i.e., v0») = 4.5“?)12 - (96], (3-1) where r is the distance between two atoms. Since the internal virbrational fre- quency of each molecule is very high due to strong bonding between the two atoms, the molecules are assumed to be rigid.33 The potential parameters used in the TMK work are appropriate for the oxygen molecule,3°'34 i.e., e = 54.34ch and 69 0' = 3.0521 with internuclear distance d = 1.20824. For monolayer density in the neighborhood of 10 molecules per 100(A)”, this system closely resembles the 6- phase of oxygen molecular monolayer adsorbed on a graphite substrate. However, more accurate modelling of this 6-phase requires the inclusion of other interactions which may be rather small. For example, the oxygen molecule has a nozero spin (5 = 1) which leads to a magnetic exchange interaction between two molecules.35 This interaction is responsible for the low temperature (below 11K ) antiferromag- netic structure in the e-phase of the 02 on graphite system.36 The corrugation of the substrate potential is nonzero but small (about 5K )5" In addition, the oxy- gen molecule has a nonzero but small electric quadrupole moment which leads to quadrupole-quadrupole interaction between two molecules. All these three small interactions have been neglected in the present study. Also the out-of—plane mo- tions of molecules have been neglected, thus all molecules are constrained to be parallel to the substrate. The ground state of this system is an IT lattice with ferroelastic (FE) molec- ular ordering.5'22 The lattice constants are a. = 3.332A and b = 8.054A; molecules are all parallel to the y-axis [see Fig. 3.1(a)]. The surface coverage at the ground state is n = 1.172po; where p0 is the coverage of the J3 x J3 struc- ture on the graphite surface which corresponds to 0.0636 molecules per (A)2. 29'” it is found that the ground state of the low-density 6-phase Experimentally, of oxygen monolayer on the graphite has an IT lattice structure with a. = 3.25.21 and b = 7.98:1. The oxygen molecules lie flat on the substrate surface and are orientationally ordered along the b-axis. In contrast to oxygen molecules, another 2D diatomic molecular solid, the ground state of the nitrogen (N2) molecules ad- sorbed on graphite,39‘41 has herringbone (HB) ordering with ET center of mass 70 lattice structure which is commensurate with the graphite hexagon.29'42 However, for the oxygen system, the energy for the HB ordering is higher than the energy of the FE—IT structure.5 In the previous constant-pressure MD simulation study by TMK,22 it was found that the oxygen system undergoes a first order phase transition from an ori- entationally ordered ferroelastic phase to an orientationally disordered paraelastic (or plastic) phase at temperature Tc 2 20.1K. This order-disorder transition was accompained by a lattice structure change where the center-of-mass of the molecules transfrom from an IT structure to a equilateral triangular (ET) struc- ture. This result is in reasonable agreement with LEED experiment.38 In the present thesis, I have reinvestigated the above phase transition and have in par- ticular explored the dynamic correlations as a function of temperature which was not done by TMK. In the course of this work, I have been able to generalize the expression for the internal stress tensor first derived by Tang to the case of polyatomic molecules and have correctly. incorporated the kinetic contribution to the stress tensor.43 The remaining part of this chapter is organized as follows. In section 3.2, I introduce a new procedure43 to obtain the internal stress tensor in the constant- pressure molecular dynamics simulations for anisotropic molecular solids. In par- ticular, I emphasize the importance of incorporating the rotational contributions to the internal stress tensor correctly. In section 3.3, I give static structure factor obtained from the MD simulations. In section 3.4, I use the self-consistent lattice dynamics to study the phonons and librons in the low temperature ferroelastic phase. In particular, I emphasize the role of the translational-rotational coupling. In sections 3.5 and 3,6, I discuss the dynamic density-density correlation function 71 and associated dynamical structure factor obtained in the MD simulation. In ad— dition, I present phonon dispersion curves for temperatures both below and above the ferroelastic transition. I will also compare these MD results with those ob- tained by the self-consistent lattice dynamics calcualation (Sec. 3.4) and address the question of phonon softening near the structural transition. In section 3.7, I give a brief discussion and a summary of the main results. (I) (b) ’ \ I \ I, \ I \ / \ I \ ’ \ I \ ’ \ I \ I: \ i b \ \ I ‘ I, \ ’ ‘ I \ I \ ~ I \ I I \ I v t—a—t L. Fig. 3.1 The miniminum energy configurations of a two-dimensional diatomic molecular solid with the ferroelastic molecular ordering (a) or the herringbone ordering (b). The ground state of the 02 /graphite system has the ferroelastic molecular ordering with an isosceles triangular lattice structure. 72 73 3.2 Internal Stress Tensor of Anisotropic Molecular Solids In recent years, the constant—pressure (in general, constant external stress) molecular dynamics has been extremely useful to study structural phase transitions in solids."4 Such studies in molecular solids have elucidated the underlying physics of ferroelastic phase transitions which are usually accompanied by an orientational 4'5'20'22'45 In constant pressure simulations,”‘48 the order—disorder transition. periodically repeating molecular dynamics cell is assigned a fictitious mass and the volume and the shape of the MD cell are allowed to change, this change being determined by the internally generated stress tensor and the externally applied pressure (stress). For our diatomic molecular system, the orientational dynamics is expected to play a crucial role in both the structural phase transition involving a shape change and the dynamic properties of the system. To properly account for the rotational motion, in this section, I develop a new procedure43 to obtain the internal stress tensor in the constant—pressure molecular dynamics simulations of molecular solids. The results presented in this section is quite general and can be applied to any polyatomic molecular system. In most of the early works on molecular solids,48 the internal stress tensor is determined by the positions and the velocities of the molecular center—of-mass. The molecular orientations only appear indirectly in the calculation of the inter- molecular forces. It was first shown by Tang49 that orientational degrees of freedom have significant contributions to the elements of the stress tensor. I will review 74 the main idea, generalize to polyatomic molecules and properly take into account the effect of rotational velocities. In this section, contributions to the stress tensor Puv which depend explicitly on the molecular orientations and angular velocities are derived and their significance in the orientationally ordered phase is discussed. In the constant-pressure (isoenthalpic—isobaric) molecular dynamics of monoatomic system,“*‘" the essential idea behind the calculation of the internal stress tensor that determines the dynamics of the MD cell is to start from the the Lagrangian [.1 for the system given by N [:1 2‘ $27111? — Z Z V(I‘,'j) , (3.2) and scale the coordinate of the particles 1‘,- by the vectors a,b in 2D, and a,b,c in 3D of the parallelogram (parallelepiped) MD cell, i.e., r, : hs, , (3.3) where h is a 2 x 2 (in 2D) or 3 x 3 (in 3D) transformation matrix given by h 2 (a,b) , or (a,b,c) . (3.4) All the particles are located inside this cell and the cell is repeated in space by periodic boundary conditions. The vectors a, b and c are allowed to change in the simulation. Both the kinetic energy and the potential energy terms in L, are now expressed in terms of the scaled coordinates and velocities s,- and s,- respectively and one usually neglects a term proportional to h in the kinetic energy part.47 The dynamics of the variables 5,- and the matrix h are determined by the following Parrinello-Rahman Lagrangian“ 75 Zm,s,Gs, Z Z V[h(s ,—s,-)] —peQ + éWTrh‘Lh , (3.5) " J'(>i) where the gauge matrix G: h+h (h+ is transpose of the matrix h),fl — —det(h) is the area of the MD cell, pc is the externally applied (constant) pressure, and W is the mass associated with the MD cell. The equations of motion (EOM’s) for h which determine the dynamics of the MD cell can be obtained from the Lagrangian equation Why” = (535:) which leads to Ith (P— peI)A, (3.6) where I is the identity matrix and A:fl&fl”. an In Eq. (3.6), P is the internally generated instantaneous stress tensor. For monoatomic systems, P is given by47 PM, =% gmrh, + Z 2 F5. (r.- 4,)" , (3.8) i J°(>i) where 11,11 = x,y,z, and F), = —8V(r,'j)/8r,-j is the force between atoms i and j. Now I consider a rigid polyatomic molecular system described by the La- grangian = 122%: Z mag. — Z Z. 2 new). (3.9) 76 where i, j are the molecular and lc,l are the atomic indices, and Pik,jl = r11 -— r5), (see Fig. 3.2). There are two ways of obtaining the internal stress tersor. In the first one, coordinates of each atom can be separated as Pu. = R.‘ + Pile, (3.10) where R, is the center-of-mass coordinate and pi), is the relative coordinate. Then the conventional procedure48 to handle the MD cell dynamics is to scale only the c.m. coordinates R, by h and treat p5,, as a constant. In this case the internal stress tensor in Eq. (3.6) is given by 1 . 'v ' y Pm, : 6 ZMRQ‘R, + Z Z)F{; (11,-— 11,-) , (3.11) i i j(>:' where III 2 Zk mi), and Fij = Z Z Fik,j1', (3.12a) l: l _ aV(r,-,.,,-,) 3.125 31‘;sz ( ) Firm = is the force between atoms ilc and jl. Note that R3 and hence Pm, do depend implicitly on the molecular orientations, although only the center—of—mass coordi- nates R,- and their velocities R; appear in the expression for Pm, explicitly. This approach has been followed by many researches in the past. A second procedure for calculating the internal stress tensor which I discuss here43 is to scale the position vectors (r,,,) of individual atoms of the molecules by h and then use the rigid molecule conditions through the introduction of fictitious time dependent forces. There will be additional terms in PM comming from the 77 kinetic energy term of the Lagrangian and these terms will depend upon the an~ gular velocities («15). The potential energy contribution to Pm, now depends upon the forces acting on the individual atoms of each molecules. In addition to the extrenal forces, one has to take into account the constraint forces acting on the atoms. The constraint forces fix the internuclear distance of the molecules. In systems where intra-atomic vibrations are important, one will have to replace the constraint forces by the actual forces acting between the atoms of a given molecule. The constraint force, which in general depends on the angular velocities, can be obtained from the equation—of—motion of the atomic coordinates. The dynamics of the MD system is described by the Lagrangian43 5M: — :2 Z mgksgthik — Z Z Z V(r,-;c j;)— pcfl +é-WT1-hl'h (3. 13) ik j( >i) I From the Lagrangian equation for h with EM and adding the constraint force, I obtain the same equation (3.6) for the dynamics of the MD cell but with the internal stress tensor given by Pun/:02: mile (1150:)” (hsikly +QZZZF“ ikjl rilc,jl+flZf ikpilc’ ilc j(>i) l (3.14) where fa, is the constraint force acting on the atom ilc with the condition 2,: fa, = 0. Through a long algebraic manipulation, Eq. (3.14) can be written in a more transparent form: Puv =P£T+ 512 mu. puma. +512”? +5.) pa. (3-15) where P3," is the stress tensor given in Eq. (3.11); F”, = 21123021 Fuhjg is the force acting on the atom ilc. The stress tensor given in Eq. (3.15) has translational 78 and rotational contributions from both the kinetic (kin) and potential (pot) energy parts of LM- Equation (3.15) is the centeral result of this section. I can formally write P 2 P"‘" + PP“, (3.16a) with Plain : Pkin,cm + Phinmot, (3160) Ppot 2 Ppot,cm + Ppotmot, (316C) and Plein,cm + Ppot,cm ___ Pcm . (316d) To illustrate the significence of the above results in MD simulations, I apply the above results to the system of homonuclear (rigid) diatomic molecules which are confined to move in a two dimensional plane with their rotational motion confined to the same plane (see Fig. 3.2).43 The rigid rotor condition yields I311. = %d define. (3-17) 1 - 1 .. . C] = — i2 = —[1Md03 + §(Fi1 — Fi2)'l‘11]1‘i1, (3-18) $111.15,- = F(0,-), (3.19) where F(0,~) = (Fgl — Fig) - {1,1 is the total force acting on the molecule i in the direction perpendicular to the internuclear axis. Eq. (3.19) is simply the equation- of—motion for the orientational variable 0;. From Eqs. (3.15), (3.17) and (3.18), I obtain Fig. 3.2 A diatomic molecule confined to the XY— plane. 79 80 cm 1 ' .. Au .. ~12 1 1 . ..., Piw 2 PM, + Il- 2;: 19311571,, " rim-1) + II 2‘: §dF(9i) "$117311 (3°20) where I 2: AI dz/ 4 is the moment of inertia of the molecule. Using 13,1 2: (c089,, sindg), and fig = (—sin0,-,cos€,-), I find that the rotational kinetic energy contri- bution to the stress tensor [the second sum on the right-hand-side of Eq. (3.20)] is - 1 - —c0820- —sin20' kin,rot _ _ 2 i i P _ fl 210‘ (—sin20,- c0820,- ) ’ (3'21) and the rotational potential energy contribation [the third sum on the right-hand- side of Eq. (3.20)] is 1 1 —sin29- cosZU- — 1 pot,rot _ __ _ . 1 g P _ Q 2: 4dF(0‘) (c0320,- +- 1 sin29; ) ' (3'22) 3 The internal pressure is defined by p‘ = TrP. From Eqs (3.21) and (3.22), it is interesting to note that p depends only on Pc'", i.e, p = TrPc'". However, the antisymmetric component A = %(P,,,, — Pyz) has important 0g-dependent terms, i.e., 1 cm cm 1 1 A —_- 5(1),, — py, )— ,2— Z EdF(0,—). (3.23) We have found5 that the last term is extremely important in MD simulation as it makes a significent contribution towards keeping the total angular momentum conserved. Next I consider the strength and the temperature dependence of the rotational contributions to the stress tensor in Eq. (3.20). To do so, I calculate the ensemble 81 average of (P —- Pc’") and for simplicity I consider the low temperature ferroelastic phase22 (see Fig. 3.2). As will be shown later in this chapter, at low temperatures, the rotational motion in this phase is deminated by the libron spectrum. The rotational harniltonian can be written as 1 ., 1 - - H, = 52:16,. + 2;D99(R‘ — R,)9,0,, (3.24) with d,- = 0,- — 90°. Then I can calculate the following averages (63) = (kBTVm. <6?) = 3(kBT)’(gz +gf). (3.25) and (160,) = ~k3T, (Iii?) = —6(kBT)zgl, (3.26) where g" = / (2%:1);[D99(q)]_n, (n = 1,2,3,...), (3.27) and 09.01) = Z enema), R. where the integral being over the first Brillioun zone. Using Eqs (3.21), (3.22), (3.25) and (3.26), I obtain the thermal-averaged values of Pkin'r” and PP°"'°‘ as (P*"""°‘) = (n)kBT[1 — 2g1k3T + 2(g2 + 990.3102 + ...] 2,, (3.28) and (PP°""°‘) = (n)chT[—1 + 49,1237” + ...] 23,, (3.29) with 2;, z ((1) 31), (3.30) 82 where (n) is the average molecular surface density. As can be seen from Eqs. (3.28) and (3.29), the off-diagonal contributions to the stress tensor coming from the rotational part vanish in the harmonic theory. This is also true if one includes the effects of linear rotation—translation coupling on the libron dynamics. How- ever, near the ferroelastic structural phase transition where the harmonic theory breaks down, these off-diagonal terms are found to be significent in our molecular— dynamics simulations.5’22 Comming back to the diagonal terms, I give in Fig. 3.3 the temperature dependence of the axe—components of the internal stress tensor. It is interesting to note that the linear term in T of (PH - P3,") is identically zero due to a cancellation between the kinetic and potential contributions and the net result is a quadratic (T2) increase. The increase in Pg" and the corresponding decrease in Pg," with increasing T make physical sense since in the ordered state the molecules are orientated along the y—direction at T 2 0. In addition to the MD cell dynamics, the MD equation for the scaled coordi- nates s,- are d , + mE(Gs,-) = h 2: F,,-, (3.3171) J¢t i.e., m§,- = h-1 2 19,,- — rue—lea. (3.31b) #1 The instantaneous temperature in defined by the kinetic energy, i.e., 3 1 N 1 __ _ . , . __ 2 '? ~2-NkBT(t) _. 2 ; m(v, v, + 4d 6,), (3.32) where v; is the molecular c.m. velocity adjusted to make the total momentum and angular momentum zero.5 In MD simulation, this instantaneous temperature 83 should then be adjusted again to the required temperature T by scaling the veloc- ities and the angular velocities by a factor [T/T(t)]l/2. In summary, I have shown that both the center—of—mass and the rotational contributions to the stress tensor should be considered in the molecular dynamics simulations of molecular systems. In particular, the latter can be quite important near a structural phase transition involving orientational degrees of freedom. F l r l I. 4 r —Kinetic part o,' .. ----Potential part 57 ....... TOtal .3/ /l/(831ll Fig. 3.3 Temperature dependence of the thermal averaged rotational contri- butions to the internal stress tensor P3" in the ferroelastic phase. T is mea- sured in units of 891k3, where g1 has the dimension of inverse energy. 84 85 3.3 Static Structure Factor Before proceeding to the dynamical study, I discuss the static structure of two- dimensional molecular solids in this section. I first give a brief description of the MD simulation precudures. Then I discuss new results beyond those obtained by TMK of the static structure factor, namely, the radial distribution fucntion (RDF) both below and above the ferroelastic phase transition. The results of dynamic structure factor will be given in section 3.5. The equations of motion for a collection of rigid rotors given by Eqs. (3.6), (3.19) and (3.31b) with the internal stress tensor given by Eqs. (3.11) and (3.28) are solved numerically by a fifth-order predictor-corrector algorithm.49 The potential parameters I use here are the same as those used in TKM’s work.22 A cutoff distance of 50' is used for the LJ potential, which is modified so that the potential and the force are continuous at the cutoff. The periodic boundary condition is applied to the MD cell. All the important parameters involved in the simulation are listed in the Table 3.1. I have carried out a series of simulationss'50 starting from the ground state and heating the system slowly. The c.m. positions, velocities, accelerations, ori- entations, and anuglar velocities etc. are updated with an integration time mesh 2 r = 0.01102 = 0.02 x 10*”.sec. In calculating the dynamic correlation functions, the data points are ”coarse- grained” by recording the data each 5 MD time stpes to avoid possibile statistical correlations between them. At each reduced temperature T’, 9000 to 10000 MD time steps are used to equilibrate the system and another 5000 to 20000 MD steps Table 3.1 Parameters used in the MD simulation Time unit MD time step cutoff N Pe parameters value Length unit a' 3.05 A Mass unit m 15.99 a.u. Temperature unit e/ka = 54.34 ,/-'"f—’ = 1.815 x 10-13 sec 1' = 0.02 x 10‘12 sec = 0.01102 50' 400 0.0 ma3 3 86 87 are monitored to obtain the dynamic variables. The temperature was adjusted by scaling both the linear and angular velocities. Near the structural phase transi- tion, averages of physical quantities are obtained from constant-volume MD after the constant-pressure MD equilibrium run. This was done to reduce the large fluctuations occuring near the structural phase transition. Other details of the MD simulation are similar to those of TMK’s work. All our MD simulations were performed on a Cyber 170 Model 750 computer to produce the raw data and on a VAX 8650 computer for extensive data analysis. From this simulation, the order- disorder transition temperature was found to be T: = 0.37 (20.1K), in agreement with TMK’s previous results.22 The melting transition from the orientationally disordered plastic phase to an isotropic 2D liquid phase at T; = 0.70 (38.0K) was also confirmed. Detail results of thermodynamic quantities have been given in Ref. 5 and will not be repeated here. The radial-distribution function g(r) is defined by (n(r2, r2 + Ar2)) “IrAr2 ’ 90‘) = (3.33) where (11(1‘2,r2 + Ar2)) is the average number of molecules in a ring with area «Ar2 and at a distance r from a given molecule. Results of g(r) obtained from the MD simulation and corresponding plots of molecular configurations are shown in figures 3.4-3.10. Some characteristic features of these figures are discussed below.50 Fig. 3.4 shows g(r) for temperature T“ = 0.36 which is just below the fer- roelastic transition. In this figure, the numbers under the smooth curve indicate the total numbers of molecules in the successive coordination shells obtained by integrating g(r) and the numbers marked with arrows give the ideal distribution of 88 molecules in the IT lattice structure. One can clearly see two rather sharp peaks at small r due to the two nearest neighbor molecules and four second neighbor molecules. The peaks associated with farther neighbors are slightly broadened but well defined. Fig. 3.5 gives a snapshot of molecular configuration obtained from the MD simulation at this temperature. One sees that the molecules are ordered more or less along the y-axis. Fig. 3.6 shows g(r) at T" = 0.40 when the system is already in the orien- tationally disordered phase. In this figure, the numbers marked with arrows give the ideal distribution of molecules in the ET lattice structure. In this case g(r) has a single peak at small r corresponding to the six nearest neighbors in the ET structure. The peaks associated with farther neighbors are considerably broadened indicating a rather short correlation length. Fig. 3.7 gives a snapshot of molec— ular configuration obtained from MD simulation at this temperature. One sees clearly that the molecules are oriented randomly but there is some indication of short range orientational correlations. Fig. 3.8 shows g(r) for three temperatures T‘ = 0.36,0.38 and 0.70. Note that there is a gap [g(r) 2 0] between the first six neighbors and the other neighbors in the solid phases. However, at T" = 0.70 there is no such gap and g(r) shows a liquid-like structure. In Fig. 3.9, I give g(r) and a configuration of the system quenched rapidly (in 500 MD time steps) from T“ = 0.38 which is just above the ferroelastic phase transition temperature to T" = 0.01. The radial distribution function shows a three-sharp-peak structure for small 7‘. These three peaks are narrow due to the absence of thermal vibrations at this low temperature. Among the first three peaks, the middle one reflects the local ET structure while the other two reflect the local IT structure. This can be seen from the insert where both local HB and FE orderings can be seen. 89 1" = 0.36 G” <-& fi-N é—N <—A é—& $10 <—Jb at?) n I l L j . 2 . . . o 1 1 1 1 1 , 0.0 2.5 5.0 7.5 10.0 12.5 ‘I 5.0 r2/a2 Fig. 3.4 Radial distribution function at T" = 0.36 obtained from the MD simulation. Tit-10.36 .. ~'"‘--'-33MMM¢MM §¢+11313¢¢¢¢¢¢t+¢¢§- MMtMM,,,¢¢¢M¢M MtMMMMMquww MMMMMMMssMH thtttttthtttttttttt MasMMMMHMMM 5+1¢¢¢t¢¢¢¢¢¢ttthw MMMMMstiMM tittttttlhtttitlswtt tMMMuMMMtMM MHMMMMMMMM MMMMMMMMsaM httttttttlttttttttht tttttlttttttttthttt tittitttitttttttrttt ihtttftfittttttttttt (HMHMMMMMM MilttMHMtnMMM 999392”: .1 ....3321211 t Fig. 3.5 A snapshot of molecules configuration at T‘ = 0.36 (below the ori- entational order-disorder transition ) obtained from the MD simulation. 90 91') n T «a «a <— ‘ .— 0 12 o L l L 0.0 2.5 5.0 7.5 10.0 12.5 15.0 Fig. 3.6 Radial distribution function at T‘ = 0.40 obtained from the MD simulation. 92 T*=o.4o Elsa-93245? ‘?1+‘S~‘9~33’3+ gaspaq颢zsaua¢+¢a¢ z¢K¢zzz9¢b~33éabbz¢1 $3g2¢¥>zz3¢2$2¢zgfla+ 1+¢f¢pz$ngcpgzsfin~agtzz ¢X¢¢R+$+¢Q¢guzzzfl¢f 39¢aqz’i’tzgfigs‘e'h‘82'2rd? ¢a¢+¢2¢¢¢‘3¢99¢z¢az ”finkfllitsfitflsit‘stazz‘e” fig$¢¢hfizqtzfi¢¢i12$2+¢ X‘a‘aszei‘Qths‘Q‘a-fis‘sl‘ ¢¢¢3¢eai$$¢s~¢zzgxlii fni¢¢z*?3‘asnszuf9$9 bb‘i‘ahq‘fi‘s‘ufiza'zdifi‘9‘ Q¢¢~al¢¢azs+ggb¢¢31+~s ¢¢tgaa¢.¢&aa.§g_,zta¢, egg¢31$+ia¢+zzztézsz tea’QN‘U’elavsi-‘evfissa +Q¢9l$b§¢‘5¢¢—?¢Qzfl Z '9-'9~.~ Fig. 3.7 A snapshot of molecules configuration at T" = 0.40 (above the ori- entational order-disorder transition) obtained from the MD simulation. 93 56¢ pay—5.39 Baggage“... our: .8.‘ 28.32:... acaaaEEmp 153% o. a firm .«8 :5: n. 2833353 Q: 22 mu ON m-. 0.. 1| . . _ . .\l. lolllli/ \I ‘\\ . \\ III ’Io/ \ x1... \ ..- . .. . \ /. /..| /. r- Ohd HIP lllll and "L. 111.1! 1 00.6 "L. i _ _ _ _ _ _ _ P . F N a uoflaund "01100111910 IBIPBH rI'TW'T I ll z :l¢’\¢“$.fi‘i4" i“ ‘ I A“ “‘o"§§i‘.b~§‘§a‘-‘.*§"m _ III...§ ‘t\§§‘°‘*?§’o“ ‘ *I’ O . 5 3‘ #’,'A¢j""a’¢:‘¢‘*°§‘ ‘ 'o-on°x “it“siw‘i‘rt’f- __ "Jffit'd‘tt' "‘ Q I t 1‘ x 6' 4" N . n I '0 D K N N’ v u x°,:‘¢¢,fl¢*°< “§ 5 1‘. 51° §$‘- ‘ j 3‘ ’ I .; ’3'\¢ «1‘0‘.°;t{‘111113 ‘¢\\\ §§§“'-\‘2l§.$‘- Radial Distribution Function 5 10 15 20 25 r2 (unit 02) Fig. 3.9 Radial configuration and the molecular configuration (insert) for a state quenched to T’ = 0.01 fi’om T‘ = 0.38 in 10 psec. 95 3.4 Self-Consistent Phonon and Translational-Rotational Coupling in the Ferroelastic Phase In molecular solids, due to the presence of the translation-rotation coupling, phonons and librons are combinations of two motions: molecular c.m. trans- lational vibration and molecular orientational libration.20 In this section I use the self-consistent lattice dynamics theory50 to explore the dynamics of the two- dimensional molecular solids at different temperatures below the ferroelastic tran- sition where the molecules undergo small amplitude translational and librational motion. I consider small displacements of molecular c.m. and small angle devia- tions of molecular orientations from their ordered state. The coupling between the c.m. translational motion and the librational motion is included in the leading or- der. As will be shown later, the theory works reasonably well up to the ferroelastic transition temperature if proper allowance of the lattice expansion is made. But in the high temperature orientationally disordered phase where molecules undergo orienational diffusion and large amplitude displacements, the anharmonicity and translational-rotational coupling are very strong. The weak coupling assumption is not adequate. A proper treatment in this case has to be done following the theory of Sahu and Mahanti developed for 3D molecular solids.1'13 In the self-consistent lattice dynamics,51‘53 the equilibrium positions of lat- tices are determined by the phonons and librons whose frequencies are in turn determined self-consistently to minimize the total free energy of a crystal. Thus thermal expansion of solids due to anharmonicity which can affect the phonon and 96 libron frequencies at nozero temperatures can be include in a SCLD theory. An- other important effect of the anharmonicity in the potential is to give the phonons and librons a finite lifetime and to cause interactions between two or more phonons. The standard approach to handle this anharmonic effect is to use phonon thermal Green’s function and Dyson equation."'51 By summing irreducible diagrams one can calculate the phonon self-energy, whose real part is related to the frequency and imaginary part to the lifetime of the phonons. However in this section I make a pseudoharmonic approximation in which the effects of thermal expansion and renormalization due to the anharmonicity are included but the lifetime of phonons and librons are assumed to be infinite. This approximation is good when the an- harmonicity is not strong. Then the phonon and libron frequencies are simply related to the eigenvalues of the dynamic matrix. I combine the SCLD formalism and the MD simulation to take into account the effects of anharmonicity induced thermal expansion. In this approach the temperature dependent equilibrium lat- tice constants entering the elements of the dynamic matrix are taken from MD simulations. The total potential energy of the two-dimensional diatomic molecular system is given by 1 U = 5 Z: Z V(|r.-,c — rm). (3.34) th 1,1 Let the lattice constants at temperature T < Tc be a1 and £12 with |a1| :- a and a temp = b [see Fig. 3.10(a)]. I expand the potential energy U around this finite temperature equilibrium position with the ferroelastic ordering. Neglecting the anharmonic terms, I obtain 97 =52 2 DM 2' ,- )(up 0M1) (3.35) iij “V where [1,1/ = :c,y,0; uz(i),uy(i) are the c.m. displacements of molecule i and u9(i) = 0.- is the angle between the y-axis and the axis of the molecule i. The coupling matrix Dw(i, j) between molecules i and j is given by the second order derivate of U with respect to up(i) and uu( j ) The equation of motion for u“ are “may“ i): Z Z” Dun/(1:3) all (j) + ZDp9(i j) 0j: (’1' _ 3:31) (3'36) 1 ”=2 and —Ié,~ = ZD99(i,j)0 + Z 2” Dg“( i ,j )uu( (j) (3.37) Jfl=3 For nearest neighbors and next nearest neighbors (see the Table 3.2), the coupling matrix D(i, j) can be calculated explicitly and are given by50 a, 0 0 D(0,nn1) = 0 61 A] , (3.38) 0 —A1 91 31 771 —¢1 D(0,nn2) = 171 71 —£1 , (3.39) 1111 51 t1 51 —7h “-11% D(0,nn3) = —771 71 {1 , (3.40) 101 —€1 t1 C11 0 0 D(0,nn4) = 0 61 —/\1 , (3.41) 0 A1 91 98 31 771 191 D(0,nn5) = 771 ’71 £1 , (3.42) —i/J1 —51 ii 31 “771 191 D(0,nn6): ~11; ‘yl —§1 , (3.43) ‘d’l 51 t1 C22 0 A2 D(0,nnn1)= 0 62 0 , (3.44) —)\2 0 92 32 ‘71—2 $2 D(0,nnn2)= —7)2 72 —£2 , (3.45) —i/J2 52 t2 32 712 492 D(O,nnn3)= 712 72 ——§2 , (3.46) 192 {2 t2 (12 0 —A2 D(0,nnn4)= 0 62 0 , (3.47) A2 0 92 H2 —772 *1/22 D(U,nnn5)= —n2 7; £2 , (3.48) 192 —52 t2 .32 772 1(‘2 D(O,nnn6): 772 72 {2 , (3.49) ”(’2 -1=i2 152 where the parameters a1 etc. for the nearest neighbor coupling are given by a, = 2a2(A1 + A2) + 2(Bl + 8;), (3.50) 99 A A A 61=a2[—21+——‘: 5] +233+B4+Bs, (3.51) 61 = 2d2A2 + 2(Bl + B2), (3.52) A A A 7, = 1,2 [3i + —‘3—:—5—] )— 612(A4 + A5) — bd(A4 — A - 5) + 233 + .134 + 35, (3.53) A A + A 1 m = a,b[—2£ + 447.3] _ 53.164. — A5), (3.54) 1 2 2 1 2 91 — 5“ d (A1 " A2) + 5d (Bl — Bz), (3.55) 1 A4 + A5 1 B4 + BS tz—zdzA————— -d’B———, 3.56 1 80 [ 3 2 ] + 2 I: 3 2 J ( ) A] Z adzAz, (357) 1 2 1 ' . (1,, : gn .1(A4 . A5) — §d(B4 — 135), (3.58) 1 1 2 £1 = gabd(A4 — As) — Zad (A4 + A5), (3.59) and 012 etc. for the next nearest neighbor coupling are given by 62 = 236 + B7 + Ba, (3.60) fiz = 202(22‘19 + A10 + A11) 4- 239 + Blo + B“, (3.61) 12 = ibz(2A9+A10+A11)+d2(Alo+A11)+bd(A11—A10)+239+Blo+Bu, (3.62) 3 3 772 '3 352(2249 + A10 + An) + 50.6“.411 — A10), (3.63) 62 = b2(2A3 + A7 + A3) + d2(A7 + A8) + 2bd(A3 — A7) + 2.33 + B7 + Ba, (3.64) 9 1 ¢2 : Za2d(Al] " A10 + 2%“le — 310), (3.65) 3 3 £2 = gabd(A11 - A10) + 10410111 + A10), (3.66) 1 A2 2 §d(B8 — B7), (3.67) 100 1 2 92 = id (236 — B'r - Ba), (3.68) 9 1 t2 = {6036129149 —- A10 - .411) + Zd2(2Bg — BIO — Bu). (3.69) The parameters A and B are calculated from . 1 1 A(0k,zl) = EV'UO) - ;g—V"(ro), (3.70) B(0k,il) : _3_V'(ro), (3.71) 7‘0 where V’, V" are the first and second order derivatives of the LJ potential, r0 is the equilibrium distance betwen atom (0k) and atom (it). The values of 1‘0 and associated indices are listed in the Table 3.3. A1 to All and B; to Bu can be obtained from Eqs. (3.70) and (3.71) and r0 values from the Table 3.3. Note that the coupling matrix has the permutation symmetry Du,,(i, j) = Dw(j,i). In addition for a given D(i, j), the elements of this D matrix are symmetric in the pure translational part but antisymmetric in the part of translation-rotational coupling?“55 Also the self-interaction term can be obtained from the sum rule55 6 201 + 431 0 0 ' D(0,0) = — ZD(0,nni) = -— o 2.51 + 471 0 , (3.72) i=1 0 0 T where 1 T = —d2(a1 + 3’6) + 2d2(Bz + B4 + B5) — bd(B4 — B5), (3.73) and only the first six nearest neighbor interactions are include because the next nearest neighbor interactions fall off by a factor of 2'“. The elements of the dynamic matrix D(q) = Z Dane-"‘1'R R. for only nearest neighbor coupling is given by 101 1 1 Du(q) : —4[(%al + 31) —— 301 cos(q,,a) — )61 cos(§q,a) cos(%qyb)J , (3.74) . 1 . 1 D,y(q) : Dyz(q) = —4m szn(—2—q,a) szn(-2-qyb), (3.75) 1 1 I 1 Dyy(q) = —4 [(561 + 71) — 561 cos(q,,a.) — 71 cos(-2—q,a) cos(-2-qyb)] , (3.76) . 1 , 1 19.961) = D;.(q) = -4. 101 608(5920) 3‘"(§¢be)a (3.77) . . . 1 1 Dyg(q) = D;y(q) = 21. [A1 sm(q,,a) — 251 sm(§qza) cos(§qyb)], (3.78) I 1 D99(q) = T + 291 cos(q,a) + 4t; cos(§q,a) cos(§qyb). (3.79) Also D(q) is hermitian. The translation-rotational couplings are given by D,9(q) and Dyg(q). Note that the elements of the translation-translation partiand the rotation-rotation part of the dynamic matrix are real. However, the element of the translation-rotation coupling part of the dynamic matrix is pure imaginary. The eigenvalues and eigenvectors of the dynamic matrix D(q) give the phonon and libron frequencies and their displacement patterns, i.e., Duh) -.mwz D..(q) 0.9m) D2110!) Din/(‘1) ‘ "“92 Dy9(¢l) 2 0° (3°80) D;0(Q) D;g(Q) D99(q) — In)2 The unit vectors of a two dimensional IT lattice are given by [see Fig. (3.10a)] a] = as, a2 : 30?: + tangb 3)), (3'81) and the corresponding reciprocal lattice vectors are bl = b‘(—12-tan¢ :i: — in), 102 where b" = a—é—‘gfi; = $5" The first Brillioun zone (BZ) is shown in Fig. 3.10(b) along with high symmetry points. The phonon and libron dispersion curves for T" = 0 and 0.12 calculated from the secular determinant Eq. (3.80) are shown in Figs. 3.11-3.13. The finite temperature equilibrium lattice parameters a. and b used in the calculation are given by the MD simulation and listed in the Table 3.4. Figure 3.11 shows the LA and TA acoustic phonons and the libron dispersions associated with the ground state without the inclusion of translation-rotational coupling. Figure 3.12 gives the same dispersions in the presence of translation—rotation coupling. There are three branches. The coupling causes a splitting of the LA and TA phonon branches, particularly near the zone boundary N point where the TA phonon frequencies go down. In the regions where the phonon and libron branches cross, the modes are admixtures of the translational and rotational motions. Also note that the libron branch is quite dispersive. In Fig. 3.13, the LA phonon frequencies at T“ = 0.12 are given as the solid lines. The dots in this figure are the phonon frequencies obtained from the MD simulation (to be discussed later) which takes into account the thermal expansion and anharmoncity automatically. One sees that at low temperatures, the SCLD results are quite good. I will come back to this point again in the next section. I have also used the SCLD results to calculate the Debye temperature which can be written as“ 00 = mom/24'?)r {2(i>}—1/2, (3.83) 2 e v“ 103 where ve is the acoustic sound velocity, 11. is the density of molecules and f is total number of modes per molecule. The angular bracket in Eq. (3.83) refer to an average over directions in q-space. Low temperature heat capacity in the Debye model is given by 2 2 C = 24ng ((3)62) = 28.8498 (03) , (3.84) D where ((23) is the Riemann zeta function. If the motion of the molecules perpen- dicular to the graphite substrate is ignored, I obtain 00 z 87K. On the other hand, if the transverse out of plane motion of the molecules are included and it is assumed to have the same sound velocity as the inplane transverse phonon, then one obtains 00 z 78K. These give low temperature heat capacities which compare favorably with the experimental values."w Table 3.2 104 Coordinates of the first six nearest neighbor molecules and the next six neigh- bor molecules of the molecule at the origin. The integrals mhmz denote the compounts of R. along the basis vectors a, and a2. pair coordinates (m1, m3) 0-nnl ((1,0) (1, 0) 0-nn2 (g-a, §b) (o, 1) 0-nn3 (mg—.1, 46) (-1, 1) 0-nn4 (—a,0) (-1, 0) 0-n115 (—-%a, —%b) (0, -1) 0-11116 Ga, g-b) (1, -1) O-nnnl (0,b) (-1, 2) 0-nnn2 (—%a, :b) (-2, 1) O-nnn3 (-§a, —§b) (-1, -1) 0-nnn4 (0, -b) (1, -2) 0-nnn5 (go, -—%b) (2, -1) 0-nun6 (ga,%b) (1, 1) Table 3.3 Equilibrium atomic distances. index r0 1 a 2 W 3 51-11 sec¢ 4 x/(taPth-d)’ 5 (6a)” + 66 + «0' 6 b 7 lb -— dl 8 b + d 9 ((6a)? + (4b)? 10 (Age): + (p — .1)2 11 \/(ga)= + (§b + d)2 105 Lattice parameters a and 6 given by the RID simulation 3.4 T‘ a/a' b/a' 45 0 1.0924 2.6406 67.52 0.12 1.1710 1.2889 66.87 0.36 1.1308 ‘ 1.3249 66.89 0.40 1.3201 1.1409 59.95 0.50 1.3441 1.1487 59.67 106 107 m (b) Fig. 3.10 (a). Two dimensional (isosceles) triangular lattice with lattice vec- tors a; and .3, (a = '81.) (b). The first Brilliouin zone corresponding to (a) and the symmetry points. 108 .4 5.:5 10152 mm 255 1333.3: 433355052: 948 4.2:. 38:: A: 5 vanguaxo 9.3 338.» 93.: 2:. $5325 123.53cuéaaozflmzah 22 305.2. 5822 GROW 2: E merino commuommmu :95: ES 5.20:9 an.” .33 Z M L <9. . m: 109 55.508 icouaicub .293: ~53 3:33.93: 50.23.? 23%.:an :5 mm 82: .3 ~23: 34‘ .?=3.:a~m:ub 39:5 8m 8428.5 «Du 33 V4 2: 3.53 .RcocSa: 59:3 mm c233 Q 2: .3 H a 3k 222 892 .‘H 5.; @209: mm ....on 35595.: éggzaawoum 23.4 4.2:. @8398 5 ammmouQHo 0.3 32...?» 93>. ~45 $575.8 356338-}:333233 of Sr: beef QQDm. 2: 5 813:9 =ommuummmw :95: was =o=04m~ Nu.» .wmnu <9. .. M L OH ma (ZH at01) (*9 110 6.8.32.5..3 Q2 be 3:58 2: P3 gen .58.: 38 6 £32 2: 6: 6:... 36m .2 u :55 ....s 28%.: .... 3863: GE: coca—UPS .338, 953 vm L ON cm on: 333.» 9:3: 048 .Nnd H :u .a mezzo nomauumamw 5533 «5 fig.» .3h 2 «dud—qqiu—q111qqdd1—qufifi qqq-qnquq—qqdq—-qu—qq 1 1“”1‘ f f V U . m; u .... In" 111 3.5 Density Correaltion Function and Dynamic Structure Factor The dynamics of molecular crystals in the orientationally disordered phase is a rather difficult problem. Furthermore, experimentally, only limited data are available. Neutron scattering which is widely used to study phonon dynamics and dispersions in 3D molecular crystals56 has not been carried out extensively for 2D systems. From the theoretical side, one can study the phonon dynamics using the SCLD. However, this method is only appropriate for the low temperature orientationally ordered phase where the phonons and librons are well defind and weakely coupled excitations. MD simulation is a direct probe of the dynamical properties and in particular; it is quite useful in understanding the dynamics in the orientationally disorded phase as the effects of strong translation-rotation coupling and large anharmonicity are automatically included in these calculations.“-58 Molecular dynamice is perhaps the only viable method available at the present time to understand strongly coupled translational-rotational dynamics.59 Dynamical properties can be studied by the time-dependent correlation functions.60 I have calculated the density-density correlation function defined by F(q,t) =< Pq(t) P-q(0) >5 (3°85) where 1 ZN —iq - Mt) pq(t) = —\/—_: e 1 , (3.86) i=1 is the space Fourier transform of the density; q is wave vector, x;(t) is the c.m. position of ith molecule at time t. The <> in Eq. (3.85) refers to an ensemble 112 average (which equals to the time average in MD). In neutron scattering experi- ments, F (q,t) is not measured directly, but its time Fourier transform, S(q,w) is releated to the neutron scattering crossection,“l S(q,w)= / 113i“)t F(q,t)dt. (3.87) Since F (q,t) is an even function of time, S(q,w) is a real quantity. Longitudinal phonons show up directly in the function F(q,t). To obtain transverse phonons from F(q,t), a reciprocal lattice vector K which is perpendicular to q should be added to q. Other types of collective excitations such as librons will affect the correlation function F (q, t) through their coupling with the phonons. Using the trajectories of molecules generated by the MD simulation, I com- puted F(q, t) by AI—l F(q.t) = Z ,1; 72.0 + temp—4km. (3.88) k=0 where M is the total number of configurations of the system used in the calculation; A is the time between two successive recorded configurations, A 2 711' where 1' = 0.02 x 10“”sec. is the (real) time between each intergation step and n = 5. The total MD steps used is n11! (corresponding to T = all! 1' in real time). To avoid rapid oscillations in S(q,w) due to finite T, I average 5 (q, w) with a Gaussian weight function.62 This is equivalent to multiplying a Gaussian smoothing function to F(q,t) while carrying out the Fourier transform, i.e., T . 2 S(q,w) 2111—120 / ezwt F(q,t) e_a(t/T) dt. (3.89) —T 113 Thus the frequency resolution Aw becomes Aw = (3%) - l/W—a (3.90) I usually choose a between 2 and 6. MD results for the time dependent correlation functions at temperature T" = 0.12 and wave vectors along the FN direction are given in Fig. 3.24. This tempera- ture is well below the ferroelastic phase transition temperature (T: = 0.37). These F (q, t)’s show well defined oscillations with periodicities corresponding to those of longitudinal acoustic (LA) phonons travelling in the PN direction. The dynamical structural factor obtained form any one of these correlation functions shows a sin- gle peak at finite frequency corresponding to the LA phonon frequency. Fig. 3.15 gives one example of such S(q,w) for q in the FN direction and magnitude equal to 0.3II‘N I. The phonon width is primarily determined by the resolution function in Eq. (3.89). In Fig. 3.16, I give the dynamic structure factors for q parallel to the FN direction (see Fig. 3.10b) obtained from the F(q,t) shown in Fig. 3.14. One can clearly see the LA phonon dispersion. Fig. 3.13 gives the LA phonon frequencies for T“ = 0.12 and q going from the N point to the zone center (1") and to the K point (solid dots). The solid lines in this figure are obtained from the SCLD calculations as discussed in the last section. The overall aggrement between MD results and the SCLD results are quite good at this low temperature. rYTIUIr'IT'UUlIIII W W W Ph.!) (Ofb- 0M“) \/\/\/\/\/\/v innilmnnLLririllrn 0 1 2 3 4 ups») Fig. 3.14 Time dependent density-density correlation function of molecular center-of-mass at T‘ = 0.12. The wave vector q increases from the zone center to the zone boundary (N) point. Curves from top to bottom are for lql = 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9 and 1.0 (in units of (PAN). 114 S (q.w) (arb. unit) Fig 1,5,-...,...-,., '1 ' () T'=0.12 J ‘ =O.31"N ‘ 1.01- 3, __,' l- : QIII'N 1 0.5F- .. L_\J K“ 1 0.0 '2 .1 .1....TA.'TT‘1—T‘ 10 5 ca (1012112) . 3.15 Dynamic structure factor for T" = 0.12 and for Iql = 0.3)I‘Nl. 15 115 116 _ Z .1 .3 2.3: .... 338.35 9:. 339....» 9:2: 04% .Nm 3.5 2: =.. >7— 3 1:3an 3 .BM 35. and .H :h .a 933$ 0.5.233 9.3354. 945 an.” .mmra a: «.3.... .6. on m. o. m o .«d \ \ \ g 2. \ \ a mdvo . . \ \ 5.00.0 . \ \ 0.0 \ \\ \ o .mdle/ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\\ \ ($11an was) (m ‘b)s 117 3.6 Orientational Disorder and Phonon Anomaly The behavior of phonons across the ferroelastic phase transition is a very inter- esting problem. For example, in molecular solid K C N , ultrasonic measurments63 have shown that the shear elastic constant (C44) softens drastically when the tem- perature approaches the ferroelastic phase transition temperature. This leads to a soft transverse acoustic (TA) phonon dispersion which has been seen in neu- tron scattering measurment.64 MD simulation in K UN has also observered this elastic softening.“ In the present MD simulation, I found that for the 2D molec- ular system, the LA phonons also soften as one approaches the phase transition temperature.“5 This result will be discussed in detail below. Fig. 3.17 shows the density-density correlation function F(q, t) at T“ = 0.36 and for q along the I‘ — N direction in the first Brillioun zone. Again F(q, t) shows well defined oscilations. This implies that the phonons are still well defined in the 1" — N direction with long life time even near the structural phase transition temperature. In Fig. 3.18 the phonon frequencires at T" = 0.12,0.36 and for the wave vector q along the FN direction obtained from the MD simulation are shown as solid dots together with the SCLD calculation results. One sees that the overall agreement between the the SCLD results and the MD results are very good. This observation suggests that a strong ferroelastic order is maintained even upto to T“ = 0.36 < Tc. This is also consistent with large order parameter and strong first order phase transition seen in this system.22 The overall phonon softening with increasing temperature can be clearly seen in Fig. 3.18. Another important feature in Fig. 3.18 is that thare is a anomalous softening of the long wavelength 118 (q _<_ IFN I) LA phonons as T —+ Tc. Although the SCLD theory can account for the overall phonon softening, the anomalous softening in the low wave vector regime is beyound the scope of pseudoharmonic approximation in the SCLD. In Fig. 3.19, I give the phonon dispersion relations at T“ = 0.36 and for q vector going from the zone center (I‘ point) to the zone boundary K point and N point. The agreement between the SCLD results and the MD results are good in the F — N direction, but not as good in the I‘ — K direction. However, the agreements at T" = 0.12 shown in Fig. 3.13 are good in both the I‘ -— K and F — N directions. This suggests that as one approaches the structural phase transition from below, the anharmonicity effects are different depending on the directions of the wave vector q. Figure 3.20 shows the density-density correlation functions at T“ = 0.40 and for q in the FN direction. Figure 3.21 shows the density correlation and the corresponding structure factors at T“ = 0.36 and T“ = 0.40 and for the wavevector at the N point of the zone boundary. Similar to Fig. 3.17, at T" = 0.36, F(q,t) is a well defined oscillating function of time and S(q,w) has a single peak. But at T" = 0.40, F(q,t) shows oscillations with rapid dacay. The structure factor at T“ = 0.40 is considerably more complicated. There is a large central peak coming from the decaying bahavior of the time correlation function and a broad peak which can be identfied as a renormalized LA phonon peak. Similary, Fig. 3.22 shows the MD results of F(q,t) and S(q,w) at T“ = 0.12,0.36 and 0.40 and for q half way from the I" point to the N point. One sees that For T" 5 Tc"', the density correlation functions have well defined oscilations with sharp phonon peaks in the corresponding dynamic structure factors. As T > Tc, the density correlation function has a decaying envelope in addition to the oscillating behavior and the dynamic structure factor develops a central peak. Figure 3.23 gives the 119 MD results of the density correlation functions at T“ = 0.40 and 0.50 and for q at the K point of the zone boundary. Compared with Fig. 3.14, it is clear that F (q,t)’s at T > Tc in the I‘ — K direction are fast decaying functions of time and there are almost no oscillations. These observations suggest that phonons are heavily damped in the high temperature orientationally disorded phase. This is due to the increased importance of the intrinsic anharmonic effect6 and the strong coupling between the translational and rotational degrees of freedom.“’*66 In. Fig. 3.24, I give the MD results for the LA phonons in the I‘ — N direction at high temperatures (T‘ = 0.40,0.50). There is a small hardening of the phonons with increasing temperature which can be traced to the effect of rotation-translation coupling.1 3 In Fig. 3.25 I show the sound velocity of the LA phonons along the I" — N direction extracted from the molecular dynamics at several temperatures, both below and above the ferroelastic phase transition temperature.4 It is evident that the elastic constant softens as the temperature approaches Tc. For T << Tc, I can use the ground state (T = 0) values of the lattice parameters to calculate the phonon velocities.:"4 The sound velocities of the LA and TA modes obtained from the SCLD calculation are (2080.9,905.6)m/3ec in the I‘ — K direction and (2179.1, 1281.8)m/3ec in the I‘ — N direction. The value 2179.1m/sec for the LA sound velocity in the F—N direction is about 10% higher than the T“ = 0.12 value, suggests that our T = 0 rsults are consistent with the T" 72’: 0 results obtained from MD and SCLD calculations. 120 ritrrlltllrrrlTlil T'=o.36 qHI'N F(q,!) (orb. units) 1 1 1 1 l 1 L; L l 1 m1 1 l L 1 1 O l 2 3 4 ”9'06, Fig. 3.17 Time dependent density-density correlation function of molecu- lar center-of-mass for T‘ = 0.36. The wave vector q increases from the zone center to the zone boundary (N) point. Curves from top to bottom are for M = o.1,o.2, 0.3, 0.4, o.5,o.e,o.7,o.8,o.9 and 1.0 (in units of erl). 121 ON out (Furs) a cm on fi 06 md v.0 md .31..." -3 qum 2: an 3:: 1.42.. 22 3.3 31.3.— 3333? 33032: .2: 93 .....3Q 45:: «ZV Defiance 0.2 3 «r: .3239 23s 0.: 59¢ 7 2% van .3.—3393:». neuuohmn 9:. 3m 0323 umfifioohfi 22 5 maowmammmfi .3523 <4 wad .mmh 3:3 e830,; ”SSE, we; U U V Vii T I rr' l dld ‘ d V O ‘..1“‘l‘1“|‘.““‘|‘ddl o\. l l l l cu ma (ZH 3,01) m 122 .8on <8 1.8 45 ha 338a QGDW 2: 0.3 8a.: 33m. 638$ 8.58.8? 9585:. 85 .3 acomgmmon «8m 2: 89: «855.88% 8.8 .83.: 31.8.. Q2 .2: 8.8 flow 33% .mumuaw==oa 8:3 5.55am 8.: .«o «nonmmoa 8.: 8.8 8:: ~53 18.38.» 2:. 6.33836 55835.? «8.35.6 .a 3 SH 3:: «.3.—93:9: 83.? 9.33888.‘ 8.: 333 8.5 cm... H :u .e 30.38%? .3.—23 «3 an.” .mmh SE: vooavmmv v z a a ma (2H3,01) (*1 rIfiTI'UVTUrIU'I'rIU 1‘30.“ qHI'N F(q,!) (orb. unite) AlllllllAILLJlLlL41 o 1 2 3 4 t (pace) Fig. 3.20 Time dependent densitrdensity correlation function of molecu- lar center-of-mass for T‘ = 0.40. The wave vector q increases from the zone center to the zone boundary (N) point. Curves from top to bottom are for |q| = o.1,o.2,o.3,o.4,o.5,o.6,0.7,o.8,o.9 and 1.0 (in units of erl). 123 124 1 4 an 1 4 u- t 1 1 1 ff'v 10 T' - I .0 w a) 101) vf'v'v ‘l5 '- .. l A A A A l A A A A1 9’ u- S(q, w) (arb. units) ‘VVIV'Vj'VTTVIVV F’ o '1“ = 0.40 2.5 F(q, t) (arb. units) L5 LO AlAAAAlA‘AALIAAAlAAALlAA (L5 JAAIAAAAILAAALLAAA (LO 5 10 15 20 t (psec) a) (1012 Hz) Fig. 3.21 Comparison ofF(q,t) and 5(q,w) obtained from the MD simula- tions at temperatures T‘ = 0.36,0.40 and for q corresponding to the point N of the Brillouin zone. 125 J O J O 4.. o 3 0 A33: .83 3.3m ‘ z o q“““J‘“q““‘-““ "bhbb’tpb’r’ppy’h ' D D - ’ h b ' -L ’ ’ by i O ’bb’bD-rl. ’ pb”.-b’. “‘q““ ’ bahfihher-F’bDv-L‘Pb-t .<.q< ‘ 1““1‘ ‘ ‘ 4 ‘ ‘ i- ‘ ‘ ‘ ‘ r-om iH.FDrD.-Dbbi ‘ ‘ ‘ 1 A 0 nm 0 4. O . A33: .553 36E 2 o 1 LAAAI 4L n- 15 20 10 a) (10 “Hz) t. (psee) Fig. 3.22 Compaison of F(q,t) and 5(q,w) obtained from the MD simula- tions at temperatures both below and above the ferroelastic phase transition and {or q halfway fi'om the F point to the N point of the Brillouin zone. 126 -mmw 5?:232880 of. 5 $3.. Ah «0 3:5 .5 85:. 2: an q 1 d d I 1 d — d 1 d d dot-2:58 Q: 2: 59¢ uni—~36 8.5; ~88qu 8:5 coafioauo... .....8:ofi-b.~a=um~ «a» an.” .mmfl Pd u.— (0 'bM/(Tbm 127 .333358 Q3 2: 39¢ 303.821. 821 1883.5 -mmt 53:33.8..8 2: E 28.32:... 55.3088 >...88v:3.88Q T: and .wmn— A... mo 33:: :3 055. can :2 , 2: S T . . . _ . . . . _ . w . . _ J . . . u.— 06 (o 'bM/(z‘bm 128 _vvvvl 1T I varfi I 'fl 15~ - T' 040 i _ o T'—0.50 : g3? 10; qllf’N f? e: ’ {I j 3 7 v I- 0 u a 5f , '. 3 —. .2 3 0 .. 1...-1....1-.- 1.4..1‘ 0 0.2 0.4 0.6 0.8 1 WAVE VECTOR (reduced unit) Fig. 3.24 LA phonon dispersion in the orientationally disordered paraelastic phase obtained fi'om molecular dynamics. SOUND VELocmuo'm/sec) 129 25“ ‘I ‘ l I "I‘ '1' ‘ 15- i J 10- J : T. m : 5 ...J....l--..1.L.‘.1....l.-.. O 0.1 0.2 0.3 0.4 0.5 0.6 Temperature Fig. 3.25 Temperature dependence of LA phonon velocity in the FN direc- tion. 130 3. 7 Summary Using the constant-pressure molecular dynamics simulation, I have investi- gated the thermodynamics and the dynamics of a two dimensional diatomic molec— ular monolayer undergoing a ferroelastic phase transition. This system closely re- sembles the 6-phase of oxygen molecule adsorbed on graphite surface. For Lennard- Jones parameters appropriate for the oxygen molecules, I find a first order tran~ sition from an orientationally ordered distorted trangular structure (ferroelastic phase) to an orientationally disordered equilateral triangular structure (paraelas- tic phase). The transition temperature is 20.1K compared to 26K for oxygen on graphite (coverage z 8 molecules per 100 4:12). There is a softening of the elastic constants near the transition, particularly in the paraelastic phase; this can be understood in terms of translation-rotation coupling. Comparsion between phonon frequencies for certain symmetry directions obtained by using SCLD the- ory and MD simulation clearly shows the presence of large anharmonicity effects in the paraelastic phase. A rapid quench from the high temperature phase to very low temperatures indicates the presence of small clusters (consistenting of 6 to 12 molecules) with both ferroelastic and herringbone ordering. In addition, I find a large density of equilateral trangular plaquettes. These give rise to a three-peak structure in the center-of-mass radial distribution function. Although the intermolecular interaction parameters are those appopriate for oxygen molecules, there are important differences between the system I have stud- ied and oxygen molecules adsorbed on graphite. Two significent physical effects not considered in the present simulation study are (1) orientations of the molecules away from the graphite plane with a concomitant motion of the center-of-mass 131 perpendicular to the graphite substrate, and (2) corrugation of the graphite sub- strate. For a detailed comparision between the MD results and experiments on Oz/graphite, the above two effects must be taken into account. Our MD simulations clearly brings out the physics of ferroelastic phase tran- sition in two dimensional molecular solids. The presence of an intermediate plastic phase is clear and the physical properties of this phase is dominated by strong rotational-translational coupling. In contrast to 3D systems, where one of the many transverse elastic constants soften as a result of the above mentioned cou- pling, there is only one transverse elastic constant (C44 = CH — 012) for a 2D triangular lattice and this softens as one approaches the ferroelastic transition temperature from above. The system therefore behaves almost like a liquid al- though the center-of—mass diffusion is absent. In fact, this is perhaps the reason t,38 it is not easy to distinguish this intermediate plastic phase why in experimen from the usual liquid phase. In addition, I find that the longitudinal acoustic phonons also soften as T approaches the ferroelastic phase transition temperature. The phonons in the plastic phase are strongly damped whereas in the low temperature ferroelastic phase they are rather well defined and can be understood within the SCLD theory. A direct experimental observation of the soft phonons and their behavior as a function of temperature will be of great help in elucidating the physics of strongly coupled rotation-translation system in 2D. The other sig- nificant aspect of this study is the possible effect of hydrostatic or uniaxial stress on the ferroelastic phase transition. As I have disscussed in section 3.2, there are important orientational coordinate and velocity dependent contributions to the internal stress tensor whose effect on the thermodynamic and dynamic properties will be of interest to explore further. List of References for Chapter 3 1 2 3 4 5 6 10 11 12 13 14 132 List of References for Chapter 3 D. Sahu, Ph. D Dissertation, Michigan State University, 1983. S.D. Mahanti and S. Tang, Superlattices and Microstructures l, 517 (1985). W. Jin and S.D. Mahanti, Bull. Am. Phys. Soc. 33, 626 (1988). W. Jin, S.D. Mahanti, and S. Tang, Solid State Commun. 66, 877 (1988). S. Tang, W. Jin, S.D. Mahanti and R.K. Kalia, Phys. Rev. B 39, 677 (1989). AD. Bruce and R.A. Cowley, Structural Phase Transition (Taylor and Francis, London, 1981); R.A. Cowley, Adv. Phys. 29, 1 (1980); AD. Bruce, Adv. Phys. 29, 111 (1980); AD. Bruce and R.A. Cowley, Adv. Phys. 29, 219 (1980). L.D. Landau and E.M. Lifshitz, Statistical Physics (Pergamon Press, Oxford, 1958). G. Shirane and Y. Yamada, Phys. Rev. 177, 858 (1969); R.A. Miiller and W. Berlinger, Phys. Rev. Lett. 26, 13 (1971); J. Feder and E. Pytte, Phys. Rev. 131, 4803 (1970); A.D. Bruce and R.A. Cowley, J. Phys. C 6, 2422 (1973). ‘ S.D. Mahanti and G. Kemeney, Phys. Rev. B 20, 2105 (1979). ML. Klein, I.R. McDonald, and Y. Dzaki, Phys. Rev. Lett. 48, 1197 (1982) F. Schwabl, Ferroelastic 24, 171 (1980). RA. Fleury, J. Acoust. Soc. Am. 49, 1041 (1971); Comments Solid State Phys. 4, 149 (1972); ibid., 4, 167 (1972); E. Pytte, Comments Solid State Phys. 5, 41 (1973); ibid., 5, 57 (1973); G.A. Samara, Comments Solid State Phys. 8, 13 (1977). D. Sahu and S.D. Mahanti, Phys. Rev. B 26, 2981 (1982); 29, 340 (1984); Solid State Commun. 47, 207 (1983); S.D. Mahanti and D. Sahu, Phys. Rev. Lett. 48,939 (1982). K.M. Rabe and J.D. Joannopoulos, Phys. Rev. B 32, 2302 (1985); Phys. Rev. Lett. 69, 570 (1987); Phys. Rev. B 36, 6631 (1987); K.M. Rabe, Bul. Am. Phys. Soc. 33, 276 (1988). 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 133 T. Schneider and E. Stoll, Phys. Rev. B 13, 1216 (1976); ibid. 17, 1302 (1978). V.K. Wadhawan, Phase Transitions, 3, 3 (1982). T.A. Scott, Phys. Reports 27C, 89 (1976). G.C. DeFotis, Phys. Rev. B 23, 4714 (1981). K.H. Michel and J. Naudts, Phys. Rev. Lett. 39, 212 (1977); E. Ohno, K. Yoshimitsu, Prog. Theo. Phys., Suppl, 80, 180 (1984). S.D. Mahanti, in Condensed Matter Theory, Vol. 2, Edited by P. Vashishta, R.K. Kalia, and R.F. Bishop (Plenum, New York, 1987). W. Dultz, M.H. Otto, H. Krause, and J.L. Buevoz, Phys. Rev. B 24, 1287 (1981). S. Tang, S.D. Mahanti, and R.K. Kalia, Phys. Rev. Lett. 56, 484 (1986). CE. Vilches, Ann. Rev. of Phys. Chem. 31, 463 (1980). J. Krim, J.P. Coulomm, and J. Bouzidi, Phys. Rev. Lett. 58, 583 (1987). RA. Heiney, P.W. Stephens, S.G.J. Mochrie, J. Akimitsu, R.J. Birgeneau, and RM. Horn, Surf. Sci. 125, 539 (1983). DD. Awschalom, G.N. Lewis, and 8. Gregory, Phys. Rev. Lett. 51, 586 (1983). J. Stoltenberg and O. Vilches, Phys. Rev. B 22. 2920 (1980); R. Marx, Phys. Rep. 125, 1 (1985); D. Kirm, B. Kuchta, and R.D. Etters, J. Chem. Phys. 87, 2332 (1987). M. Nielsen and J.P. Mctague, Phys. Rev. B 19, 3096 (1979). M.F. Toney, R.D. Diehl, and 5.0. Fain, Phys. Rev. B 27, 6413 (1983). R.D. Etters, R.P. Pan, and V. Chandrasekharan, Phys. Rev. Lett. 45, 645 (1980). ' R.P. Pan, R.D. Etters, K. Kobashi, and V. Chandrasekharan, J. Chem. Phys. 77, 1035 (1982); O.M.B. Duparc and R.D. Etters, J. Chem. Phys, 86, 1020 (1987). M. Parrinello, in Molecular-Dynamics Simulation of Statistical-Mechanical Systems, ed. by G. Ciccotti and WC. Hoover (North-Holland, Amsterdam, 1986), p. 204. 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 134 G. Herzberg, Molecular Spectra and Molecular Structure, I. Spectra of Di- atomic Molecules, 2nd ed. (D. Van Nostrand, Princeton, 1965). C.A. English and J.A. Venables, Proc. R. Soc. Lond. A 340 , 57 (1974). C.A. English, J.A. Venables, and D.R. Salahub, Proc. R. Soc. Lond. A 340, 81 (1974). Y.B. Gaididei and V.M. Loktev, Phys. Lett. A 123, 85 (1987). W.A. Steele, Surf. Sci. 36, 317 (1973). M.F. Toney and 5.0. Fain, Phys. Rev. B 36, 1248 (1987). J.K. Kjems, L. Passell, H. Taub, and J.C. Dash, Phys. Rev. Lett. 32, 724 (1974). J. Eckert, W.B. Ellerson, J.B. Hastings, and L. Passell, Phys. Rev. Lett. 43, 1329 (1979). A.D. Migone, H.K. Kim, M.H.W. Chan, J. Talbot, D.J. Tildseley, and W.A. Steele, Phys. Rev. Lett. 51, 192 (1983). R.D. Diehl, M.F. Toney, and 5.0. Fain, Phys. Rev. Lett. 48, 177 (1982). W. Jin, S.D. Mahanti and S.Y. Tang, Phys. Rev. B 39, 11928 (1989). M. Parrinello, in Molecular-Dynamics Simulation of Statistical-Mechanical Systems, edited by G. Ciccotti and W.G. Hoover (North-Holland, Amster- dam, 1986), p. 204. M.L. Klein, Ann. Rev. Phys. Chem. 36, 525 (1985). H.C. Anderson, J. Chem. Phys. 72, 2384 (1980). M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980); J. Appl. Phys. 52, 7182 (1981); J. Chem. Phys. 76, 2662 (1982); in Melting, Localization and Chaos, ed. by R.K. Kalia and P. Vashishta (Elsever Sci. Pub. Co., 1982), p. 97; M. Parrinello, A. Rahman and P. Vashishta, Phys. Rev. Lett. 50,1073 (1983). S. Nosé and M.L. Klein, J. Chem. Phys. 78, 6928 (1983); Mol. Phys. 50, 1055 (1983); Phys. Rev. Lett. 50, 1207 (1983); S. Nosé, J. Chem. Phys. 81, 511 (1984); S. Nosé and F. Yonezawa, J. Chem. Phys. 84, 1803 (1986); M. Sprik, R.W. Impey, and M.L. Klein, Phys. Rev. B 29, 4368 (1984). 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 135 S. Tang, Ph.D dissertation, Michigan State University, 1985. W. Jin and S.D. Mahanti, Phys. Rev. B (to be published). A.A. Maradudin and V. Ambegaokar, Phys. Rev. 128, 2589 (1962); Phys. Rev. 135, A1071 (1964); A.A. Maradudin and A.E. Fein, Phys. Rev. 128, 2589 (1962); V. Ambegaokar, J. Conway, and G. Baym, in Lattice Dynamics, edited by R.F. Wallis (Pergamon, New York, 1965), p. 261. T.M. Hakim and H.R. Clyde, Phys. Rev. B 37, 984 (1988); J.P. Hansen and M.L. Klein, Phys. Rev. B 13, 878 (1976). S. Califano, V. Schettino, and N. Neto, Lattice Dynamics of Molecular Crys- tals, Vol. 26 of Lecture Notes in Chemistry (Springer, Heidelberg, 1981). G. Venkataraman and V.C. Sahni, Rev. Mod. Phys. 42, 409 (1970). O. Schneep and N. Jacobi, in Dynamical Properties of Solids, edited by G.K. Horton and A.A. Maradudin (North-Holland, New York, 1975), chapter 3, page 151-204. Condensed Matter Research Using Neutrons (Plenum Press, New York, 1984), edited by S.W. Lovesey and R. Scherm. M.L. Klein, Ann. Rev. Phys. Chem. 36, 525 (1985). M. L. Klein, in Molecular-Dynamics Simulation of Statistical-Mechanical Sys- tems, edited by G. Ciccotti and W.G. Hoover, (North-Holland, Amsterdam, 1986), p. 424. S.D. Mahanti, G. Zhang, and G. Kemeney, Phys. Rev. B 35, 8551 (1987). S.W. Lovesey, Condensed Matter Physics: Dynamic Correlation, 2nd. ed, (Benjamin/Cummings, 1986). L. Van Hove, Phys. Rev. 95, 249 (1954); 93, 1374 (1954); L. Van Hove, N.M. Hugenholtz, and LP. Howland, Quantum Theory of Many-Particle Sys- tems (Benjamin, New York, 1961). T.R. Koehler and RA. Lee, J. Comput. Phys. 22, 319 (1976). S. Haussiihl, Solid State Commun. 13, 147 (1973). A. Loidl, J. Knorr, J. Daubert, W. Daltz, and W.J. Fitzgerald, Z. Phys. B 38, 153 (1980). M.L. Klein and LR. McDonald, J. Chem. Phys. 79, 2333 (1983). 136 66 K.H. Michel, Z. Phys. B 54, 129 (1984); K. Zabiriska and P. Zieliriski, Z. Phys. B 58, 173 (1985). 137 Chapter 4 Theory of Neutron Scattering From Paired Fermion Superconductors 4.1 Introduction A large number of metals and alloys undergo a transition to a phase called the superconducting state below a certain critical temperature Tc. The spectacu- lar physical properties of this phase, such as the zero electrical resistance, perfect diamagnetism and magnetic flux quantization, have made it one of the most fasci- 1 nating areaof condensed matter physics. The basic theory of superconductivity is the Barden-Cooper-Schriefi'er (BCS)2'3 theory and it’s strong coupling version - Eliashberg theory.4"l° Recently, a great flurry of activity has been stirred up by the discovery of superconductors with critical temperature near 100K in the 11—15 layered oxide perovskites. In this chapter, I will develop a theory of inelas- tic neutron scattering from superconductors and discuss the applicability of this 14,15 theory to the real oxide superconductors as well as layered superconductors in general. 138 In the remaining part of this section, I will briefly review currently known structural and physical properties of the high-Tc oxide superconductors, introduce different theoretical models which have been proposed for these and present my motivations to study the problem of inelastic neutron scattering from the gap excitations in these superconductors. Structure and Properties The new high-Tc oxide materials include hole-doped single-layered Lag_,R,Cu04 (R = Ba, Sr, Ca) (also called 214) compound,11 double-layered YBagCu307_,(123) compound,12 multi-layered BiCaSrCuO and TlCaBaCuO compounds,16 and cubic Ba1_,K,BiO3 compound.17 In addition, electron-doped N d2-,Ce,Cu04 compound has been recently discovered.18 The 123 and 214 compounds are the better investigated materials. The common structural feature of these high-Tc superconductors (with the exception of the cubic compound) is the presence of one or more planes of copper atoms with four strongly bonded oxygen atoms in a square planar arrangement (Cqu planes, also called the ab—plane). In some materials these planes occur in groups with individual planes inside a group being separated by one Y or C'a layer and the groups being intercalated by a variable number of CuO (or BiO, etc.) layers. There is a slight orthorhombic dis- tortion of the tetragonal lattice at low temperatures in both La2_,Sr,Cu04 and YBag Cu307_,. The quasi-two-dimensional layered character of these compounds is studied by varieties of experiments such as X-ray and neutron diffraction,19 and is also supported by band-structure calculations” which reveal an antibond- ing copper 3dz'-v' and oxygen 2pm, band near the Fermi level with very little dispersion along the c-axis. 139 The most important element in these oxide superconductors is the quasi-two- dimensional 0120; plane. In the undoped material, each copper site is associated with one hole. These holes do not move freely but are localized near the individ- ual copper sites by strong onsite Coulomb repulsion. The spins of the holes are ordered antiferromagnetically (AF) due to the superexchange interaction between them,21 and their magnetic properties can be understood very well using the two- dimensional spin 1 / 2 Heisenberg model.22'23 Doping creates additional holes going into predominatly oxygen 2p orbitals. The local moments are disordered as seen in Raman scattering, infrared optical24 and neutron scattering25 experiments and the system looses the AF order and becomes superconducting at low temperature. The NMR results26 show that the nuclear relaxation rate of the copper nuclear spins varies drastically across Tc, indicating that the holes (or electrons) created by dopping in the Cqu layers are the bearers of the supercurrent. This. picture is also supported by high-energy spectroscopy such as X—ray absorption near-edge spectroscopy” and electron energy loss spectroscopy?"B Experimental observations show that the unit of the superconducting cur- rent in the oxides is twice the electronic charge. These expreriments include flux 9 in units of 5%, the AC Josephson effect30 showing superconductive quantization2 tunneling of charge 2e, and Andreev reflection.31 These observations indicate that Cooper pairs are still responsible for superconductivity in the new oxide supercon- ductors. In addition, the observation32 of a DC Josephson tunneling through a junction between a singlet state superconductor PbSn and YBag Cu; 07 indicates that the pair state in YBazCu307 has ”s-wave” symmetry, since a Josephson supercurrent would be forbidden between a singlet s-state superconductor and a superconductor of different pairing symmetry.32 140 All these layered superconductors are highly anisotropic in their physical properties33 with London penetration depth in the ab—plane about 1500.4 which is much smaller than that in the c-direction. The superconductivity is strongly type-II,“ with He; much higher than Hcl(-Hc2 > 100T). The coherence length 6, the size of a Cooper pair, is short - about 10.21 in the ab—plane for YBazCu307, which is less than twice the average distance between oxygen holes. The coherence length is also much shorter than the magnetic field penetration depth A. These last two properties have lead to interesting fluctuation phenomena in these oxide superconductors, like melting of magnetic flux lattice.35‘37 Theoretical Models In the past, a number of possibile theories have been suggested,”*38 aimed to explain, in addition to the high-T.2 itself, the aforementioned properties of these layered oxides, both in the normal and the superconducting states. These theo- retical models can be generally divided into two main catagoreies according to the statistics of the underlying quasiparticles (QP) which form Cooper pairs.15 One is the model of superconductivity with paired fermions2 although the attractive 39—50 mechanism is still uncertain. The other is that of paired charged spinless 51—53 bosons. Some essential characteristics of these two types of models are briefly discussed below. In the paired fermion models, the carriers in the normal state form a normal Fermi liquid54 and condensation of these charged carriers near the Fermi level gives rise to the superconductivity. The quasi-two-dimensional character of the elec- tronic states with ”nesting” of the Fermi surface plays an important role40 in the 141 observed large Tc, although the superconductivity in these layered oxides is strictly three-dimensional in nature because of weak but finite interlayer coupling."’°’55’57 The source of the electron-electron attraction which results in Cooper pairing can involve lattice degrees of freedom (phonons) and other attractive forces of electronic 20,39 origin. However, estimation of Tc using strong coupling theory indicates that the electron-phonon attraction alone is not strong enough to give the observed high Tc. Many theoretical arguments38 and experimental spectroscopic results58 such as the existence of satellites in the resonant photoemission spectra59 have demon- strated strong electronic correlation effects in these oxide superconductors. N on- conventional pairing mechanisms which are purely electronic in origin based on Coulomb correlations between copper 3d electrons (and oxygen 2P=m electrons) have also been proposed.” Some of these are base on spin fluctuations, like 43 antiferromagnetic spin fluctuations,“"42 ferromagnetic spin fluctuations, and local-moment fluctuations.“ These have been studied using a Hubbard 60—62 Hamiltonian, usually the single band Hubbard model on a square lattics, i.e., H = —t 26")“, c$cj,+ U 2, nngil, [where Cluck) is the creation (annhi- lation operator of an electron (or hole) with spin 0‘(=T, l) at site i, and ng, = citcg, J with nearest neighbor p — d hopping integral t and on-site Coulomb repulsion U 51 between electrons with opposite spins. In the strong correlation limit (U >> t), this one band model can be expanded as a series in t/ U by use of canonical 40'“ and can be written as an effective Hamiltonian for mobile va- transformation cancies moving with a hopping integral t through an antiferromagnetic backgroung with exchange interaction J 2 4t2 / U between spins S; [where S.- = %c;;&'agc;g]. 142 For exact half filling, this is an antiferromagnetic Mott insulator?“62 The pa- rameters (t a: 0.1eV, U z 5eV, etc. ) entering into the Hubbard Hamiltonian can be calculated from the electronic structures using different approximations."3 The other paired fermion models are based on charge fluctuation effects, such 8 49 as excitons" or plasmons. In the usual ”excitonic” superconductivity,“4 the conduction electrons in conducting plane induce virtual electronic transitions on nearby easily polarizable molecules or atomic complexes, which result in an ef- fective attractive interaction between conduction electrons. In a different version of the ”excitonic” model,48 the charge transfer excitations between copper and oxygen atoms in the conducting Cqu planes (Cuz'l'Oz‘ ——» Cu‘l'O') serve as the pairing exchange force. This has be studied48 using an extended Hubbard Hamiltonian containing copper (122-1,: orbitals and oxygen Pam orbitals. Next, I will introduce the boson models, where the underlying quasiparticles are spinless charged bosons. As discussed above, the 3d,2_y2 copper electrons and 21),,” oxygen electrons in the 011.02 planes are strongly correlated and form a system which is close to the metal- insulator Mott transition limit.51 Numeri- cal studies"5 find that near the half-filled limit, the one-band Hubbard model can have antiferromagnetism for a wide range of t / U- as a result of the Coulomb in- teraction. Increasing U over a critical value U c causes a band spliting and results in the metal-insulator transition. At exactly half filling, the system is described by the spin 1/2 antiferromagnetic quantum Heisenberg model. This model has long range ordered antiferromagnetic ground state in two dimensions. Doping by a small number of holes introduces frustration into the quantum antiferromag- net and destabilizes the antiferromagnetic order. The new ground state of the 143 doped system was hypothesized“ to be a quantum-spin—liquid-state (or resonat- ing valence bond state) which can be approximated by an appropriate quantum mechanical superposition of states where two spins are paired into singlets. On the basis of the topological considerations, the quasiparticles in the spin-liquid-state have been conjectured" to be neutral spin 1/2 fermionic excitations and charged spinless bosonic excitations. The neutral spin 1/ 2 fermions are the quasiparticles of the undoped system. They consist of unpaired electrons moving in the spin- liquid background with properties analogous to solitons in polyacetylene.‘m The effective mass (m) of these fermions is of the order of hz/(azl) where a is the lattice constant constant and J = 47‘]: is the exchange interaction between nearest neighbor spins. Doping generates holes (removes electron of charge —e and spin 1 / 2) in the Fermi sea of these neutral fermions. The hole can bind to a neutral fermion and form soliton-like excitations. These are spinless charged quasiparticles and obey Bose statistics and have an effective mass of the order of the electron mass. Thus in this picture, a real electron or a hole (a fermion with charge :l:e . and spin 1 / 2) is a composite particle consisting of a neutral spin 1/2 fermion and a charged spinless boson. These bosons are the charge carriers of the doped spin- liquid-state. In addition, any effective attractive boson-boson interaction, which can be mediated for example through an interlayer Josephson coupling,” or any other mechanism will produce pairing of these bosons in quasi two-dimensions. Condensation of these paired charged bosons results in superconductivity.53 However, neither the paired fermion nor the boson models have been com- pletely successful in explaining all the experimental observations. 144 Neutron Scattering From Excitations in Superconductor From the above discussions, we see that the different theories proposed for high-Tc oxide superconductors up to now are only partially successful. In princi- ple, a proper theory should not only give the high Tc, but also explain all other 69 This requirement places properties and phenomenon of these superconductors. a strong constraint on the possible theoretical models on high-Tc oxide supercon- ductors. Besides, accurate determination of the superconducting parameters such as the superconducting order parameter A, symmetry of the order parameter, the coherence length 5 etc., can provide important information about the possi- ble microscopic mechanisms. Many spectroscopic and scattering experiments have become more and more useful in this regard. Some photoelectronic probes,58 like the infrared and optical absorption, Ra- 7 and man scattering, photoemission, inverse photoemission, x-ray absorption,2 other probes such as positron annihilation and electron energy loss spectroscopy“ can give information about the temperature dependent spectrum of elementary excitations in a superconductor, starting from meV to eV energy scales.” Neu- trons, due to their magnetic moment, can directly couples to the magnetic fields produced by these elementary excitations. Thus neutron scattering is an excellent probe of magnetic correlations.” Neutron scattering has been used in studying the 11 magnetic excitations and magnetic orderings in heavy-electron superconductors, antiferromagnetic superconductors," and the high-Tc oxide superconductors.25 There have been primarily three experimental techniques which have been use- ful in probing the superconducting gap excitations directly. These are tunneling, 73 far—infrared absorption and Raman scattering measurements. In this chapter 145 and the following one, I will discuss how neutron scattering can directly probe the gap excitations and other superconducting correlations. In particular, inelastic and quasielastic neutron scattering from superconductors can directly probe the order parameter A, the coherence length 6, and most importantly, the statistics of the quasiparticles.1“*l5v74v75 The remaining part of this chapter is organized as follows. In section 4.2, I summarize some most important general formulas about inelastic neutron scatter— ing from superconductor. In particular, the relation between the inelastic scatter- ing cross section and transverse current-current response function will be discussed. In section 4.3, I investigate the dynamic response functions from the paired fermion superconductor and discuss some important features of inelastic and quasi elastic neutron scattering from the gap excitations. Analytical and numerical results will be given in this section. Application of this theory to the high-Tc oxide supercon- ductors will be discussed. Finally in section 4.3 I'give a summary. The case of paired spinless cahrged boson superconductor will be discussed in next chapter. 146 4.2 Inelastic Scattering Cross Section of Neutrons From Superconductors In this section, the essentials of neutron scattering by supercondutors will be discussed in order to bring out certain salient features. When a neutron is scattered by a superconductor, the magnetic moment of the neutron interacts with the effective local space and time dependent fluctuating magnetic field B that is intrinsic to the quasiparticle spin, its itinerancy, and also any local magnet moments in the solid. The interaction can be written as” Ti" 2 —'yufo - B(r,t), (4.1) where 6" denotes the Pauli matrices, 7 = —1.91 is the gyromagnetic ratio of the neutron, and [LN = eh/Zmpc is the nuclear magneton (mp is the mass of proton). Note that the circumflex above a letter denotes an operator. The differential scattering cross section dza/dfldE per unit outgoing solid angle 0 and outgoing energy E within the Born approximation can be written as 70'" A 2 < A'k'a'lé-BIARU > 6(hw+EA—E;‘.), dza = ( mp )2(7uN)2E Z PAP, dfldE 21rh2 lc MWU' (4.2) where k(k') and 0(0' ) are initial (final) neutron wave vector and spin, MN) spec- ifies the initial ( final) state of the QP in the solid, hw = 7120c” — lcz)/2mN is the energy transfer (mN is the mass of neutron), PA and P, are two probability weights of the initial QP and neutron states. For unpolarized neutrons, the sum 147 over the initial neutron state can be carried out. Then by using an integral rep- resentation of the energy conservation 6-function in Eq. (4.2), the scattering cross section can be written in terms of a correlation function of the magnetic field." The effective microscopic magnetic field B is self-consistently generated ac- cording to Maxwell’s electrodynamic equations by virtue of the currents produced in the solid. This current is composed of an orbital paramagnetic part jo, due to the motion of QP, an intrinstic spin part j., due to the spin of the moving quasi- particles, and the orbital diamagnetic current jg, due to its interaction with the electromagnetic fields. Using the microscopic Maxwell equations including both the instantaneous Coulomb field and the effective fluctuating magnetic field,77 and the fluctuation-dissipation theorem," the magnetic field correlation function can be related to the transverse current-current response function with a diamagnetic screening factor. Thus in general, the total differential scattering cross section con- sists of three parts; one, arising from the paramagnetic current; two, a diamagnetic contribution due to the intrinsic electromagnetic fields in the system and is usu- ally combined with the first for reasons to be explained later; three, spin current density contribution due to the intrinsic fluctuating spin moments of the itinerant QP. In addition, localized magnetic moment, for example those arising from the Cu‘l‘+ ions in the copper-oxygen complex in the high Tc oxide superconductors, also contribute to the scattering. Here I will only focus on the first three parts, leading to a possible determination of the superconducting characteristics of the system over and above the contribution due to the magnetic moments of the fourth part which may or may not be directly related to the superconducting aspects of the system and may even be absent in copperless oxide superconductors such as Ba1_,K,Bi03. 148 These contributions together make up the expression for the inelastic differen- tial scattering cross section in the form of the transverse current-current response function and can be written as" d2 2 2 k! dmiE = (21:21:?) (7:2,) f(41r)2[1 + "MM” 1 n(q,w)’ x [—ImP1(q.w)/(czqz)] - (4.3) where q is the momentum transfer k - k', and n3(hw) is the usual Bose factor [n3(:c) = (emB — 1)‘1,B = Ei—T' is the inverse temperature]. Pi(q,w) is the trans- verse part of the wave vector and frequency dependent current-current response function, i.e., P-L(Q:w) = :05!» " qu‘lu)PMV(q,w)2 (4-4) pu where p, V(= 2:, y, z) are cartesian indices and Pm,(q, w) is the spatial and temporal Fourier transform of the total current-current response function Pm,(:c, t; m', t' ), i.e., ddq dw i-x—x' —iw -' Ppp(a:,t;z',t')=/(27r)d 2; e‘“ )e (t t)P,,,,,('q,w). (4.5) In the above d(= 2,3) is the spatial dimension. The denominator K.(q,w) in Eq. (4.3) is due to the self-consistent electrodynamic screening effect discussed above and is given by77'" + [flL—Impumw] 2. (4.6) 149 2 p = 41me2 /m is the plasma frequency, 77. is the average density of the where w quasiparticles and m is the mass of the quasiparticles. The electrodynamic screen- t" only at very small q and in our numerical calculations, I omit ing is importan this screening effect by putting n(q,w) % 1. Eq. (4.3) is the basic equation for studying the inelastic neutron scattering from QP in solids in general and superconductors in particular. The current-current response function represents the response of the total currents to an applied magnetic field in superconductor. Let 7:! be the halimitonian for the superconductor without the pertubating field and 72' = f dzC(:c,t)é(:c,t) be the perturbation due to an applied time-dependent field 3(3, t) (such as electric field or magnetic field), where 0 is a physical operator (such as particle density or electric current operator). Then the lchange in the diagonal matrix element of the operator C is given by 6(C(z,t)) = 71‘- } dt' fdrc'DR(a:t,:c't') E(:c',t'), where the time retarded correlation function is-doefined byn'so DR(m,t;z',t') = —i < [Ox(z,t),Cx(m',t')] > 9(t — t'). Here [ ] denotes the commutator, <> stands for a thermal average in the grand canonical ensemble over the unperturbed total grand canonical Hamiltonian If: 2 7f — uN (p. is the chemical potential, N is the total particle number operator and usually [77, IV] = 0), and Cx(r,t) : emf/”C(zk'm‘fl‘. For superconductor in a magnetic field of vector potential A(:c,t), the induced current is given by 6j,,(:c,t) = 7,17: fddz' [dt' hezn m Pm,(:c, t; :c' , t') 6(2 — z')6(t — t')6,,,, A,(a:’,t’) . (4.7) The diamagnetic current contribution to the linear response function is given by (hezn/m)6(:c -:c')6(t —t')6,w which is longitudinal to A. The paramagnetic orbital 150 current and spin currnt contribution is contained in the following current-current response function Pm,(:c,t; m',t') = —i < [jp(:c,t) , ju(:c',t')] > 0(t —— t'). (4.8) where the current operator J for itinerant quasiparticles has contributions from jo and 1,, i.e, J 2J0 +j,. (4.9) The orbital paramagnetic current operator is given by 34m) = Z 8"" Human) — (v¢: 9(t — t') , (p,v = :c,y,z, +, —). (4.17) If spin-orbit interaction is neglected from the total Hamiltonian 7'2 of the su- perconductor, the orbital and spin components can be separated. Then P _L(q,w) is the sum of PO,J_(q,w), the orbital contribution and P,,J_(q,w), the spin contri- bution (no cross terms), i.e., PJ-(q’w) : PO.—L(q:w) + P0.-L(wa) = Po..L(q,w) + (CQMB)242X+—(Q:w)- (4-18) 152 However in the presence of the spin-orbit coupling, spin is not a good quantum number and the current-current response function will get contribution from the mixing of spin and orbital currents.81 It can be shown that in the appropriate localized spin approximation, the above result reduces to the familiar spin correlation function70 usually employed in neutron scattering experiments in localized magnetic systems.25'72 The differential scattering cross section integrated over all positive energy transfer can be measured in quasi elastic neutron scattering. From Eq. (4.3), I can define the quasi elastic scattering cross section by 1 ~(q.w) 9m) = F / d(hw)[1 + n3(hw>1[—ImP1(q.w)/q2] - = @001) + .(Q) , (4.19) where F is the constant prefactor in Eq. (4.3) (independent of w), and (I’D, ‘1), are orbital and spin contributions, respectively. The results mentioned above are completely general and independent of the statistics of the (quasi) particles. Formally they have the same appearance for both the paired fermion and the paired spinless charged boson superconductors. In the later case, where the quasiparticles have no spin, the fermion operator #3, should be replaced by the boson operator ,3 without the spin indix 0'. In this case there is no spin current, only the orbital current contributes to the neutron scattering. 153 4.3 Dynamic Response in Paired Fermion Superconductor Form the previous section, we know that the important quantity characteriz- ing the neutron scattering from a superconductor is the current-current response function defined in Eq. (4.8). In this section, I calculate this quantity by using the singlet pairing (s-wave) BCS scheme2 and a standard self-consistent Gor’kov factorization method.°°’°"“"4 For simplicity, I will only consider the leading order Feynman diagrams self-consistently to include the Cooper pair contribution. Thus the ladder diagrams or vertex corrections are ignored.83 This kind of framework has been proved to be quite reliable and has been applied in calculating, for ex- ample, the Raman scattering from the gap excitations in superconductor.85 As already discussed in section 4.1, within the fermion QP models, condensation of paired fermions is still the basis of the superconductivity in these high-Tc ox- ide layered compounds, and the BCS model is the most common paired fermion model appropriate for a diverse variety of pairing forces. This justifis our using BCS model for the purpose of exhibiting the important characteristics of neutron scattering from paired fermion models for these high-Tc oxide superconductors. Consider a many fermion (charge —e, mass m) system in the presence of an applied magnetic field specified by the vector potential A(:c). The total grand canonical Hamiltonian without spin-orbital coupling is given by 15 = 7? - 9N = 5.3/443 1l3.‘f(=v){ [-i-(-ihV + Elf-42 ~u}1/3«(m) 2m +% Effddzddr' v(:c — x') $103) A:.($')tlia'($')§5a(x)- (4°20) 0,0' 154 The pairing interaction of two fermions near the Fermi surface is given by 12, whose microscopic source need not to be specified here. The pairing mechanism may be one or several combinations of those already mentioned in section 4.1. However, I will assume that v is a constant interaction, i.e., v(:c — :c') = —V6(:c — :c')(V > 0) with an energy cutoff hwc. For fermionic superconductors, the chemical potential p is approximately independent of temperature2 and equals to the Fermi energy 5;: = hzkfy/Zm. As will be discussed later, this is one major diference between the fermionic system and bosonic system where the chemical potential depends strongly on temperature and interactions.15'52153 In order to study the finite temperature response functions, I use the Mat- subara Green’s-function formalism.°°’“'” I first briefly introduce some basic formulations of superconductivity to establish our notations. Then I will proceed to calculate the dynamic response functions. To do so I introduce an imaginary time operator by substituting t —» —ir, i.e., 1,5,,(231') = e’ic"'/"zlv°.(:c)e"’tf/h and 1i»:(z‘r) = ek’/"1/3j(m)e’k’/". Note that d:(::) is hermitian conjugate to 190(3), however 1,13: (31') is not hermitian conjugate to 1]},(31'). For any operator C(z), it’s Heisenberg representation is given by 0x037) = emf/”C(x)e‘xr/h, and the ensem- ble average of O is defined by < C >= Tr[e-5kO]/Tre'fik. The single particle Green’s-function is defined by”0 9(3T,:c"r') = — < T, [1,9,(1:T)1lr:(a:'r')] >, (4.21) and two anomalous Green’s functions for superconductor are defined by 77(231', :c'r') = — < T.[,£,(m),£,(z'r')] >, (4.22) and mew = — < 919199919291 >, (4.23) 155 where TT is Wick chronlogical operator, i.e., r is decreasing from left to right. T,- also gives a change of sign (-1) for every permutation of fermion operator. The pairing amplitude (or the ordre parameter) for superconductor is defined by 24(4) = W(m +0;:c,1') = v < 431(4),;44) > . (4.24) These Green’s functions satisfy the well know Gor’kov equations.°°'3""'84 For a uniform system one can introduce the Fourier transforms gm, 2: r) ——Ze -..... (r *'> f—— d 1‘ aim-7") g(k ....) (4.25) :Ii—h (21r)de ’ ’ and f+(23'r, :D T) 1:6 '4“ (T T)/_ dd k ea“: 3) f+(k w”) (426) Zfih (2a)“ ’ where w" = %(2n + 1) (with integer n) is the Matsubara frequency for a fermion. Then 9(k,wn), .7:(k,w,,) and .7-‘+(q,w,.) can be solved from the Gor’kov equations. If A is taken to be real and independent of 2:, the solutions for A = 0 are well known"0 and are given by 9?. vi 9(k w"):z wn—Ek/h+iwn+Ek/h, (4.27) 1 1 _ + .. __ .7:(k,wn) — .7: (k,wn)— — —ukvk wn—Ek/h —iwn+Ek/h . (4.28) In Eqs (4.27) and (4.28), E). = (M: + A2 (£3, = 5;, — p) is the energy of new quasiparticle excitations (bogolonsa7) in the superconductong state, where ck is the band energy; A is also the energy gap. The coefficients upc and v). are the amplitudes of the bogolons and are given by ui=1(1+€—"), v§=1(1—§5—), (4.29) 156 which satisfy A 6!: ukv;B = -2--E—;, u: + v: = 1, u: — v} = 25;. (4.30) g and .7: , f1" can be represented diagrammatically as shown in Fig. 4.1(a) and 4.1(b). The order parameter A is determined by the BCS gap equation anh[ E $th if "lo 0(hwc—lékl)- (4.31) To simplify numerical calculations, in the following I will neglect the detailed band structure of the CuO; planes and use an effective-mass approximation, i.e., 5;. = link2 / 2m which has been found to be quite reasonable. Then in quasi two- dimensions, Eq. (4.31) can be written as ..., tanhkflm] 2 - d ‘ ‘l 5" W ’ (4.32) Q where A = N2(cF)V 2 fig;- is the dimensionless coupling constant. In gereral, the Fermi energy EF is of the order of leV, while the gap A(O) at T ~ 0 is about 10meV in the high-Tc oxides.73 If I take a strong coupling limit A z 1 and require the zero temperature gap A(0) = 0.05cp, then A(T) has to be solved numerically from Eq. (4.32) as a function of T/Tc. This numerical solution will be used in the calculation of the response functions. In order to calculate any correlation function involving two particle coordi- nates, like the current—current response function, I first consider a function defined by Puy(m,r;:c',r') = -— < Tr[ju(mr)jy(:c'r')] > . (4.33) 157 It can be easily shown that the grand canonical average of r-ordered product of an even number of fermion operators or boson operators such as in 73,“, defined above is periodic in each 1' variable with period [371 in the range 0 < 1' < flh. Thus the Fourier representation of 17”.,(21', 2:1") can be introduced as 1 —iu r—r' dd 1' - x—x' ’Puu(x,r;:c',r') = BE 28 "( ) (27731; 8 q( ’Puu(q,vn). (434) where 11,, = 2n1r/flh with integer 11.. Using the Gor’kov factorization method,3°'82 or Feynman diagrams,“ I can calculate 73“,,(sr, :c' 1" ) or 77,“,(q, u") directly. Then through analytic continuation (iun —r w + in), the retarded correlation function P".,(q,w) is related to the thermal function 17”,,(q, 11,.) according to the following prescription” PIW(Q:“’) : puu(QaVn) _ _ ' (4'35) nun—ow-Hn In the remaining part of this section, I will use the above formalism to calculate the orbital current, the spin current and the density response functions of the paired fermion superconductor, and study the neutron scattering characteristics by analysing the behavior of these response functions. Orbital Current Response Function From Eq. (4.10b) and (4.33), the orbital current contribution to Puu(mt,m’t') can be brought into the form (424111 0,0' 158 ><(T.[1133“(2)1/3a(1)ti$(2')1l’av(1')[>} , (4-36) 2=1+O 3’=1'+0 where 1 stands for a space-time point (211,71 ), etc. Using the Gor’kov’s factoriza- tion procedure,84 the above average of r-ordered product of four fermion operators can be factorized into a sum of the products of average of pairs of the field op- erators and in turn can be written in terms of the Green’s functions 9,]: and 1“”. After substituting the Furier representations of 9,]: and .7?+ and doing the differential operations, one readily obtains W.) :42)” / (33,4. 9.9+ 3). x Iii Z [90. + (1,11,, + um)g(k, um) + f-(k + (1,11,, + um).’F+(k, u...)]. (4.37) The factor 2 comes from the spin summation. This can be represented by Feynman diagrams as shown in Fig. 4.2. Using the Green’s-functions given in Eqs (4.27) 3'me and (4.28), the m-sumation can be carried out readily, i.e., Em , e : ,Bhf(:c) zw m—z and f(:c)+f(—:c) = 1, where f(:c) = (613: +1)"1 is the Fermi distribution function. Then (1) shifting the k variable to (k — 52‘) and symmetrizing the new k-integral by adding another integrand of the same form but with k replaced by —k, (2) using the symmetry 6.], 2: 53, (3) regrouping the various terms in the integrand and (4) using Eq. (4.35), I obtain eh 2 am. P04111312”) — h(;) [W kukv 1 §+€-+IAIZ X{§(1+ E+E- )[f(E—)—f(E+)[ 1 1 .(hw+in—E++E_ _hw+in+E+—E_) 159 1 €+€-+|A|2 “ii—(1- E+E_ )[1’f(E-)‘f(E+)] 1 1 .(hw+i77—E+—E__hw+in+E++E_)}a (4.38) where the subindex i means k d: q/Z. The factors $(1 :l: gigfg) in Eq. (4.38) are the superconducting coherence factors. The first term in the big bracket in Eq. (4.38) is the contribution arising from scattering of thermally excitated quasiparticles, whereas the second term is the contribution due to creation or destruction of two quasiparticles. Note that both the normal and the anomalous polarization diagrams in Fig 4.2 contribute to each of these two processes. At zero temperature and A # O, f(E) : 0, thus the first term vanishes, since there are no quasiparticles excitated at zero temperature for A > 0. To get the transverse part, I use the identity Z(6,,., — q“q")k,,k., = M. (4.39) W q2 q2 Thus the orbital contribution to P i(q,w) is given by Eq. (4.38) with the factor Ionic, replaced by gq—Efl, and the real and imaginary parts of P _1_(q, w) are obtained by using the 6-function formula (2: + in)"1 = P(%) —i7r6(:c) where P stands for the principal part. The diamagnetic orbital current does not contribute to the neutron scattering directly, but indirectly through the dynamic screening and has no effect when K. z 1. 160 §pin Current I_Lesponse Function To calculate the transverse spin-spin susceptibility X+-: I first consider a function defined by X+_(x1',z'7") = — < T,[3+(z1')$_(x'r')] >. Using Eq. (4.16) and the Gor’kov factorization, I have x+_(1,1') = —(T.[11(2)2/31(1)¢r(2')¢31(1')]) 1=1+0 2'=l’+0 = [g(1.2')g(1',2> + f+(2',2)f(1.1'>] (4-40) 2=1+o ’ 2'=1'+o where 1 = (231,11). Using the same procedures as in the calculation of the orbital part, I get the final result: d 2 X+_(q,w) = h (31% x -12-- {$(1 + “:2“ )[f(E_) — 1(a)] 1 1 .(hw+in—E++E_ _hw+in+E+-E_) 1 {+5- + IAI2 [ — 1 -— 1 -— E- — E ] +2( E+E.. f( ) f( +) 1 1 . — 4.41 (hw+in—E+—E_ hw+in+E++E_)}’ ( a) and P.,1(q,w) = (cgu3)292x+—(q,w)- (4-41b) Note that for an isotropic s-wave superconductor, x“ = ny = x" = %X+-- Eq. (4.41) together with Eq (4.38) completes our theoretical formula for the transverse current-current response function whose imaginary part gives the inelastic neutron scattering from the gap excitations of the paired fermion superconductors. 161 Density-Density Response Function The external electric field (or scalar potential) couples to superconductor via the local charge density. The induced local charge density is contained in the retarded density-density response function86 D(mt,:c't') = —i < [6fi(a:t),6fi(:c't')] > 9(t — t'), (4.42) where 613(3) 2 11(3) — n(:c) is the density fluctuation operator, fi(:c) = Z,1&j(m)1/i,(z) is the density operator, n(:c) =< 71(2) >= —2§(:c,‘r;a:,‘r + 0) is the average density. The imaginary part of the Fourier transform of the density- density response function I mD(q,w) is related to the inelastic electron scattering cross section from a superconductor.80 Let us first consider the thermal conterpart D(l, 1' ) [here 1 = ($1,T1)] and use the Gor’kov factorization to get D(m') = - EXT. Momma(2016.41)]> - n<1)n(1') ma’ 3=1+0 2' =1'+0 = 2[g(1,2')g(1',2) — J-‘+(2',2)}‘(1,1')] (4.43) 2=1+o ‘ 3'21, +0 Fllowing the same procedures used in obtaining PM", and PM 1, I obtain ' d _ 2 D(q,w) =1: (37}; {$(1+5+§+E'_A' )[f(E—)-f(E+)] 1 1 .(hw+in—E++E_ —hw+in+E+—E_) 1 5+§-—|A|2 +§(1—— E+E_ )[1—f(E_)—i(E+)] 162 1 1 ' " . 4.44 (W+in-E+—E_ hw+in+E++E_)} ( ) Note that this result corresponds to the two ring diagrams in Fig. 4.2 with ap- propriate vertex functions. The coherence factors In Eq. (4.44) are different from those in the current-current response function [Eq. (4.38) and (4.41)] by a sign in front of |A|2. In addition, correction to Eq. (4.44) due to the Coloumb interaction between the fermions can be obtained by using Dyson equation.80 This interaction also produces collective modes propagating in the superconductor. 163 I ’ g($171,$272) f+(3317'1,$27'2) H f(317'1,$27'2) k,w,, ' _’— 9(k,wn) kawn I 7+(k,wn) kawn ’ ‘— 70902..) Fig. 4.1 (a) Diagrams representing the thermal Green’s functions in paired fermion superconductor. The inward arrow at one end represents destruction (1b) of electron while the outward arrow represents creation ((1;+ ) of electron. Thus f+ represents creation of electrons at the two ends. (b) Diagrams of Green’s functions in Fourier space for paired fermion superconductor where each line is labeled by a momentum k and frequency w". .7:+(k,wn) and F(k,w,,) represent creation and destration of Cooper pairs. 164 k+q Vn‘l'Vm k+q Vn+Vm Fig. 4.2 Feynman diagrams that contribute to the response functions within the Gor’kov factorization scheme. 165 Norrrgl State Limit; In order to distinguish between the characteristics of scatterings from the normal and the superconducting states, I first give the results for the normal state (A = 0). In this limit, u), = 0(6),),121, = 0(—£k). The coherence factor -1- 1 + 6+6- = 1 if k and k- are on the same side of the Fermi surface and 2 I£+£ l + {+5-— 1 zero otherwise, and — (1 —- 2 |€+€—| the Fermi surface and zero otherwise. After considerably long algebra, from Eqs ) = 1 if It... and k- are on the different sides of (4.38) and (4.41), I obtain Po,uu(chw)} _ h ddk (512)20‘ + 3),.(1‘ + 21);; _ 1 X+-(Q1w) (2W)d 5 X2 “5km “- f(€k+q)] 1 1 (hw+iTI—Ek+q+£k _hw+in+ek+q_5k)’ (4'45) and D(q,w) = 4xu(q,w). We see that when A = 0, Eqs (4.38) and (4.41) properly go over to the normal fermi liquid RPA (random-phase-approximation) response functions corresponding to creation of electron hole pairs.5“'80 At zero temperature, f(£,.) = 0(kF — Ikl), the k-integral in Eq. (4.45) can be carried out exactly. In polar coordinates, the angular integral can be done first, leaving the IkI-integral from km," to kp. Then by considering positions of the two Fermi circles (d = 2) or spheres (d = 3) with their centers separated by q and relations of the two vectors k and k + q with respect to the two Fermi surfaces, together with the 6-function 6(hw — ek+q + 5),), the lower limit of the k-integral can be figured out easily for several regions of (q,w) values and are listed in the Table 166 4.1. The calculated imaginary parts of the response functions for the normal fermi liquid are listed in the Table 4.2 for d = 2 and in the Table 4.3 for d = 3 in units of h(cgp.3)2Nd(5p)lc%.. We see that at small q, the spin and orbital contributions 1 3 to the inelastic scattering cross section are proportional to q“ and q“ , respec- tively. A comparison of the numerical values of I mPo, _(_(q,w), I mp“ _L(q,w) and I mD(q,w) for the normal fermi liquid are given in figures 4.2-4.9. We see that when w —> 0, they all go to zero. There is a q-dependent energy gap for q > 2kp and no gap for q < 2kp. Also there is a q-dependent maximum (at hw ~ 5p) in both the spin and orbital contributions. At higher energies (hw ~ 105p), all the responses go to zero. These features can show up in the inelastic neutron scatter- ing from normal metals. In the superconducting state, since A is usually much smaller than 613', at high energy the responses from the superconducting state will be very close to their corresponding normal state values for hw >> A. Thus any significant diference between the A = O and the A 74 0 results will appear only in the low energy and/ or small q regimes. I will discuss this feature in the later part of this section again. Using expressions in the Table 4.2 and 4.3, the quasi elastic scattering [Eq. (4.19)] can aso be calculated exactly. For three-dimension normal fermi liquid at zero temperature, I obtain aways)? + 2(4)? — 4], for q 3 2km «10(q) = 50 x (4.46) 9.45:)”. for q 2 21w ; and ”(Mil _ 2-(.3_)2], for q 3 21cm 2 1 (4.47) 2n, for q 2 WW; 167 and the total quasi elastic scattering cross section is given by §(§;)[§+2(£§)2—%(f;)2]. forqSZkF; M®=&x ' ma) 1.1% + §(*—:—)’]. for q 2 21m; where m mlc p 50 Z Fh(Ch/m)2Nd(€F)€F; N2(EF) : m, N3(€F) = m. Eq. (4.48) agrees with the results given by Lovesey70 and is a check on our A —+ 0 limiting behavious. As q —i O, the spin part goes to zero, but the orbital part goes to infinity, i.e., : LF ‘I’o(.(q) 5(5) Thus for a normal fermi liquid, the orbital response displays divergence for small wave vectors. This divergence is removed upon the proper use of the electrody- namic screening.77 For electron system, the electron-electron Coulomb interaction will renormal- ize the above results. The spin part acquires the well-known enhancement factor whereas the orbital part is suppressed-’9'81 This may be understood on physical grounds by noting that the mutual electron interactions tend to enhance the spin correlations while suppressing the orbital motions of the itinerant electrons, i.e., tend to ”localize” the electrons. However this is not necessarily true if one takes account of other effects such as band structure and correlation effects beyond the Coulomb interactions.“ 168 Table 4.1 Regions of q,w and limits of the k-integral in calculating the normal fermi liquid response functions. Here q and km... are in units of lap, 11 = 5%.. q V kmin I 52 05usl§q’-ql \/1-2V 11 _<_? l%q’-q|SVS§-q’+q V/q--§-q III >2 4q2-qs»s;q2+q u/ -1. 169 Table 4.2 Imaginary parts of the response functions for a three-dimensional isotropic normal fermi liquid at zero temperature. Here q is in unit of hp and u = 2t—i' _ImPO.L(Q1w) -ImP,,i(q,w) I 1%[1— V " (3" " fly] §wqu ‘l’ u 2 3 it v I" 111 Ell—(E‘iql l Ell“(;-%Q)2] 170 Table 4.3 Imaginary parts of the response functions for a two dimensional isotropic normal fermi liquid at zero temperature. Here q is in unit of kp and V = %. -ImPo,J_(q,w) I %{[1—(§‘ iflzlm- [1-(§-+-§-q)']m} 3/2 11,111 _3_[1_(z__1q)2] —ImP,HL(q,w) I q{[l’(i‘i")’lm-[1—(%+%q>’]m} [/3 11,111 q[1_(%_%q)2] 171 1.5 T U I I I 1 T r I I I i r U I I I l I . A = 0.0 (1) q/k, = 0.5 (2) Q/kr " ~ (3) q/k, = 1.0 ~ 1.0 P' (4) Q/kr = 1.2 _ l (5) q/k, = 1.5 . 2 . .. g ' . . . 0.5 r- "\ 5 ~ I P a: —ImPo(q, w) 0.0111111mLJ .11 1;; O 1 2 3 4 hes/6F Fig. 4.3 Imaginary parts of the transverse orbital current-current response functions for a two dimensional normal fermi liquid at zero temperature in units of h(cgp3)2N2(ep)lc%.. 172 .r..l....l.r..l....l... . d=2 . _ A=0.0 (1)Q/k'=0.5 * (2)q/k,=o.a r- (4)Q/k'=1.2 (5)Q/k'=1.5 A .. 3 6* ” 5 1 ‘27, . 0.. E 1r- 4 '1 ' ' 3 F 2 ' 1 r/ I 0 ..él ...! . . hes/6F Fig. 4.4 Imaginary parts of the transverse spin current response functions for a two dimensional normal fermi liquid at zero temperature in units of h(cgp3)3Nz(ep)lc}. 173 )- L I d A 8 Q ’7 " l VjIWYI‘l’jll’Ur'rl d = 3 A = 0.0 (1) Q/kr = 0-5 (2) cm... = 0.3 "' 1.0 — 1.2 2 (5) Q/kr = 1-5 A 4': Q R“ '1 l "“ . 1.5 . 1.0 - ’5 6* ‘6 ii T 0.5 - 0.0 L 0 Tim/6F Fig. 4.5 Imaginary parts of the transverse orbital current-current response functions for a three-dimensional normal fermi liquid at zero temperature in units of h(cgp3)2N3(ep)lc%.. 174 1.5 ' r ' ' I I ' U l I I I r I " d = 3 .. . A — 0.0 , P (1) hU/EF = 0.5 I (2) hrs/e,- = 1.0 l b (3) rim/E’- = 1.5 l 1.0 '_ ._( 3 _ 2 ‘ d" P 1 ‘6 '- .1 0... g - . I 0.5 r- LLL 0.01 *1 Q/kr Fig. 4.6 Imaginary parts of the transverse current-current response functions for a three-dimensional normal fermi liquid at zero temperature in units of h(691113):1‘73(¢‘=F)ki=-- 175 .. .1....[. .rr. .1.. " d - 3 2.0 "‘ A = 0.0 _ : (1) Q/kr = 0-5 i . (3) Q/kr =1-0 4 1.5 L. (4) q/kp = 1.2 _g - (5) Q/kr =1-5 3 i (5* * 1 ‘7: 1.0 ‘- _ D-a ~ . i i 5 . .. , 4 . 005 i—- -—-J .. . 3 q r 2 . l .( 0.0 A J. 1 LL 1 1 1 hos/eF Fig. 4.7 Imaginary parts of the transverse spin current response functions for a three-dimensional normal fermi liquid at zero temperature in units of h(C9l‘B)2N3(5F)ki='° 176 P T l ' FT . ’ (a) = 3 ‘ 0.9 " A a 0.0 — : q/k, . 2.5 < A - '1 a 04 L- _ c7. " .4 ‘5 4 0.. . E l 0.2 (b) o 1 - n l A - . . l o s 10 15 flu/e, Fig. 4.8 (a). Imaginary parts of the transverse orbital current-current re- sponse functions for a three-dimensional normal fermi liquid at q = 2.51:; and zero temperature in units of h(cgu3)’N3(ep-)k%~. (b). Imaginary parts of the transverse spin current response functions for a three-dimensional normal fermi liquid at q = 2.5kp and zero temperature in units of h(¢9#s)’N3(¢r)ki=~ Fig. 4.9 Imaginary parts of the density-density response function for a two- 1 l L 4 1 .J. 1 had/6F 1.5 177 dimensional normal fermi liquid at q = 0.5kp and at zero temperature in units of th(ep). 178 3fi..rr. ‘I""l“ - d=3 L- A-0.0 (1) Q/kr‘0-5 . (2) Cl/kI-‘=O'8 ’5 U.“ '1 ., .... 4 I 1 .“ .1 5 1 g 4 0 1 11J 1 1 1 1 - 2 3 4 hes/6F Fig. 4.10 Imaginary parts of the density-density response functions for a three-dimensional normal fermi liquid at zero temperature in units of hN3(€F). Results for Superconducting State 179 When T < T6 and A 75 O, the integrations required in Eqs (4.38) and (4.41) are quite complicated. In order to investigate the effects of superconducting cor- relation on the neutron scattering, I first give the approximate results for the dynamic response functions at small q. Then I give results for arbitary q and zero temperatures from an extensive series of numerical calculations. Use the following expressions hz £k+q :£k+5q + —k'q, m h” 5,. h’ A” VkaEk = —— + (—-)2 —kk, m m E13. I expand Ek+q and f(E+), etc. to second order in q. For example, Bk 2 h2 {k 1 ’12 2A2 Ek+q : Ek + [Eq 'l' :n—k "1] —“ + —(';n-) EEO‘ ' (1)2- Then to second order in q, I have #5,. E —E_=——- + THE], q, 1 h’ 2A2 E+ + E- = ZEk + Zeq/z—é—k- + _<_) __(k 'Q)2, Bk 4 m E13: 2 2 2 2 flzhl E. A_(k.q)2, E+E_ 2 m E: and 32f(Ek) 3E; f(E+)-f(E-) _ 61w.) +1(§:)’A_2 E+—E_ " 8B,. 3 m E}; (k - <1)" (4.50) (4.51) (4.52) (4.53) (4.54) 180 83f(E,,) 6E: (k ' (1)2: (4'55) 1.1111: 24 m E: —f(E+)- f—(E )=g(E.)+— ;(’,"—,,-)1"‘2E,i 39$) (1..)2 +35%) 5:"; 3;?) (1.12, <4“) where 9(5),.) = 1 — 2f(Ek) = tanh(%flEk) . (4.57) These relations enable me to write the response functions as ( X 1‘)2 {P0,1(Qaw)} : h(f§)2 ddk 371-2—— 4 P:,J_(q1w) m (27") g-qz ><§{( )flf(Eu [1-f(Ee)] :(k ”(11.14... _2E++E_) +102) «Eu-E734 q)’(n..+.-.,__1E1_E )}. (4.58) To first order in q, the first term inside the big bracket vanishes because (k - q) or (k . q)(k x q)2 integrated over angle between R and q give zero. The first term will contribute to second order in q. For T ;E 0, this contribution is due to thermally excited quasiparticles. At very low temperatures, the derivative of the fermi distribution function vanishes, the first term in the big bracket is close to zero as T —> 0 and only the second term mainly contributes with g(Ek) = 1 — 2e‘3Eh z 181 1 — 26’5““. I can put E+ = E- = E]. in the denominator of the second term, then the required k-integral can be proformed. It can be shown that for T << Tc and as q —+ 0, P0,J_(q,w) and P.,i(q,w) are functions of the variables qé and (hw — 2A) / A only, where 6 = Fwy / A is the coherence length2 which represents the size of a Cooper pair (up = hkp/m). Also ImP_L(q,w) = 0 for flu: < 2A, i.e., one has a threshold behavior at twice the gap energy. However, the gap will be filled by thermal excitation of the quasiparticles at T > 0 coming primarily from the first term. I have q k 2 { ‘ImPOvi(q’w)} 1 d‘k (‘1; x 12;) : N(e )e (21r)d 1 q 2 —I P, ,w d F F _ _ m .1(q ) 2(kF) A2 3 k 2 Xg(E],) EiF( k2q) 6(hw - 2E1.) (4.59) k F where I mPo and I mP, are expressed in units of h(cgp3)2Nd(ep)lc§.-.. Note that the coherence length 5 defined here is wfgcs, where {gas is the coherence length defined in the BCS theory.2 For general wave vector and temperature, I can calculate Po,_)_(q,w) and P,,_L(q,w) from Eqs. (4.38) and (4.41) by doing numerical intergations. The nu- merically calculated results in the limit T = 0 are shown in figures 4.11-4.15. Figure 4.11 shows the imaginary parts of .Po,J_(q,¢.v)/q2 and P,,i(q,w)/q2 for a quasi two- dimensional superconductor [in units of h(gn3)2N2(5p)lc%~] with A(0) 2 0.0555» and for q/kp = 0.5. Figure 4.12 shows the same quantity but with q/lcp = 1.5. For a three-dimensional nearely isotropic system like Ba1-,K,Bi03, the calcu- lated results are given in figures 4.13 and 4.14. From these figures one can see the threshold behavior for inelastic scattering at hw 2: 2A. Another important 182 feature of these numerical results is that there is a large reduction in the inelastic scattering at small w in the superconducting state compared to the normal state, in both the orbital part and the spin part. For example, large reduction of the scattering in the energy range 0 < hw < 6A is clearly seen from Figs 4.11 and 4.12. This is quite smilar to that observed in the Raman scattering experiments in the high Tc superconductors YBazCU307 by Cooper et. (11.73 For q/lcp < 1, the orbital contribution dominates the spin contribution by orders of magnitude and should be observable in experiments.“'“ The quasielastic scattering cross section also gets drastically modified in the superconducting state in the small q regime. Since A provides a lower energy cutoff, the corresponding q cutoff would be given by hqu ~ A. Clearly the coherence length 5 = 7112;» / A provides a lower cutoff for q”. Thus the superconducting correlation removes the divergence in the normal state orbital scattering. In fact, from Eq. (4.59) I have xk2 (Mg) _ 5'0 _d‘_k (144—)— x (Effigylbql T49) _Nd(5F)5F (27!)! 1 g k 2 The above integrals can be calculated exactly at T = 0. I obtain for a three- dimensional isotropic system like Ba1_,K,Bi03 , 443—» 5365;) (51.1), (440), (4.61) and 1r 4 4.(q) a 6—0-(441)[1 + (cwl’ (q —4 0), (4.62) (in units of So). In particular, as q —+ 0, the normal state spin contribution (1r/2)(q/lcp) becomes (7r/24)(q/kp)2(§lcp). The divergence of the normal state 183 orbital contribution (7r/2)(kF/q) is removed and becomes (1r/60)(§kp). These results are summarised in the Table 4.4. In Fig. 4.16, I plot the quasi elastic orbital scattering as function of (q/lcp). In Fig. 4.17, I plot in as functions of (qfi). Again one sees the orbital response dominates for small q and it should be possible to see the effect of the superconducting gap in the neutron scattering experiments. Also one should be able in principle to obtain the coherence length from the small q limit of the orbital response. Finally I plot in Fig. (4.18) (In, and ‘1’, for a quasi two-dimensional system. The reduction in the total scattering cross section in the superconducting state is proportional to the factor N (cp)W where W is the width of the band over which the attractive interaction is operative. Typically for a traditional low temperature superconductor such as lead, N (e F) z a few states / eV, W x .10 meV (phonon energy), the gap A z 0.8 meV, and the coherence length 6 z 500021. In contrast, for the high Tc system, N (e p) is about the same, W is perhaps electronic in origin and is of order 100 meV, A is estimated to be of order 8 — 32 meV, and 6 z 10 — 15A. Thus the scattering strength from the superconducting pairs is likely to be larger by an order of magnitude in the oxide superconductors. From these estimates, it seems that the neutron scattering experiments may yield valuable infromation about A and C for the new high T; system as the energy and momentum parameters are in the accessible experimental range. 184 N —ImPo(q. w)/q2 ...a L L L 0 O 2 4 6 8 10 hw/A Fig. 4.11 Orbital contribution to the inelastic neutron scattering from a quasi two-dimensional paired fermion superconductor (solid line) at zero tem- perature with A(O) = 0.05:5- and q = 0.51:;- as a function of flu: / A in units of Mega)” N3(ep). The dashed line represents the normal state result. 185 08 I I I I D (lOSew 0.6- q=o.5k,. J ~1mPs(q. w)/ <12 00 l l l 0 2 4 6 8 10 hw/A Fig. 4.12 Spin contribution to the inelastic neutron scattering from a quasi two-dimensional paired fermion superconductor (solid line) with A(O) = 0.05”», q = 0.51:1: and at zero temperature as a function of flu: / A in units of Mcgug)2 N3(ep). The dashed line represents the normal state result. 186 0.10 .- d H N I 0.00 l ' . 0 5 10 hw/A Fig. 4.13 (a) Orbital contribution to the inelastic neutron scattering from a quasi two-dimensional paired fermion superconductor (solid lines) at zero temperature with A(O) = 0.05“: and q = 1.5kp as a function of hw/ A for different q/lcp. The dashed line represents the normal state result. 187 00 0.1 *- 5 10 hw/A Fig. 4.13 (b) Spin contribution to the inelastic neutron scattering from a quasi two—dimensional paired fermion superconductor (solid lines) at zero tem- perature with A(O) = 0.05cp and q = 1.5kp as a function of hw/A for differ- ent q/kp. The dashed line represents the normal state result. 188 1.0 0.5 mot °'°2 A =0.05 2,, ‘ 0.01 _Impo(qa Gil/(12 0.00L fun/2A Fig. 4.14 Orbital contribution to the inelastic neutron scattering from the three-dimensional s-wave paired fermion superconductor (solid lines) at zero temperature with A(O) = 0.05ep as functions of Fun/2A for different q/kp, in units of Magma)2 N3(ep). The dashed lines represent the normal state results. 0.15 0.10 “ImPJq. 0)/ q2 0.05 fun/2A Fig. 4.15 Spin contribution to the inelastic neutron scattering from a three- dimensional s-wave paired fermion superconductor (solid lines) at zero tem- perature with A(O) = 0.05cp as functions of hw/ZA for different q/kp. The dashed lines represent the normal state results. 189 190 d: 3 - Normal State l't A A = 0.05 5' A=O.I Er 4’0“!) Q/kr Fig. 4.16 Orbital contributions to the quasi elastic neutron scattering cross section for a three-dimensional s-wave paired fermion superconductor at zero temperature with A(O) = 0.055;: and 0.1:; as functions of q/lcp in units of 1rSo. A(O) = 0 corresponds to the normal state. ¢S(Q) 0.10 0.00 191 I r I I, A = 0.0], 5' [_z‘ 2, A = 0.05 5', ,~’. 3, A =3 0.1 6’. .’.\/ .4 d :3 3x .’.’. ’\2 i / ’. ...................... V...- Fig. 4.17 Spin current contributions to the quasi elastic neutron scattering cross section for a three-dimensional s-wave paired fermion superconductor at zero temperature with A(O) = 0.05p and 0.15p as functions of {q in units of «So. A(0) 2 0 corresponds to the normal state. d : 2 1, Normal State 2, A = 0.05 e;- 3, A = 0.1 65- 10 - '. a «2’2 95001) O l l l _ 0 0.2 0 4 0.6 0.8 Q/kF Fig. 4.18 (a) Orbital contributions to the quasi elastic neutron scattering cross section for a quasi two-dimensional paired fermion superconductor at zero temperature in units of 1rSo. 192 1, A L. 0.4 2, A d: ¢s(Q) 0.0 L Fig. 4.18 (b) Spin contributions to the quasi elastic neutron scattering cross section for a quasi two-dimensional paired fermion superconductor at zero temperature in units of «So. qé 193 194 Table 4.4 Comparison of quasi elastic scattering cross section at small q for a three dimension normal Fermi liquid and a three-dimensional s-wave paired fermion superconductor at zero temperature. (In, and 4’, are expressed in units of So, 5 is the choerence length. Normal state Superconducting state '9' o A Q v an A e I? V 5% {ks as)?» '6' . A .Q V up. A ~71» v 195 4.4 Summary In summary, I have studied the dynamic response functions and neutron scat- tering characteristics of paired fermion superconductors. The effects of supercon- ducting correlations appear both in the orbital and spin parts. Both orbital and spin currents associated with the Cooper pairs make equally important contribu- tions to the scattering, the former dominating the latter in the regime of small wave vector transfer. The calculation presented here is within the BCS scheme but appropriate estimation of the effects are made for the new high-Tc oxide su- perconductors using the accepted features of the new superconductors. I have shown that a direct experimental estimation of the superconducting gap A and the coherence length E at low temperatures can be made by inelastic and quasi elastic neutron scattering experiments in the new high Tc superconductors. In the above study, I have assumed that the pairing state has the s-wave sym- metry and is a spin singlet. As discussed in section 4.1, experimental observations31 and theoretical arguments“9 suggest that the superconductivity in the high-Tc ox- ide superconductors is in the s-state. However, the calculation procedures and results given in this chapter can be generalized to superconductors with other pairing symmetry or anisotropic states [where A(k) depends on the direction of k] 9 or spin-triplet p-wave state as in superfluid like heavy-electron superconductors,” 3He (Ref. 90). Recently Joynt and Rice91 examined theoretically the neutron scattering from anisotropic heavy fermion superconductors using only the spin sus- ceptibility. They have found that for a singlet state, the spin susceptibilities are suppressed from the corresponding normal state values. Their analysis is a special case of our general expression when a spin-only approximation is employed. In this case the neutron scattering cross section is directly related to the imaginary part 196 of the frequency and wave-vector dependent spin susceptibility [see Eq. (4.3)]. As I have discussed in this chapter, the complete expression for neutron scattering has orbital contribution which dominates the spin contribution in the regime of small momentum transfer. The expressions given in Eqs (4.38), (4.41) and (4.44) can be generalized to include several bands. In YBantt307, where one has chains and planes contributing to the superconducting property, a corresponding two—band scheme may be appropriate.50 In that case, one obtains intra- and inter-band contri- butions even at T = OK. These would involve quasi-particle excitations above 2A1,2A2,A1 + A2, and IA] — Azl, where A1 and A2' are the relevent gaps as- sociated with the chains and planes. For T yé 0K, there is the usual thermal smearing of the effects mentioned above but also other processes contribute. How- ever the gap structure will still be exhibited. It is therefore preferable to perform the experiment at low enough temperatures to probe at these superconducting parameters. 197 Chapter 5 Response Function Characteristics of Paired Charged Boson Superconductor 5.1 Introduction As I have already discussed in chapter 4, section 4.1, one of the theoretical models suggested to explain the unusual properties of the high-Tc oxide super- conductors is the paired charged spinless boson model.5l"53 In this model the individual quasiparticles which attract to form the pairs are spinless bosons. In contrast, in the fermion models, the charge carriers which form pairs are either 1 holes or electrons. These are fermions with charge iIeI and spin 5. In chapter 4, I have already investigated the dynamic and quasi elastic neu- tron scattering properties of a fermion superconductor to see how these response functions can provide a direct experimental estimate of the superconducting gap A and the coherence length £. There I showed that the scattering from the orbital currents can be quite important and can dominate the spin current contribution in the small momentum transfer regime. In the boson models, where the elementary 198 excitations do not carry spin degree of freedom, only the orbital currents associ- ated with these excitations can couple to the neutrons. It is therefore important to explore how the neutron scattering from these orbital curents change when the bosons pair in the superconducting state. Furthermore, the nature of pair correlation and the dynamic response of the two types of superconductors discussed above are likely to be vastly different. In this chapter, a careful study of the dynamic response of a charged boson supercon- ductor is carried out. In particular, I have studied the dynamic current—current and density—density response functions since the imaginary parts of them are re- lated to inelastic neutron and electron scattering cross sections, respectively. In addition, I will compare these results with those for the fermion system given in the previous chapter. The remaining part of this chapter is organized as follows. In section 5.2, I introduce the spinless charged boson model and discuss the nature of quasiparticles and their energy spectrum. In section 5.3, I study the dynamic current-current and density-density response functions and discuss the characteristics of inelastic neutron scattering in this model. In section 5.4, I discuss the quasi elastic neutron scattering and application to the high temperature oxide superconductors. Finally in section 5.5, I give a summary. 199 5.2. Paired Spinless Charged Boson Superconductivity Consider a spinless charged boson gas with attractive interaction. Further- more, the bosons move only in two spatial dimensions.51 The Hamiltonian for this system is similar to Eq. (4.20) with the fermion field operator 1/3, replaced by the boson field operator cf) without the spin indix. In momentum space, this Hamiltonian can be written as . 1 _. 2 : + E : r + + ’C — k (6k — [1.)bkbk - 5 k k, l’ bk b_kb_krbkr, (5.1) where b: (bk) are the boson creation (destration) operators and V(> 0) is accractive (assumed to be a constant) pair interaction. In Eq. (5.1), ck = £225.: is the boson kinetic energy and m is the effective mass of the bosons. The thermal Green’s functions for the boson system are similar to those de- fined in Eqs. (4.21)-(4.23). From the corresponding Gor’kov equations, I obtain the solutions for these thermal Green’s function as "3. vi 93(k,wn) 2 icon — Ek/h — iwn + Ek/h’ (5'2) 733(k w") = .7:+(k 0),.) = —ukvk 1 — 1 . (5.3) ’ B ’ iwn—Ek/h iwn+Ek/h where w” = 2n1r/ fih. The new quasiparticles in the superconducting phase have Ek = V512. — |A|2, (5-4) energy 200 where £1, = 6;. — p, u(< 0) is the chemical potential and A is the superconducting order parameter. In Eqs. (5.2) and (5.3), It]. and v]. are the amplitudes of the new quasiparticles in the superconducting state and are given by 1 5k) 1( £1.) 1"“ 2<1+Ek’v'° 2 +E,.’ ( ) which satisfy A El. 11.51)]; = m, It: - 1)::— 1, 11.: + I); = E. (5.6) The temperature-dependent parameters A and p. are determined from the following two equations, mamas-:1. k and n = Z Elk—{£k[n3(Ek) + g] - $19k}, (5.8) k where n is the density of the bosons and-113(E) = (653 — 1)‘l . The R summation in Eq. (5.7) is cutoff at an upper limit A (which also defines an energy 6A = h2A2 / 2m). One of the fundamental differences between the fermion and boson systems is the strong temperature and V dependence of the chemical potential It. For the boson system, one can define a wave-vector scale k1, = (41I'n)l/2 and an energy scale kBTo = hzkg / 2m. In addition, one can define a dimensionless coupling constant /\ = fiyV. In two dimensions, when lul > [A], the k—integral 2 in Eq. (5.8) can be carried out explicitly to give ZkBTo — |p| = —%ln [23inh('—g- #2 — A2)], (5.9) and Eq. (5.7) can be written as ‘A 2 1 ,6 — = d —- th(—E ). 5.10 A / 51‘ ER co 2 k ( ) 0 However, if |A| = lul, the divergent parts in the above integrals must be separated out. In particular I have 2’i’BTo — IMI = —\/M2 - A2 + dy 8&1, (5.11) / W and 2 _ —1 lfll+5A —1 WI A—cosh [ A ] cosh [A 3\/(|#|+¢A)'-A’ 2 + / d2: (6’ _ ”“32 + (HAW . (5.12) (5 u’-A’ —-1 It can be shown52 that at zero temperature, when A < Ac 2 [ain’t-“fin z 1.03, the gap E9 = VIM? — |A|2 is zero whereas when A > AC, 13,, > 0. At nozero temperatures, there is always a gap in the excitation spectrum.53 Figure 5.1 shows the zero temperature chemical potential [1(0) and the superconducting order parameter A(O) as functions of A. A comparison of the basic properties of ‘ paired fermion and paired charged boson superconductors is given in the table 5.1. A fundamental difference75 between the two models is that for fermions, the gap E9 in the quasiparticle spectrum is identical to IAI, the superconducting order parameter whereas for the bosons, the gap depends both on the chemical potential and the superconducting order parameter. 202 Table 5.1 Comparison of properties of the paired fermion and the paired charged boson superconductors. fermion boson Chemical potential p 1:: at: p S 0 Quasiparticle energy (M; + A2 {i — A2 Gap Eg A Vfl! — Ai u: we) Mus) m. we) %(-1+‘:) ui+vg 1 g: u,—v,g a. 1 203 A (Cl/kaTo Fig. 5.1 The chemical potential [.1 and the superconducting order parame- ter A at zero temperature for a paired boson superconductor as functions of coupling constant A. 204 5.3 Dynamic Response Functions To study the dynamic response of the paired boson supercOnductor, I have cal- culated both the current-current response function Pm,(:ct; z't') and the density— density response function D(zt;:c't') as already defined in chapter 4. These re- sponse functions can be obtained by using the corresponding thermal Green’s functions and the Gor’kov factorization method.80 I obtain the spatial and tem- poral Fourier transform of the current-current response function as 7’uu(q.vn) = (33)“ 20‘ + 52!)"(1‘ + 3),, k Egg—9.0. + q,Vn+Vm)QB(kan)+-7:B(k + q.u.+vm)f; 0, there is a q-dependent threshold at energy 2Eq/2 = 2\/(5q/2 _ “)2 _ IA|2 : 2\/(5q/2 )2 + Zluleq/Za (519) - due to the 6-function 6(hw—E... —E_) in the imaginary part of Pi(q,w). Note that in the paired fermion superconductor, the threshold is always at twice the energy gap 2A for all q < 2165'. For the paired boson superconductor, for energies below 2E9”, I mP i(q,w) = 0. Just above this threshold, there is a large and broad peak in I mP _L(q,w). Note that the q x R factor in the transverse response function cancels the coherence factor singularity (i.e., when Ekiq/z = 0). The imaginary part of the response function decays rather rapidly with w and goes to zero at high energies. Also the peak position moves to the higher energies and the peak broad- ens as q increases. These features will show up in the inelastic neutron scattering experiment from superconducting boson systems. At finite temperatures, the gap in I mP _L(q,w) will be filled up by thermally excited quasiparticles. However, at low temperatures, well below the Superconducting transition temperature Tc, the results will not be quite different from those given above. For I mD(q,w), three kinds of factors namely the coherence factors, the Bose distribution functions and the 6—functions control the q and or dependence. Let us first look at the Bose functions. In contrast to the superconducting fermion system, all the terms involving Bose functions at Ekiq/Z = 0 have to be seperated 209 out. The value of this condensate contribution should be obtained through the self-consistent particle density equation [Eq. (5.8)] and the gap equation [Eq. (5.7)]. Secondly, the singularity of the coherence factors will produce a logarithmic singularity in D(q,w) at energy hw = Eq. Note that Eq is large than the threshold energy 2Eq/2 in the gapless regime. In Fig. 5.4, I give the imaginary part of the density-density response function D(q,w) for A = 0.5. at zero temperature. As expected from the above disscussions, in addition to the q—dependent threshold at 2E9”, there is a singular peak at a higher energy value Eq. These features will show up in the inelastic electron scattering from the gap excitations of the charged boson superconductors.15 0.015 , I d: A=O.5 1,q=0.1kb 2,q=0.5 kb 3,q=1.0 kb ’ 4,q=1.5 kb 0.0101- _ on 6' 2 3 é a. .3. I 0.005- 4 - 000i L I 1‘ 0' O 4 6 8 fun/A 10 210 Fig. 5.3 Imaginary parts of the zero temperature transverse current-current response functions divided by q2 in the superconducting state for a paired bo- son superconductor in units of h( gym 10 —ImD(q, w) ..l A 1q=0.5kb 2q=1.0kb I 3,q=1.5kb 4q d 2.0 k, Tim/A 211 Fig. 5.4 Imaginary parts of the zero temperature transverse density-density response functions divided by q2 in the superconducting state for a paired bo- son superconductor in units of hNg(1/ 7r) . 212 5.4 Quasielastic Scattering Using the results obtained in the last section, I can also obtain the quasi elas- tic neutron scattering cross section defined by Eq. (4.19). Figure 5.5 gives the quasielastic scattering cross section at zero temperature as functions of q in both the gapless and the finite gap regime. We see that @(q) decreases rather slowly with increasing q. This is quite different from the paired fermion superconductors where the orbital scattering dominated at the small q regime. The supercon- ducting correlations rcmoved the divergence of (q) at small q occuring for the normal fermion state and lead tD(q) to decay rapidly with increasing q. These differences provide important experimentally observable signatures to differentiate paired fermion from the paired boson models. In addition, for q ~ 0, (q) goes to a A—dependent constant. Its value can be obtained from a small-q expansion of the coherence factor in Eq. (5.9). Using the identities Eqs (4.50)-(4.54) with _the substitutiOn A —+ iA, I obtain the q = 0 limit of the boson quasi elastic response <1>(0) = I. l’” {1 .. _———"‘l2 _ Wimp“ ‘ 'A' } (5.19) ZlflllAl lul + IAI [in units of Fh(%)2NngTo]. Thus in the gapless regime where Ipl = IAI, the zero momentum transfer quasi elastic scattering is directly proportional to IA]. The inset in Fig. 5.5 gives (0) as a function of the coupling constant A. One sees that it increases rapidly with A and saturates at A ~ 2Ac. In the new high temperature oxide superconductors, the important electrons are those belonging to the 0110; plane. The holes or electrons created by doping the 0110; plane are the charge carriers. From various experimental measurments such as AC conductivity and specific heat measurements, the density of holes is 213 estimated33 to be m, z 8 x 1014cm‘2. For the boson model, if I approximate the density of bosons (n) as nearly the density of holes and the mass of bosons as roughly 1 ~ 5 electron mass, then k5 = (4"Ir'rz.)1/2 a: 131‘1 and kBTo z 0.1eV which are not very different from the fermi momentum (hp) and fermi energy (sip) in the paired fermion model. Thus physically relevant regime of A is 0.1 -- 0.5. From the above discussion and the results given in chapter 4, I now compare the neutron scattering from the paired fermion and paired boson superconductors. In general, the q—dependence of the inelastic and quasi elastic scatterings are quite different for the two superconductors. In Fig. 5.6 I have plotted the absolute value of the difference of the normal and superconducting state inelastic scattering cross sections as a function of hw/A for q/kp = 0.5 and 1.5 for the fermi and q/kb = 0.5 and 1.5 for the hose cases. The peak at 50: = 2A for all q in the fermi case is a consequence of the threshold for scattering of a pair of quasiparticles with total energy 2A. In contrast, in the hose case, there is a low energy threshold not directly related to A and a peak structure which are q-dependent moving to the right as q increases. The values of the cross section difference decrease as q increase in both cases. For paired fermion superconductors, the chemical potential p.(> 0) is always near the fermi energy 51:». Both the orbital and spin degree of freedoms con— tribute to the scattering. The threshold at low temperature inelastic scattering is always at twice the gap energy, 2A, as long as q < ZkF, from which the su- perconducting gap parameter can be extracted. The scattering decays to zero at high energy (above 5F). The orbital part dominates the quasi elastic scatter- ing at q < 56 "1. @(q) falls off rather rapidly over a characteristic wave vector ...-ii 214 of ~ 6“. For the charged spinless paired boson superconductor, the chemi- cal potential p(< 0) and the order parameter A (which is not the gap for the boson system) are crucially dependent on the temperature and the coupling con- stant. Only orbital current contributes to the scattering. The low temperature inelastic scattering threshold is q—dependent (ZEq/z). ‘I’(q) displays a much slower fall-off with q on the scale of kb and @(0) is proportional to A in the gapless regime. This feature is intrinsic to the boson system. A basic comparison” of the neutron scattering from the paired fermion superconductor and the paired charged boson superconductor is given in the Table 5.2. 215 (344 (1J2 ¢ (q) 0.0— ' I' 0 1 2 3 q/kb Fig. 5.5 Quasi elastic neutron scattering cross section in the superconducting state in units of F(eh/m)2 N2 kgTo. Inset gives the zero momentum transfer quasielastic scattering amplitude as a function of the coupling constant A. 216 1.0 I r I I 0.03 FEWHWI A = 0.05 e, 0.02 0.01 0.0 l l ‘ 1 0.00 0 2 4 6 8 10 hw/A Fig. 5.6 The absolute value of the difference of the normal and the super- conucting state inelastic scattering cross sections as a function of hm / A for the paired fermion superconductor. 217 Table 5.2 Comparison of the neutron scattering characteristics from the pair-ed fermion and paired boson superconductors. fermion boson Threshold at 2A for T = 0 at 2Eq/2 for T = 0 smeared for T 75 0 smeared for T ,4: 0 q-independent q-dependent Type of Response Orbital and Spin Orbital only Small q orbital dominates spin Peak at ~ §sp at ~ 2A q-dependent ‘ q-dependent Scattering -& 0 for ha: .z c; for hw z 10A Quasi elastic orbital response falls of rapidly with q falls off slowly with q 218 5.5 Summary In summary, I have made the first detailed theoretical investigation of the dynamic current-current and density-density response functions in the charged spinless paired boson superconductor. I have pointed out in this chapter that the elastic and quasi elastic neutron scattering from the charged spinless paired boson superconductor is quite different from a s-wave paired fermion superconductor, and therefore serve as diagnostic tools for the possibility of characterizing the nature of the quasiparticles involved in the mechanism of the high-Tc superconductors as well as determining their gap energy and coherence characteristics. List of References for Chapters 4 and 5 219 List of References for Chapters 4 and 5 1 Superconductivity, Vol. I and Vol. II, edited by R.D. Parks (Marcel Dekker, New York, 1969). 2 J. Bardeen, LN. Cooper, and J.R. Schrieffer, Phys. Rev. 108, 1175 (1957). 3 J .R. Schrieffer, Theory of Superconductivity, (Benjamin, Reading, Mass., 1964); I.M. Khalatnikov and A.A. Abrikosov, Adv. Phys. 8, 45 (1959); G. Rickayzen, in Ref. 1, Vol. I, chapter 2, page 51-116; J. Bardeen, in Carge’se Lectures in Theoretical Physics, edited by M. Lévy (Benjamin, New York, 1963), chapter 6; J. Bardeen, LN. Cooper, and J .R. Schrieffer, in Phys. Today, July, 1973, page 23-45. - 4 G.M. Eliashberg, Sov. Phys. JETP 11, 696 (1960); Sov. Phys. JETP, 12, 1000 (1961). 5 D.J. Scalapino, J .R. Schrieffer, and J .W. Wilkins, Phys. Rev. 148, 263 (1966); D.J. Scalapino, in Superconductivity, edited by R.D. Parks (Marcel Dekker, New York, 1969), Vol. I, page 449-560. 6 V. Ambegaokar, in Many-Boday Physics, edited by C. DeWitt and R. Balian (Cordon and Breach, New York, 1968), page 300-350. 7 W.L. McMillan, Phys. Rev. 167, 331 (1968). 8 RB. Allen and RC. Dynes, Phys. Rev. B 12, 905 (1975). 9 RB. Allen and B. Mitrovié, in Solid State Physics, edited by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic Press, New York, 1987), Vol. 37, page 1-92; P.B. Allen, ih Dynamic Properties of Solids, edited by G.K. Horton and A.A. Maradudin (North Holland, Amsterdam, 1980), Vol. 3, p. 95. 10 D. Rainer, in Progress in Low Temperature Physics, edited by D.P. Brewer (Elsevier, North-Holland, 1986), Vol. 10, chapter 4, page 371; G. Bergmann and D. Rainer, Z. Phys. 263, 59 (1973); D. Rainer and G. Bergmann, J. Low. Temp. Phys. 14, 50 (1974); V. Ambegaokar, P.G. deGennes, and D. Rainer, Phys. Rev A9, 2676 (1974); J.M. Daams and J.P. Carbotte, J. Low Temp. Phys. 43, 263 (1981). 11 J.C. Bednorz and K.A. Miiller, Rev. Mod, Phys. 60 585 (1988); J.C. Bed- norz and K.A. Miiller, Z. Phys. B 64, 189 (1986). 12 13 14 15 16 17 18 19 20 21 22 220 M.K. Wu et. al., Phys. Rev. Lett. 58 908 (1987); C.W. Chu et. al., Phys. Rev. Lett. 58, 405 (1987); C.W. Chu, Physica Scripta T27, 11 (1989); R.J. Cava et. al., Phys. Rev. Lett. 58, 1676 (1987); see also Phys. Today 40(4), 17 (1987). See articles in: Novel Superconductivity, edited by S.A. Wolf and V.Z. Kresin (Plenum, New Nork, 1987). A.K. Rajagopal, S.D. Mahanti, and W. Jin, Solid State Commun. 69, 847 (1989). W. Jin, S.D. Mahanti, and A.K. Rajagopal, Solid State Commun. (in press, 1989). H. Maeda et. al., Jap. J. Appl. Phys. 27, L209 (1988); see also Phys. Today, April, 1988, page 21-25. R.J. Cava et. al., Nature 332, 814 (1988); D.G. Hinks et. al., Nature 333, 836 (1988); C.K. Loong et. al., Phys. Rev. Lett. 62, 2628 (1989). Y. Tokura et. al., Nature 337 , 345 (1989); J .M. Tranquada et. al., Natrue 337, 720 (1989); V.J. Emery, Nature 337, 306 (1989). J.D. Jorgensen et. al., Phys. Rev. Lett. 58, 1024 (1987); For a review of the structure of the high-Tc oxide superConductors, see: R. Beyers and T.M. Shaw, in Solid State Physics, edited by H. Ehrenreich and D. Turnbull (Academic Press, New York, 1989), Vol. 42, page 135-212. L.M. Mattheiss, Phys. Rev. Lett. 58, 1028 (1987); L.F. Mattheiss and D.R. Haman, Solid State Commun. 63, 395 (1987); J. Yu and A.J. Freeman, and J.H. Xu, Phys. Rev. Lett. 58, 1035 (1987); A.J. Freeman, J. Yu, and CL. Fu, Phys. Rev. B 36, 7111 (1987); WE. Pickett, H. Krakauer, D.A. Papaconstantopoulos, and LL. Bogen, Phys. Rev. B 35, 7152 (1987); W.B. Pickett, Rev. Mod. Phys. 61, 433 (1989). D. Vaknin, S.K. Sinha et. al., Phys. Rev. Lett. 58, 2802 (1987); G. Shirane et. al., Phys. Rev. Lett. 59, 2802 (1987); K. Yamada et. al., Phys. Rev. B (to be published, 1989). S. Chakravarty, B.I. Halperin, and D. Nelson, Phys. Rev. Lett. 60, 1057 (1988); Phys. Rev. B 39, 2344 (1989); S. Ty‘c‘, B.I. Halperin, and S. Chakravarty, Phys. Rev. Lett. 62, 835 (1989); E. Manousakis and R. Salvador, Phys. Rev. Lett. 60, 840 (1988); ibid. 62, 1310 (1989); Phys. Rev. B 39,575 (1989). 23 24 25 26 27 28 29 30 31 32 33 221 A. Auerbach and D.P. Arovas, Phys. Rev. Lett. 61, 617 (1988); D.P. Arovas and A. Auerbach, Phys. Rev. B 38, 316 (1988); Y. Okabe, M. Kikuchi, and A.D.S. Nagi, Phys. Rev. Lett. 61, 2971 (1988); C.L. Kane, P.A. Lee, and N. Read, Phys. Rev. B 39, 6880 (1989). KB. Lyons et. al., Phys. Rev. Lett. 60, 732 (1988); R.J. Thomas et. al., Phys. Rev. Lett. 61, 1313 (1988). R.J.Birgeneau et. al., Phys. Rev. B 38, 6614 (1988).; R.J.Birgeneau et. al., Phys. Rev. B 39, 2868 (1989); R.J.Birgeneau et. al., Physica C 153-155, 515 (1988); M.A. Kastner it. al., Phys. Rev. B 38, 6636 (1988); Y. Endoh et. al., Phys. Rev. B 37, 7443 (1988); G. Aeppliet. al., Physica C 153-155, 1111 (1988); G. Aeppli et. al., Phys. Rev. Lett. 62, 2052 (1989); G. Shirane et. al., Phys. Rev. Lett. 63, 330 (1989). W.W. Warren et. al, Phys. Rev. Lett. 59, 1860 (1987); RE. Walstedt et. al., Phys. Rev. B 36, 5727 (1987); RE. Walstedt et. al., Phys. Rev. B 38, 9299 (1988); M. Mali et. al. ,.Phys lett. A124, 112 (1987); Y. Kitaoka et. al., J. Phys. Soc. Jpn. 57,31 (1988); K. Ishida et. al. J. Phys. Soc. Jpn. 57, 2897(1988). J.M. Tranquada et. al., Phys. Rev. B 35, 1613 (1987); J.M. Tranquada et. al., Phys. Rev. B 36, 5263 (1987). N. Niicker et. al., Z. Phys. B 67, 9 (1987);. N. Niicker et. al., Phys. Rev. B 37,5158 (1988). CE. Cough, et. al., Nature 326, 855 (1987); D. Esteve et. al., Europhys. Lett. 3, 1237 (1987); P.L. Gammel, D.J. Bishop, et. al., Phys. Rev. Lett. 59,2592 (1987). J.S. Tsai, et. al., Phys. Rev. Lett. 58, 979 (1987); T. Yamashita, et. al., Jpn. J. Appl. Phys. 26, L635 (1987). A.F. Andreev, Sov. Phys. JETP 19, 1228 (1964); V. Bentum et. al., Physica C 153-155, 1379 (1988). J. Niemeyer, M. R. Dietrich, and C. Politis, Z. Phys. B 67, 155 (1987); 0.5. Akhtyamov, Sov. Phys. JETP Lett. 3, 183 (1966). For a review on physical properties of the high-Tc oxide superconductors, see M. Tinkham and C.J. Lobb, in Solid State Physics, edited by H. Ehrenreich and D. Turnbull (Academic Press, New York, 1989), Vol. 42, page 91-134; see also H.E. Fischer, S.K. Watson, and D.C. Cahill, Comments on Condensed Matter Physics 14, 64 (1988). 34 35 36 37 38 39 222 A.L. Fetter and RC. Hohenberg, in Ref. 1, Vol. II, page 817-923; A.A. Abrikosov, Sov. Phys. JETP, 5, 1174 (1957). P.L. Gammel et. al., Phys. Rev. Lett. 61, 1666 (1988); P.L. Gammel et. al., Phys. Rev. Lett. 60, 144 (1987); P.L. Gammel et. al., Phys. Rev. Lett. 59, 2592 (1987); RS. Markiewicz, J. Phys. C 21, L1173 (1988); see also Phys. Today, March, 1989, page 17-21. D.S. Fisher, Phys. Rev. B 22, 1190 (1980); B.I. Halperin and D.R. Nelson, J. Low. Temp. Phys. 36, 599 (1979); S. Doniach and B.A. Huberman, Phys. Rev. Lett. 42, 1169 (1979); V. Ambegaokar, B.I. Halperin, D.R. Nelson, and E.D. Siggia, Phys. Rev. B 21, 1806 (1980); D.R. Nelson, in Fundamental Problems in Statistical Mechanics, edited by E.G.D. Cohen, (North-Holland, New York, 1980), Vol. V, page 53-108; D.R. Nelson, in Phase Transition and Critical Phenomena, edited by C. Domb and J .L. Lebowitz (Academic, New York, 1983), Vol. 7, page 1-99; S.N. Coppersmith, D.S. Fisher, B.I. Halperin, P.A. Lee, and W.F. Brinkman, Phys. Rev. Lett. 46, 549 (1981); 46, 869(E) (1981); Phys. Rev. B 25, 349 (1982). . D.R. Nelson, Phys. Rev. Lett. 60, 1973 (1988); D.R. Nelson and HS. Se- ung, Phys. Rev. B 39, 9153 (1989); T. N attermann and R. Lipowsky, Phys. Rev. Lett. 61,2508 (1988). For reviews on theoretical models of the high-Tc oxide superconductors, see: T.M. Rice, Z. Phys. B 67, 141 (1987); Nature, 332, 780 (1988); ibid. 337, 686 (1989); J. Mag. and Mag. Mater. 76-77, 542 (1988); N.F. Mott, Nature 327, 185 (1987); L. Garwin and P. Campbell, Nature, 330, 611 (1987); P. Fulde, Physica C 153-155, 1769 (1988); Physica Scripta, T23, 101 (1988); V.J. Emery, MRS Bulletin 14, 67 (1989); V.J. Emery, in Phys. Today, January 1989; V.L. Ginzburg, Phys. Today, March 1989, page 9; C.M. Varma, in High Tc Superconductors, edited by H.W. Weber (Plenum, New York, 1988), page 13-25. W. Weber, Phys. Rev. Lett. 58, 1371 (1987); W. Weber and L.F. Mattheiss, Physica B 148, 271 (1987); Phys. Rev. B 37, 599 (1988); E. We- ber, in High Tc Superconductors, edited by H.W. Weber (Plenum, New York, 1988), page 153-161; W. Weber, in Electronic Structure of Complex Systems, edited by P. Phariseau and W. Temmerman (Plenum, New York, 1984), page 345; RE. Cohen, W.P. Pickett, L.L. Boyer, and H. Krakauer, Phys. Rev. Lett. 60, 817 (1988); RE. Cohen, W.P. Pickett, and H. Krakauer, Phys. Rev. Lett. 62, 831 (1989). 40 41 42 43 44 45 46 47 48 223 J.E. Hirsch, Phys. Rev. Lett. 59, 228 (1987); ibid. 54, 1317 (1985); D.J. Scalapino, F. Loh, and J.E. Hirsch, Phys. Rev. B34, 8190 (1986); B 35, 6697 (1987); J.E. Hirsch, Phys. Lett. A 136, 163 (1989). RA. Lee and N. Read, Phys. Rev. lett. 58, 2691 (1987); RA. Lee, G. Kotliar, and N. Read, Physica B 148, 274 (1987); G. Kotliar, P.A. Lee, and N. Read, Physica C 153-155, 538 (1988); G. Kotliar, in Field Theo- ries in Condensed matter Physics, edited by Z. Tesanovic (Addison-Wesley, New York, 1989); C. Castellani and G. Kotliar, Phys. Rev. B 39, 2876 (1989); D.M. Newns, M. Rasolt, and RC. Pattnaik, Phys. Rev. B 38, 6513 (1988); D.M. Newns, Physica Scripta T23, 113 (1988). G. Kotliar and S. Liu, Phys. Rev. B 38, 5142 (1988); G. Kotliar and S. Liu, Phys. Rev. Lett. 61, 1784 (1988). K.S. Bedell and D. Pines, Phys. Rev. B 37, 3730 (1988). V.J. Emery, Phys. Rev. Lett. 59, 2794 (1987); V.J. Emery and G. Reiter, Phys. Rev. B 38, 4547 (1988); J. Zaanen and A.M. Oles’, Phys. Rev. B 37, 9428 (1988); P. Prelovsek, Phys. Lett. A 126, 287 (1988); F.C. Zhang and T.M. Rice, Phys. Rev. B 37, 3795 (1988). J .R. Schrieffer et. al., Phys. Rev. Lett. 60, 944 (1988); J.R. Schrief- fer et. al., Physica C 153-155, 21 (1988); J.R. Schrieffer et. al., Phys. Rev. Lett. 61, 2814 (1988); J.R. Schrieffer et. al., Phys. Rev. B 39, 11663 (1989); G. Vignale and K.S. Singwi, Phys. Rev. B39, 2956 (1989); J .R. Schrieffer et. al., Physica Scripta T27, 99 (1989). Y. Guo, J .M. Langlois, and W.A. Goddard III, Science 239,- 896 (1988); G. Chen and W.A. Goddard III, Science 239, 899 (1988); M.L. Cohen and L.M. Falicov, Science 243, 547 (1989). P. Prelovsek, T.M. Rice and F.C. Zhang, J. Phys. C 20, L229 (1987); P. Prelovsek et. al., Physica, 3148, 268 (1987); T.M. Rice, Nature 332, 780 (1988); K. Nasu, Phys. Rev. B 35, 1748 (1987); R.J. Birgeneau, M.A. Kastner, and A. Aharony, Z. Phys. B 71, 57 (1988); D. Emin, Phys. Rev. Lett. 62, 1544 (1989); D. Emin and M.S. Hillery, Phys. Rev. B 39, 6575(1989) C.M. Varma, S.Schmitt-Rink, and E. Abrahams, Solid State Commun. 62, 681 (1987); Physica B 148, 257 (1987); see also Ref. 13, page 355; E. Abra- hams, Physica C 153-155, 1622 (1988); A.J. Freeman et. al., J. Appl. Phys. 63, 4220 (1988); J. Ashkenazi and C.G. Kuper, in Studies of High Temper- ature Superconductors, edited by A.V. Narlikar (NAVO Science Publishers 49 50 51 52 53 54 55 56 57 58 224 Inc., New York, 1989), Vol. 3; UL. Cox et. al., Phys. Rev. Lett. 62, 2188 (1989); A.M. Oles and J. Zaanen, Phys. Rev. B 39, 9175 (1989). V.Z. Kresin, Phys. Rev. B 35, 8716 (1987); A. Griffin, Phys. Rev. B 37, 5943 (1987); ibid. 38, 8900 (1988); J. Ruvalds, Phys. Rev. B 35, 8869 (1987); A. Griffin and A.J. Pindor, Phys. Rev. B 39, 11503 (1989). A.K. Rajagopal, R.V. Kasowski, and S.D. Mahanti, Solid State Commun. 67, 883(1988) P.W. Anderson, Science 235, 1196 (1987); P.W. Anderson, in Fron- ties and Borderlines in Many Particle Physics, edited by J .R. Schrieffer and R.A. Broglia (North-Holland, Amsterdam, 1988); P.W. Anderson et. al., Physica C 153-155, 527 (1988); P.W. Anderson, Physica Scripta T27 , 60 (1989). M.J. Rice and Y. Wang, Phys. Rev. B 37, 5893 (1988). J.M. Wheatly, T.C. Hsu, and P.W. Anderson, Phys. Rev. B 37, 5897 (1988); Nature 333, 121 (1989). D. Pines and P. Noziéres, Theory of Quantum Liquids (Benjamin, New York, 1966); P. Noziéres, in Cargése Lectures in Theoretical Physics, edited by M. Lévy (Benjamin, New York, 1963), chapter 4; P. Noziéres, Theory of Inter- acting Fermi Systems (Benjamin, New York, 1964); G. Baym and C. Pethick, in The Physics of Liquid and Solid Helium, edited by K.H. Bennemann and J .B. Ketterson (John Wuley, New York, 1978), Vol. II, chapter 1, page 1- 122; J .W. Negele and H. Orland, Quantum Many-Prarticle Systems, (Addison-Wesley, Reading, Mass., 1988), Chapter 6, page 296-331. 5.8. Jha, J. Phys. C 29, L615, (1987) ; S.S. Jha, in Studies of High Tem- perature Superconductors, edited by A.V. Narlikar (Nova, New York, 1989), Vol. 1; 3.5. Jha and A.K. Ragajopal (to be published). A.K. Rajagopal, W. Jin, S.D. Mahanti, and 8.5. Jha, in Highlights in Can- densed Matter Theory, edited by P. Jena, R. Kalia, M. Tosi, and P. Vashishta (to be published by World Scientific, 1990). R.A. Klemm, A. Luther, and M.R. Beasley, Phys. Rev B 12, 877 (1975); S. Takahashi and M. Tachiki, Phys. Rev. B 33, 4620 (1986); Y. Suwa, Y. Tanaka, and M. Tsukada, Phys. Rev. B 39, 9113 (1989); J. Ihm and RD. Yu, Phys. Rev. B 39, 4760 (1989). G. Wendin, Physica Scripta T27 , 31 (1989) and references therein. 59 60 61 62 63 64 65 66 67 225 R.L. Kurtz et. al., Phys. Rev. B 35, 8818 (1987). J. Hubbard, Proc. Roy. Soc. London, Ser. A 276, 138 (1963); 277,237 (1964); 281,401 (1968). M.C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963); Phys. Rev. 134, A923 (1964); 137, A1726 (1965); T.M. Kaplan, P. Horsch, and P. Fulde, Phys. Rev. Lett. 49, 889 (1982); P. Horsch and T.A. Kaplan, J. Phys. C16, L1203 (1983). C. Gros, R. Joynt, and T.M. Rice, Phys. Rev. B 36, 381 (1987); C. Gros, R. Joynt, and T.M. Rice, Z. Phys. B68, 425 (1987) ; W.F. Brinkman and T.M. Rice, Phys. Rev. B 2, 4302 (1970); A.E. Harris and R.V. Lange, Phys. Rev. 157, 295 (1967); K.A. Chao, J. Spalek, and A.M. Oles’, J. Phys. C10, L271 (1977); Phys. Rev. B18, 3453 (1978); K. Huang and E. Manousakis, Phys. Rev. B 36, 8392 (1987); K. Huang, in Studies of High Temperature Superconductors, edited by A.V. N arlikar ( Nova Science Pub. Inc., Commack, New York, 1989), Vol 2; A.H. MacDonald, S.M. Girvin, and D. Yoshioka, Phys. Rev. B 37,9753 (1988). J. Hybertsen, M. Schliiter, and N .E. Christensen, Phys. Rev. B 39, 9028 (1989); M. Schliiter, J. Hygertsen, and N.E. Christensen, Physica C 153- 155, 1217 (1988); A.K. MacMahan, R.M. Martin and S. Satpathy, Phys. Rev. B 38, 6650 (1988); N.M. News at. al., Pys. Rev. B 38, 7033 (1988). W.A. Little, Phys. Rev. 134, A1416 (1964); V.L. Ginzburg, Sov. Phys. JETP 40, 397 (1964); Phys. Lett. 13, 101 (1964); Physica C 153-155, 1617 (1988); Physica Scripta, T27, 76 (1989); D. Allender, J. Bray, and J. Bardeen, Phys. Rev. B 7, 1020 (1973); H. Gutfreund and W.A. Little, in Highly Conducting One-Dimensional Solids, edited by J .T. Devreese et. al. (Plenum, New York, 1979). J.E. Hirsch, Phys. Rev. B 31, 4403 (1985); J. Stat. Phys. 43, 841 (1986). P.W. Anderson, Mater. Res. Bull. 8, 153 (1973); P. Fazekas and P.W. An- derson, Phil. Mag. 30, 423 (1974). S.A. Kivelson, D.S. Rokhsar, and J .P. Sethna, Phys. Rev. B 35, 8865 (1987); S. Kivelson, Phys. Rev. B 36, 7237 (1987); S. Kivelson and D.S. Rokhsar, Physica C 153-155, 531 (1988) ; Dzyaloshinski, A.M. Polyakov, and RB. Weigmann, Phys. Lett. A127, 112 (1988); RB. Weigmann, Phys. Rev. Lett. 60, 821 (1988); Physica C 153-155, 103 (1988); see also Phys. Today, February 1988, page 19-23. 68 69 7O 71 72 73 74 75 76 77 78 79 226 W.P. Su, J .R. Schrieffer, and A.J. Heeger, Phys. Rev. Lett. 42, 1698 (1979); Phys. Rev. B 22, 2099 (1980); J.R. Schrieffer, in Highlights of Condensed Matter Theory, edited by F. Bassani (Elsevier, New York, 1983), page 300-348; A.J. Heeger, S. Kivelson, J .R. Schrieffer, and W.P. Su, Rev. Mod. Phys. 60, 781 (1988). W.A. Little, Science 242, 1390 (1988); M.R. Beasley and J.W. Geballe, Phys. Today 37(10), 60 (1984). S.W. Lovesey, Theory of Neutron Scattering From Condensed Matter, Vol. 2 (Clarendon Press, Oxford, 1984); For a general reference on neutron scat- tering, see Physics Today, January, 1985. G. Aeppli et. al., Phys. Rev. Lett. 58, 808 (1987); ibid. 60, 615 (1988); G. Aeppli et. al., J. Mag. Mag. Mat‘er. 76-77, 385 (1988). S.K. Sinha, G.W. Crabtree, D.G. Hinks, and H.A. Morok, Phys. Rev. Lett. 48, 950 (1982); S.K. Sinha, H.A. Mook, O.A. Pringle, and D.G. Hinks, in Superconductivity in Magnetic and Erotic Materials, edited by T. Matsub- ara and A. Kotani (Springer-Verlag, Berlin, 1984), page 14-28; S.K. Sinha, in Condensed matter Research Using Neutrons, edited by S.W. Lovesey and R. Scherm (Plenum, New York, 1984), page 249-276; see also Physics Today, March 1982, page 17. BL. Cooper et. al., Phys. Rev. B 37, 5920 (1988) and references therein. W. Jin, S.D. Mahanti, and A.K. Rajagopal, Phys. Rev. B (to be pub- lished). A.K. Rajagopal, W. Jin, and S.D. Mahanti, Physica C, in press. For the magnetic neutron scattering from solid, see R.M. White, Quantum Theory of Magnetism, (Springer-Verlag, Berlin, 1983), chapter 8, page 248- 270; T. Izuyama, D.J. Kim, and R. Kubo, J. Phys. Soc. Japan, 18, 1025 (1963); S. Doniach, Proc. Roy. Soc. 91, 86 (1967); J.E. Hebborn and N.H. March, Adv. Phys. 19, 175 (1970). K. Sasaki and Y. Obata, Prog. Theor. Phys. (Supplement) 69, 406 (1980). R. Kubo, J. Phys. Soc. Japan. 12, 570 (1957). A.K. Rajagopal and K.P. Jain, Phys. Rev. A 5, 1475 (1972); Intern. J. Magnetism, 2, 183 (1972). 80 81 82 83 84 85 86 87 88 89 227 A.L. Fetter and J .D. Walecka, Quantum Theory of Many-Particle Systems (McGraw-Hill, New York, 1971). S.K. Misra, P.K. Misra, and S.D. Mahanti, Phys. Rev. B 26, 1903 (1982); H. Fukuyama, Prog. Theor. Phys. 45, 704 (1971); H. Fukuyama and J.W. McClure, Phys. Rev. B 9, 975 (1974); D.J. Kim, Phys. Report, 171, 129 (1988). L.P. Gor’kov, Sov. Phys. JETP 34, 505 (1958); ibid. 9, 1364 (1959); V. Ambegaokar, in Superconductivity, edited by R.D. Parks (Marcel Dekker, New York, 1969), Vol. 1, page 259-319. A.A. Abrikosov and L.P. Gor’kov, Sov. Phys. JETP 35, 1558 (1959) ; ibid. 35, 1090 (1959); A.A. Abrikosov and V.M. Genkin, Sov. Phys. JETP 38, 417 (1974). A.A. Abrikosov, L.P. Gor’kov, and LE. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics, (Prentice-Hall. Inc., Englewood Cliffs, N.J., 1963), Chapter 7. ‘ A.A. Abrikosov and L.A. Falkovsky, Sov. Phys. JETP 13, 179 (1961); A.A. Abrikosov and L.A. Falkovsky, Physica C 156, 1 (1988); Physica Scripta T27, 96 (1989); K.P. Jain and C.C. Chancey, Phys. Rev. B 39, 9049(1989). J .W. Negele and H. Orland, Quantum Many-Particle Systems (Addison-Weslwey, Reading, Mass, 1988). N.N. Bogoliubov, Sov. Phys. JETP, 7, 41 (1958); J .G. Valatin, Nuovo Cimento, 7, 843 (1958); S.T. Beliaev, in The Many Body Problem, edited by C. DeWitt (John Wiley, New York, 1959), page 343. G. Baym and ND. Mermin, J. Math. Phys. 2, 232 (1961); RC. Mar- tin and J. Schwinger, Phys. Rev. 115, 1342 (1959); L.P. Kadanoff and RC. Martin, Phys. Rev. 124, 670 (1961); L.P. Kadanoff and G. Baym, Quantum Statistical Mechanis, (Benjamin, New York, 1963); L.P. Kadanoff, in Lectures on the Many-Body Problem, edited by ER. Caianiello (Academic, New York, 1964), Vol. 2, page 77-112; L.P. Kadanoff, in Statistical Physics, Phase Transitions and Superfluidity edited by M. Chrétien et. al. (Gordon and Breach, New York, 1968), V0. 1, page 109-178. P.A. Lee, T.M. Rice, J.W. Serene, L.J. Sham, and J.W. Wilkins, Comments on Condensed Matter Phys. 123, 99 (1986); C.M. Varma, Comments Solid State Phys. 11, 221 (1985); G.R. Stewart, Rev. Mod. Phys. 56, 755 (1984); H.R. Ott, in Progress in Low Temperature Physics, edited by D.F. Brewer (North-Holland, New York, 1987), Vol. 11, p. 215; Z. Fisk, 90 91 92 93 228 D.W. Hess, C.J. Pethick, D. Pines, J.L. Smith, J.D. Thompson, and J.O. Willis, Science 239, 33 (1988). A.J. Leggett, Rev. Mod. Phys. 47, 331 (1975); P.W. Anderson and P. Morel, Phys. Rev. 123, 1911 (1961); R. Balian, L. Nosanow, and N.R. Werthamer, Phys. Rev. Lett. 8, 372 (1962); R. Balian and N.R. Werthamer, Phys. Rev. 131, 1553 (1963); R. Balian, in Lectures on the Many-Body Problem, edited by ER. Caianiello (Academic, New York, 1964), Vol. 2, page 147-158; P.W. Anderson and W.F. Brinkman, Phys. Rev. Lett. 22, 1108 (1973); W.F. Brinkman, J.W. Serene, and P.W. Anderson, Phys. Rev. B 10, 2386 (1974); P.W. Anderson and W.F. Brinkman, in The Physics of Liquid and Solid Helium, edited by K.H. Bennemann and J .B. Ketterson (Wiley, New York, 1978), Vol. 2. R. Joynt and T.M. Rice, Phys. Rev. B 38, 2345 (1988). J. Gavoret and P. Noziéres, Ann. Phys. (N.Y.) 28, 349 (1964); RC. Ho- henberg and RC. Martin, Ann. Phys. (N.Y.) 34, 291 (1965); W.A.B. Evans and Y. Imry, Nuovo Cimento 63B, 155 (1969); Y. Imry, in Quantum Liq- uids, edited by N. Wiser and D. Amit (Gordon and Beach, New York, 1970), p. 603; A. Griffin, Phys. Rev. B 19, 5946 (1979); A. Griffin and E. Talbot, Phys. Rev. B 24, 5075 (1981); P. Noziérex ans D.S. James, J. Physique 43, 1133 (1982). W. Jin, A.K. Rajagopal, and S.D. Mahanti, in Studies of High Temperature Superconductors, edited by A.V. Narlikar (to be published by NOVA Science Pub. Inc., Commack, New York, 1990), Vol 6. “10900“