.14]? . 4. .,4 .au'4. u «4' J. 4 4 ":l‘i - '.I\‘ \'~ '.4LIL.\’L n .mr. ”.42. $4444.14.“ L". \24 1.4.. u .4 4...... .41 J l 4 4. 4.1.4.1. .- ..L . 14.4“. "4.4.4. I “A. 4.1.4 .4..4 vim. .1} .4. » . .3. :.:i.L.A.4u4.¢u.:__f “mi-1. . . 4.4..lr44....4u 44-..“:- . 4 4‘. . 4. 41.4 “4.. mug-4 . .4. -. 4 ‘l: .414. ..4. A.'. ~. 4 ‘ {up 2:...1'; ' .L':.\l;l\l| .: > 4 ....44.; 1.. . 4-. quEE .1; III-M"! .4..44..4-.4~ 4 A u... 4.- 4;, .. '1,.7|;.. 4. 1.-.44::L&:44.:l ¢ 444 ...4....-."L".4'.4 . ‘ ' {'3 raw n1. 1 'r 4. ”4.4.: .4 4-4:.4.L.4~ :1. .4 4L 4 . 1)...“ L154, :1 gilu...n'J ... h ;.\ ..4....n.. . " .545.....‘...4 :..;.:..14. a... ‘II'JLxull H15“ . 'Q $1.1»... ‘11... a. . 4; _, .. VF‘Z-m - . . . "4- In a}; 'u. L. 4; “4}! u. .:1.- 3.4;,” - 4.44 .414.n4.#.'1. :.I1l.1L’-.~lu 4 144.4 4.4.114- .4) . 4 .. ..4 “n u" I: - 4 4 44.4. rut-1 44 4 .,, 1 , {$1.41. 4......14 14 :....'u:' -.. 4. . ‘k l 4.24:; .‘U'm‘ Il 4 .m.. 4.31 .4... (.1 LI' I u 4l.'uu- .. .... L. 44. ‘2... m - :4 't. 44.4w“, '1} -.-.,.. .4 1441.414 .. I I ' l 1 \‘.l» In" 4 N4“; 4 , . 1.x: .wdu..\.w.. if... :4 va- .14.“; 441mm ;—14 44m 1...4 :L ,4., “4. -L ..\" g- .7; (".4 . . n4 .. ,. 4 hum..-‘ . . 4.- .. .4...» r». .m v.4 mm ..4-...- 4. 3...:yi z. 441 “14.4. n4'L..\ .2“ .4. , (1. huh“. 14 ‘ if" H CC . 2312...: ,. .3...u;......, . ‘ ".}..4f:......4: -. . . .’. n m «NV-mu wait ‘4' M.-. Lt}... .. , «4 4.4.:«4'C‘Li.4.a :AI‘IIJU 14-.--4LL- «:1... ‘5“; .. .‘..r<-»1A:~l'h 4.4 3. .\:L- Myaw' -..4\. 4.1.54. .. ..L: 4- L..4Jf :43'..-L.. ML" III n-leu LW “14.4... 4...... rm... «1.: u 4 I. \.K-_ -~,-u-...im.h 1:44,. 4 4... 44: - 4 v_ .. ,., A ,. . .4 -4m. .44 .w -.4 x- 14.! 4.: 4 (W 'l 3:1,. 4.. .~ . .. .-.. vu. m . 4 .v u: 4'; «9-4 4» 4~r .44 :w-z v. n- 4"“ «1"!4......... 41:f_"1""":"‘ ..f l {'I I' .4..4.4.44~',.1.m nu. ,. 4.. 4.. .u _...... «WilliWillilililillllifllllfill 3 1293 00903 4970 THESIS This is to certify that the dissertation entitled Structural Phase Transition in Two Dimensional Layered Systems presented by Zhi—Xiong Cai has been accepted towards fulfillment of the requirements for Ph.D. degree in Physics th / Major professor Z Date 9 LI [770 Msu;.....1m...m.;...,1 ‘ r1 m" ~ , .- .- 042771 LIBRARY Mlchlgan State Unlverslty PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE MSU Is An Affirmative Action/Equal Opportunity Institution c:\clvc\datndw.un3—or STRUCTURAL PHASE TRANSITION IN TWO DIMENSIONAL LAYERED SYSTEMS By Zhi-Xiong C ai 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 1990 ABSTRACT STRUCTURAL PHASE TRANSITION IN TWO DIMENSIONAL LAYERED SYSTEMS By Zhi-Xiong Cai In this thesis we perform Monte Carlo simulations on four types of layered systems to study the structural phase transitions in three and diffusion in the fourth of these systems: (1).The orientational order-disorder transition of N 3 adsorbed on graphite mono- layers, for which we have established a. 2-dimensional anisotropic planar rotor model by calculating the quadrupole-quadrupole interaction of N 3. Using the recently devel- oped histogram method combined with analysis of finite size scaling and the fourth order energy cumulant, we have determined unambiguously that the herringbone orientational transition of N2 molecules adsorbed on graphite (001) surface is a con- tinuous transition contrary to the RG results. The critical exponents are determined to be those of q=3 Potts model in two dimension, in good agreement with the recent experimental data on a system of the same universality class. (2).The melting transition and superlattice ordering of ternary graphite interca- lation compounds 11.81403 is studied by performing Monte Carlo simulation of constant concentration (1:) on a spin-1 Ising (lattice gas) model appropriate for this system. The concentration (2:) dependence of the lattice-gas melting temperature Tm is studied for different strengths of interaction parameters. We find that the melting temperature Tm of the ternaries is always greater than that predicted by Vegard’s law. A simple physical argument for this behavior within a Landau expansion is also given. (3). Oxygen ordering in the YBagCuaO-M; basal plane is described by a 2-dimensional lattice gas model with anisotropic long range interaction. Using Monte Carlo sim- ulations, we studied the orthorhombic to tetragonal phase transition, the effect of thermal quenching and the effect of substituting Cu by trivalent impurities on the oxygen ordering. (4). Finally, we have constructed a lattice model to probe the 2-dimensional mi- croporosity of intercalated layered systems of the type A1_,B.L where A and B are two types of intercalants sandwiched between host layers-L. Percolation of a probe particle in this 2-d microporous medium has been studied analytically and by Monte Carlo simulation. Depending on the transverse layer rigidity which is characterized by a healing length A, we find two “normal” percolation thresholds which coalesce for a particular value of /\ = A1. Percolation at this point belongs to a different universality class. Anomalous diffusion has also been investigated for these models. ACKNOWLEDGEMENTS I am deeply indebted to my advisor, Professor 3. D. Mahanti, for his constant guidance and encouragement in this work. During the last three years I have been benefitted immensely from his broad knowledge on various areas of physics and chem- istry and his approach to research will always been an inspiration for me. I would like to thank Dr. S. A. Solin and Professor T. J. Pinnavaia of Chemistry Department for their advice and collaboration on the problem of intercalated layered systems. I would like to thank Professor T. A. Kaplan for his advice on the problem of Y-Ba—Cu-O systems. I would also like to thank Dr. Alan Ferrenburg of University of Georgia for his advice on the histogram method. I would express my appreciation to Dr. Xi-Yiao Chen of East China Normal University, my B.Sc. thesis advisor, for teaching me the Monte Carlo simulation method. I am grateful to Professor M. F. Thorpe, Professor J. Cowen, Professor J. Kuhn and Professor D. Stump for their service on my Ph. D guidance committee, and I would like to express my special thanks to Professor J. Kovacs for his patient assistance and Mrs. Janet King for her help. i I gratefully acknowledge financial support from Michigan State University, the Na; tional Science Foundation and the Center for Fundamental Material Research. I dedicate this work to the memory of my father. iii Contents 4 Introduction Application of Monte Carlo Method to Study Phase Transitions 2.1 Introduction ................................ 2.2 Order Parameters, Critical Exponents, Scaling, and Universality . . . 2.3 Monte Carlo Technique ...................... . . . 2.4 Finite-Size Effects ............................. 2.5 Optimized Treatment of Data: Histogram Method ........... N3 Adsorbed on Graphite: Anisotropic Planar Rotor Model 3.1 Introduction ................................ 3.2 The Anisotropic Planar Rotor Model .................. 3.3 Previous Work Using Mean Field Theory, RG, MC, MD simulation 3.4 Present Monte Carlo Studies ....................... 3.5 Summary and Conclusion ........................ tempura 11 13 17 17 19 22 24 32 Lattice-Gas Melting of Two-Dimensional Alloys: Application to Ternary Graphite Intercalation Compounds 4.1 Introduction ................................ 4.2 Review of Previous Work on Ternary GIC’s .............. 4.3 Lattice Gas Model And Its Parameters ................. 4.4 Mean Field Approach ........................... iv 38 38 42 47 4.5 Monte Carlo Simulation ......................... 51 4.5.1 Melting Transition ........................ 51 4.5.2 Order-Disorder Transition on the 2 x 2 Superlattice ...... 58 4.6 Summary and Conclusion ........................ 60 5 Oxygen Ordering in YBEQCU3O7_5 and Related Systems 65 5.1 Introduction ................................ 65 5.2 Oxygen Ordering and Superconductivity ................ 66 5.3 Orthorhambic to Tetragonal Phase Transition ............. 67 5.4 The Lattice-gas Model and its Parameters ............... 72 5.4.1 Earlier Theoretical Work ..................... 74 5.4.2 Monte Carlo Simulation ..................... 76 5.4.3 J ustification of the Parameter Values .............. I 78 5.5 Effect of Thermal Quenching on Oxygen Ordering ........... 81 5.6 Effect of Substituting Cu(1) by a Trivalent Impurity on Oxygen Ordering 84 5.7 Summary and Conclusion ........................ 86 6 Dual Percolation Threshold in 2-d Microporous Media 110 6. 1 Introduction ................................ 110 6.2 The Model System ..................... . ....... 1 11 6.3 Solution for the Case of A = A1 ..................... 118 6.4 Solution for the Case of A > A1 ..................... 126 6.5 Summary ................................. 127 7 Summary and Discussion 133 List of Tables 3.1 Comparison of theoretical and experimental critical exponents for the model which belongs to the universality class of Heisenberg model with face type cubic anisotropy compared with the critical exponents of q=3 Potts model in two dimension. ..................... vi List of Figures 3.1 3.2 3.3 3.4 3.5 3.6 3.7 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 5.1 5.2 5.3 The herringbone structure ........................ 21 Temperature dependence of VL ..................... 26 temperature dependence of specific heat ................ 28 Lattice size dependence of specific heat peak .............. 29 Temperature dependence of susceptibility ................ 30 Lattice size dependence of susceptibility peak ............. 31 Lattice size dependence of Tc ....................... 33 Graphite and its intercalation compounds. ............... 40 A graphite layer and intercalants .................... 44 Various superlattice structures ...................... 48 Two sublattice model. .......................... 49 Mean field results for Al/aBz/3V3 .................... 52 MC simulation results for 141/38“an .................. 55 Melting temperature versus concentration :1: .............. 56 Configuration of A1/332/3V3 ....................... 57 Correlation function 0,, .......................... 59 Ground State configuration of A1/3Bg/3V3 ................ 61 The ground state of YBEzCU3O7. .................... 68 Effect of oxygen content on T6 ....................... 69 Fractional site occupancy of oxygen: experimental data ........ 72 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 5.19 5.20 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 The ground state structure of YBagCugO-r basal plane ......... 73 The order parameter and susceptibility ................. 77 Oxygen concentration and fractional site occupancy: theory ..... 79 Oxygen concentration isotherm ..................... 89 Temperature dependence of Cu-O chain length ............. 90 The Phase Diagram of YBagCu307_5 .................. 91 Oxygen configuration of the annealed YBagCu3O” .......... 92 Oxygen configuration of the quenched YBagCugOM . . . . . . . . . 93 Oxygen configuration of the annealed YBagCu3Ou .......... 94 Oxygen configuration of the quenched YBagCu303_7 ......... 95 Oxygen structure factor of the annealed YBagCuaOM ........ 96 Oxygen structure factor of the quenched YBagCuaou ........ 97 Variation of T6 of the YBagCu3_.,Ga,O-, .......... . ...... 98 Structure of YBagCu3_,Ga,O7 basal plane with single impurity M. . 99 Order parameter of YBagCu3_.,Ga,O-,. ................. 100 Oxygen configuration for a: = 0.02 .................... 101 Oxygen configuration for a: = 0.12 .......... . .......... 102 The unit cell of pillared clay ....................... 112 Top view of the gallery layer ....................... 113 Lateral constraints ............................ 115 Vertical constraints ............................ 116 Possible configuration for A = A1 ..................... 119 N,~sforA=A1 ............................. 123 Sm“ ~ L for A = A1 ........................... 124 ~tforA=A1 ........................... 125 W(z) ~ :2: for different healing length .................. 128 The phase diagram ............................ 129 viii Chapter 1 Introduction One of the fast growing areas of condensed matter physics is to understand and predict structural and transport properties of layered systems starting from a microscopic basis. There are large number of solids such as graphite, chalcogenides, pillared complex layered oxides, and the new high-Tc oxide superconductors all of which have layered structure. A common feature of these systems is that the interaction 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 systems even though the interlayer interaction cannot be completely ignored. The primary objective of this thesis is to study the structure and dynamic properties of different types of layered systems using computer simulation of lattice models. This dissertation may be divided into three parts. The first part consists of chapter 2 where the general theory of phase transition and critical phenomena are discussed. Of special importance is the application of Monte Carlo simulation and finite size scal- ing theory to understand the critical phenomena. We will also introduce the newly developed energy histogram method and discuss its implementation. This method is very useful in calculating the thermodynamic quantities of model systems within rea- sonable computer time. In the second part, we study the structural phase transition of various layered systems using the method outlined in the chapter 2. In chapter 3 we will study a simple 2-dimensional lattice system: the orientational phase tran- sition of N2 molecules adsorbed as a monolayer on a graphite substrate. Extensive experimental and theoretical work have been done on this system. It is well estab- lished that an anisotropic planar rotor model describes adequately the orientational phase transition observed in this system. We have studied this model using Monte Carlo simulation. In particular, we give a clear indication that this model has a con- tinuous transition, contrary to the prediction of renormalization group calculation, but in good agreement with the experimental results. In chapter 4 we will study a more complicated layered system: the ternary graphite intercalation compound. We will show, via mean field theory and computer simulation how the lattice-gas melting temperature changes as a function of the relative concentration of the two intercalants. Superlattice ordering for certain alloy concentrations at low tempera- ture will also be discussed using the same model. In chapter 5 we will study the oxygen ordering and the orthorhombic to tetragonal phase transition in the newly discovered high temperature superconducting oxide YBa20u307 and related oxygen deficient compounds. A simple lattice gas model with short range anisotropic inter- action is established based on the phenomenology of screened Columb interaction. Using this model, we are able to simulate the oxygen ordering in'the basal plane of YBagCu307_5 under various preparation conditions and our results are in excellent agreement with experiments. A slightly modified model based on essentially simi- lar arguments is proposed to study the effect of substituting Cu atom in the basal plane by trivalent nonmagnetic impurities. The physical meaning of this model will be discussed. The relationship between structure and electronic and superconducting properties of these systems will be briefly reviewed. The third part of this thesis con- sists of chapter 6 in which we construct a lattice model to probe the 2-dimensional microporosity of mixed ion pillared clay systems, a new class of ternary intercala- tion compounds. Percolation and diffusion of rare gas atom in this system is studied both analytically and by Monte Carlo simulation. Depending on the transverse layer rigidity of the intercalating medium, characterized by a healing length A, one finds two “normal” percolation thresholds.These two thresholds coalesce for a particular value of A, where the percolation and diffusion belongs to the universality class of hull percolation which is different from the usual 2—d site percolation. Chapter 2 Application of Monte Carlo Method to Study Phase Transitions 2.1 Introduction It is well established that phase transition, which is characterized by abrupt changes, discontinuities, and strong fluctuations in thermodynamic quantities, is a consequence of a cooperative phenomenon and thus is intimately related to the interaction between the microscopic constituents of matter. The number of constituents in a macroscopic system is typically of the order of 1023 and each constituent may carry several degrees of freedom. For such a large system statistical mechanics provides the theoretical basis for the relationship between the microscopic and macroscopic descriptions of a physical system. In classical statistical mechanics, a configuration of a system is described by a set of mechanical variables, I‘. I‘ contains the value of all possible degrees of freedom for each particle of the system, e.g. spatial position, velocity, and magnetic moment. The phase space, {I‘}, is the parameter space spanned by all possible configurations of the system. Properties of the system are determined by a Hamiltonian function, H (1"), defined on the space of mechanical variables. In statistical mechanics, a probability is associated with each configuration e-H(I‘)/kBT p(F) = _—Z—_’ I (2'1) where Z is a normalization factor (the partition function) Z = e‘Hirl/kaTdI‘. (2.2) {1‘ } T is the absolute temperature and k3 is Boltzmann’s constant. Given the probability distribution of the configurations, the thermodynamic value of a measurable physical quantity, A(f), is obtained in the canonical ensemble as < A(T) >.—. z-1 [{r}A(r)p(r)dr. (2.3) Equation 2.3 constitutes the formal connection between the microscopic and macro- scopic physical worlds. It is obvious that this equation involve multi-dimensional integrals. Only in a very few special cases can these integrals be evaluated exactly. Approximation schemes have to be introduced for most of the physical systems of interests. It is at this point that the modern computer enters as an indispensable tool to make progress. By using a fast computer, the integrals of statistical mechanics may be evaluated with an accuracy which is only limited by the computer power available. Computer is used to simulate directly the behavior of a physical system taking as its starting point the fundamental equations of statistical mechanics. It is like an exper- iment: it is built on as little bias as possible. Only the fundamental physical laws governing the interaction between the microscopic constituents are invoked. Contrary to a real experiment, the simulation is carried out on a well-defined system and there is full control over every “experimental parameter”. However, simulation shares with the real experiment an important potential, namely that it allows for possible new discoveries which could not have been inferred from the basic physical laws of inter- action. The outcome of a simulation may be thought of both as new experimental data and as a theoretical result which can be used to assess the validity of basic assumptions and predictions of analytical theories. Therefore, computer simulation interpolates between theory and experiment. By this unique ability, computer simu- lation may serve to illustrate and illuminate subtle and, unfortunately, not commonly recognized basic conceptual relations in scientific reasoning. There exists two major types of computer simulation techniques which have proved exceedingly successful in the study of phase transitions, namely molecular dynamics method and Monte Carlo method. In a molecular dynamics simulation, the deterministic time-evolution of a classic many body system is calculated by a numerical integration of Newtonian equations of motion. In a Monte Carlo simulation, certain stochastic elements are introduced which facilitate the evaluation of statistical mechanical averages. The underlying principles of applying Monte Carlo simulation techniques to study phase transitions will be described in Sec. 2.3 of this chapter. Before we proceed with a detailed description of the Monte Carlo method (Secs.2.3 and 2.5), we will discuss briefly the theoretical treatments of phase transitions and critical phenomena in general (Sec.2.2). These remarks serve to lodge computer studies within the overall framework of theoretical approaches. 2.2 Order Parameters, Critical Exponents, Scaling, and Uni- versality Yang and Lee were the first to associate phase transitions with the singularities of partition function Z [1] Using the definition of the free energy F = —k3Tln Z, (2.4) we shall classify phase transitions as continuous or first-order according to whether the first derivatives of F are continuous or not at the singular point, i.e. at the phase transition temperature[2]. A phase transition is characterized by a spontaneously broken symmetry [3]. The phase with the broken symmetry has a symmetry which is lower than that of the Hamiltonian. Symmetry-breaking is conveniently described in terms of an order pa- rameter, Q. For a given problem, there is no unique way of choosing Q. Its magnitude measures the degree of long-range order present in the system. It has the symmetry of the ordered phase and it may have several components. The order parameter can be associated with a thermodynamic conjugate field, ha, which couples directly to Q. The corresponding term in the free energy is —hq.Q. Consequently, the order parameter is defined by Q = —(8F/8h§ )T. According to the Fisher’s classification of phase transitions[2], we then see that Q is discontinuous at first-order transition and goes to zero continuously at a continuous transition (critical point). In experiments as well as in computer simulations, the behaviour of the order parameter at the phase transition is usually an important criteria to determine the nature of the transition. It should be noted, however, that there can be phase transitions without long range order. The most famous example is the Kosterlitz-Thouless transition of XY model in two dimension[4, 5]. The susceptibility associated with the order parameter is an important quantity which is expected to be influenced by critical fluctuations particularly for a continuous transition. The isothermal ordering susceptibility, X9, is defined by X9 = (BQ/Bhsh- (2.5) The fluctuation-dissipation theorem[6] relates the fluctuations in the order parameter to the corresponding susceptibility. Introducing the order parameter “operator” M, we therefore have X4, = (kBT)—1(< f”: > -T2), (2.6) where Q =< AI >. This is an important result for computer simulations because it implies that response functions may be evaluated directly from the thermal fluctua- tions and a perturbing ordering field need not be introduced. Another useful version of the fluctuation-dissipation theorem relates the fluctuations in internal energy to the specific heat 313 '2 -1 a 2 Ch 2 (Eh; 22' (kBT) (< H > — < H > ). (2.7) Widom was the first to advance the scaling hypothesis for the static critical phenom- ena [7]. According to this hypothesis, the singular part, F, of the free energy is a generalized homogeneous function, i.e., 171m, Abh.) = A°+"17’(t, 11.), ' (2.8) where t E (T — T,_.)/Tc measures the relative distance from the critical temperature, Tc. Taking the derivative of F with respect to h; to obtain the corresponding singular behaviour of the order parameter, we find emhzoyvpu” (am with the critical exponent ,6 = (1 - b)/a.. Similarly, power-law singularities may be derived for other thermodynamic quantities, e.g. Ch(t) ~ t“°‘,T > Tc ~(4rflT Tc ~(4rflT" news—anal") where E” represent the energy eigenvalues. Since only a few states a which are close (2.14) to the most probable values contribute significantly to the average, it is sufficient to generate only the important states by assigning enhanced probabilities to them. This 10 is known as the “important sampling” and the averages are then defined as : 2 22:1 AuPp‘lezp(—Eu/kBT) < A > 2:1 P;133P(—Eu/kBT) (2.15) where the Pu represent the probabilities assigned and the summations are over the n configurations generated. The Metropolis method is to choose P” as the Boltzmann probability itself: emp(-E,,/kBT) P = . 2.16 “ ezp(—E../kaT) ( ) With this choice, the estimate for < A > reduces to the simple expression < A >,,= n‘1 2 A". (2.17) u=l More over, the relative probability of states [1. and V is given by Pu f = F = erp{—(E,, -— E'Q/kgT}. (2.18) V This eliminates the need to calculate Pu explicitly and only the simpler quantity f is required. A possible implementation of the Monte Carlo important-sampling method is to start with some initial configuration, and then repeat again and again the following six steps which simulate the real thermal fluctuations: 1. Choose an arbitrary (e.g. random) trial configuration, 2. Compute energy difference AE connected to the initial and trial configurations. 3. Calculate the “transition probability” f = e$p(AE/k3T) for that change. 4. Generate a random number R between zero and unity. 5. If f > R, the trial configuration is accepted as the new configuration of the system, otherwise the system remains in the old state. 6. Analyze the resulting configuration as desired, store its properties to calculate the necessary averages, etc. 11 The trial configuration is in most cases chosen so as to make the sampling procedure as fast as possible. A popular choice for lattice models corresponds to a sequential or random visit of lattice sites and to the simplest possible single-site excitation. For a model with q discrete single-site states (Ising and Potts models), this choice may be a uniformly random choice among the q possible states. For a model with a continuous degree of freedom, e.g. characterized by a unit vector, the simplest possible excitation corresponds to choosing a new uniformly random direction on the unit hypersphere. All these types of excitation may be called Glauber-like excitation since they resemble the dynamics of the kinetic Ising model [15]. Use of Glauber—like excitations leads to dynamics with non-conserved order parameters. Another possible mechanism of excitation follows Kawasaki-like dynamics[16] which conserve the order parameter. A simple Kawasaki mechanism is two-site exchange of single-site properties. With these two very simple methods of excitation, it is often possible to construct, for a given problem, a mechanism which to a very high degree mimics the real physical excitations. 2.4 Finite-Size Effects The most serious drawback of a computer simulation approach to the study of phase transition is that one must, out of necessity, deal with finite systems. No finite system with a non-singular Hamiltonian can exhibit a true phase transition[1]. Nevertheless, finite systems have reminiscence of phase transitions, and systematic studies of these pseudo-transitions as functions of system size can reveal information about the phase transition in the infinite system. A finite system gives an accurate description of the infinite system as long as the correlation length, 5, does not exceed the linear extension, L of the system[17]. When 6 2 L (e.g. near a critical point), the properties of the finite system will reflect the nature of the boundary conditions. For systems with periodic boundary conditions, the fluctuations will be over-correlated and the various properties will be “rounded”. 12 As an example, the long range order will persist above the phase transition and singularities in the specific heat and the ordering susceptibility will be rounded and shifted in temperature. Fisher et al.[17, 18] have developed a finite-size scaling theory for critical phenom- ena. This theory provides a useful guide for the extrapolation of Monte Carlo data of finite-system properties to the thermodynamic limit. According to this theory, the free energy of a finite system of size L and containing N interacting entities is given by a homogeneous function F(N,T) = L‘(""‘)/".7-'(tL1/"), (2.19) where or and u are the critical exponents pertaining to the specific heat and to the correlation length, respectively; .7" is a scaling function involving the scaled variable tLl/V only, and t = (T — Tc)/Tc with Tc = T¢(L = 00). .Fisher[17] has suggested the position of the maximum of the specific heat ( or alternatively the ordering susceptibility) as an appropriate definition of the “transition temperature” of the finite system, Tc(L). According to the finite-size scaling theory, the shift in the critical temperature then scales as 6T. = T.(L = 00) — TC(L) ~ L‘l/V. (2.20) From the free energy in Eqn. 2.19, the scaling properties of other thermodynamic functions may be derived, e.g. the order parameter Q is given by = L‘5/”M(tL1/"). (2.21) In the limit, t —-> 0 and L —+ 00, Eqn. 2.21 has to reduce to the infinite system singular behavior, Eqn. 2.9, and therefore, in this limit, the order parameter scaling function is given by M(tL1/”) ~ (flay/")3. (2.22) The various scaling functions and amplitudes are non-universal properties which de- pend on the details of the system under consideration. In particular, these nonuni- versal properties are functions of the boundary conditions chosen. 13 The validity of the finite-size scalingtheory has been demonstrated through ex- tensive Monte Carlo simulations on two- and three-dimensional Ising models with periodic boundary conditions as well as free surfaces[19, 20]. It appears from these calculations that systems with L 2 10 are well inside the asymptotic region described by Eqn. 2.19. This is a very important result because it anticipates that Monte Carlo simulations of critical phenomena are feasible using system sizes and statistics compatible with current computer capacities. 2.5 Optimized Treatment of Data: Histogram Method Good and reliable Monte Carlo work requires large amounts of computer time, often of the order of 10 — 100 cpu hours on standard mainframe computers. Therefore highly optimized and careful treatment of the data obtained is needed. The requirement of minimizing the statistical error is to some extent in conflict with the need to deal with a large number of particles. The desire to study larger and more complicated systems has been the impetus behind many advances in computation algorithms and computer hardware, aimed at increasing the speed at which a simulation is performed. Advances such as vectorization and parallel processing have made it possible to study systems which would have been impossible to examine only a few years ago. A different approach to improve the efficiency is to optimize the treatment of the information obtained from a simulation. Recently, F errenberg and Swends.en showed [21, 22] that histograms can be used to greatly increase the amount of information obtained from a single computer simulation in the neighborhood of a critical point. In this section we briefly introduce the method and discuss its applications in the Monte Carlo simulation study of phase transitions. Let us consider a Monte Carlo simulation of some physical model. Each configura- tion is generated with its proper thermal weight. Eqn. 2.17 then gives the equilibrium averages of any quantities of interest. These averages are the usual output of Monte Carlo simulations. However, because the form of the probability distribution is known, 14 it is possible to extract even more information from the simulation. To see this, consider a model with Hamiltonian 'H. The probability distribution of energy E and order parameter Q at temperature T and external field h = h, is given by Pumas) = fiMEflkmd-BE + M). (2.23) where H = Iii—1" N (E,Q) is the number of configurations at the point (E, Q) in the phase space, and Z (T, h) is the canonical partition function given by Z(T,h) = Z N(E, Q)e:cp(—3E' + hQ). 3,9 The histogram of values of (E, Q) generated by the Monte Carlo simulation is propor- tional to P(T,h)(E, Q). By storing this histogram on the computer, it is easy for one to generate the normalized probability distribution. The histogram can then be used to generate data for different parameters. The normalized probability distribution with new parameters (T',h’) can be expressed in terms of the distribution with (T, h) in the following way: P(r,h)(Ea ‘I’lezl’lW' - (”E + (h, - hl‘I’l Zane P(T,h)(Ev (Plexplw' — (”E + (h' " hl‘I’l. Since T’ and h’ are continuous variable, any quantities of interest can be calculated P(Tr.hr)(E,(P) = (2.24) as continuous functions of these parameters. The denominator in Eqn. 2.24 serves as an estimate for the partition function[21, 22, 23, 24, 25]. The thermodynamic average of any variable A(E, Q) at (T’, h') can then be calculated using the following formula: < A(T,h) >= 2 A(E,)P(T,,,,(E,). (2.25) 3,4: This method is most useful in simulations in the critical region. Since a single simulation provides a continuous function of A near critical point, critical exponents can be determined very accurately. This, combined with the finite size scaling theory, can determine unambigously the nature of a phase transition, as we will see in the next chapter. Bibliography [1] C. N. Yang and T. D. Lee, Phys. Rev. 87, 404 (1952). [2] M. E. Fisher, Rep. Prog. Phys. 30, 615 (1960). [3] P. W. Anderson, Basic Notions of Condensed Matter Physics, (The Ben- jamin/ Cummings Publishing, 1984). [4] H. E. Stanley and T. A. Kaplan, Phys. Rev. Lett. 17, 913 (1966). [5] J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 (1973). [6] L. E. Reichl,A Modem Course in Statistical Physics, (University of Texas Press, 1980). [7] B. Widom, J. Chem. Phys. 43, 3898 (1965). [8] H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena, (Clarendon Press, 1971). [9] H. A. Meyer (ed.), Symposium on Monte Carlo Methods, (Wiley, 1953). [10] J. M. Hammersley and D. C. Handscomb, Proc. Camb. Phil. Soc. 60, 115 (1967). [11] K. Binder (ed.), [Monte Carlo [Methods in Statistical Physics, Topics in Current Physics, Vol.7 (Springer-Verlag, 1979). [12] K. Binder (ed.), Applications of the Monte Carlo Methods in Statistical Physics, Topics in Current Physics, Vol.36 (Springer-Verlag, 1984). 15 16 [13] O. G. Mouritsen, Computer Studies of Phase Transitions and Critical Phenom- ena, (Springer-Verlag, 1984). [14] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953). [15] R. J. Glauber, J. Math. Phys.,4, 294 (1963). [16] K. Kawasaki, In Phase Transitions and Critical Phenomena. (C. Domb and M. S. Green, eds.) Vol II. p.443 (academic Press, 1972). [17] A. E. Ferdinand and M. E. Fisher, Phys. Rev. 185, 832 (1969). [18] M. E. Fisher, In Proc. Int. School of Physics Enrico Fermi. (M. S. Green, ed.) Course LI p. 1. (Academic Press, 1971). [19] D. P. Landau, Phys. Rev. B 13, 2997 (1976). [20] D. P. Landau, Phys. Rev. B. 14, 255 (1976). [21] Alan M. Ferrenberg and Robert H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988). [22] Alan M. Ferrenberg and Robert H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989). [23] G. Bhanot, s. Black, P. Carter, and R. Salvador, Phys. Lett. B 183, 331 (1987). [24] G. Bhanot, K. M. Bitar, S. Black, P. Carter, and R. Salvador, Phys. Lett. B 187, 381 (1987). [25] G. Bhanot, K. M. Bitar, and R. Salvador, Phys. Lett. B 188, 246 (1987). Chapter 3 N2 Adsorbed on Graphite: Anisotropic Planar Rotor Model 3. 1 Introduction In this chapter we investigate the nature of the orientational order-disorder phase transition in a monolayer of N3 molecules adsorbed on graphite surface. The prop- erties of such overlayers of both spherical and nonspherical molecules adsorbed on solid surface have received enormous attention among experimentalists as well as the- orists during the last decade. In particular, the structure and the phase transitions of such systems have been of major interest. The main reason for this interest has been that these layers are realizations of two-dimensional condensed systems. From a theoretical point of view, two-dimensional systems are often more tractable than three-dimensional ones and a large body of theoretical study is available for two- dimensional models [1]. A clean surface of a bulk crystalline material presents itself to exterior gas molecules as a regular two-dimensional array of adsorption sites. Low energy electron diffraction (LEED) and neutron diffraction studies have shown that monolayers at appropriate coverages and under suitable thermodynamic conditions may form condensed phases 17 18 characterized by superstructures which may be in registry (commensurate) or out of registry (incommensurate) with the substrate. These phases may undergo phase transitions induced by changes in temperature or coverage[2]. Here we shall only be concerned with commensurate phases and order-disorder transitions in these phases. There is a vast number of adsorbed systems undergoing order-disorder transitions. However, the concept of universality (see Sec. 2.2) imposes a simplifying constraint by hypothesizing that many different transitions belong to only a small number of uniVersality classes. The most important ones are the universality classes of the ferromagnetic Ising model, the three and four state Potts models, the X Y model with or without cubic anisotropy, and the Heisenberg model with cubic anisotropy [3, 4, 5, 6]. Universality then implies that if the Landau-Ginzburg-Wilson Hamiltonian of the physical system under consideration is the same as that of one of the models listed above, the system then exhibits a phase transition of the same nature as the model. Specifically, if the transition is continuous, the critical exponents will be the same. Statistical surface physics and order-disorder phase transitions of commensurate adsorbed monolayers are conveniently studied by means of lattice models. The lattice sites represent preferred adsorption sites. Since we are not interested in properties of the substrate itself, this may simply be thought of as a symmetry-breaking field which, in combination with the interactions within the monolayer, gives rise to the lattice and an effective anisotropic interaction between the adsorbate molecules. Computer simulation has been an important theoretical technique to unravel the phase behavior of models for adsorbed monolayers. Since many powerful theoretical tools, such as mean field theory and renormalization theory, break down in two- dimensional systems, Computer simulation sometimes is the only reliable method to investigate the phase transition in these systems. A realistic microscopic interaction model of N3 adsorbed on graphite will be de- scribed in sec. 3.2. In Sec. 3.3 we review some of the earlier work on this model. 19 Our Monte Carlo studies will be presented in Sec. 3.4. 3.2 The Anisotropic Planar Rotor Model At temperatures below 47K and at coverages less than § of the number of hexagonal sites of the graphite substrate, N2 molecules form a commensurate \/§ x \/5 overlayer structure. The centers of mass of the molecules constitute a triangular lattice with the lattice points lying at the centers of the honeycombs of the substrate mesh [7]. Thus, the translational degrees of freedom are effectively frozen at low temperatures and the transition into an orientationally ordered phase is separated from two-dimensional melting. The triangular lattice of adsorption sites has a lattice spacing of a = 4.26A, i.e. V3 times larger than that of the graphite mesh. Obviously, the x/3 x s/3 structure may be placed on three equivalent positions on the honeycomb substrate. As the temperature decreased below 30K, a further transition occurs which was first detected by specific heat measurements [8]. Later, neutron scattering studies [9] showed that this transition involves the formation of a superlattice structure with a doubling of the unit cell along one of the crystal axes. Low-energy electron-diffraction (LEED) measurements [10] led to the unambiguous identification of the low-temperature phase as a 2 X 1 herringbone structure as shown schematically in Fig. 3.1. The transition around 30K thus proceeds via the orientational degrees of freedom. The herringbone structure and its thermodynamic properties are determined by the anisotropic interactions of the N2 molecules with the substrate and with each other. In the case of N3, the substrate potential is very large and negative [11]. Therefore, the interatomic axes are effectively confined to the graphite plane. The intermolec- ular interactions are, as in solid N2, well described by classical electric quadrupole- quadrupole interactions [12]. That these interactions indeed have the herringbone structure as the ground state has been demonstrated by Monte Carlo calculations [13] and by energy minimization calculations based on the electric quadrupole-quadrupole interaction and atom-atom potentials [14]. For planar quadrupoles on a triangular 20 lattice, this interaction takes the form of an anisotropic planar rotor Hamiltonian [15] H = J 2 cos(2<,o,- —- 290,-) + K 2 608(290,‘ + 2% — 40,-,), (3.1) where the summations are over the six nearest-neighbor pairs only. The angle 30.- measures the orientation of the ith molecule and 0,3 specifies the direction of the line connecting sites i and j of the lattice (cf. Fig. 3.1). For pure electric quadrupole- quadrupole interactions, the ratio of the coupling constants is K / J = 35/3 and the isotropic X Y-like term in eqn.(3.1) can therefore be neglected. It is the remaining anisotropic term which favors the herringbone structure. The coupling constant K is related to the nearest-neighbor electric quadrupole-quadrupole coupling constant I‘o through the relation 875 875 6(eQ)2 . K=—I‘=—————=33i5Kl 3.2 64 ° (64) 25.13 ( ) 8W“ ( l where we have used that the electric quadrupole moment of N3 is eQ = (1.4 :i: 0.1) x 10‘2°esu cm2 [16]. With the value of K given in Eqn.(3.2) and with J : 0, the anisotropic planar rotor Hamiltonian is believed to be a highly realistic lattice model (with no adjustable parameters) to describe the orientational properties of N2 physisorbed on graphite. The order parameter of the 2 x 1 herringbone structure has three real components corresponding to the three-fold symmetry axis of the triangular lattice. Therefore there are six thermodynamically equivalent ground states. The components of the order parameter may be written as N a(T) = N“ < |Zsin(2cp,- — 20a)e‘q°"‘| >,a = 1,2,3, (3.3) 121 where (1,, are the wave vectors describing the herringbone ordering and Oa’s are the directional angles of the triangular lattice 217 —. 0, —— 0 = 0 q‘ ( .fi)’ 1 1r 7r 211' ‘12 — (—;’—a_\/—§)’02_3 (3.4) 21 Figure 3.1: Schematic representation of the herringbone structure on a triangular lattice. The planar rotors represent the interatomic axes of diatomic molecules N; physisorbed in a commensurate \/3 x J3 (referred to the graphite unit cell) over- layer on the hexagonal substrate of the graphite. The angles «,0,- and 49,-, enter the Hamiltonian in Eqn. (3.1). (3.5) The nature of the phase transition from the herringbone phase to the orientationally disordered phase has been a subject of great theoretical interest and will be discussed in the following sections. 3.3 Previous Work Using Mean Field Theory, RG, MC, MD simulation Since the transition from the (2 x 1) herringbone ordering can be described by a system with a 3-component order parameter, to facilitate the discussion of the phase transi- tions of the anisotropic planar rotor model, it is useful to introduce the Heisenburg model with cubic anisotropy. Through this model, two universality classes related to this phase transitions may be defined. The Heisenberg model with cubic anisotropy is characterized by the Landau-Ginzburg-Wilson Hamiltonian [6] 3 [1(a) = r 2 31+ a=l p "M“ ... won)” + 11(sz 18,)“ + v i: o; (3.6) with a = 1, 2, 3 referring to the three components of the ordering field Q = ((1)1, (1’3, (1)3). The last term in Eq. 3.6, which is invariant only under the symmetry operations of the cube, breaks the rotational symmetry of ‘5. Under a renormalization group trans- formation, this term Will become dominant implying that only the stationary points of the cube will play a role [17, 18, 19]. Thus Q will point either to the faces or the corners of the cube. In these two cases, the fourth order coefficients will take on the values (a + 1))]§|4 and (u + §)|§|‘ Both these fourth order term will have to be positive in order for the Hamiltonian to stay thermodynamically stable under the renormalization group transformation. For 12 < 0, the faces are favored and for v > 0, the corners are favored. This leads to the notion of face and corner type cubic anisotropy. 23 The two types of cubic anisotropy give rise to two universality classes. Renormal- ization group calculations of (4 - 6) type predict that for 2-d systems, in the case of face-type cubic anisotropy, the phase transition is always of first order[19, 20, 21]. In the case of corner type anisotropy, the transition may either be continuous or of first order. If it is continuous, the critical behavior is predicted to be that of the two-dimensional Ising model [6, 22]. It can be shown that the Landau-Ginzburg-Wilson Hamiltonian, which may be con— structed from the Herringbone order parameters in Eqn. 3.6, is that of the Heisenberg model with face type cubic anisotropy [23]. Therefore by universality, the renormal- ization group argument predicts a first order transition for this planar rotor model. However the mean-field approximation [23] and recently, the calculation using cluster variation method[24], predict a continuous transition. Recently aniir and Piercy[25] studied the phase transition of p(2 X 1) structure of oxygen on Ru(001) surface us- ing LEED technique. This system has a three-component order parameter and is believed to be in the universality class of the (n=3)-c0mponent Heisenberg model with face type cubic anisotropy. Contrary to the prediction of RG calculation, the transition was found to be continuous with critical exponents close to the three-state Potts model. Both mean-field theory and renormalization group calculations tend to give inac- curate results in 2-dimensional systems. Therefore computer simulations have been particularly useful in resolving the nature of the phase transition in this situation. The first extensive Monte Carlo simulation on the anisotropic planar rotor model was performed by Mouritsen and Berlinsky[15]. They claimed to have found a fluctuation driven first order phase transition at T/K = 0.775, based on an apparent disconti- nuity in the order parameter at Tc for large systems (L ~ 100). However the short simulation time (500 ~ 5000 Monte Carlo steps/site) they used may have caused large statistical error and could have given inaccurate results. For example, the tempera- ture dependence of the specific heat they calculated showed a broad peak and was 24 found to be insensitive to the system size for L 2 20, which actually indicated that the specific heat does not diverge at the “critical temperature”. Subsequent Monte Carlo work performed by Allen and O’Shea[26] use a larger system and longer Monte Carlo sweeps on the model system using the full quadrupole-quadrupole interaction in Eqn. 3.1. However no conclusive evidence about the nature of the transition was given, although their study of the fractal properties of the cluster distributions give some weak evidence in favor of a continuous change while RG calculation predict a first order transition for their model[27]. Cai and Mahanti [45] performed Monte Carlo simulations using a Gaussian ensamble technique introduced by Challa and Hetherington[29] (which is particularly suitable to detect first order phase transition) using small lattices and one million MC sweeps. Again the results were inconclusive but give some weak evidence in favor of a continuous transition. 3.4 Present Monte Carlo Studies To discuss the issue of the precise nature of the phase transition in the anisotropic planar rotor model, we present an extensive simulation study of the this model supple- mented by finite size scaling ideas. There are two issues here: first, the identification of a continuous or first order transition and second, if the transition is continuous, to evaluate the critical exponents in order to determine which universality class it belongs to. The method we have used to determine the nature of phase transition depends on the identification of the size dependence of the forth order cumulant at Tc. Binder has argued that this quantity is very sensitive to the nature a phase transition[30]. For the energy, the fourth order cumulant V1,, is defined by V =1————-—. L 32 (3.7) For a continuous transitions, VL —-> 32; for all temperatures as L —> 00, while for a first-order transition VL —> g for high and low temperatures but tends towards 25 a nontrivial value at the transition temperature [31]. For q=10 Potts model, for example, VL = 0.51. The problem of the precise determination of the peak position and the height of V1, for large system , where the peaks are very narrow, is solved by a technical advance due to Ferrenberg and Swendsen [32, 33]. The histogram method they have recently developed enables us to obtain complete thermodynamic information over the entire scaling region near the phase transition using MC simulation at a single temperature point. Therefore we can determine the position and height of various quantities precisely with reasonable computation effort. For a detailed discussion of the histogram method, see Chapter 2 of this thesis. Monte Carlo simulations are performed on a triangular lattice with periodic bound- ary condition, the system size ranging from 12 x 12 to 120 x 120. 1—5 million Monte Carlo sweeps per site are used to get accurate energy histogram after discarding 5 X 10‘ Monte Carlo sweeps per site to get the system into equilibrium state. In Fig. 3.2, VL calculated from simulation at T = 0.76 for each lattice size is shown plotted versus temperature. It indicates clearly that in the thermodynamic limit, VL —v g strongly suggesting a continuous transition. It should be pointed out, however, that the theory of the fourth order cumulant is phenomenological and has only been tested for Monte Carlo simulations of q-state Potts model [31]. Thus it is important to calculate the critical exponents to determine the nature of the phase transition and the universality class. Fig. 3.3 shows the temperature variation of the specific heat for different system sizes. Contrary to the simulation results of Mouritsen et al.[34], we find that the specific heat peak has a strong size dependence. As the lattice size increases, the peak grows while the width narrows. The difference between the present results and the earlier ones by Mouritsen et al. is mainly due to large statistical errors in Mouritsen et al.’s simulation because of short Monte Carlo runs (~ 5000MCsteps/site) compared to ours (~ 5 x 10°MCstep3/site). Our results are also consistent with the existence o.seses7.ref,....I [fl 0.6666 I—J > _ -\ , / 0.6665 - \ l. — . \ / 4 ,. \ / . ( / . '\ / 0.6664 I l I A l 1": l l A I I l L L L l 0.7 0.75 0.8 0.55 0.9 Figure 3.2: Plot of Binder’s forth order cumulant V1, versus temperature for the two dimensional anisotropic planar rotor model of L x L lattices while L = 60 (dash-dot); L = 90 (dashes); L = 120 (solid). The results are from single simu- lations at T/K = 0.76. 27 of a finite temperature phase transition in this 2-dimensional anisotropic planar rotor model because they indicate that the specific heat of the system diverges when L ——> 00. Fig. 3.4 shows the lattice size variation of the specific heat peak. It shows that the finite size scaling theory is valid and C ~ La/V. By least square fitting of this log-log plot and using the hyperscaling law [35]: ud=2—a, (am where d = 2 is the spatial dimension, we find that a = 0.335 i 0.002;V = 0.832 :f: 0.002, (3.9) while for the first order transition we would have C ~ L2. Thus the above value of a and V give a clear evidence that this system undergoes a continuous transition. Fig. 3.5 shows the temperature variation of the susceptibility of the order parameter (as defined in Eqn. 2.6) for different lattice sizes and Fig. 3.6 shows the lattice size dependence of the susceptibility peak. According to finite size scaling assumption, the susceptibility peak diverges in the thermodynamic limit as x ~ LV". From a least square fitting and the value of exponent u from specific heat simulation, we find 7=L6i02 (am) Using the scaling law described in Chapter 2, we can thus find [3. Table. 3.1 lists the critical exponents of the two dimensional anisotropic planar rotor model, the experimental data on the system with the same universality class as this model [25] and that for a 2-d 3-state Potts model. These results unambiguously indicate that the phase transition of the two dimensional planar rotor model, or, in a wider context, the models which belong to the universality class of Heisenberg model with cubic anisotropy of the face type, can have a continuous transition which belongs to the universality class of q=3 Potts model in two dimension. The transition temperature can be determined by extrapolating TC(L) as defined by the peak position of either the specific heat or the susceptibility as a function of 28 4 r fr r l v w. 3 "' d ' / . / ' . II' 2- I \‘_ - 1—1 _ \ D I \' a ' . \\ I -’ \ \ / . 1.. . . \ \. _ ' l / \ \ -. I \ \ 1 _’ / \ \ / l \ \I\ i 0 . l ‘ 5.]. __ _ 0.7 0.75 0.8 I 0.85 0.9 T/K Figure 3.3: Plot of specific heat 0;, versus temperature for the two dimensional anisotropic planar rotor model of L x L lattices while L = 60 (dash-dot); L = 90 (dashes); L = 120 (solid). The results are from single simulations at T/K = 0.76. 29 Cmax(L) Figure 3.4: Plot of specific heat peak Cm(L) versus L for the two dimensional anisotropic planar rotor model of L x L lattices at the critical temperature Tc( L). The solid lines are the least square fit. 30 400 . . . , l . . , r l J 300 _ 200 _ u-J X 4 100 _. 0 . '7 0 . 9 Figure 3.5: Plot of susceptibility XL versus temperature for the two dimensional anisotropic planar rotor model of L x L lattices while L = 60 (dash-dot); L = 90 (dashes); L = 120 (solid). The results are from single simulations at T/K = 0.76. 31 1000 _ I f l r I r l I '-' 1 '- 4 A 100 :- :4 A - 4 v * d N " .4 g " -+ X - '1 10 l l L l l i l l l 10 100 Figure 3.6: Plot of susceptibility peak xm(L) versus L for the two dimensional anisotropic planar rotor model of L x L lattices at the critical temperature T¢(L). The solid line is the least square fit. 32 Table 3.1: Comparison of theoretical and experimental critical exponents for the model which belongs to the universality class of Heisenberg model with face type cubic anisotropy compared with the critical exponents of q=3 Potts model in two dimension. Exponent Planar rotor model(theory) Experiment q=3 Potts model a 0.335 4: 0.002 0.30 i: 0.06 0.333 (g) ,3 0.03 a: 0.05 0.13 i 0.02 0.111 (g) 7 1.6 :1: 0.2 1.35 i 0.15 1.444 (35‘!) V 0.832 d: 0.002 0.74 :i: 0.08 0.833 (g) the lattice size L to the L —> oo limit. From finite size scaling theory we have TC(L) ~ Tc(oo) + L-l/V. (3.11) F ig.3.7 shows the lattice size dependence of Tc as defined by the specific heat and susceptibility peak positions. We can see that the peak positions of specific heat and susceptibility are quite different for any finite size system. However they approach to the same value Tc(oo) = 0.76 in the thermodynamic limit. 3.5 Summary and Conclusion In summary, We have carefully analysed the herringbone orientational transition of N 2 molecules adsorbed on graphite (001) surface by Monte Carlo simulation. Using the recently developed histogram method combined with analysis of finite size scaling and the fourth order cumulant of energy, we have identified unambiguously a continuous transition in this model contrary to the real-space RC results. The critical exponents have been determined to be that of q=3 Potts model in two dimension and are shown to be in good agreement with the recent LEED experiment on a system which 33 0.78 ....l 0.77 0.76 fl I I I \‘l V U U I r I I V I 0.75 ..‘.1-...1....1....1.... 0 0.002 0.004 0.006 0.008 0.01 L-l/V Figure 3.7: Plot of critical temperature T¢(L) determined by the peak position of susceptibility (x) or the peak position of specific heat (0) versus L. The lines are the least square fit. 34 belongs to the same universality class, i.e. the Heisenberg model with face type cubic anisotropy. These results indicate that although the Ginzburg-Landau-Wilson analysis of the universality may still be valid in two dimension, the 4 — 6 expansion may give unphysical results in certain types of models. There still some points, however, which may cast doubts to our results. From the above simulation we still cannot exclude the possibility that 3 = 0, which may indicate a discontinuity of the order parameter at the transition temperature. More precise simulation is needed to clarify this issue. It is still no plausible physical arguments such as why the anisotropic planar rotor model has a continuous transition while the 2-d q = 6 Potts model, which also belongs to the universality class of the Heisenberg model with face type cubic anisotropy, has a first order transition. One can argue that although these two models both have six degenerate ground states and the GLW Hamiltonians have the same form, the energy needed to create domains are quite different. One can also argue that this difference may be caused by the frustration of bonds in the anisotropic planar rotor model in triangular lattice. The mechanism which drives this transition to second order can only be revealed by careful analysis which will be the aim of the future studies. Bibliography [1] S. K. Sinha (ed.), Ordering in Two Dimensions, (North-Holland, 1980). [2] P. Bak, Rep. Prog. Phys. 45, 587 (1982). [3] E. Domany and E. K. Riedel, Phys. Rev. Lett 40, 561 (1978). [4] E. Domany and M. Schick, Phys. Rev. B 20, 3828 (1979). [5] M. Schick, Physica 109-1103, 1811 (1982). [6] M. Schick, Surf. Sci. 125,94 (1983). [7] J. K. Kjems, L. Passell, H. Taub, J. G. Dash, and A. D. Novaco, Phys. Rev. B 13, 1446 (1976). ‘ [8] T. T. Chung and J. G. Dash, Surf. Sci. 66, 559 (1977). [9] J. Eckert, W. D. Ellenson, J. B. Hastings, and L. Passell, Phys. Rev. Lett. 43, 1329 (1979). [10] R. D. Diehl, M. F. Toney, S. C. Fain, Phys. Rev. Lett. 48, 177 (1982). [11] W. A. Steele, J. Phys. (Paris) Colloq. 38, 04-61 (1977). [12] T. A. Scott, Phys. Rep. 27, 89 (1976). [13] S. F. O’Shea and M. L. Klein, Chem. Phys. Lett. 66, 381 (1979). [14] C. R. Fuselier, N. S. Gillis, and J. C. Raich, Solid St. Commun. 25, 747 (1978). 35 36 [15] O. G. Mouritsen and A. J. Berlinsky, Phys. Rev. Lett. 48, 181 (1982). [16] A. D. Buckingham, R. L. Disch, and D. A. Dunmur, J. Am. Chem Soc. 90, 3104 (1968). [17] E. Domany, E. K. Riedel, Phys. Rev. B 19, 5817 (1979). [18] B. Nienhuis, R. K. Riedel, and M. Schick, Phys. Rev. B 27, 5625 (1983). [19] P. Pfeuty and G. Toulouse, Introduction to the Renormalization Group and Crit- ical Phenomena, Chap. 9, (Wiley, 1977). [20] A. Aharony, Phys. Rev. B 8, 4270 (1973). [21] I. J. Ketley and D. J. Wallace, J. Phys. A 6, 1667 (1973). [22] G. s. Crest and M. Widom, Phys. Rev. B 24,6508 (1981). [23] A. B. Harris and A. J. Berlinsky, Can. J. Phys. 57, 1852 (1979). [24] E. Chacon and P. Tarazona, Phys. Rev. B 39, 7111 (1989). [25] H. aniir and P. Piecy, Phys. Rev. B 4, 582 (1990). [26] M. P. Allen and S. F. O’Shea, Mol. Sim. l, 47 (1987). [27] P. Tarazona and E Chacon, Phys. Rev. B 39, 7157 (1989). [28] Zhi-Xiong Cai and S. D. Mahanti, unpublished. [29] Murty S. S. Challa and J. H. Hetherington, Phys. Rev. Lett. 69, 77 (1988). [30] K. Binder, Phys. Rev. Lett. 47, 693 (1981). [31] Murty S. S. Challa, D. P. Landau, and K. Binder, Phys. Rev. B 34, 1841 (1986). [32] Alan M. Ferrenberg and Robert H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988). 37 [33] Alan M. Ferrenberg and Robert H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989). [34] O. G. Mouritsen, Computer Studies of Phase Transitions and Critical Phenom- ena, (Springer-Verlag, 1984). [35] H. E. Stanley, Introduction to Phase Transitions and Critical Phenomena, (Ox- ford University Press, 1971). Chapter 4 Lattice-Gas Melting of Two-Dimensional Alloys: Application to Ternary Graphite Intercalation Compounds 4.1 Introduction The weak interlayer coupling in graphite gives rise to an extensive class of materials known as graphite intercalation compound (GIC’s), in which foreign species maybe introduced between the layer planes of the host matrix [1, 2]. This subject spans the fields of physics, chemistry and materials science [3]. A great rebirth of interest in GIC research can in part be traced to the availability after 1960 of a synthetic carbon commonly referred to as highly oriented pyrolytic graphite (HOPG) [4]. Exceptional attention from the chemists has also been paid in recent years to GIC’s, the reasons of which are manifold. Of practical consideration was the possibility of synthesizing compounds from plentiful materials (e.g. carbon) which were lightweight and had very high electrical conductivity of an order of magnitude or greater than that of 38 39 copper. Of more fundamental concern is the quasi-two-dimensional character of graphite and its intercalation compound. They may therefore serve as a test ground for phys- ical concepts in quasi-one and two-dimensions [6]. The basic structure of the graphite intercalation compounds is discussed in several review articles [5, 4, 6, 7] and is shown in Fig. 4.1. For pure graphite the hexag- onal carbon layers are stacked parallel to each other in the staggered configuration ABAB of Fig. 4.1(a). The honeycombclike in-plane carbon-carbon bond length is 1.42 Awhile the layer separation is 3.35 A. Intercalating graphite leads to an ordered array of 71. carbon layers (C) followed by an intercalated layer (I). This is called a stage-n compound, Fig. 4.1(b). Here the stage index, n, denotes the number of C- layers between two adjacent I-layers. For example, a stage-3 GIC has the sequence, ICCCICCCICCC, Fig. 4.1(b). For an alkali-metal stage-1 GIC, except for lithium, the general chemical formula of MOS reflects the nominal in-plane layer density of eight carbon atoms per alkali metal atom M. In the higher stage GIC’s (n _>_ 2) the correlation between the I-layers and there bonding to the C-layers become sufficiently weak that the intercalated layers in graphite resemble physisorbed monolayers on graphite as we discussed in the previous chapter. Therefore, intercalation physics is in a certain sense equivalent to surface physics with the important advantage of dealing with multiple “surfaces” within the bulk[6]. Over the years, substantial understanding of these materials have been accom- plished. Never-the-less, only recently the structural properties of ternary GIC’s which have two foreign species within the same graphite layer has been studied. In this chapter we discuss the results of a systematic investigation of a spin-1 Ising (lattice—gas) model on a 2-dimensional triangular lattice with nonzero interactions up to third nearest neighbors. Our motivation for studying this model is threefold. First of all, there has been considerable interest recently in understanding the structural properties of ternary graphite intercalation compounds (GIC’s) of the type A,B1_,CB, 40 ut sue: mine: in sac: 4th snot O... O... O...‘ O...‘ O... . aa— 3— —_ ‘ O... A ‘ O... 3 J 0000 0— ‘— 0000.] O O a— ° ' c—— 3— C . CARBON HEXAGON LAY“ e o e a LAYER 0F INTEKCALATE (b) Figure 4.1: Structure of prystine graphite and its stacking sequence; (b) the stacking sequence of carbon and intercalate layers in stage-1 to stage-4 compounds. [6] 41 where A and B are alkali-metal atoms[4, 8]. These systems appear to be well described by a lattice-gas picture in which the intercalant atoms occupy specific sites between the layers of the graphite substrate[8]. This lattice gas model can be mapped into a spin-1 Ising model. There is a charge transfer from alkali-metal atoms to the graphite layers and one expects the effective interaction between pairs of intercalant alkali-metal ions to be of the screened Columb type. While there is some interaction between intercalant atoms in different layers, this interaction is much smaller than the dominant intralayer interaction. Therefore for an initial attempt to understand the lattice gas melting problem, we restrict ourselves to ternary systems with intralayer interaction only. It is then possible to investigate some of the properties of such systems on the basis of a two-dimensional spin-1 Ising model[lO]. Our second motivation for this study was to investigate the usefulness of employing Monte Carlo simulation techniques to study the lattice gas melting problem and the nature of different ground state structures in ternary GIC’s. Recent experiments on ternary GIC’s by Solin et al.[9]have indicated phonon softening of K1_,,Rb,,Cg for 3: ~ 3. They suggest that clusters of 2\/3 x 2\/3 ordered structure might be present at low temperature. By studying this model in detail, we wanted to understand not only the concentration (3:) dependence of the lattice gas melting temperature Tm(a:), but also address the issue of the'existence of above clusters at low temperatures. and provide useful directions for future experimental work. Finally, the model that we have studied is a generalization of the well known Blum- Emery-Griffiths (BEG) model[lO]. Properties of this model have not been studied previously for a triangular lattice and/ or for nonzero magnetization and its critical properties are of basic theoretical interest. The structure of this chapter is as follows. In Sec.4.2 we describe some of the previous work in this field. In Sec.4.3 we describe the model and discuss the values of its parameters appropriate for the ternary GIC’s for a fixed :0, rather than for fixed chemical potential. Sec.4.4 presents the results for the concentration dependence of 42 Tm obtained within a mean field approach. In Sec.4.5 we present the Monte Carlo (spin-exchange method) simulation (the so called Kawasaki dynamics[15]) results for the melting transition and the order-disorder phase transition on the 2 x 2 superlattice at temperature below Tm, the latter transition has been studied for one value of a: _ 1 (a: — 3) only. In Sec. 4.6 we summarize this work and present our conclusions. 4.2 Review of Previous Work on Ternary GIC’s There has been considerable work on alkali-metal GIC’s in the last decade [4]. Most of these concentrate on binary GIC’s and only recently some experimental work has been done on ternary GIC systems. Solin et al [9] found anomaly in Raman spectrum for K1_3Rb¢Cg when a: = 32,- and they ascribed the origin of this anomaly to the presence of ordered clusters in this system. Similar anomalies have also been seen in other properties of this ternary system [11, 12]. Recently, Carnerio et al. [13] have proposed a three-state Potts lattice-gas (PLG) model to study the order-disorder phase transition of mixtures of 3He-‘iHe adsorbed on graphite. These systems can be represented by (A3B1_¢)V2 where V is a vacancy site, A and B refer to 3He-4He respectively. On the other hand, the Blume-Emery-Grifiiths model [10] has attracted much attention for its ability to describe ternary systems. Particularly, Lee and Landau [14] proposed a modified BEG model to study the phase diagram of model ternary alloy (AB),,V1_,, system on a square lattice. The reason that the BEG model has not been applied to a triangular lattice and for general (AcB1_,)yV1_y systems is that the spin Hamiltonian describing such a system becomes complicated and it is difficult in practice to determine all the inter- action parameters in a reasonable region so that the results can be compared with experiment. Furthermore, the chemical potential of the system usually changes with temperature in experiment, thus assuming a constant chemical potential in this model is unreasonable. We shall show in Sec. 4.3 that these difficulties can be overcome by doing the Kawasaki dynamics[15] instead of spin-flip Monte Carlo simulation [16] 43 (the so called “Glauber dynamics”). 4.3 Lattice Gas Model And Its Parameters From X-ray structure measurements [1, 9] we know that the alkali-metal ions interca- lated into graphite layers form a stage-1 ordered 2 x 2 superlattice at low temperature (see Fig. 4.2). At high temperature, these systems undergo an order-disorder'phase transition and become liquid like. The graphite substrate on which alkali-metal ions sit defines a periodic potential with well-defined minima at the center of the hon- eycomb lattice for possible adatom occupation. We shall restrict ourselves here to the case of very deep potential wells which effectively force the intercalant atoms to occupy positions on a rigid lattice formed by this substrate potential. Two types of intercalant atoms are present and hence each triangular lattice site formed by the center of the hexagon of Fig. 4.2 may be filled with either intercalant species or remain empty. The general expression for such an alkali-metal GIC system can be written as A¢B1_2V3. A simple Hamiltonian may then be written in terms of site occupational operators Pi", where P} = 1 if site i is occupied by species A and zero otherwise. The generalized Hamiltonian is that of the form[14] H = Z Z 1*”(r,,r,)P3P;’ + Z: Vin-010,.A + H0, (4.1) (I‘VWJ) A 5 where the parameter I A” describes the interaction between species A and u on sites i and j respectively, and the second term in Eq. 4.1 describes the binding of the intercalants to the substrate. Other thermodynamically “irrelevant” terms are in- cluded in Ho. This model can be related to a spin-1 Ising model by the following transformations: A 1 2 Pi = 5(S£+Si)3 (4-2) 1 P? = E(53—5.), (4.3) Figure 4.2: (2 x 2)R0° triangular lattice on a graphite substrate. The solid circle represents the occupied sites when temperature T < Tm(z:), the melting temperature. 45 pi” = l—Sf, (4.4) I» where S; = 1,0, —1 is the spin-1 Ising variable associated with the ith lattice site, A and B are the two intercalant species and V refers to the vacant sites. The resultant Hamiltonian which represents a generalized BEG model is H : Z: K134522533 + Z JngiSj + Z LijSESj (i,j) (57.7,) (‘91.) +ZA45? +2145; + H6, (4-5) where K,, 2 £03.44 + 1,33 + 213.3), (4.6) J,-,- = Eat“ + 158 — 213.3), (4.7) Lij = $15,314 - 158), (4-8) A.- = $04" + V.” — 2W“). (4.9) m = $04" —- v.3). (4.10) The concentration a: of A atoms is given by _ZiSi(1+5i) _ 225512 ' (4.11) The ground-state properties of H obviously depend upon the values of the interac- tion parameters and we can alter the region of stability of different ordered phases or even eliminate some of them completely by choosing the parameters suitably. Since the interactions between intercalant atoms are probably for the most part indirect and therefore difficult to calculate theoretically, the estimates of the actual interaction parameter appropriate for the system should come from a comparison between calcu- lated and observed structures. In this paper we restrict ourselves to nearest-neighbor (NN), next-nearest-neighbor (NNN), and third-nearest-neighbor (3NN) interaction parameters. We choose KNN = 1, KNNN = 0.3, J3NN, and L3NN about one-tenth of K NNN, and all others are set equal to zero. This is because we intend to apply this 46 model to ternary systems such as K1_,Rb.,Cg, where we expect the interaction to fall off rather rapidly with interparticle separation. In addition as we will see, this choice of parameter values leads to a ground state with 2 x 2 superlattice structure stable over a considerable range of concentration a: as seen experimentally. This indicates that if a generalized BEG model is indeed appropriate for these systems, then the nearest- and the next-nearest-neighbor coupling (NN and N NN) are predominantly repulsive but an attractive third-nearest-neighbor coupling between unlike species is present. We have set L nonzero since in experiments [4] KC; and Rng have different melting transition temperatures. It is easy to see that the L term distinguishes AV3 from B V3 and acts like an effective field. We assume that the concentration a: is unchanged during the melting process and use canonical ensemble so the last three terms of the Hamiltonian become constant and can be ignored. The resulting Hamiltonian which we the finally study is H = K1 2 535,? + K2 2 535,? + J 2 5.5,- + L E 535,-, (4.12) NN NNN 3NN 3NN Where K1 = KNN, K; = KNNN: J = J3NN, and L = LBNN- The energy (per site) of different ordered structures (not including the single-site contributions) shown in Fig. 4.3 (only the sites on the 2 x 2 superlattice are shown) and for a: = 0, 1 (AV3, BV3) are as follows: 19..., = i(J+L), (4.13) E”, = $(J—L), (4.14) EM”, = —11—2(J+L), (4.15) 13A,”, = —11§(J—L), (4.16) E43,“, = —gL, (4.17) EAJBVn = 3L. (4.18) We can see that the ground-state energy of the system is independent of K1, K2, Slnce no two intercalant atoms are either nearest or next-nearest neighbors of one 47 another. Since in the model we described above there is only “nearest-neighbor” antiferromagnetic coupling on the 2 x 2 superlattice, the ordered structures shown in Figs. 4.3(c) and 4.3(d) can only exist at zero temperature. 4.4 Mean Field Approach In order to understand qualitatively the a: dependence of melting transition within this model we have developed a mean-field theory based on the two sublattice model (see Fig. 4.4). As we are interested in the relation between the melting transition temperature and the concentration :0, we use the Bragg-Williams approximation [17] so that the concentration (magnetization in the Ising language) is fixed. Since long-range order on the 2 x 2 superlattice has not been found in experiments for any 2:, it seems reasonable to assume that the A, B order-disorder phase transition within the 2 x 2 superlattice occurs at much lower temperature than the melting tran- sition. For certain concentrations, this transition may be absent altogether. Monte Carlo simulation work (see Sec. 4.5) confirmed this assumption. Therefore, near the melting transition temperature the occupation probability of A and B ions within the 2 x 2 superlattice can be regarded as completely disordered. Thus the two sublattice picture (see Fig. 4.4) is assumed to be valid for all values of a: (0 < a: < 1). Let us define the probability of finding an atom A in the sublattice V as P)“, (A = A, B, V; V =I,II). Then we have in terms of W, the long range order parameter, PA, = 2(1— W)a: + Wc, (4.19) P3, = in — W)(1 — 2:) + W(1— 3:), (4.20) PW = 2(1- W), (4.21) PA” = 3(1— W), (4.22) P3” = in — 3:)(1— W), (4.23) PV” 2 1(3 + W). (4.24) 4 48 Figure 4.3: Various ordered superlattice structures for ternary systems. The lattice points represent the solid circles in Fig. 4.2. (a) 14.83%; (b) AgBVg; (c)ABaVu; (d) A33 V". Unoccupied sites of Fig. 4.2 are denoted as vacancy sites V and are not shown in the above figure. AVAVAYAVAVAVA AVAYAYAYAYAVAV . AVAVAYAVAYAVAV AVAVAYAYAYAYAV AVAVAVA'AYAVAV ‘ i AVAVAVAVAVAVAV VAVAVAVAVAVAV 3:13:33: :. Figure 4.4: The two sublattice model. In the above figure representing the ordered 8 ure of stage-l Avg, 8V3, A.Bl-.V3 systems, A and B atoms always occupy the sublattice I sites (a) and vacancies occupy the sublattice II sites (0). The ground state for either AV3 or BV; is a (2 x 2)R0° structure. 50 The free energy per particle in the Bragg-Williams approximation (mean-field theory) is given by F = gag + K2)(PAI + PB;)+;(K1 + K2)(PAII + PIE!!!)2 +%J[(PA, — P13!)2 + 3(PAII - P3102] +2L[(P§I — P§,)2 + 3(P211 _ 131911)] +iT[(ZPz1n P.) + 3(2 P,- 111 PM, (4°25) where i indicates AI, BI, VI and j indicates AII, BII, VII. To locate the transition temperature[9] we expand the free energy given by Eqn. 4.25 up to order W4 as 3 F = C + Iii-[3.194; —- 1)2 + 3L(2:r — 1) — (K1 + K.) + 8T]W2 +;—:—TW3 4. 0(W“), (4.26) where C is a constant independent of W. Because of the third-order term, the transition is first order. However, to under- stand qualitatively the a: dependence of the melting transition, we will first analyze the behaviour of T"(:c) where the second Landau coefficient (the coefficient of the W2 term in Eqn. 4.26) goes to zero. T‘(:c) is given by K T‘ = 42832 _ §[J(2:c — 1)2 + 13(22: — 1)]. (4.27) Since 82T —— = —12J .28 6232 9 (4 ) and for J > 0, g < 0, T‘ is higher than that predicted by the Vegard’s law according to which T’(:r) = (1 — 3:)T'(0) + xT'(1). (4.29) In Fig. 4.5, the order parameter W versus the reduced temperature T (measured in units of Icy/Kl, where 103 is the Boltzmann constant) is plotted for J = ngv = 0.05, 51 L = LsNN = 0.01, and a: = §. As discussed above, the third-order term in Eqn. 4.26 makes the melting transition first order for all values of 2:. The phase boundary separating the “liquid” and “solid” phase Tm(a:) in the mean-field approximation is given in Fig. 4.7 along with the Monte Carlo simulation results (to be discussed in Sec. 4.5). 4.5 Monte Carlo Simulation 4.5.1 Melting Transition We have performed Monte Carlo simulation of the generalized BEG model with Hamiltonian Eq. 4.12 for a triangular lattice. In order to study the relationship between the concentration :1: and the melting transition temperature Tm, the average concentration :0 has been kept constant. Accordingly, the Monte Carlo simulation has been performed by trials of exchange of two spins on the nearest-neighbor sites (Kawasaki dynamics) [15]. Compared with the spin-flip Monte Carlo (MC) method, the Kawasaki dynamics operate with a much simplified Hamiltonian since many fac— tors in the general Hamiltonian remain constant thus enabling the interaction param— eters to be determined with relative ease. It is also a more appropriate way to do the simulation since in experiments the concentration a: usually does not change much at the transition temperature, rather the chemical potential changes with temperature. The order parameter is defined by the summation of correlation function of 52(r,-) in the q space[16]. Since at the present stage we are mainly interested in the melting transition, the order parameter is expressed in terms of 5,2 rather than 5,, thus it does not distinguish whether an A or B ion occupies the ith site. In terms of the two quantities 5., and Cq given by 5,, = 252(r,)eiq"‘, (4.30) i 80 and C, :< .25-, >, (4.31) 52 . v v I I I f U V U I I I r T r V U I r l I T ' .' I ' r r I '-°° Janos 1 . Lao.“ . - X31]! 4 I d D 1 h 1 D . 1 3 ' l o.” [- '1 E d b 0 A A A L l L l A LL m A_ A L l A A_L LLA L L A l a L m 0.10 0.20 0.30. 0.40 0. 50 0. 00 0. 70 T Figure 4.5: Temperature (T) dependence of order parameter W' obtained in mean field theory indicating the lattice-gas melting for parameter values J = 0.05, L = 0.01, and a: = 3‘, for the ternary A,B1_.V3. 53 the order parameter mg can be defined as 3 mG = [Z Cc,/3C,=o]%, (4.32) i=1 where quo = N/16. The size dependence of C, may signal the difference in the ordered and the disor— dered structures: Cq ~ 0(1) except at q=0 in the disordered phase, and C,, ~ O(N) for the commensurate wave vectors q 2 G1, G2, and G3 in the ordered phase[16]. Above the melting temperature Tm, the correlation function 0,, has peaks at q=a', q=b‘, andq:a‘+b‘, where the reciprocal-lattice vectors corresponding to the trian- gular lattices of Fig. 4.2 are a’ = 23(1, —1/\/§),b* = mom/3). (4.33) Where the lattice constants has been chosen as 1. Below Tm, additional peaks appear at the 2 x 2 commensurate wave vectors within the first Brillouin zone. These are G1 = aye. = 92—h. : :44: + b‘). (4.34) Simulation has been performed on a 36 x 36 lattice with periodic boundary condi- tions. The two dominant interaction parameters are chosen at K1 : 1.0 and K2 = 0.3. For J and L we choose the sets (J 20.05; L =0.01,0.02,0.03,0.05,0.08, and 0.10) and (J 20.0; L 20.05 and 0.10). 5000 MC steps per NN bond are performed. Since it is difficult to find the ground state of these ternary systems for arbitrary 2:, we have performed Monte Carlo simulation by heating the system to a very high tem- perature first to get a completely disordered state (liquid state) and then by cooling down the system very slowly. This gives us both the melting transition and the low- temperature configuration. We can then use this low-temperature configuration as the ground state and gradually increase the temperature until the system “melts”. We have found the cooling and the heating curves to agree very well (see Fig. 4.6). This confirms that we are indeed going through a reversible melting transition and also strongly suggests that this transition is a second-order phase transition for all 54 values of 3:, and for the two sets of interaction parameters that we have chosen. This is in contrast to molecular—field theory (MFT) results. This is understandable since the phase transition from a 2 x 2 ordered structure to a disordered one is predicted to be in the same universality class as the four-state Potts model in 2-dimension which is known to have a continuous transition[18]. For binary systems (AI/3 or BV3), our results are consistent with the work of Glosli and Plischke[19] (in their case K2 : 0.1). It should however be pointed out that in some experiments[20], the melting transition of binary alkali-metal GIC’s have been found to be weakly first order. We suggest that this is either due to interlayer coupling not included in the present calculation or the inadequacy of the lattice gas picture. In Figs.4.6(a)-(c) the order parameter IVI, energy E, and specific heat C of the ternary A1/3B2/3V3 are given as a function of the reduced temperature T for J = 0.05 and L = 0.01. These figures strongly suggest that the melting transition is second order in contrast to the MFT results. In Figs. 4.7(a) and 4.7(b), the melting transition temperature Tm versus concentration a: is plotted as a function of a: for J = 0.05 and L = 0.01 and for J = 0.05 and L = 0.08 obtained both from the Monte Carlo simulation and MFT calculation. We find from the simulation data that the vaersus a: curve tends to agree with the Vegard’s law when J << L and for J/L : 0, we find that Tm obeys Vegard’s law. The low-temperature configuration shows that the system tends to phase separate when J << L. Since phase separation has not been observed in alkali-metal GIC’s, we believe that J is larger than L and therefore we suggest that melting transition temperature for the ternary GIC’s such as in K1_,,Rb,Ca (0 < :c < 1) may be higher than that for the corresponding binary systems (K08 and Rsz). In Fig.4.8, a typical low temperature configuration of the ternary A1/3Bz/3I/E; is given for J = 0.05 and L = 0.01. Since we only allow the nearest-neighbor spin exchange, and since the order-disorder phase transition on the 2 x 2 superlattice (see Sec. 4.5.2) occurs at a temperature well below the melting transition temperature, 55 '_'_f T b t.os.—._ (a) g : 4 i a : o ‘F culls. ’ ° Mum x i ‘ O.” . < b . o t . § . Q , P b t . A A A, l 3 D O.” o ‘ 1' I ' (b) i b I 0.0! P . ’ 1 > g 4 U f O 1 . . + 6.1“ l 0.04 I- O o noun. . ’ . 4 t t / 1 4 l a.” 0.19 0 so 1' I E (c) , -t 2.8 I: 9 . : o _ 4- null" p ‘60110' L: E .° 1., .P Q % q - b O - o i J a : “6‘ .3 Q t 8.30 0 1° 0 40 1' Figure 4.6: Thermodynamic properties of Al/ng/3V3 ternary alloy as a function of temperature (T) near the melting transition temperature for K1 = 1.0, K; = 0.3, J = 0.05, and L = 0.01: (a) T dependence of order parameter mK(= M) [defined in Eq. 4.32]; (b) T dependence of energy per site E; (c) T dependence of specific heat per site C. 0.0 0.0 0 0.4 [.- 0.3 0.2 0.1 56 (O-l AAA‘AALA‘ 'VV'VVij'Vj'VTr1VVT"f" AHLLLLLILALALLLLLIA A1 A 'vvvvr'frVI"VV'TV—ffrffvvf (N LAAAIAAAA‘AALA‘ W N, \ r t D D b I D p t p D p p p D I t i D D a D D D D 0.0 0.0 9 a 9 a .0 so \ LLLLLLLALAIILLLIILLAIAALAm 0 0.2 0.4 0.0 0.01.00 0.2 0.4 0.0 0.0 1.0 X X 0.1 Figure 4.7: Concentration dependence of the melting temperature Tm in molecu- lar-field theory and Monte Carlo simulation: (a) for J = 0.05 and L = 0.01; (b) for J = 0.05 and L = 0.08. 57 erAAoeoeAAoAA ‘ 0‘. 0 x as us ‘ ‘ . ‘ J i 0.05 0.0.: e O.A 0.0.0": L- om 0AA. OOAOOAAO :3133; erA A A ooAAoooA Figure 4.8: A typical configuration of A1/3Bg/3V3 ternary alloy first being heated to a temperature above T6 = 0.7 (above the melting transition temperature Tm), then slowly cooled to a temperature T = 0.1 (below Tm) for J = 0.05 and L = 0.01. 58 the energy barrier prevents the system from going into thermal equilibrium within our limited number of MC runs. So far, all the experimental samples available are made well above the temperature of the A and B order-disorder transition. Thus when cooled, the system actually freezes into a kind of “glass” phase in which different ordered “islands” float in the “sea” of disordered phase. Upon cooling, the disordered system transforms to an ordered 2 x 2 superlattice structure at Tm, the correlation function C, has additional Bragg peaks at the wave vectors Gg, while above Tm there are only peaks at the wave vectors corresponding to the reciprocal of the graphite lattice (see Fig.4.9).The appearance of new Bragg peaks for in the ordered phase indicating new superlattice ordering of A and B atoms on the 2 x 2 superlattice is not possible if we allow the spin exchange to take place only between nearest neighbors of the original triangular lattice. We must instead introduce an exchange between A and B atoms directly. This is discussed below. 4.5.2 Order-Disorder Transition on the 2 X 2 Superlattice We now consider the order-disorder phase transition in the solid alloy phase which. takes place within the 2 x 2 superlattice at temperatures well below Tm. In general the possibility and the realization of the ordered phase depends on the concentration cc. Let us look at the ternary A1/3B3/3VE; system. Since for T <‘ Tm almost all the atoms are on the 2 x 2 superlattice site, we can define a new lattice-gas model within this 2 x 2 superlattice with a Hamiltonian given by H = J 2 0,173+ L 2 030,-, (4.35) NN NN where NN indicates the nearest-neighbor interaction on the 2 x 2 superlattice sites. Since the spin variable 0',- can only have two values (:izl) because these lattice sites are occupied by either an A or a B atom, the second term of Hamiltonian (4.35) is a constant if we fix the concentration a: at 51,-. Thus Hamiltonian (4.35) becomes H = J 2 0,17,. (4.36) NN 59 9.0 vavrrvI—VfivTIvfrrlvvrrTrtTYIVY'1r"'VY (a) v r l L 4.0 [ 5(a) 0.0 L 1 L l 1. l L L A I lllLllllllILJLLLLAIlLLJlllllllLlLllllLL -2.0 -1.5 -1.0 -0.5 O 0.5 1.0 1.5 2.0 0/00 9.0 erfrv—rrIvrrvvv'rrvlrrrTrTTTt‘TUVV;']V'a .4 .7 _] J .] 4 d d v U r r I 4.0 5(0) I I r T T l' T T V 0.0 '- h I‘lllllellIllllllj_l'lLlllJmlLLllLlLLLJ—IL -2.0 -1.0 0 ' 1.0 2.0 moo Figure 4.9: Correlation function 0,, in the direction of n: (3?, §) for J = 0.05, L = 0.01, and :1: = g at'temperatures: (a) T = 1.0 (T > Tm); (b) T = 0.4 (T < Tm). 60 This is actually the antiferromagnetic Ising model with the nearest-neighbor in- teraction only which has been studied extensively in the past[20]. Monte Carlo sim- ulation is performed by exchanging the N N spins (on the 2 x 2 superlattice, which correspond to the third nearest-neighbor on the original lattice). The system size is 18 x 18, the interaction parameter considered is J = 0.05, and 5000 MC steps per NN bonds are performed. At low temperature, the correlation function 0,, has commensurate wave-vector peaks at G1 = (27r/3)(1,\/3), G2 = (27r/3)(1,——\/3), and G3 = (47r/3)(1',0), and the corresponding order parameter is given by mK = [15 23: C'c:.-/C'q=0]l/2 (4.37) i=1 We find that the system undergoes a second-order phase transition to an ordered (2\/3 x 2x/3)R30' structure (see Fig. 4.10 and 4.3a). The transition temperature that we find is Tc/ K 1 = 0.06, which is indeed much lower than the melting transition temperature Tm/ K1 = 0.285. At temperatures higher than Tc, we find that the configuration of the system is similar to that shown in Fig. 4.8. This indicates that the long equilibration time needed for the spin exchange (equivalent to the dissimilar atom exchange) may be the reason why long-range order has not yet been seen in experiments. Up to now all the samples are made at high temperature from where the system is cooled. The system possibly freezes into a glass like state so that long- range order cannot be achieved, even when the ground state can in principle have a long-range-ordered structure. 4.6 Summary and Conclusion Using Monte Carlo simulation we have determined the melting curve Tm for a spin-1 Ising (lattice-gas) model appropriate for the stage-1 alkali-metal ternary GIC’s and have compared our results with the MF T results. The physical origin of the :1: depen~ dence of melting transition temperature Tm is discussed within a Landau expansion. 61 0A AO0A0 AO0AO0AO0AO0AO0A00 AAATOMS 0AO0AO0AO0AO0AO0A0 Figure 4.10: Typical configuration of A1/3Bg/3V3 ternary alloy first being heated to the temperature above Tc, then slowly cooling to the temperature lower than Tc. 62 A brief discussion of the order-disorder phase transition on a. 2 x 2 superlattice is given using MC simulation. Substantially more experimental data are needed, however, to determine the physical nature of this melting process. This work has been published in Physical Review B[22]. Bibliography [1] G. S. Parry, Mat. Sci. Eng. 31, 99 (1977). [2] J. E. Fischer, Comments Sol. St. Phys. 8, 153 (1978). [3] A. Weiss, Ange. Chem. 75, 755 (1961). [4] S. A. Solin, Adv. Chem. Phys. 49, 455 (1982). [5] M. S. Dresselhaus and G. Dresselhaus, Adv. Phys. 30, 139 (1981). [6] H. Zabel and P. C. Chow, Comments Cond. Mat. Phys. 12, 225 (1986). [7] H. Zabel and S. A. Solin (ed.), Graphite Intercalation Compounds, Vol. I, Struc- ture and Dynamics, (Springer—Verlag, 1988). [8] Per Bak and E. Domany, Phys. Rev. B 20, 2818 (1979). [9] S. A. Solin, P. Chow, and H. Zabel, Phys. Rev. Lett. 53, 1927 (1984). [10] M. Blume, V. J. Emery, and R. B. Griffiths, Phys. Rev. A 4, 1071(1971). [11] P. Chow and H. Zabel, Synth. Met. 7, 243 (1983). [12] B. R. York, S. K. Hark, and S. A. Solin, Solid State Commun. 54, 475 (1985). [13] G. M. Carneiro and E. V. L. de Mello, Phys. Rev. B 35, 358 (1987). [14] H. H. Lee and D. P. Landau, Phys. Rev. B 20, 2893 (1979). 63 64 [15] K. Kawasaki, in Phase Transition and Critical Phenomena, edited by C. Domb and M. S. Green (Academic, New York, 1972), Vol.2, p443. [16] Y. Saito, Phys. Rev. B 24, 6652 (1981). [17] L. E. Reichl, A Modern Course in Statistical Physics, (University of Texas Press, 1980). [18] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982). [19] James Glosli and Michael Plischke, Can. J. Phys. 61, 1515 (1983). [20] N. C. Bartelt, T. L. Einstein, and L. D. Roelofs, Phys. Rev. B 35, 1776 (1987). [21] Zhi-Xiong Cai and S. D. Mahanti, Phys. Rev. B 36, 6928 (1987). Chapter 5 Oxygen Ordering in YBazCu3O7_5 and Related Systems 5.1 Introduction In 1986, Bednorz and Miiller discovered superconductivity above 30K in the La-Ba- Cu-O system[1, 2]. Several groups quickly reproduced their results and identified the superconducting phase as Lag-,Ba,CuO4. Many researchers then began to de- termine the effects of elemental substitutions and different processing conditions on the structure and superconducting properties of this oxide. Following this approach, in 1987, several groups discovered superconductivity above 90K in the Y-Ba-Cu-O system[3, 4, 5]. As before, the discovery was quickly reproduced, the 90K phase was identified as YBa3Cu30-r, and studies of the effects of substitutions and processing conditions were initiated. This cycle has been repeated twice in 1988 with the inde- pendent discoveries of superconductivity above 100K in the Bi-Ca-Sr-Cu-O system by Maeda et al.[6] and in the Tl-Ca-Ba-Cu-O system by Sheng and Hermann[7, 8]. In this thesis, we are concerned exclusively with the structural phase transitions observed in the YBagCu307_5 and related systems since these systems have been most extensively studied experimentally as well as theoretically. The observed struc- 65 66 tural phase transition can be directly-linked to the ordering of oxygen atoms in one of the three Cu-oxygen planes which is oxygen deficient. In Sec. 5.2, we will describe the basic atomic arrangements in YBa2CuaO-r and how these arrangements affect the superconducting transition temperature. In Sec. 5.3 we will discuss experimental and theoretical studies of the orthorhombic to tetragonal structural phase transition in YBagCu307 systems. We then establish in Sec. 5.4 an isotropic lattice gas model using a simple screened Columb repulsion and a Cu2+ mediated attractive interaction between oxygen atoms. In Sec. 5.5, we use this lattice gas model to study the effect of different thermal cooling conditions on the oxygen ordering in oxygen deficient systems. In sec. 5.6, we modify this lattice gas model to study the effect of substi- tuting trivalent impurities for Cu2+ ions on the structural properties. These studies have both theoretical and practical interests because they provide clues about the effect of different processing procedures on the structural features of the YBagCuaO7 systems. In the last section, we will summarize our results and discuss possible di- rections for future theoretical study on the structural aspects of these high-Tc oxide superconductors. 5.2 Oxygen Ordering and Superconductivity Soon after the discovery of high-temperature superconductivity in the Y-Ba-Cu-O system, the phase responsible for the 90 K transition was found to have a cation stoichiometry of 1Y:2Ba:3Cu. The unit cell dimensions determined by electron and x-ray diffraction identified the structure as being related to a cubic perovskite with one of the cube axes tripled [9, 10, 11, 12, 13, 14, 15, 16, 17, 18] and there are three copper-oxygen planes in one unit cell. The neutron diffraction studies using Rietveld refinement technique showed unambiguously that the oxygen atoms in the oxygen deficient copper plane (basal plane with layer stoichiometry CuO) of the unit cell were ordered in the orthorhombic structure as shown in Fig. 5.1 and the 90 K material contained nearly seven oxygen atoms per unit cell[19, 20, 21, 22, 23, 24, 25, 67 26, 27]. The plane containing Cul in this figure is the oxygen deficient basal plane. In the perfectly ordered situation, this plane consists of Cu-O chains parallel to the crystallographic b-axis. Superconducting properties of this system is closely related with the oxygen con- tent and oxygen ordering. It has been found [28] that in the thermally quenched samples, the superconducting transition temperature changes from 92 to OK as the oxygen content is varied from 7.0 to 6.5. On the other hand, transition tempera- ture measurements on samples obtained by the gettered annealing method at low temperatures[29] show a clear plateau in Tc versus the oxygen concentration curve as shown in upper panel of Fig. 5.2. In this case Tc —» 0 somewhere between 6.4 and 6.3 for the oxygen concentration. In the lower panel the width of the transition which nar- rows in the “plateau” region is shown. This plateau at Tc z 60 K is a strong evidence for the existence of two superconducting phases in the YBagCu3O7_5 system and the tWO transition temperature 92K and 60K have been related to different hole concen- trations in the CuO; plane [30] which in turn depends on the oxygen structure in the basal plane. The plateau is thought to be related to the existence of a double-cell orthorhombic structure, which has been seen in several experiments[31, 32, 33, 34]. Recent experiments where samples have been annealed for a long time with fixed oxygen content shows an increase in superconducting transition temperature with in- crease in the orthorhombicity of the system [35]. It has been shown that the increased orthorhombicity is caused by increased degree of oxygen ordering in the basal plane. These experiments clearly indicate that there is a direct link between oxygen order- ing and oxygen content in the CuO basal plane and the superconducting transition temperature. 5.3 Orthorhambic to Tetragonal Phase Transition The observed structure and properties of YBagCuaO7_5 (we refer to this as 123 com- pound) below room temperature depend critically on how the material is processed 68 Figure 5.1: The ground state of YBa;Cu307. TlKl K/(decode M) 69 100 , , , TRANSITION 3 - - WIDTH 2 .— 1 ' l o 1 1 1 l 7 X In BOZYCIJgO; ‘ Figure 5.2: Effect of oxygen content on Tc. 70 at higher temperatures. This situation arises because the superconducting properties are largely controlled by the oxygen content and ordering in YBazCugO7-5, which in turn, are controlled by the annealing times and temperatures,the oxygen partial pres- sures, and the quench rates used in preparing the material. In this section we review briefly the experimental facts about the orthorhombic to tetragonal phase transition in 123 compound in order to lay the foundation for the theoretical model we will establish in the following sections. The orthorhombic to tetragonal phase transition in 123 compound was initially identified by electron beam heating in a transmission electron microscope (TEM) [36] and was subsequently studied in situ by numerous techniques. Hot-stage X~ray diffraction[36, 37] showed a sudden change of lattice parameters above 500°C and the a and b axes became equal. The orthorhombic 123 compound turns tetragonal at a temperature which depends on the oxygen partial pressure [37, 38]. Thermogravi- metric studies [38, 39, 40, 41, 42] have shown that 07 starting materials begin to lose oxygen reversibly above ~350e400°C. The orthorhombic to tetragonal transition is assumed to occur when there is a discOntinuity in the curve of weight loss versus temperature plots[39]. . In situ neutron diffraction[43] showed that the observed temperature dependent changes both in the lattice parameters are caused by changes in the oxygen content and oxygen order in the basal copper plane as shown in Fig. 5.3. The oxygen that is lost above ~ 400°C comes primarily from the O(4)[=(0,1/2,0)] site and goes to the normally vacant (1/2,0,0) site, resulting in the contraction of the b axis and expansion of the a axis just below the transition. These results show that the oxygen structure changes from fully ordered at low temperature, to partially ordered at higher temperatures in the orthorhombic phase, to completely disordered at still higher temperatures in the tetragonal phase. Jorgensen et al. [43] also find that although the orthorhombic to tetragonal phase transition temperature depends on the oxygen partial pressure, the oxygen concentration at the transition temperature is relatively 71 insensitive to the oxygen partial pressure and the transition always occurs near the oxygen concentration of 6.5. These experimental data clearly indicate that the orthorhombic to tetragonal struc- tural transition in Y832011307_5 is an oxygen order-disorder transition occurring in the CuO basal plane. Thus a two-dimensional lattice gas model is a reasonably good starting point to study this phenomena, as we shall see in the next section. 5.4 The Lattice-gas Model and its Parameters To develop a lattice gas model for the order-disorder transition, the structure of CuO basal plane can be described by a 2-dimensional square lattice consisting of Cu atoms, oxygen atoms, and oxygen vacancies (V). In the ordered (orthorhombic) phase oxygen atoms and vacancies are ordered so that there are linear Cu-O chains parallel to the b-axis and Cu-V chains parallel to the a-axis, as show in Fig. 5.4. The stoichiometry of the arrangement is YBanll307. If the oxygen atoms are absent from this plane, then one has YBagCuaOo and the structure of the compound is tetragonal. To obtain a lattice gas model to describe the oxygen order-disorder transition let us assume that the interaction between oxygen ions in this system is predominately Columb interaction with metallic screening for 6.5 < 7—6 < 7. The metallic screening is due to the holes in the oxygen band. For 5 > 0.5, it is believed that the oxygen holes are localized and do not take part in metallic screening. The lattice gas Hamiltonian used in our study has the form: HLG = V1 2 ngnj + V: (Z): ninj + V}, gain,- + (6(a) — a) Z 72,-, (5.1) ij Cu ij 2' Where the single site energy 6(2)) can be defined as C(23) 2 E5 + Ed/Z + 012:, (5.2) where E5 is the binding energy of the oxygen atom in the basal plane and its origin is the Madelung energy; Ed is the disassociation energy of an oxygen molecule and a: 72 1.0 (a) 100% oxygen 0.0 r 00. ‘ 0.7 ‘ (0.1/2.0) slte >0 2 as g as- 0.5- 2 a 0.4- 3 0.34 .9. g 02- 1.1: 0.1- (1/2,0.0) site 0 0.0 - '9 -o.1 - r . . - - e - o 200 400 800 300 1000 Temperature (‘0) Figure 5.3: Fractional site occupancy of oxygen in the CuO basal plane versus tem- perature for a sample heated in one atmosphere oxygen. (Jorgensen et al.[43]) 73 1| 13u O Oars-n 5 D Oxygen vacancy Figure 5.4: The ground state structure of YBa3Cu307 basal plane. The back dots indicate the positions of Cu ions, the open circles indicate the positions of the oxygen atoms and the open square indicate oxygen vacancy positions. 74 is the average oxygen site occupancy which is given by: 1 :13 2 N27“. (5.3) The physical origin of the as: term in Eqn. 5.2 can be ascribed either to a change in the copper valency with increased oxygen concentration a: [44] or to a mean field representation of the longer range part of the oxygen-oxygen repulsion[45]. The chemical potential p. can be calculated using the equation of state of the diatomic ideal gas[46]: Pl— —h u=ksTlnf ( if“ “11“. (5.4) 0 where P is the oxygen partial pressure, T5 is the characteristic temperature defining the vibrational energy of the oxygen molecule. The characteristic pressure P0 is given W 271'ng 4W21k3 P0 = (—'h,—‘)3/2(T—WWBTT/2 (5-5) where the moment of initiaI of the 02 molecule is given by the relation #1 = 2.1K; M is the mass of the 02 molecule. In Eqn. 5.1, V1 is the interaction between two nearest neighbor oxygen atoms. The interaction between two next nearest neighbor oxygen atoms is either V; or V3 depending on whether a copper atom bridges these two or not, respectively. de Fontaine et.al[47] were the first to realize that V; and V3 could be different and even be of opposite sign. The physical meaning of these parameter will be discussed later. 5.4.1 Earlier Theoretical Work Before we explore this lattice gas model in detail, let us briefly review some of the other works related to the theoretical study of oxygen ordering in this system. The first lattice gas model calculation was performed by Bakker et al.[48] who assumed that the interaction between oxygen atoms in the basal plane of YB82011307-5 can be modeled by an isotropic lattice gas Hamiltonian with only nearest neighbor interactions (V, = 75 V; = 0) which is equivalent to a 2-dimensional antiferromagnetic Ising model. Using quasi-chemical approximation, they found that for P = 1atm., the system showed an order-disorder transition which agreed with the experimental observation. However later Monte Carlo simulation results[45] showed that the oxygen concentration at the critical temperature was always near a: = 0.37 and was relatively insensitive to the oxygen partial pressure. This result is in disagreement with the experimental results of Jorgensen et al.[43], Which shows that the concentration of oxygen at the critical temperature is near a: = 0.25. The value of the critical concentration (2,, = 0.37) for the nearest neighbor re- pulsion model was first discussed by Binder and Landau[49] who argued that it is the percolation threshold of a system consisting of a filled (fl X J2) lattice (the open circles in Fig. 5.4). This gives :1:c ~ 0.25 + 0.25 x 0.59 = 0.37, where 0.59 is the site percolation threshold for a square lattice. Therefore Bakker et al.’s model with isotropic repulsive interaction is not adequate to describe the orthorhombic to tetragonal phase transition of the YBagCu307_5 system. In addition to the problem of) me, this model also gives a much stronger temperature dependence of a: near the critical temperature T' which is not seen experimentally. This latter problem was addressed by Salomons et al.[50] who introduced a concentration dependent binding energy term (E, + as: of Eqn. 5.2) as discussed earlier in this section. The isotropic repulsive model has another drawback, namely it does not favor chain formation for oxygen deficient systems as seen in experiments. The fact that one finds Cu-O chains in oxygen deficient systems strongly suggest that V3 < V3 or that V; be negative. Bell[51] introduced a model with isotropic attractive interaction for next nearest neighbor oxygen atoms. Although Bell’s model gives the temperature dependence of the oxygen fractional site occupancy in good agreement with exper- iments, there is a phase separation at low temperature and Cu-O chain formation is not energetically favored. On the other hand, de Fontaine[47] and Wills et al.’s model[52] which has V; < 0 and V2; > 0 removes this difficulty. They have obtained the 76 temperature-concentration phase diagram using a cluster variational method which agrees very well with our Monte Carlo simulation results[45]. However the cluster variation method does not work well at low temperature. 5.4.2 Monte Carlo Simulation We will now discuss our Monte Carlo simulations of the lattice gas model with Hamil- tonian defined in Eqn. 5.1. These simulations were carried out on a 40 x 40 square lattice as shown in Fig. 5.4 using periodic boundary condition and 10,000 ~ 20,000 Monte Carlo sweeps per site. We have done both heating and cooling runs and have used different shapes of the boundaries to perform the simulation. These studies give results which agree very well with each other and the studies of the size dependence of the critical behavior gives evidence that for the parameters we have chosen, the transition is continuous. _ First, the Monte Carlo simulations were carried out for a number of different sets of interaction parameters. We find that the oxygen concentration (23¢) at the transi- tion temperature T‘ depends on the relative strengths of the next nearest neighbor interaction parameter V3 and V}. For V2», 2 [Va], 2,, ~ 0.37. However for V3 < IVzl, we can have values ranging from 0.25 ~ 0.33 depending on their relative strengths. To fit the general features of the experimental data shown in Fig. 5.4, we have made the following choice for the values of the parameters in Eqn. 5.1: V1 = 0.25eV, V, = —0.2eV, V; = 0.1eV, E(:c) = —1.45 +1.53. In Fig. 5.5 we show the temperature dependence of the order parameter 17 (sub- lattice magnetization in the Ising language) and the susceptibility x. The transition temperature T' is identified with the peak position of the X ~ T curve. Its values are 700°C for P = 1bar, 660°C for P = 0.2bar and 620°C for P = 0.02bar; These values agree very well with the experimental data (700, 670, and 620°C for the three values of pressure respectively). With the above values of the parameters and for the pressure given above, the low-T phase has a: = 0.5 which corresponds to the 77 stoichiometry of YBaQCuaO-r. In Fig. 5.6, we give the oxygen concentration a: and the fractional site occupancies 1:1, 2:2 as a function of T for different values of P. It can be seen that the occupation of the “wrong” sublattice (:03) decreases drastically below T‘. The value of the oxygen concentration at T‘ is close to 0.25 in good agreement with the experimental results. Also the T-dependence of a: for T > T’ is very well reproduced by our Monte Carlo simulation. To see if the model and the parameter values we choose accurately represent the physical properties of the system we study, we show in Fig. 5.7 the pressure dependence of oxygen concentration isotherm compared with the experimental data obtained by Salomons et al. [50]. The simulation results agree very well with the experimental data. To further quantify the oxygen ordering , we have calculated the average and max- imum chain lengths along a and b axes respectively. These are shown in Fig. 5.8. We can see that even at temperature well above T’, the average and the maximum Cu-O chain length still have very high value, indicating strong short range order in the tetragonal phase. 5.4.3 Justification of the Parameter Values Since Monte Carlo simulations of this lattice gas model with the interaction parameter chosen as above give results which are in very good agreement with experimental data, it is important to discuss in some detail the physical meaning and the assumptions made in establishing the values of these parameters. In particular, we would like to discuss why a model with only short range interaction is a good approximation to the system YBagcu307_5 which is apparently dominated by ionic bonding. First, let us consider the system YBazcuaOg. This is an antiferromagnetic insulator With no oxygen atoms in the basal plane. The charge state of Cu(1) is 1+. When Oxygen is intercalated into the system, the number of Cu” ions per unit cell (pCuH') decreases linearly as the oxygen concentration 2: increases, as shown by the Xoray 78 .0-P31bn 1 +F=IOJ bar —O- P8 0.02 bar : 0-5P . q 1 Figure 5.5: Monte Carlo simulation results of the temperature dependence of (a) the order parameter 77; and (b) the susceptibility for different oxygen partial pressures. 79 l 1070 1170 770 970 1’ (In K) Carlo simulation studies of the temperature dependence of (a) the total oxygen concentration; and (b) the fractional site occupancies for different OXYSen partial pressure. lfigure 5.6: Monte 80 absorption measurements of the Cu K -edge by Tranquada et al.[53]. The slope of the p0,,” versus a: curve is between 1 and 2 depending on whether the sample is in the tetragonal or the orthorhombic phase. Thus our assumption of the binding energy of oxygen as a linear function of z: with the slope of 1.5 is quite reasonable for the whole range of .oxygen concentration of YBagCu3Og+h system (0 < a: < 0.5). For the YBagCu307 system, the nominal charge state of the basal plane Cu is 3+. However, due to the strong on site Columb repulsion, d8 state (two holes in the d- shell) is highly improbable theoretically and also has not been seen experimentally. XP S experiments indicate that the extra “hole” goes to the oxygen site and is mobile. Thus when the concentration of oxygen in the basal plane a: Z 0.25, the system goes through an metal-insulator transition and becomes metallic. We suggest that it is the metallic screening by these mobile holes which makes the interaction between oxygen ions in the basal plane short ranged so that the lattice gas model with short range interaction is a good approximation to the real system. If we assume (1) that there is one mobile hole per unit cell and the distribution of the “hole gas” is such that there is no mobile hole near the Y layer and (2) that holes are uniformly distributed between two nearest CuO; layers (see Fig. 5.1), then including both hole-induced Thomas-Fermi and ion dielectric screening (6 z 7), we estimate that VI z 0.3eV and V2, z 0.1eV, which are consistent with the choice of the parameters given in Eqn. 5.1. The interaction between two oxygen atoms bridged by a Cu atom, Vg, is negative primarily due to the indirect interaction caused by the transfer of oxygen holes via Cu. Unfortunately there is no satisfactory method to calculate V; because it depends not only on the oxygen hole density but also on the interaction between the mobile holes with the localized spins at the Cu sites. To estimate the order of magnitude of V2, we consider the CuO chains only and construct an one dimensional tight binding model to describe the motion of the oxygen hole. We further assume that the localized spins on the Cu sites are either parallel or antiparallel to each other. Using the values of the parameters derived from the Emery’s model[54], we find that V; is negative and 81 its magnitude is of the order of 0.leV. Thus we believe that this lattice gas model is physically reasonable and the oxygen-oxygen interaction is basically anisotropic and of short range nature in the “metallic” phase. It is known that the Thomas-Fermi screening length increases with the decreasing of oxygen content which decreases the hole density, and the lattice gas model with only short range interaction fails to be a. good approximation for a: < 0.25 when the system becomes insulating and the interaction becomes long ranged. Although we have not investigated the phase diagram in the The T ~ :1: plane, such a phase diagram has been studied by various groups using cluster variation method, transfer matrix method, and by Monte Carlo simulation[55]. For similar sets of parameters, these calculations show that in addition to the orthorhombic phase of YBagCu3O7 (OI) and the tetragonal phase of YBagCuaOg, there is a thermody- namically stable phase of YBEzCIl30gj with alternate full and empty Cu-O chains as shown in Fig. 5.9. The nature of the thermodynamic transition between the tetrago- nal and the above two orthorhombic phases is found to be second order. In the next section, we will discuss the nature of the oxygen ordering under different processing conditions. 5.5 Effect of Thermal Quenching on Oxygen Ordering Soon after the orthorhombic YB83C11307 and tetragonal YBagCuaoe phases were identified, researchers began to investigate the structures and properties of YBAzCu3O7_5 samples with intermediate oxygen contents (0 < 6 < 1). In a large number of such studies [56, 15, 57, 57, 58, 59, 60], the samples were prepared by rapidly quench- ing the high-temperature, low oxygen content tetragonal phase from 600—1000°C to room temperature or liquid nitrogen temperature. These studies usually found a Smooth variation of T6 with oxygen content with the samples becoming insulators near 035, corresponding to an orthorhombic-to-tetragonal phase transition[56]. Al- ternatively, intermediate oxygen content samples were prepared by removing oxygen 82 at lower temperatures (usually 400—500°C) from 07 starting material. These studies employed a variety of methods to remove the oxygen, including annealing in reducing atmospheres[60, 61], equilibrating 06 and O7 123 mixtures in sealed quartz tubes[29], gettering annealing with zirconium in a sealed tube[5], and titrating electrochemi- cally in a solid-state cell[62]. The variation of T5 with oxygen content in samples prepared at lower temperatures is different from that observed in samples quenched from the tetragonal phase. The superconducting transition temperature does not vary smoothly with oxygen content in these samples. Instead, Tc remains at 91K between 07 and 0”, then falls to a ~ 60K ”plateau” between 05.7 and OM, and finally drops to zero (i.e. not superconducting) between 05.5 and O“. In their attempt to explain the plateau in the oxygen concentration dependence of Tc, Zaanen et al.[63] devel- oped a simple model by assuming that the superconducting transition temperature is directly linked to the hole concentration in the Cunglane. The hole concentration in turn depends on the oxygen concentration and oxygen ordering in the CuO plane through the electronic structure of the CuO chains formed in the CuO basal plane. To apply this model to the real system, however, it is important to know the CuO chain length distribution under various procession methods. In order to see how the CuO chain length distribution depends on the oxygen con- centration and quench rates, we have carried out a series of Monte Carlo simulations for constant oxygen concentration (Kawasaki dynamics). We will explore the low temperature structure of the CuO” (0 < y < 1) plane for different oxygen concentra- tions and thermal cooling rates. Structure factors associated with oxygen ordering along the a and b axes are calculated to monitor the growth of the orthorhombic domains. A 40 x 40 lattice was used to calculate the structure factor of oxygen con- figuration. The high temperature configuration was obtained by using Monte Carlo simulation at 2300K (well above the order-disorder phase transition temperature for O7) and the system was equilibrated at this temperature with 10,000MCS / site (MCS = Monte Carlo steps). The low temperature data were obtained by either slowly 83 cooling (annealing) the system with a temperature interval of 100K or rapid cool- ing (quenching) the system to 300K. From the oxygen configuration generated from Monte Carlo simulation, the oxygen structure factor was calculated through equation iniea: i -r 2 Sq:f(q) where N = Z, n,- and f(q) is the atomic scattering factor of oxygen. 10,000MCS/site (5.6) were used to obtain the thermodynamic averages for each temperature. Figs. 5.10 and 5.11 show the oxygen configurations of the annealed and the quenched systems at 300K for the YBagCuaO“ system. We find that the struc- ture for the annealed samples is a double-cell orthorhombic (OII) structure with half of the CuO chains intact and half empty. The quenched system, on the other hand, shows very small domain of orthorhombic structure, with CuO chains oriented along both a and b directions. The overall symmetry is, however, tetragonal as indicated by the calculated structure factor along the a and b directions. The difference between the structure of the quenched and the annealed systems has also been studied for the oxygen concentration 06,7. Since this oxygen con- centration is not appropriate for either the OI structure or the 011 structure, new structural features arise. Typical oxygen configuration and the structure factor for the YBflzCll30gJ compound at 300K are shown in Fig. 5.12 and 5.13 and Fig. 5.14 and 5.15, respectively. From these simulations we find that the low-temperature structure of the annealed system (Fig. 5.12) is a highly degenerate one. Instead of a single domain of ordered structure or a two-domain structure which is phase sep- arated into OI and OII domains, we find long Cu-O chains “intercalated” into each other so that one sees a mixture of structure with OI and OII symmetry as indicated by the calculated structure factor (Fig. 5.14) as well as the “snapshots” of the oxygen configurations (Fig. 5.12). From the width of the diffuse peak at (0, 11', 0) along 0. axis [Fig. 5.14(b)] we estimate the typical OII domain size along the 0. axis to be about 5 - 6 unit cells corresponding to about 20-25A, which agrees rather well with the experimental results[33, 34]. It is also shown, both by the oxygen structure factor 84 calculation and the “snapshots” of the oxygen configuration, that the average Cu-O chain length depends sensitively on the thermal history of the sample. This may in turn, as indicated by Zaanen et al.[63], affect the hole density in the CuO; plane, and the superconducting transition temperature Tc. The quenched YBaQCU30gJ system shows short chains in both a and b directions (see Fig. 5.13) and the overall symmetry is tetragonal, as seen in the structure factor for both a: and y (i.e., b and a ) directions (Fig. 5.15). 5.6 Effect of Substituting Cu(1) by a Trivalent Impurity on Oxygen Ordering In this section the results of Monte Carlo simulations to understand the effects of substituting basal plane Cu by trivalent impurities on the oxygen ordering and struc- ture are discussed[64]. The motivation for both experimental and theoretical studies for these substittuted systems has been to elucidate the role of Cu in the CuO chains and CuO; planes in the observed supercOnducting properties. Most of the work has focused on elements in the series from Fe to Ga, although some work has been done on the Al doped systems. It is generally agreed that the trivalent ions Fe, Co, Al and Ga preferentially substitute into the Cu(1) sites, the so called chain sites[65, 66, 67, 68]. In this thesis, we concentrate only on the nonmagnetic impurities such as Ga and Al. X-ray and neutron diffraction studies indicate that the structure of YBagCu3_¢(Ga,Al),O7 becomes tetragonal for a: > 0.18. Unlike the tetragonal phase of YBagCu307_5 (6 > 0.5), the tetragonal YBagCu3-.(Ga,Al),O-; system is a superconductor with Tc as high as 80K[65], as shown in the Fig. 5.16. The Cu(1) site can also be preferen- tially substituted by magnetic elements such as Fe and Co [66, 67, 68]. The structure of these systems becomes tetragonal at a much smaller value of :1: compared to the Ga or Al doped systems. These impurities introduce additional complexity into the Systems due to their magnetic moments. We will briefly discuss how to model these 85 systems to study the oxygen ordering problem by developing a generalized lattice gas model with the assumption of metallic screening of oxygen interaction similar to the case of “pure” system as we discussed in the last section. When a divalent or monovalent Cu ion is substituted by a trivalent M (M =Ga or Al) ion, holes associated with the oxygen sites will tend to localize at the M sites shown in Fig. 5.17. This will have two major effects which are (a) the binding energy of the oxygen ion is increased near the impurity site and (b) the oxygen-oxygen repulsion is increased near the impurity site either due to the increase of the negative charge of oxygen ions or equivalently due to the decrease in local screening because of the decrease in the hole density. Based on this physical picture, we modify the lattice-gas model of Eqn. 5.1 to incorporate the effects of the impurities. The Hamiltonian for the modified lattice gas model is now given by H = Z Vujngnj + V: Z nin,‘ + Z ngjngnj + 2 E5724, (5.7) (‘1') (ii)Cu (ii) i where V; : —0.2eV, ~ (5.8) f/lij = 0.256V, V3,} = 0.18V, (5.9) if neither site i nor site j are the nearest neighbor of M; V2 = —0.2eV, (5.10) Vh'j 2' 0.56V, Vggj = 0.28V, (5.11) if either site i or site j is the nearest neighbor of It/I; V2 = —0.2eV, (5-12) Vuj = 1.06V, I/3ij = 0.46V, (5.13) if both sites i and j are the nearest neighbors of M. In addition, the single site binding energy E, = -4eV if the site i is the nearest neighbor of M and E; = 0 otherwise. 86 The above choice of parameters is a simple interpolative approach to take into account the inhomogenuous charge distribution and screening in these substituted systems. Since the oxygen concentration in YBagCu3_,(Ga,Al),O7_5 is very close to 7, we assume that the oxygen concentration is independent of a: and remains to be 7. We also assume that the impurity M is randomly distributed in the CuO plane. To study the effect of impurities on oxygen ordering, we have performed Monte Carlo simulation with fixed oxygen concentration and for several values of :c on a 60 x 60 square lattice with periodic boundary conditions. For a given configuration of impurity M, 10,000 MCS / site are used to calculate the order parameter ‘1’, which is given by the difference in the oxygen fractional site occupancy 2:1 and 1:; defined in Sec. 5.4 (for the YB&2011307-5 system this order parameter was defined as 17): 231—232 $1+132 (I): (5.14) For each value of 2:, (I’ is averaged over 20 configurations of the impurities. The calculated order parameter Q at T = 300K is shown as squares in Fig. 5.18. The order L-s. b+a parameter deduced from the experiment of Xiao et.al[64] [assuming that (I) = is also shown in Fig. 5.18 as crosses. The general features agree rather well, thus providing a justification for the basic assumptions behind our modified lattice gas model. The discrepancy between theory and experiment can be either due to the long- range effect of the impurity on local oxygen ordering or the fact that the distribution of impurity atoms is not completely random. From a “snapshot” of the oxygen configuration of the Ill-substituted system we find that local “cross-links” are formed around the M sites at the low impurity limit (see Fig. 5.19). These cross-links destroy the Cu-O chains along the b axis and form short Cu-O chains along the a. axis. As the concentration of the impurity increases, the number of the short Cu-O chains along the a axis also increases. At high impurity density, there are long Cu-O chains in both directions thus overall the system shows tetragonal symmetry (see Fig. 5.20). However, unlike the tetragonal phase of the quenched YBagCuaO-r-5 systems, the length of the Cu-O chains are rather long (due 87 to large oxygen concentration). We hypothize that the chain lengths maybe sufficiently large to produce sufficiently large density of conducting holes in the Cqu plane. This may be the undelying reason why the superconducting transition temperature T, is still quite high in these Ga or Al substituted compounds even though the structure is tetragonal. Careful electronic structure calculations following the precedures of Zannen et. al should be carried out to test this hypothesis. 5.7 Summary and Conclusion In this chapter we have established a simple lattice gas model with anisotropic interac- tion to explore the structure of YBagCu307_5 and related systems using the argument of metallic screening Columb interaction. It is shown, via Monte Carlo simulation, that this model represents very well the structural aspect of the YB32011307_5 sys- tem. In particular, it describes the effect of temperature, oxygen partial pressure, and cooling rate on the oxygen ordering in the CuO basal plane very well. The parame- ters of this model have been obtained recently via a first-principles linear muffin-tin orbital (LMTO) calculation by Sterne and Wille[68]. Their results agrees very well with the parameter we derived in this chapter. A major finding from our simulation studies is that the oxygen oi'dering is very sen- sitive to the thermal history of the samples. The cooling rate should be slow enough to get long Cu-O chains. The low-temperature structure of the oxygen deficient sys- tem obtained by annealing is shown to be a mixture of thermodynamically stable OI and 011 structures. The effect of substituting Cu by trivalent impurities on the oxygen ordering has been investigated using a modified lattice gas model established under the assump- tion of screened Columb interaction. The orthorhombic to tetragonal transition in Ca and Al doped system is found to be due to the formation of short cross-linked CuO chains near these trivalent impurities. Recent transmission electron diffraction 88 measurements[69] shows clearly the cross-link (“tweed”) structure as shown in Fig. 5.20. It is important, however, to realize that the lattice gas model we have established is only a very crude model of the effective interaction in YBagCu3O7_5 system. In this model we have completely ignored the long-range interaction due to the lattice distortion and have treated the long range part of the Columb interaction in a mean field fashion which clearly is not applicable for the case of 6 > 0.5. Khachaturyan et al. [70, 71, 72] have argued that, to correctly describe the ordering process, the full three-dimensional nature of the materials must be used. Although we have shown that this lattice gas model gives a good description of the oxygen ordering under various processing conditions, it obviously cannot explain the stability of the twinn boundary and some of the superlattice features found in oxygen deficient compounds. More realistic Hamiltonian is needed to describe these phenomena. -3 10 c 89 " —T=030K ‘o-‘L —-—T3MK . . ----T= 1001K ’ d t ’ d i e --1'=eeex ,4 -.-1'assex ! Figure 5.7: Pressure dependence of the oxygen concentration at different temperatures compared to the experimental data of Salomons et al.[50]. The dots in the figure are experimental data and different lines are Monte Carlo simulation results. (a) 1 ‘0—Patbar $980.20" , ‘Q— P a 0.02 bar-l 1 l 1 L 970 T (In K) Figure 5.8: The temperature dependence of (a) average Cu-O “chain” length; and (1)) maximum Cu-O “chain” length for different values of oxygen partial pressure. 91 0.. : f F v I l r r . l r v I r r I I ' 2 mos 0.0 r 4 - Tetra-0 : E 0.4 '- 1 _ _ - I“ s -'- . .. ; : .- °*’ .- Orthe-l/I . .3. .}—.——. ._.__: rose 0 0.2 b—o—d ._..._.,‘.. i *1 i—H 1—-o—-1‘1 E b—--o see 9 0 1 ih—I—q i—‘J '- ' : ones-V4 ”"1 I 1 0.0 A A l L A L L l A A A A l A A A L L A L A A o 0 0.1 0.2 0.3 0.6 0.5 Basal Plane Oxygen Concentration. 2 Figure 5.9: The phase diagram in the 2:, T plane for interactions V1 = 1, V; = 0.5, and V3 = -0.5 [55]. Three second order lines, separating one disordered phase (Tetrago- nal) and two ordered phases (CI and 011), meet at a multicritical point. Also shown are experimental points due to Specht et al.[73] (filled circles). Monte Carlo results are presented for simulations using Glauber dynamics (x) and Kawasaki dynamics (0). The location (it) and half width at half maximum for the non-scaling specific heat peak is shown, giving the location of the disorder lines. The solid line is the transfer-matrix results (see ref. [55] for details). ()1) uni-undue; Palm-1103 92 Figure 5.10: Oxygen configuration of the annealed YBagCu305_5 ....... -'-v'—'— . __-,..__._...-.__——...eene _ 0...»... .. W m — ‘e 0’... .l.. M..— . e o e I -.va—-— a 94 # ...... e e on eeeeeee* ccccccccc ate 0 ee- 00 e _— oooooooo e ...... eeeeeeee— — e at e. e e 0 save.“ ——-—.eeeee e e we. eeee—_ ssssssssssss e —eeee e eeeeeeefleee — .......... e eee eee—eeeeeee ————eeee ee eeeee— e ...... c—eeeeeeeeee—eee e eeeeeeee efileeOOOQQO—OOOOIIOC‘ e eeeeeeeeeeee ee—eeeeeeeeeeeeeeveeeeeee e ssssss eeeIeeeeetoot—OOOOeeeOIOOO'O'OO uuuuuu -_e see———— sees—ease no 0 ensee-_——- e ooooooooooo e to. eeeee e e—n——_eeeeeee as h— e e eeee‘eee—_eeeeo eeee v—_ —_e ...... oeeee e e as .o—-e-—_- -__. eeeeee e—eI e eee e— eeeeeeee ——-——eee e Moe-e “1.0—“e eeeeeeee see ee -* -- 0e... ecu-o Figure 5.12: Oxygen configuration of the annealed YBa;Cu30u 95 aloe. e e Eeeellelele elf-Ila I e h 111111111 lit—Ill]. 05 e v e 0 0000001! 1111111 (It 5 e tttttttttttttttt if, I. O i It 1 [Oil III. e O... ee.1 11 l i e eeee e O I 000000 eeeeeeeeeeee e I 0 e e I. l l e eeeell e e 1’0 0 i e life ee 1111111111111111111 114 E 0 i ll. eeeeeeeeee f _ e esi1ll1111 111111111 1i 1 e Eeeee e ee ef‘ eeeUeee e e eel1l 11.1. ee-ee l e l1 1] .e .010, a 0 100m 0_m 000 0 m e000, 0 I e v 11111 1i 0 e ee - 1111111 l _ ee -ee 00000000 0 i yum ll 0 0 i t [11 O 0 i 0 tttttttttttt ill 0 e eeeee of! ee 0 e 0 11111114— 1] ttttttttttttt 0 0 I 000000 0 0 0 fl 11111111111 1114 111111111 e ”is”; e e O. efee If. eel-e 111111 1 e 111111 l 111 oe e e l eeeeifi eel 0 CODE 0 E e 1 e. _ 0 00000 0 00000” 0 E 0 e a.” e e 1111111111111 O Oi 0 — O 1 11 [[4 ”FE efie [1 e e e ems 1l11 1 [1 e.eeee e e e ell-Ill OOOOOOOOOOOOOOOO I Ell I10 1 l f! Figure 5.13: Oxygen configuration of the quenched YBagCuaO” 96 e.ee - ,f‘—b.r,!~.. L h e.4e- ' 1 020- ' ..._- -J --J 101 Figure 5.14: Oxygen structure factor of the annealed YBagCugou 97 0.0a 111.r-1-fi,e- 0.0a-...,.. a 1 t b t b p 4 b D ‘ b ...e.eie]- '1 1.0.010- » t a b 1 . p 4 t b ‘ . 11.1-L Mm 00 s 10 °o" KX Figure 5.15: Oxygen structure factor of the quenched YBagCuaou 3.92 3.9 3.88 3.88 3.84 3.82 Lattice parameters (A) Tc (K) hi I W I I l 7M7 T I I l l' r I V l [T r I - c/S - _ _ _ _ - b . '- _ I 1'- - - l - fl. - q - d - I F" — P 1 I- ‘ u p < L— ( _ h C l A l I l L I L I L l I I L l I I l L l I. I ‘ I I' U l r I V r l I I Y I I I I I I I - q I G h 1 t- - l l LLLLL L L L L L ALL I L L L L LL o _ 5 10 15 20 Gallium content (at. 7;) Figure 5.16: Variation of T, of the YBagCu3_,Ga,OT 99 II P! 0 Cu 0 Oxygen E] Oxygen vacancy Figure 5.17: Structure of YBagCu3_,Ga.O-r basal plane with single impurity M. 100 1.0‘ . r . , . _+- Xiao'l data, + node! 0.5 -,' . o a L in .7. 1 i 1 L 1 > 0 0.10 0.20 Ga concentration U. EVIILCHUI .IIUIL Figure 5.18: Oxygen order parameters ¢I> of the YBagCu3_,_.M,07 system as a function of the impurity concentration a: at 300K from Monte Carlo simulations and from data of Xiao et al.(Ref.[64]). 101 o.o....o....o.c......-...........-.-. L. .. . o .o-. o- o.-.- . .. . ...... o.. h‘. I-o 0-o-I' ....-. . .... o . c.....- - ....... -.. . ..-.-.- -o......-..4 0 0 0 ........ ..... . . ...... ..... .o-........ > a 0 sec .. . ........ o o . OII‘DII-.¢ - . . O. t ~0t0-I0' e ' u .-.>¢-o‘o~0 ' .0.....-.-.-.OO . .I.C.v...l.'... ...n-ou-....u..nun....-.u.......o.o-...-‘... ouc-oooonnooouonuuo-o-o.0---c-o-coooaooooooocoo-o .I...C.I.I.I.....D...I.O.’.....I.I...I.I...I.l.‘. ell.I.IIODOIIIO.IIv.-.0.IOIOIOII-QUOIQIOIIOOOIC0v. O .'.'.'...‘.I.....ODO I OOIUOOIOIIO.OOIIIODCOCe‘I- . 00.0.c00coon-IOUIOOI’QDOOOIOIIIoaeochIInOIO-o-o ... . . .uc...-o-o.....o D I OIOOIOOIOOIIOI’ on 00.....og-u-o-o .Cl08......-C|...'.."O'.‘OIO..'I C'....".'I‘.II coeca'o-guoo.oo-o-o..n'.-.......-...........-.-.-. I I 00.00.00.00'00 ' II-OIIC O OOIIUIIOOOIOIOOOOUOGII . - -.....-...-.. ...u...........u....-. ...-....... O O C u-o-euc-o-o-n... o c o - no.0-nuo-o-o .‘I-ottll o . . -¢-a»o- . .. . . . -... e .-........o-.-o~- o o ......--,o............ ....»o......-.-o-....-.. .e-o-n- . ............ e . ... ............... o - ...o--....-u.......... .o . . a a”... . ...- o . ......... ....... .. . . . . .....-.‘.....‘ o . --.o...-o- . -...... . . ... .. .... . ... o s o . ...- - o-..-. .......... . ....... .-. ...... . . . o o .o............ ... .......... ..... . ...... o . . . - .uu-o... . .. . . o . .. .... . . . . . .........‘... . . . . .0 . ..... .-..- ..... . . ...... . ..... . . Figure 5.19: Oxygen configuration of the YBa;Cu3_,M,O-r system for x=0.02 ob- tained from Monte Carlo simulation at 300K. 102 Figure 5.20: Oxygen configuration of the YBa3Cu3-.M,O-, system for x=0.12 ob- tained from Monte Carlo simulation at 300K. Bibliography [1] J. G. Bednorz and K. A. Miiller, Z. Phys. B 64, 189 (1986). [2] J. G. Bednorz, M. Takashjge, and K. A. Miiller, Euophys. Lett. 3, 379 (1987). [3] W. K. Wu, J. R. Ashburn, C. J. Torng, P. H. Hor, R. L. Meng, L. Gao, Z. J. Huang, Y. Q. Wang, and C. W. Chu, Phys. Rev. Lett. 58, 908 (1987). [4] Z. Zhao, L. Chen, Q. Yang, Y. Huang, G. chen, R. Tang, C. Liu, C. Cui, L. Chen, L. Wang, S. Guo, S. Li, and J. Bi, Kexue Tongbao, 32, 1098 (1987). [5] R. J. Cava, B. Batlogg, R. B. van Dover, D. W. Murphy, S. Sunshine, T. Siegrist, J. P. Remeika, E. A. Rietman, S. Zhurak, and G. P. Espinosa, Phys. Rev. Lett. 58, 1676 (1987). [6] M. Maeda, Y. Tanaka, M. Fukutomi, and Asano, Jpn. J. Appl. Phys. 27, L209 (1988). [7] Zhengzhi Sheng and A. M. Hermann, Nature 332, 55 (1988). [8] Zhengzhi Sheng and A. M. Hermann, Nature 332, 138 (1988). [9} W. R. McKinnon, J. M. Tarascon, L. H. Greene, G. W. Hull, and D. A. Hwang, Phys. Rev. B. 35, 7245 (1987). [10] W. J. Gallagher, R. L. Sanstrom, T. R. Dinger, T. M. Shaw, and D. A. Chance, Solid State Commun. 63, 147 (1987). 103 104 [11] D. G. Hinks, L. Soderholm, D. W. Capone II, J. D. Jorgensen, and Ivan K. Schuller, Apll. Phys. Lett. 50, 1688 (1987). [12] T. Hatano, A. Matsushita, K. Nakamura, K. Honda, T. Matsumoto, and K. Ogawa, Jpn. J. Appl. Phys. 26, L374 (1987). [13] Y. Kitano, K. Kifune, I. Mukouda, H. Kamimura, J. Sakurai, Y. Komura, K. Hoshino, M. Suzuki, A. Minami,Y. Maeno, M. Kato, and T. Fujita, Jpn. J. Appl. Phys. 26, L476 (1987). [14] M. Hirabayashi, H. Ihara, N. Tereda, K. Senzaki, K. Hayashi, S. Waki, K. Mu- rata, M. Tomumoto, and Y. Kimura, Jpn. J. Appl. Phys. 26, L498 (1987). [15] E. Takayama-Muromachi, Y. Uchida, Y. Matsui, and K. Kato, Jpn. J. Appl. Phys. 26, L476 (1987). l [16] Y. Syono, M. Kikuchi, K. Ohishi, K. Hiraga, H. Arai, Y. Matsui, N. Kobayashi, T. Sasaoka, and Y. Muto, Jpn. J. Appl. Phys. 26, L498 (1987). [17] K. Semba, S. Tsurumi, M. Hikita, T. Iwata, J. Noda, and S. Kurihara, Jpn. J. Appl. Phys. 26, L429 (1987). [18] S. B. Qadri, L. E. Toth, M. Osofsky, 5. Lawrence, D. U. Gubser, and S. A. Wolf, Phys. Rev. B. 35, 7235 (1987). [19] J. J. Capponi, C. Chaillout, A. W. Hewat, P. Lejay, M. Marezio, N. Nguyen, B. Raveau, J. L. Soubeyrux, J. L . Tholence, and R. Tournier, Europhys. Lett. 3, 1301 (1987). [20] M. A. Beno, L. Soderholm, D. W. Capone, D. G. Hinks, J. D. Jorgensen, I. K. Schuller, C. U. Segre, K. Zhang, and J. D. Grace, Appl. Phys. Lett. 51, 57 (1987). l2 1] F. Beech, S. Miraglia, A. Santoro, and R. S. Roth, Phys. Rev. B 35,8778 (1987). 105 [22] J. E. Greedan, A. O’Reilly, and C. V. Stager, Phys. Rev. B 35, 8770 (1987). [23] F. Izumi, H. Asano, T. Ishigaki, E. Takayama—muromachi, Y. Uchida, N. Watan- abe, and T. Nishikawa, Jpn. J. Appl. Phys. 26, L649 (1987). [24] M. Francois, E. Walker, J. L. Jorda, and Y. Yvon, Solid State Commun. 63, 1149 (1987). [25] T. Kajitani, K. Oh-Ishi, M. Kikuchi, Y. Syono, and M. Hirabayashi, Jpn. J. Appl. Phys. 26, L1144 (1987). [26] Q. W. Yan, P. L. Zhang, Z. G. Shen, J. K. Zhao, Y. Ren, Y. N. Wei, T. D. Mao, C. X. Liu, T. S. Ning, K. Sun, and Q. S. Yang, Phys. Rev. B 36, 5599 (1987). [27] W. I. F. David, W. T. A. Harrison, J. M. Gunn, O. Moze, A. K. Soper, P. Day, J. D. Jorgensen, D. Hinks, M. A. Beno, L. Soderholm, D. W. Capone II, I. K. Schuller, C. U. Segres, K. Zhang, and J. D. Grace, Nature 327, 310 (1987). [28] W. K. Kwok, G. W. Crabtre’e, A. Umezawa, B. W. Veal, J. D. Jorgensen,, S. K. Malik, L. J. Nowiciki, A. P.‘Pa.ulikas, and L. Nunez, Phys. Rev. B 37,106 (1988). [29] R. J. Cava, B. Batlogg, C. H. Chen, E. A. Rietman, S. M. Zahurak, and D. Werder, Phys. Rev. B 36, 5719 (1987). [30] R. J. Cava et al., Phys. Rev. B, to be published. [31] C. Chaillout, M. A. Alario-Franco, J. J. Capponi, J. Chenavas, P. Strobel, and M. Marezio, Solid State Commun. 65, 283 (1988). [32] Y. Kubo, T. Ichihashi, T. Manako, K. Baba, J. Tabuch, and H. Igarashi, Phys. Rev. B 37, 7858 (1988). [33] R. M. Fleming, L. F. Schneemeyer, P. K. Gallagher, B. Batlogg, L. W. Rupp, and J. V. Wazczak, Phys. Rev. B 37, 7921 (1988). 106 [34] C. H. Chen, D. J. Werder, L. F. Schneemeyer, P. K. Gallagher, and J. V. Wazczak, Phs. Rev. B 38, 2888 (1988). [35] Veal et al., Phys. Rev. B, to be published. [36] R. Beyers, G. Lim, E. M. Engler, R. Savoy, T. M. Shaw, T. R. Dinger, W. J. Gallagher, and R. L. Sandstrom, Appl. Phys. Lett. 50, 1918 (1987). [37] I. K. Schuller, D. G. Hinks, M. A. Beno, D. W. Capone II, L. Soderholm, J. P. Locquet, Y. Ruynserade, C. U. Segre, and K. Zhang, Solid State Commun. 63, 385 (1987). [38] R. Beyers et al., Appl. Phys. Lett. 51, 614 (1987). [39] P. K. Gallagher, H. M. O’Bryan, S. A. Sunshine, and D. W. Murphy, Mat. Res. Bull. 22,995 (1987). [40] K. Kishio, J. Shimoyama, T. Hasegawa, K. Kitazawa, and K. Fueki, Jpn. J. Appl. Phys. 26, L1228 (1987). [41] S. Yamaguchi, K. Terabe, A. Saito, S. Yahagi, and Iguchi, Jpn. J. Appl. Phys. 27, L179 (1988). [42] E. D. Specht, C. J. Sparks, A. G. Dhere, J. Brynestad, O. B. Cavin, D. M. Kroeger, and H. A. Oye, Phys. Rev. B 36, 5723 (1987). [43] J. D. Jorgensen, M. A. Beno, D. G. Hinks, L. Soderholm, K. J. Volin, R. L. Hitterman, J. D. Grace, I. K. Schuller, C. U. Segre, K. Zhang, and M. S. Kleefisch, Phys. Rev. B 36, 3608 (1987). [44] T. Iwazumi, I. Nakai, M. Izumi, H. Oyanagi, H. Sawada, H. Ikeda, Y. Saito, Y. Abe, K. Takita and R. Yoshizaki, Solid State Commun. 65, 213 (1988). [45] Zhi-Xiong Cai and S. D. Mahanti, Solid State Commun. 67, 287 (1988). 107 [46] L. D. Landau and E. M. Lifshitz, Course of theoretical Physics, Vol. 5. (Pergamon Press, 1980). [47] D. de Fontaine, L. T. Wille and S. C. Moss, Phys. Rev. B36, 5709 (1987). [48] H. Bakker, D. O. Welch and O. W. Lazarethe, Solid State Commun. 64, 237 (1987). [49] K. Binder, D. P. Landau, Phys. Rev. B 21, 1941 (1980). [50] E. Salomons, N. Koeman, R. Brouwer, D. G. de Groot, and R. Griessen, Solid State Commun. 64, 1141 (1987). [51] J. M. Bell, Phys. Rev. B 37, 541 (1988). [52] L. T. Wille, A. Berera, and D. de Fontaine, Phys. Rev. Lett. 60, 1065 (1988). [53] J. M. Tranquada, S. M. Heald, A. R. Moodenbaugh, and Youwen Xu, Phys. Rev. B 38, 8893 (1988). [54] V. J. Emery, Phys. Rev. Lett. 58, 2794 (1987). [55] T. Aukrust, M. A. Novotny, P. A. Rikvold and D. P. Landau, Phys. Rev. B (to be published). [56] J. D. Jorgensen, B. W. Veal, W. K. Kwok, G. W. Grabtree, A. Umezawa, L. J. Nowicki, and A. P. Paulikas, Phys. Rev. B 36, 5731 (1987). [57] H. Nozaki, I. Ishizawa, O. Fukunaga, and H. Wada, Jpn. J. Appl. Phys. 26, L1180 (1987). [58] M. Tokumoto, H. Ihara, T. Matsubara, M. Hirabayashi, N. Terada, H. Oyanagi, K. Murata, and Y. Kimura, Jpn. J. Appl. Phys. 26, L1565 (1987). [59] S. Nakanishi, M. Kogachi, H. Sasakura, N. Fukuoka, S. Minamigawa, K. Nakahi- gashi, and A. Yanase, Jpn. J. Appl. Phys. 27, L329 (1988). 108 [60] W. E. Farneth, R. K. Borida, E. M. McCarron III, M. K. Crawford, and R. B. Flippen, Solid State Commun. 66, 953 (1988). [61] P. Monod, M. Ribault, F. D’Yvoire, J. Jegoudez, G. Collin, and A. Revcolevschi, J. de Physique 48, 1369 (1987). [62] R. Beyers, G. Gorman, P. M. Grant, V. Y. Lee, R. M. Macfarlane, S. S. P. Parkin, S. JlaPlaca, B. T Ahn, T. M. Gur, and R. A. Huggins (unpubliShed). [63] J. Zaanen, A. J. Paxton, O. J. Jepsen, and O. K. Anderson, Phys. Rev. Lett. 60, 2865 (1988). [64] Zhi-Xiong Cai and S. D. Mahanti, Phys. Rev. B 40, 6558 (1989). [65] Gang Xiao, M. Z. Cieplak, D. Musser, A Gravrin, F. H. Streitz, and C. L. Chien, Phys. Rev. Lett.60, 1446 (1988). [66] C. L. Chien, Gang Xiao, M. Z. Cieplak, D. Musser, J. J. Ryne, and J. Golaas (unpublished. [67] Yasukage Oda, Hiroshi Fujita, Hauruhisa Toyoda, Tetsuyuki Kaneko, Takao Ko- hara, Ichiroch Nakada, and Kum'suke Asayama, Jpn. J. Appl. Phys. 26, L1660 (1987). ‘ [68] P. F. Micheli, J. M. Tarascon, L. H. Greene, P. Barboux, F. J. Rotella, and J. D. Jorgensen, Phys. Rev. B 37, 5932 (1988). [69] P. A. Sterne and L. T. Wille, Physica C162-164, 223 (1989). [70] Y. M. Zhu et.al., (unpublished). [71] A. G. Khachaturyan and J. W. Morris, Jr., Phys. Rev. Lett. 59, 2776 (1987). [72] A. G. Khachaturyan and J. W. Morris, Jr., Phys. Rev. Lett. 61, 215 (1988). [73] A. G. Khachaturyan and J. W. Morris, Jr., Phys. Rev. B 37, 2243 (1988). 109 [74] E. D. Specht, C. J. Sparks, A. G. Dhere, J. Brynestad, O. B. Vavin, D. M. Kroeger, and H. A. Oye, Phys. Rev. B 37, 7426 (1988). Chapter 6 Dual Percolation Threshold in 2-d Microporous Media 6.1 Introduction Until now we have been focusing on the structural properties of various two—dimensional systems. Once a structure is given, one is interested in studying the connectivity in these lattice model, or physical problem that depends sensitively on this connectivity such as particle diffusion in 2-dimensional systems. Percolation and diffusion in two—dimensional (2d) microporous systems are of con- siderable physical interest [1, 2]. In recent years, synthesis of ternary intercalation compounds [3, 4], e.g. mixed ion clay intercalation compounds [4, 5, 6] of the form A1_,,B,L (where A and B are two types of intercalants with different sizes and shapes sandwiched between host layers-L of finite transverse layer rigidity) has made it pos- sible to study the adsorption/diffusion of probe particles C as a function of the 2d microporosity which is controlled by a: and the layer rigidity. Such a study is prereq- uisite to an understanding of the chemical, catalytic and molecular sieving properties of this important class of materials. Yet, the quantitative theoretical aspects of the interplay between topological and elastic properties of layered system have yet to be 110 111 addressed. Accordingly, in this chapter, we derive the functional dependence of the microporosity on the composition :17, the relative sizes of the species A, B, C, and on a characteristic length A, which we refer to as the healing length of the layers (A is a measure of the transverse layer rigidity). We also establish the characteristics of the probe particle percolation process for different microporosities. The ternary systems discussed here are excellent models for 2d pillared complex layered oxides which are potentially useful for microporous adsorbents and shape se- lective catalysts [5]. These layered oxides contain two types of pillar ion, A and B, which occupy the gallery space between the two oxides layers (see Fig.6.1). A chem- ically active species C can diffuse through the percolation network defined by A and B and interact with other similar species. Thus a study of both the adsorption and diffusion of C in this particular 2—d medium serve the dual purpose of clarifying the general fundamental aspects of microporosity in reduced dimension while ellucidating the corresponding properties of a novel subclass of layered solids[5, 6, 7, 8]. The arrangement of this chapter is as follows. In section 6.2 we define the porosity parameters and healing length for model 2d layered systems and give simple lattice models for three different healing lengths. In sections 6.3 and 6.4 we present analytic solutions and numerical simulation results for these models. Finally in section 6.5 we summarize our results. 6.2 The Model System The physical systems we are interested in are ternary layered intercalation compounds of the type A1_,B¢L. We assume that the lateral distribution of A and B ions is frozen, i.e. the system is in the quenched limit. The A and B ions randomly occupy the sites of a triangular lattice with probability 1 — an and a: respectively (Fig.6.2). Assume that the ions A and B and the diffusing particle C can all be represented by ellipsoids of revolution with diameters (2rA,hA),(2r3,hB) and (2rc,hc), such that 7‘A < 7'3, and hA < hg. 112 Figure 6.1: Schematic illustration of the tetrahedral and octahedral sites in a 2:1 layered silicate. Open circles are oxygen, closed small circles are cations in tetrahedral (Si, Al) and large circles in octahedral (Al, Fe, Mg, Li) positions. The intercalants such as Co(en)3 occupy the interlayer galleries. 114 We also assume that each of these species is oriented in the host gallery with the h axis perpendicular to the layers. The topological constraints on the access of 0 into A1_,B¢L depend upon the lateral porosity parameters: a-—-2rA ): (6.1) .91 2 sign( C a-—-rA -T£3 .32 2 sign( 2rc ), (6.2) . a — Zr .93 = szgn(—2;£), _ (6.3) where the function sign(a:) = 1 if a: > 1, sign(a:) = 0 if a: < 1; and a. is the distance between the two nearest neighbor intercalants (see Fig.6.3). In addition, these topological constraints depend upon two vertical porosity parameters: .34 = sign(hA/hc), (6.4) 35 = SignUlB/hc), (6.5) and the length A[9] which is a characteristic distance over which a local deformation caused by the insertion of a single large ion B into AL relaxes to he (see F ig.6.4). For soft (“floppy”) layers A = a/Z and for rigid layers A >> a. For a given geometric condition, the adsorption/ diffusion of C in this A1_,B.,L is determined by the set of porosity parameters {.9} = {abs-2,33,34,95} and A. Fur- thermore, we assume that the sites which the species C can occupy are the ver- tices of a honeycomb lattice which is dual to the original triangular lattice decorated by A and B. Diffusion of 0 takes place on this “dual honeycomb lattice” (DHL) whose bonds are open or closed depending on the local distribution of A and B ions as shown in Fig.6.3. We define a saturation value of the mass uptake function ”7(a) = M(z:)/MT(2:), where MT(2:) is the total mass of adsorbed 0 species when the entire available interlayer space is occupied by C'. M (1:) is the actual amount of C’ intercalated into the gallery. Clearly, 111(2) S Alfie) because the former depends on both availability and accessibility of a given DHL-site whereas the latter depends only on the availability. 115 Figure 6.3: Lateral constraints of C particle percolation on the triangular lattice of intercalants A, B. The dashed and heavy lines are, respectively, closed and open bonds on the dual honeycomb lattice (DIIL) when A = A1. In this case the BB bonds are closed due to lateral restriction and the AA bonds are closed due to vertical (gallery height) restrictions. Only AB bonds are open for the C particle to diffuse. 116 k A a. I O B U Figure 6.4: Vertical constraints of C particle percolation on the A1-,B,L system. 117 It is interesting to note that there are two percolation processes in the above model. The porosity subset (.91, .92, .93) determines the laterally controlled percolation process whereas the subset (s4, .95) and A determine the vertically controlled percolation pro- cess. The interplay between these two processes can give rise to 0, 1, or 2 percolation thresholds. An example of a physical system that can be modelled in the above fash- ion is the adsorption/diffusion of Argon, Krypton and Nitrogen in a clay intercalation compound [10, 11]. For Argon (3.6Ain diameter), {.9} = {1,0, 0,0, 1}, and the healing length is such that there is one percolation threshold. For Krypton and Nitrogen (3.8 Ain diameter) {.9} : {1,0,0,0,0}, and the healing length is such that there is no percolation region. In the present work we consider the case where {s} = {1, 1, 0,0, 1}. Physically this corresponds to the case where a C atom is laterally constrained to pass through only AA and AB bonds but not through BB bonds. The vertical constraints in general depend on the healing length A and in this paper we study four distinct cases: 1_ A = A1 = %, where C can only pass through an AB bond; 2, A = A2 : age, where C can pass through an AB bond and an AA bond which has‘at least one B atom at a nearest neighbor site; 3, A = A3 : 4‘1, where C can pass through an AB bond and an AA bond which has at least one B atom at a nearest neighbor or next-nearest neighbor site; 4_ A = A4 —» 00, where C can pass through both an AB bond and an AA bond. Experimental realizations of systems with different healing lengths are intercalation compounds of graphite (short), layered dichalcogenides (intermediate) and complex layered oxides such as pillared clay (long) [3, 9, 12]. Potential candidates are stage-2 graphite intercalation compounds such as Rb1_:Cs,C24. 118 6.3 Solution for the Case of A = A1 For A = A1, the possible topological configurations of the A1_,,B:L system within a triangular plaquette are shown in Fig.6.5. It can be seen from this figure that the open bonds on the DHL appear in connected pairs. There are no single open bonds (dangling ends) or triple open bonds (branches). Thus the clusters on the DHL consist of nonintersecting closed “rings” which are in 1:1 correspondence with the external perimeter (hull) of the cluster of species A and B on the original triangular lattice (see Fig.6.3). As a consequence, there is one a percolation threshold at :v = we, above or below which there is no infinite percolation cluster and we is equal to the site percolation threshold on a triangular lattice {i.e. arc : 0.5). Proof: For a: > 0.5 (:c < 0.5), there exists an infinite cluster of B (or A), which has no external perimeter. The external perimeter of the finite cluster is of course finite. For a: = 0.5, there exists an infinite cluster of B (or A) on the triangular lattice with fractal dimension D = 3% < 2 [1]. Thus, only at this point, the external perimeter of the infinite cluster is also infinite and so is the corresponding cluster on the DHL. Therefore there exists at and only at ac : 0.5 an infinite cluster on the DHL. Thus for A = A1 there exists a single percolation threshold me = 0.5 and there are no percolation clusters present on the DHL either above or below arc. Percolation is one type of phase transition. Therefore it will be interesting to compare its behavior near the percolation threshold to the critical behavior near the critical temperature in thermal phase transition as we discussed in chapter 2. We can define critical exponents and universality class associated with this phase transition just like what we have done for the thermal phase transition [1]. The critical exponent a of this model when A = A1 can be obtained from the scaling law.[1] On the triangular lattice, there is only one external boundary for a finite cluster of A or B. Thus the 119 as. as. Figure 6.5: Possible configurations of A1_,B,L in a triangular plaquette. Dashed and heavy lines are, respectively, closed and open bonds on the DHL when A = A1 (see text). In this case the BB bonds are closed due to lateral restriction and the AA bonds are closed due to vertical (gallery height) restrictions. Only AB bonds are open for the C particle to diffuse through. 120 total number of clusters on the DHL is given by EMA-’0) = 2713(3) + 271336”), (6-6) 0 0A (B where an(:c) is the number of clusters of connected A sites on the triangular lattice with size 3A and nfB(a:) is the number of clusters of connected B sites on the triangular lattice with size .93. When a: —» 0.5, 2 7134(3) or [2: — 0.5]2_°‘°, (6.7) 5A 2 naAB(.r) or [93 — 0.5]2_°‘°, (6.8) 53 where ao(= —2/3) is the known exponent for the 2-d site percolation.[1] For the case A = Al, the concentration of the open bonds on the DHL is p and it is given by the relation p = 2:6(1 — 2:); therefore the total number of clusters on the DHL, Z, n.(p) as a function of p is given by Z n.(p) 0< I? ~ 10.;I1‘“°/2 = I? - Pclz‘“, (6-9) where p6 2 p(1: = :1:6 2 0.5) = 0.5. Thus we have 2 — or = 1 — 010/2 or a = 2/3. The fractal dimension D of the percolation cluster on the DHL is equal to the fractal dimension of the hull of the percolation cluster of the triangular lattice at the site percolation threshold. An exact calculation of fractal dimension of this problem was given by Saleur and Duplantier [13], who found that D = 4:. Using these results we can obtain all the geometric exponents by using the scaling law[1]: _2 _6 _15 _1 V—§:U—?aT"'.—7_17— ‘ (6.10) Having studied the static (geometric) properties of this model, we now study to its dynamic (transport) properties. Diffusion and transport properties on percolation and fractals have been extensively studied. The analysis of transport in disordered media by means of diffusion on percolation clusters was suggested by de Gennes [15]. The idea was to perform random walks on a percolation system, for which de Gennes 121 coined the term “ant in the labyrinth”. This can be used to measure the diffusion constant of the system. The theory of random walks has been applied in many areas of science, espe- cially as a model for transport phenomena. In uniform Euclidean systems, the mean- square displacement of a random walker, < R2(t) >, is proportional to the time t, < R2(t) >0: t, for any number of spatial dimensions d (Fick’s law). However, in dis- ordered systems, this law is not valid in general. Rather, the diffusion law becomes anomalous [14]: < R2(t) >~ W“, (6.11) with dw > 2. This slowing down of the transport is caused by the delay of the diffusing particles in the dangling ends, bottlenecks and backends existing in the disordered structure. In our case, the anomalous-diffusion exponent dw for A = A1 can be calculated by noting that clusters on the DHL are topologically l-d rings. Following the calculation by Alexander and Orbach [14], we compute the spectral dimension d, = 2D/dw = 1, where dw is the anomalous diffusion exponent defined by the Eqn. 6.11. Thus we get cl", = 3.5. Therefore when A = A1, this model belongs to the universality class of hull percolation of a 2-d system. Similar critical behavior has also been found in some continuous system models, e.g. the zero threshold (pc 2 0) Trugman and Weinrib model [16], although it belongs to yet another universality class due to different scaling laws. It is however interesting to note that the fractal dimension of the percolation cluster in Trugman et.al’s model is the same as that of our model for the A = A1 case. To test the scaling theory, we have performed Monte Carlo simulation on a 20, 000 x 20,000 lattice using the Hoshen and Kopelman algorithm [17] to verify the scaling law for the cluster size distribution: n, or s‘T at the percolation threshold. We have used free boundary conditions and ten samples were used to obtain an estimate of the statistical error in the value of the r exponent. Fig.6.6 shows the result for the case of A = A1. From a least square fitting of n, or s" we find 7’ = 2.130 :l: 0.009, in 122 excellent agreement with the analytic solution (see Eqn. 6.10). We also performed a series of simulations to establish the fractal dimension D, of the percolation cluster for lattice size L varying from 100 x 100 to 2000 x 2000. Fig.6.7 shows that the size of the largest cluster Sm“, or LB, with D = 1.757 :1: 0.008. We have noticed that T and D satisfy the hyperscaling relation [1, 18]: = (6.12) where d is the spatial dimension (d = 2). These results agree well with the previous simulation studies of hull percolation by Voss [19]. As for the anomalous diffusion exponents dw which is defined by the scaling relation: < R2(t) >oc ti, where < R’(t) > is the average displacement of a random walker from the origin after t time steps when the system is at the percolation threshold, there are two different ensambles which give different values of dw.[2] 1 In one ensamble the random walker is allowed to be in any cluster for which we define the exponent as d1”. In the second ensamble the random walker is only allowed to move in the infinite cluster, and for this we define the exponent as dw. Obviously (12,, > dw. In fact they satisfy the scaling law:[2] 9.. ' D if = 1 + —2—-(2 — r). (6.13) These exponents were calculated using the exact enumeration method [2]. A 1, 000 X 1,000 lattice and 100 different configurations of cluster were used. Fig.6.8(a), (b) give dw and (11,, as a function of t. The error bar, if plotted, is smaller than the size of the dots shown in these diagrams. For long t, dw —> 3.6, whereas dz” —+ 3.97. The discrepancy of the simulation results with the above scaling law (eqn. 12) and the analytical results presented before is probably due to the finite size of the simulated system. Simulation in larger systems and longer time needs to be done to clarify this issue. 123 s_-.'T"""'-.fi- O ' 0. X8805 . ‘ ' o faadm.“' - Q < I . a O o - . - A . 1 d " o z p . ‘ a . -, . 3. n . ‘ 5 " ‘ 1 . . l! j I l I L. l A | L L I I-l l 40 2.0 4,0 5.0. 3.0 Logts) Figure 6.6: The cluster size distribution for A = A1 and for a 20, 000 x 20, 000 lattice; N, is the number of clusters with size s. 124 r f r V I I fi' fi .. 5.0 - x g . a I o A I 3.0 - . Xc-0.5 I I L I I l l l A ”1.0 2.0 3 0 Figure 6.7: The system size dependence of the size of the largest cluster 5,“, at a: = 2‘, when A = A1. 125 ‘.5 r'r‘l'r"rr1r'rrrrr|uvu - ‘3’ I 1 4.0 - _: . 3 ' 'l, ‘0 " g 3.5 - '_j - '1 P d 3.0 1. . I 1 0‘ 2000 4000 ‘.0 r r r I l I l I r l u y . I I I r F1— » (bl ‘ I. e a—q fi 2.5_ L1 1 A LL 1 L r l J L_J 1 L 1 1 a 1 A 0 1000 2000 t Figure,6.8: Plot of d", and d1” ble in which diffusing particle can be in any cluster; (b) the ensamb 61188111 the diffusing particle can only be in the largest cluster. against t for the case A = A, when a: = 0.5: (a) the le in which 126 6.4 Solution for the Case of A > A1 When A > A1, in addition to the pairs of open bonds for the configuration ABB or AAB given in Fig.6.5, a triplet of open bonds may appear in the configuration AAB under certain conditions and in the case of A = A3 or A4, may even appear in the configuration AAA if there are B atoms nearby to open up the gallery height. Therefore an infinite cluster of open bonds exists if there exits a percolation cluster of A (i.e. if a: < 0.5) on the triangular lattice. However, The cluster of B is still “impenetrable” due to the lateral constraint. Therefore there will be no infinite cluster on the DHL if there is a percolation cluster of B on the original triangular lattice (i.e. if a: > 0.5), so mcz = 0.5 is the percolation threshold for A 2 A1. This we refer to as the laterally controlled percolation threshold. It is obvious from the above discussions that when A > A1, it is essential to have a minimum concentration of B atoms (real) to open up the gallery so that enough triplet bonds can be formed to have a percolation cluster on the DHL. There are no perco- lation clusters for :c < :ccl since the diffusing particle will be vertically constrained in that region. Therefore in addition to the threshold 2:3(2 0.5) associated with the lat- erally controlled percolation process, there is another threshold, 13.:1, associated with the vertically controlled percolation process. Thus two percolation thresholds flier and 3:62 exist when A > A1 and the infinite percolation cluster on the DHL exists when :vcl < :c < 2:52. The value of (rd depends on the healing length A. When A = A; —> oo, i.e. when the L-layer is infinitely rigid, 33a —i 0, thus recovering the ordinary single threshold percolation at mc : 0.5. To find the saturated weight uptake function ”7(a) defined in section II, we per— formed Monte Carlo simulations on 1000 X 1000 lattices, averaging over 100 different configurations for each value of m. Fig.6.9 gives W(:c) for the cases in which A = A2, A3 and A4. By studying the size dependence of the percolation threshold (cc (deter- mined as the peak position of the average cluster size S' = 2, 32m) for system sizes between 100 X 100 to 2000 x 2000, we find that all three have the same upper threshold 127 {Caz = 0.500 :l: 0.002 whereas the lower threshold :0“ = 0.215, 0.117 and 0 for the three cases respectively. These simulation results confirm our qualitative discussions given in the last paragraph. Fig.6.10 gives a phase diagram in the a/A vs. a: plane separat- ing the percolating (P) and the nonpercolating (NP) regions for the 2-d microporous medium under investigation. We also did simulations to determine the critical exponents, again on a 20,000 X 20,000 lattice and averaged over 10 different samples to obtain an estimate of the statistical error in the exponent 7'. For A = A2, we find that r(:ccg) = 2.053i0.009 and r(:z:c1) = 2.03 :l: 0.01 which agrees well with the value for the ordinary 2—d percolation m 91.[1] Similarly for the case of A = A3 and A = A, we found that r = To. exponent r0 = In all of the above cases the fractal dimension D of the percolation cluster was found to be 1.889:t0.008, very close to the exact 2—d value for ordinary percolation D0 = 2715. Therefore for A > Al the model belongs to the same universality class as the ordinary 2-d percolation system. The exponents of anomalous diffusion at the percolation threshold for the two different ensambles are found to be, respectively, dw = 2.86, d; : 3.00, which are consistent with the results obtained from the simulation of a 2-d percolation system on a triangular lattice by Ben-Avraham et.al [20] and by Pandey et.al [21]. 6.5 Summary In summary, we have constructed a simple lattice model to probe the 2-dimensional microporosity of layered intercalated systems A1_,B,,L where A and B are two types of intercalants sandwiched between host layers—L. Percolation of a probe particle in this 2-d microporous medium has been studied both analytically and by Monte Carlo simulation. Depending on the transverse layer rigidity which is characterized by a healing length A , we find two percolation thresholds which coalesce for a given A . Percolation at this point belongs to a different universality class from the ordinary 2-d site percolation although the correlations are of short range. Anomalous diffusion 128 rvvl'rTrUrrrrrrVIvu't'ttrrI' 1.0 wooooygeaaaggagg..... a A a A ++++ N (X) 0 Haaagq-«H-i Ll.LLllllLlLlIIILILLLLIIILLL 0 0.10 0.20 0.30 0.40 _.0.50 0.60 X ++ O. O. 95". :0. . Figure 6.9: The order parameter (the saturated weight uptake function W(z) defined in the text) for the different values of the healing length: +, A = A2; A A = A3; 0, A = A4. 129 Figure 6.10: The phase diagram in the a/A vs. a: plane separating percolating (P) and nonpercolating (NP) regions for the A1-.B.L system. the “+” symbols are the percolation threshold obtained from Monte Carlo simulation for four different values of the healing length A and the dashed lines are simple linear extrapolations. 130 on these systems have also been investigated. On the theoretical side, we find that for certain types of microporosity the 2d ternary system possesses two percolation thresholds as a function of as; one (2261) controlled by /\ and the other (3.9) independent of x\. In the limit of infinite x\, the first threshold goes to 0 thus recovering the usual single percolation threshold response. In the limit of short x\, the two thresholds coalesce, i.e. 13d 2 ma 2 sec. Whereas for the ma 7e we; case, the critical exponents are the same as those for a normal 2d percolation, the critical exponents for the me] = 13a case belong to a different universality class which turns out to be the “hull” percolation in 2d systems. A paper based on this investigation is going to be published in the Physical Review [22]. Experiments to probe the percolation and diffusion properties of a diffusing par- ticle C in real A1_¢B,L type systems to study the different characteristic responses predicted by this model are underway [23, 24]. An excellent candidate is argon atom diffusing through a fluerohectorite clay intercalation compound where A and B are [Co(en)2]3+ and [Cr(en)3]3+, respectively, and en is an organic ethylene diamine ligand.[23] Bibliography [1] Dietrich Staufer, Introduction to Percolation Theory (Taylor & Francis Inc., 1985). [2] S. Havlin and D. Ben—Avraham, Adv. in Physics 36, 695 (1987). [3] S. A. Solin and H. Zabel, Adv. in Physics 37, 87 (1988). [4] R. M. Barrer, Zeolites and Clay Minerals as Sorbents and Molecular Sieves, (Academic Press, 1978). [5] T. J. Pinnavaia, Science 220, 365 (1983). [6] P. Laszlo, Science 235, 1473 (1987). [7] M. H. Kopelman and J. G. Dillard, Clays Clay ll/Iiner. 28, 211 (1980). [8] P. A. Plitowicz and J. J. Kozak, J. of Chem. Phys. 92, 6078 (1988). [9] S. Lee, H. Miyazaki, S. D. Mahanti and S. A. Solin, Phys. Rev. Lett. 62, 3066 (1989). [10] H. Kim, T. J. Pinnavaia, Zhi-Xiong Cai, S. D. Mahanti and S. A. Solin, Bull. Am. Phys. Soc.35, 443 (1990). [11] Ping Zhou, S. Lee, S. A. Solin, H. Kim and T. J. Pinnavaia, Bull. Am. Phys. $00.35, 626 (1990). 131 132 [12] S. A. Solin, Intercalation in Layered Materials, ed. by M. S. Dresselhaus, Plenum, New York, 1987, p.145. [13] H. Saleur and B. Duplantier, Phys. Rev. Lett. 58, 2325 (1987). [14] S. Alexander and R. Orbach, J. Phys. Lett. (Paris) 43, L625 (1982). [15] P. G. de Gennes, La Recherche, 7, 919 (1979). [16] S. A. Trugman and A. Weinrib, Phys. Rev. B 31, 2974 (1985). [17] J. Hoshen and R. Kopelman, Phys. Rev. B 14, 3428 (1976). [18] R. Zallen, The Physics of Amorphous Solids, Chapter 4, (A. Wiley-Interscience Publication, 1983). [19] R. F. Voss, J. Phys. A 17, L373 (1984). [20] Daniel Ben-Avraham and Shlomo Havlin, J. Phys. A 15, L691 (1982). [21] R. B. Pandey, D. Stauffer, A.Margolina, and J. G. Zabolitzky, J. Stat. Phys. 34, 427 (1984). 22 Zhi-Xion Cai, S. D. Mahanti, S. A. Solin and T. J. Pinnanvaia, Phys. Rev. B, S (to be published). [23] H. Kim and T. J. Pinnavaia (unpublished). [24] S. Lee, Ping Zhou, S. A. Solin, H. Kim and T. J. Pinnavaia, (unpublished). Chapter 7 Summary and Discussion In this thesis we have used Monte Carlo simulation technique to study various as- pects of structural properties of different two-dimensional model systems. From these studies we can conclude several important features of computer simulation. First of all, the simulation itself is merely a mathematical method to evaluate mul- tidimensional integral as described in Chapter 2. It will be meaningless without the properly established model which describe at least some aspects of the system we are studying. The phrase often heard in the community of computational physicists “Garbage in Garbage out” described vividly the danger to over-emphasis and over- generalize the simulation results without justification. The model hamiltom'an should be established based on the physical law of the microscopic interactions in the system, not based on the “expected” results. For example, in chapter 5, the lattice gas model describing oxygen ordering in Y-Ba-Cu-O system was derived based on the screened Columb interaction and the short range covalent bonding. From this derivation we know that this lattice gas model is only valid for oxygen concentration greater than 6.5. Therefore simulation study of this model at low oxygen concentration will not be meaningful. Moreover we can modify this model to study the effect of impurity. If we establish a lattice gas model for the sole purpose of reproducing the experimen- tal phase diagram, we cannot gain much insight about the microscopic interaction. 133 134 “Parameter fitting” should not be the main purpose of computer simulation. Secondly, we can obtain precise results of the multidimensional integral only if we scan the entire phase space of the integration. In other words, the Monte Carlo steps should be sufficiently large in order to minimize the statistical error. Exactly how long the simulation should run will depend on the quantity we are interested in and also depend on the phase space. In chapter 4,5, and 6 we are dealing with problems in the phase space with discrete degree of freedom, and we are only interested in the qualitative understanding of the phase diagram. Therefore a relatively short Monte Carlo simulation time is sufficient to obtain the results we need. In chapter 3, however, we are interested in more quantitative aspects of the phase diagram, namely the determination of the critical exponents. Also the model we studied, the 2-dimensional planar rotor model, has continuous degree of freedom which made it a rather difficult problem to tackle. We saw from the earlier simulation works that short simulation time sometimes can give qualitatively incorrect results. Since there are no other reliable methods to study the model we described in chapter 3, the incorrect results due to short computation time created some confusion over the years. Finally, thermodynamic quantities have definite meaningonly for infinite systems, therefore it is important to perform simulations on sufficiently large system. Again this depends on the model we study, in particular, it depends on the correlation length. The finite size scaling theory gives us a powerful tool to approximate the properties of infinite system by doing finite system simulations. Thus for a quan- titative understanding of the model we are interested in, it is always necessary to study the system size dependence of the quantity we are calculating. This, however, has to be compromised sometimes to allow sufficient simulation time to minimize the statistical error as we discussed above. We see from the discussion in this thesis that the study of structural phase transi- tions in two dimensional layered systems is a broad and exciting field. New materials and new phenomena offer continued challenge to the theorists. On the other hand, 135 more powerful computers and new methods of simulation make it possible for us to study many problems which previously could not have been solved. The method described in this thesis can also be used in other systems, such as commensurate- incommensurate transitions, kinetics of domain growth and low temperature excita- tion of disordered system such as spin glass, etc. NICHIGQN STQTE Hl W l . 2930 903 - ”3‘2...”- . I it" 51.: :r' . .. \r‘ .‘, 3:...1‘.