MSU RETURNING MATERIALS: P1ace in book drop to LIBRARIES remove this checkout from _;-—. your record. FINES win be charged if book is returned after the date stamped be10w. ON THE TRANSPORT PROPERTIES AND DYNAMICS OF DISORDERED SYSTEMS By Neizhu Xia A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirments for the degree of DOCTOR OF PHILOSOPHY Department of Physics and Astronomy 1988 ABSTRACT ON THE TRANSPORT PROPERTIES AND DYNAMICS OF DISORDERED SYSTEMS BY Heizhu Xia This thesis presents recent studies of the disordered systems in condensed matter physics. It addresses the problems of effective steady-state transport properties, such as electric conductivity and elastic moduli, and dynamic responses of strongly disordered systems. In Part I, the continuum percolation of a system containing random el- lipses is studied. The percolation thresholds are obtained, for various elliptical aspect ratios, from computer simulations. The macroscopic effective steady-state conductivity for this system is studied by incor- porating the properties of effective conductivity at low and critical concentrations of inclusions. A set of semi-phenomenological interpola- tion formulas is derived and agrees very well with the experimental data over the whole range of concentrations. In Part II, the elastic per- colation problem of a stretched spring model on the dilute honeycomb lattice is studied. While we find the similar second order rigid + floppy phase transition studied before on the triangular lattice, an ad- ditional first order rigid + floppy phase transition is found above a tricritical point. The bulk modulus, as well as some other elastic con- stants, behave differently when the two different phase transitions occur. A Landau type phase transition theory is applied to draw an analogy between the two types of phase transitions. A self-consistant effective medium theory is also developed for the phase boundary and the tricritical point observed in the computer simulations. In Part III, the density of states (DOS) of vibrational excitation spectra of per- colation networks at thresholds are obtained by computer simulations using the equation of motion method. A careful study of the low fre- quency part of the density of states shows that the spectral dimensionalities, extracted from DOS, agree well with the predictions of the scaling hypothesis which takes into account the critical scaling be- haviour of both mass and elastic moduli. This direct method of (finaining spectral dimensionalities is superior to the random walk method in the superconducting-normal network. Part I and Part III are the expanded versions of the following pub- lished papers: I. H. Xia and M. F. Thorpe, "Continuum Percolation of Ellipses", Phys. Rev. A. (1988) (in press) 111. R. A. Day, H. Xia and M. F. Thorpe, "Spectral Dimensionality of Random Superconducting-normal network", Phys. Rev. B, 37, 1339, (1988). This is contained in the appendix. iv ACKNOWLEDGEMENTS First of all I would like to thank my thesis advisor Professor M. F. Thorpe for his careful and patient guidance over the last four years. His knowledge of solid state physics and his approach to re- search have been a real inspiration to me. I would also like to express my appreciation to Professor S. D. Mahanti for valuable lectures in statistical mechanics, solid state physics and many helpful conversa- tions, to Professor H. Repko for interesting lectures in quantum field theory and high energy physics, to all of the three above and Professor J. Cowen for serving on my guidance committee. I gratefully acknowledge Dr. A. R. Day for his enlightening collaboration on the work described in Part II and III and for stimulating discussions. I also benefit a great deal from useful conversations with Dr. S. de Leeuw and Dr. P. Sen. Dr. E. Garboczi, Dr. H. He, Dr. N. Tang, Mr. Y. Li, J. Gales, N. Jin andli.Yan have been good friends and a source of much help and relaxed discussions throughout my graduate career. I would like to thank Dr. J. S. Kovacs and Mrs. J. King for additional assistance. I would also like to express my gratitude to my dear wife Lihua for all her patience, encouragement and assistance over last three years. Financial support from the Office of Naval Research, the Michigan State University Center for Fundamental Materials Research and the National Science Foundation is gratefully acknowledged. TABLE OF CONTENTS List of Tables List of Figures General Introduction Part I. Continuum Percolation and Macroscopic Conductivity of Random Composite Systems I. Introduction II. Continuum Percolation III. Computer Simulations and Results IV. Critique of Effective Medium Theories V. Interpolation Formulae VI. Summary of Results VII. References Patr II. Elastic Percolation of Stretched Spring Model on Honeycomb Lattice I. Introduction 11. Elasticity of Pure Honeycomb Lattice and Percolation Thresholds III. Computer Simulations and Effective Medium Theory Approach IV. Tricritical Point and Effective Medium Theory V. Conclusions vi viii ix 1“ 15 17 29 36 A0 50 51 5A 55 60 81 101 120 VI. References Part III. The Vibrational Spectra of Random Superconducting- Normal network I. Introduction 11. Density of States III. Scaling Arguments IV. One Dimensional Case V. Conclusions VI. References Appendix I. Elasticity of Pure Honeycomb Lattice Appendix II. Percolation Thresholds of Stretched Spring Model on Honeycomb Lattice Appendix III. Published Paper vii 121 123 1211 129 11111 152 158 159 162 168 17a Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. 1.2 1.3 1.“ 1.5 1.6 1.7 2.1 2.2 2.3 2.“ 2.5 2.6 2.7 2.8 2.9 LIST OF FIGURES An example of randomly oriented ellipses with aspect ratio b/a = 0.5. (a) A prolate with a>b and b=c; (b) A oblate with a>b and a=c. Percolation thresholds pc for various b/a. f2: 1 - exp(1abnc) and f1: pc's from EMT and computer simulations. Iabnc vs. b/a. Conductance from the interpolation formula. A comparison between interpolation formulae and Experiments. Bulk moduli vs. p for central force model. Honeycomb lattice. Extrapolated intercepts. Extrapolated intercepts P(E), p(T) and p(B) vs. LO Equilibrium positions and the corresponding elastic potential. Simulation results of elastic potentials. b Keff and Keff vs. LO/L. E0, T and B are plotted vs. p for various LO/L. Equivalence relations are checked. Quantities with sudden discontinuity at percolation threshold. viii /L for honeycomb lattice and triangular lattice. 19 26 32 35 39 A“ A9 57 65 7O 75 78 79 80 82 8A 86 Table 1.1 Table 1.2 Table 3.1 LIST OF TABLES Values of pa, n k, and (Aex) for randomly c, oriented ellipses for various aspect ratios b/a. 33 Values of pc, no, k, and (Aex) for ellipses with two allowed orientations. 34 Comparing the spectral dimensionality 8 obtained directly in this work with the scaling relation. 150 ix Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. 2.13 2.1” 2.15 2.16 2.17 .20 wwwm N wwww O‘U‘lt Phase diagram. Configuration with zero B for LO/L = 1. 0 eff/Leff from simulations. Varoius elastic constants from simulations. L CAN/C11 are plotted vs. p for various LO/L. Critical values of Gun/C11 vs. LO/L. Phase daigram and the tricritical point. Phase boundary and the tricritical point. nu plotted against p. E, T and B plotted against p. Diagrams of monomer, dimers, trimers. A bond percolation network on square lattice. The DOS for a random square lattice of Nl—a superconducting and normal bonds at pc = The DOS for a simple cubic lattice at pc Diagrams of 'pockets' in 2D. Diagrams of 'pockets' in 3D. r(E/Eco) vs. E/Eco for p = O.A75, O.A80, 0.u85, 0.H90 and 0.A95. The DOS for a random linear chain. The linear chain results for p = 0.999 and the two limiting forms. .249. 87 89 9A 96 99 100 106 117 118 119 133 136 137 138 1A1 1H3 151 156 157 GENERAL INTRODUCTION The study of macroscopic transport properties of inhomogeneous sys- tems has been a subject of great interest for more than a century. A good historical review on the subject is given by Landauer‘. When studying the macroscopic properties, solids are classified, at a rela- tively coarse level, into two categories: homogeneous materials and inhomogeneous ones. Inhomogeneous materials (particularly composite materials) can be treated as a mixture of different consitituent grains. The consitituents at some typical scale, for example less than grain size, can be treated as homogeneous. Mineral rocks, sandstones, con- crete, cast-iron and conductor-insulator alloys are Just a few examples. These composite materials are characterized by the so called effective transport properties such as electrical and thermal conductivities, which reflect the average behavior of the bulk materials. As the ap- plications of composite materials become more and more specialized, accurate predictions of the overall bulk transport properties are re- quired in order to make effective use of these materials. While the macroscopic transport properties such as the electrical and thermal conductivity of homogeneous materials are well understood and can be found in many handbooks, their counterparts in the in- homogeneous case are not. As a matter of fact, these transport properties are very difficult to calculate except for some special cases where the microscopic geometric structures of composite systems are simple. Maxwell, Clausius, Mossotti, Lorenz and Lorentz are the names 1 Modern theories associated with early attempts to solve this problem. started when Bruggeman published an extensive treatment of dielectric properties of two phase materials in 1935.1 Since then, there have been many theoretical approches to the problems1'u. Among the various ap- proaches, Effective Medium Theory (EMT)1’5’9 and absolute bounds6’7’8 (i.e. upper and lower bounds) are two major ones. The ab- solute bounds are derived from the variational methods and are useful when a real calculation is impossible. There are many versions of EMT which can be regarded as perturbation expansions where calculations can only be made for relative simple weak inhomogeneous systems (also called weakly disordered systems). Simple refers to the shape of the geometric boundaries between the phases and 219.115. means that the transport properties, (e.g. electric conductivities) of the different phases are similar. The effective conductivity is calculated by expanding the dif— ference of conductivities (a small parameter) in a power seriesg. The results agree well with experiments in most cases. For strong in- homogeneous systems (also called strongly disordered systems) in which the ratios of conductivites of the phases are infinite (i.e. mixture of conductors-insulators or superconductors-conductors), however, the above perturbation theories fail because the expansion parameters are no longer small. Absolute bounds (i.e. upper and lower bounds) in this case turn out to be very large and practically useless. Computer model- ing of the transport properties has also been impossible, because even the largest available machines cannot store enough information to meaningfully discretize the composite systems for complicated geometric structures. Therefore a semi-phenomenological description for transport properties of strongly disordered continuum systems is needed in order to make overall predictions for systems with certain geometric distributions. To show the underlying physics and simplicity of our ap- proach, we only consider two phase composite systems in this thesis. The equations describing the transport properties of processes such as electric conductivity, thermal conductivity, dielectric displacement, magnetic induction and diffusion all have the same mathematical structure. Therefore one only needs to study one of these and then generalize the results to all cases. We choose electric conductivity for convenience. The bulk effective conductivity, denoted by Z, is in general a mnction of both the conductivity and geometric distribution of the constituent phases. It can be written as <3) = = zeFf (1) where (J) and (E) are macroscopic current and applied electric field which can be measured experimentally. 0 denotes the average over samples with the same constituents and statistical geometric distributions. In a continuum conductor-insulator composite system, when the fraction of conductor is below a certain threshold, no current is able to flow across the system. This threshold can be well described by a continuum percolation theory‘o. Unlike the lattice percolation 11,12 problem , there is no underlying lattice in continuum percolation problem. In the case of the conductor-insulator disk continuum percola— tion problem, one can imagine that on an uniform background of conducting matrix (with finite conductivity 00) identical circular holes are punched out randomly (overlaps are allowed). As more and more holes are punched out, the system will not carry any current below a critical threshold (called percolation threshold) even if a voltage is applied. Although there is no underlying lattice in continuum percolation problem concepts such as percolation threshold, correlation length and cluster size etc. are still similar to those of lattice percolation13. The con- ductivity vanishes at percolation threshold according to an exponent depending on the dimensionality of the problem”. Of course the deter- mination of the percolation threshold is more difficult in continuum percolation problems than in lattice percolation especially when the ob- jects are irregular. As a matter of fact, the continuum percolation thresholds even for many regular objects were still unknown. Prior to this work the percolation thresholds are known only for circles and parallel ellipses in 2015. In Part I, we consider the continuum per- colation of identical random elliptical holes with various aspect ratios ranging from disk-like objects to needle-like ones. We obtain the per- colation threshold from computer simulations. Then we develop a set of semi-phenomenological interpolation formulas for overall bulk behaviors for composite systems consisting of random identical elliptical holes. The interpolation formulas are based on the percolation threshold and the single defect effective medium theory which is exact in the low con- centration limit. The transport properties such as electric conductivity, thermal conductivity, dielectric constant, magnetic permeability and diffusion coefficient are scalar quantities. Materials also possess elastic properties. Elastic constants, however, are fourth rank tensors and therefore elastic composite systems in general are more complicated and difficult to treat. Due to the nature of the problems we only consider solids with covalent chemical bondings. In studying elastic properties, chemical bonds can be thought of as springs with certain elastic energy potentials depending on the nature of the problems and the atoms serve as nodes linking up the springs. de Gennes and Stauffer were the first16 to study the sol + gel phase transition defined below using the site elastic percolation model on a lattice. In the process of the sol ¢ gel phase transition some monomer molecules are initially dissolved in a liquid solution. As the chemical reactions continue, monomers form finite polymer molecules (sol molecules) through covalent chemical bondings. In the sol phase the sol molecules (finite molecules) are separated by liquid which does not resist any shear deformation, so that the system in sol phase has zero shear modulus. As the sol molecules become larger and larger, a micromolecule crossing the whole system (gel molecule) will eventually appear. Because the gel molecule can resist the shear defamation, the shear modulus for the system is no longer zero. Of course the shear modulus will increase as more and more cross links in gel molecule are built up. The following process Na SiO 2 3 + 3H20 . 2NaOH + HuSiO u is an example of above mentioned sol + gel processm. The molecules of ”1181011 are monomers and will stick together to form sol molecules and eventually gel molecule through covalent chemical bondings. By using an idealized elastic percolation model, de Gennes and Stauffer were able to capture the major features of the above sol + gel phase transition 16. For example, adding sites and connecting springs correspond to sticking together sol molecules; the percolation threshold corresponds to the sol+gel transition point; elastic shear moduli vanishes below percola- tion threshold and sol + gel transition point etc. Other elastic percolation models can also be used to study the elastic behavior of compositional chalcogenide glass such as Se x and y17. 1-x-yASxGey as a function of While the model proposed by de Gennes (called isotropic force Mel) can be used to explain the sol + gel phase transition, it also serves as a primary model among various elastic percolation models. de Gennes pointed out that the elastic percolation problem of this model is identical to the conductivity percolation problem which had been studied “’12. This is reflected by the fact16 that the two in great detail problems can be mapped into each other and therefore (a) the percolation thresholds of elastic and conductivity percolation problems are the same; (b) the critical exponents describing the vanishing bulk moduli and conductivity near percolation thresholds are identical. Feng and Sen18 proposed the central force model and pointed out that the elastic percolation described by central force model is different from that of isotropic force model. In fact they belong to different class of the problems18’19. For example, in two dimensions the percolation thresholds, according to Thorpe's constraint counting methodzo, are pc = g for isotropic force model and pcen = 121 for central force model. Here 2 is the number of nearest neighbors of a site. In the above two models, the elastic energies are expressed in quadratic forms of small displacements. Tang and Thorpe recently” introduced a rotationally invariant stretched spring model which leads to a more natural way of expressing the elastic energy for Hooke springs. It turns out that the isotropic force model and the central force model are the two extreme limits of the stretched spring model. Tang and Thorpe pointed out that, by continuously changing a parameter, the stretched spring model can serve as a bridging model between the isotropic force model and central force model. In general there is a lot of stress associated with stretched or compressed springs in this model. Most computer simula- tions have been done on triangular lattice in the stretched region where initially every spring is stretched. In second part of this thesis we study the elastic percolation of stretched spring model on honeycomb lattice. The motivation of this study is as follows. 17 20 Both the constraint counting method and effective medium theory predict that the percolation threshold for the central force model on 8 honeycomb lattice is pcen= 13‘. This result is unphysical because p must always equal or less than one. This tells us that the central force model on a honeycomb lattice is not stable. One may find that the shear nndulus is zero for the honeycomb lattice with the central forces. There is no such instability on the triangular lattice. Therefore, by studying the honeycomb lattice with the stretched spring model, we ex- pect to see some new phenomena associated with the instability not observed on the triangular lattice. Indeed we have found that the usual second order rigid » floppy phase transition becomes a first order phase transition above a tricritical point. We also want to determine the percolation threshold for the stretched spring model in the central force limit. The properties we have just discussed above are the static properties of the disordered systems. It is also very interesting to study the dynamic properties of the disordered systems particularly in the low frequency (or energy) limits. The quantity that concerns us here is the density of states of vibrational excitations which gives the number of excitational modes per unit frequency interval. In a homogeneous system, in the very long wavelength or low frequency limit, the density of states denoted by g(E) has the Debye f'orm22 g(m) dw ~ to do) (2) where m is the frequency of the excitation modes and d is the dimen- tionality of the problem. The simplest way to obtain the density of states g(m) is to calculate the excitational dispersion relation (phonon dispersion relation) on a perfect crystal lattice and consequently ob- tain g(m). In the long wavelength limit, i.e., the wavelength of excitation is much larger than the lattice spacing, the lattice is regarded as a homogeneous continuum and (2) is obtainedzz. One can ask the question of how the density of states changes when disorder is introduced. A good review article concerning the general aspects of this question is given by Elliott et al.23 Here we only discuss the low frequency limit. From one's physical intuition, it can be expected that when the wavelength of excitation is comparable to the length scale on which the system is disordered the density of states should have dif- ferent form than (2) because the system is no longer homogeneous. On the other hand, however, (2) should still be observed in the very long wavelenLth limit because then even disordered systems are homogeneous . Alexander et al.2u pointed out that in dilute disordered system (or in- homogeneous system) (i) the density of states should have the following form 5-1 g(w)dm ~ w dw (3) where a is called the spectral dimensionality which is different from d for the homogeneous case; (ii) the density of states has a cross-over 10 from the homogeneous region to the inhomogeneous region at frequency woo; and (iii) 8 is an universal quantity for d greater than two. While it is still controversial whether the results should be applied to some 29 real amorphous systems at low temperature , there have been con- siderable theoretical studiesz‘j"?8 concerning the spectral dimensionality d. The studies are done on fractal systemszS-za. A random fractal network is inhomogeneous on a length scale less than a certain length and homogeneous at a scale larger than the length. Therefore it is an ideal system to observe the cross-over30. Of course there are many ran- dom fractal networks". The fractal network chosen in this thesis, as well as in many other studies, is a percolation network. The length scale on which the cross-over happens is now simply the percolation cor- relatitn1.length. Many other methods can also be used to obtain a. The 31 and the random walker method32 are the ones transfer matrix method used most frequently. It should be pointed out that, when considering the density of states, the random walker (RW) method is an indirect one to obtain 21'. In dilute systems, the RW method gives results as good as any other method. But this method is not useful for the model33 dis- cussed in this thesis. In part III of this thesis, we study the low frequency vibrational density of states of a superconducting-normal networks. In this model the density of states, in the low frequency limit, has the same form as (3) but with a different spectral dimen- tionality a due to the different scaling relations. These scaling relations are supported by the computer simulations. 11 REFERENCES 1. R. Landauer in AIP Conference Proceedings No.90, Ed. by J. Garland and D. Tanner (1978). 2. G. Deutscher, A. Kapitulnik and M. Rappaport, Ann. Israel Phys. Soc. 5, 207 (1983). 3. L. K. H. van Beek in Progress in Dielectrics 7, Ed. by J. B. Birks (1967). A. .L C. Garland and D. B. Tanner, Electrical Transport and Optical Properties of Inhomogeneous Media, AIP Conference Proceedings No 90 (1978). 5. R. Landauer, J. Appl. Phys. g}, 779 (1952). 6. Z. Hashin and S. Shtrikman, Phys. Rev. 139, 129 (1963). G. W. Milton, Appl. Phys. Lett. 31, 300 (1980). D. J. Bergman, Annals of Phys. 138, 78 (1982). 9. S. Torquato, J. Appl. Phys. 58(10), 3790 (1985). 10. R. Zallen, The Physics of Amorphous Solids, (John Wiley & Sons, New York, 1983). 11. S. Kirkpatrick in Ill-Condensed Matter-Les Houches 1978, Ed. by R. Balian, R. Maynard, and G. Toulouse,(North-Holland, Amsterdam, (1979). 12. D. Stauffer, Introduction to Percolation Theory, (Taylor and Francis, London, 1985). 13. I. Balberg, Phys. Rev. B 31, 2391 (1988). 1‘1. 15. 16. 17. 18. 19 20. 21. 22. 23. 2‘1. 25. 26 12 B. I. Halperin, S. Feng, and P. N. Sen, Phys. Rev. Lett. 51, 2391,(1985). I. Balberg, Phys. Rev. B Q, 11053 (1985). P. G. de Gennes, J. Phys. (Paris) fl, L-1,(1976); D. J. Stauffer, J. Chem. Soc. Faraday Trans. 11 lg, 1359 (1976). M. F. Thorpe, J. Non-Crys. Sol. 51, 355 (1983); H. He and M. F. Thorpe, Phys. Rev. Lett. 3, 2107 (1985). S. Feng and P. N. Sen, Phys. Rev. Lett. :2, 216 (19811). A. R. Day, R. R. Tremblay and A. -M. S. Tremblay, Phys. Rev. Lett. fl, 2501 (1986). S. Feng, M. F. Thorpe and E. Garboczi, Phys. Rev. B 31, 276 (1985). W. Tang and M. F. Thorpe, Phys. Rev. B 31, 5539 (1988); W. Tang, Ph. D. Thesis, Michigan State Univ (unpublished). C. Kittel, Introduction to Solid State Physics (J. Wiley 81 Sons, Inc. New York, 9th Edition, 1971). R. J. Elliott, J. A. Krumhansl and P. L. Leath, Rev. Mod. Phys. E, 1465 (19711). 3. Alexander and R. Orbach, J. de Physique Letters 33, L625, (1982). R. Ramal and G. Toulouse, J. de Physique Lettres, fl, L13 (1983). Y. Gefen, A. Aharony and S. Alexander, Phys. Rev. Lett. 59, 77 (1983). 27. 28. 29. 30. 31. 32. 33. 13 3. Alexander, C. Laermans, R. Orbach and H.M. Rosenberg, Phys. Rev. B _2_8_, 11615 (1983). A. Aharony, S. Alexander, 0. Entin-Wohlman and R. Orbach, Phys. Rev. B 3_1, 2565 (1985). R. Orbach, Science, Vol. 231, 819 (1986). B. B. Mandelbrot, The Fractal Geometry of Nature, (Freeman, New York, 1983); Ann. Israel Phys. Soc., _5_, 59 (1983); J. Stat. Phys., 33, 895 (19811). M. A. Lemieux, P. Breton, A.-M. S. Tremblay, J. Phys. (Paris) Lett 36 L1 (1985). D. C. Hong, H. E. Stanley, A. Coniglio and A. Bunde, Phys. Rev B 33, 4561-1 (1986). A. R. Day, W. Xia and M. F. Thorpe, Phys. Rev. B 31, 1339 (1988). Part I Continuum Percolation and Effective Macroscopic Conductivity of Random Composite Systems 14 15 I Introduction The universal behaviour of the critical exponents that describe transport quantities such as electrical conductivity, thermal conduc- tivity, the diffusion constant and elastic moduli of a composite system near the percolation threshold can be understood by continuum percola- tion theory.1 Experimental results agree well with the theoretical predictions of such quantities.2 In designing composite materials it is more important to know the overall behaviour of the properties of these materials which are governed by non—universal quantities away from the critical region. The location of the critical point is also non— universal. When the concentration of one of the components (for simplicity, we will consider only two component composite systems) is extremely low, the behaviour of quantities like the electrical conduc- tance can be adequately described by the Clausius-Mossotti equation.3 Inbetween these two extremes, an exact microscopic theory or detailed computer modeling of the transport properties would be very difficult and neither is currently available. Some progress has been made recently in studying two elliptical holes in a homogeneous medium." However, even here there are still unresolved problems associated with overlapping inclusions that have prevented a useful generalisation of the Clausius-Mosotti equations. Computer simulations have also not been possible because even the largest available machines cannot store enough information to meaningfully discretise such continuum systems. Thus a major reasearch tool, that has led to so much of our understanding of 16 the response of discrete lattice systemss, has not been available or ex- ploited yet in continuum systems. The purpose of this paper is to develop a semi-phenomenological description of the behaviour of the transport properties of these systems. Many effective medium theories (EMT) have been developed3'6'11 but are of dubious validity away from the dilute limit, where all these theories agree. In order to ascertain which of these theories are good for all concentrations, we have located the critical concentration of ellipses at percolation. These results are new except for the special case of circles and provide a most stringent test of effective medium theories. We find that there are E reasonable EMT for electrical conductance; all fail to predict the correct critical concentration for circles. On the other hand, we find one such existing theory to be clearly superior and adequate for elastic inclusions. This is the asym— metric reinforced model (also called SCA-A, reinforced problem in 7). section II B of Thorpe and Sen This was originally derived for cir- cular inclusions by Hill, Budiansky, Wu and Berryman 8‘10 using different self-consistant methods. In the circular limit symmetric and asymmetric theories are identical. Those results were generalised by Thorpe and Sen 7 to ellipses for which symmetric and asynmetric theories are no longer identical. Note however that this theory applies to the case when the circular inclusions are infinitely hard so that the elas- tic compliance and not the elastic modulus vanishes at the critical concentration. We would expect that this effective medium theory should also be superior for mixtures where the elliptical inclusions are hard. 17 The SCA-A for 92.135. does not give a good value for the critical concentration. The interpolation formula we develop incorporates the behaviour at: the two extremes (i.e. low concentrations and near the percolation threshold) for a system containing randomly distributed insulating el- liptical laminae (i.e. holes) embedded in the uniform background of a conducting matrix. The layout of this paper is at follows. In section II we discuss the geometric aspects of continuum percolation for a system containing random elliptical laminae in 20 and spheroids in 3D and in section III we present our computer simulations of the percolation thresholds in 20 and compare with the results of previous work. In section IV we use our results to critique effective medium theories by examining their predic- tions for the percolation concentrations. In section V an interpolation formula for the electrical conductance is developed and in section VI a comparison is made between the effective conductivity predicted from our interpolation formulae and that of experiments. II Continuum Percolation (1) Continuum Percolation of Elliptical Laminea in 2D. In the continuum percolation problem, as well as in percolation on a lattice, an important quantity to describe the onset of percolation is 12 13 the percolation threshold pc. There exist extensive studies in the literature of pc for percolation on various lattices, where pc is the 18 fraction of bonds or sites remaining, depending on the type of percola- tion being studied. In the continuum percolation problem, pc is defined to be the fractional area occupied by one phase which, in our case is the area remaining after the elliptical holes are removed. Fig.1.1 shows an example of the system under study and the area covered by el- lipses has fractional area 1 - p (at the percolation threshold p = pc). Imagine that a constant voltage is applied across a conducting sheet and randomly oriented elliptical holes with random centers are punched out. As more and more material is removed, electric current flow through the sheet is restricted and vanishes at pc. We are interested in how pc changes as the geometry (i.e. aspect ratio) of one phase changes”’15 or more presicely how pc changes as the inclusions change from circles to needles. 20 We use the aspect ratio b/a to describe the asymmetry of the ellipse where a and b (with a > b) are the major and minor semi-axes respectively. Note that the eccentricity of the ellipse is given by e = [ 1—(h/a)2 11/2. In the following discussion, identical, but randomly centered holes each with area A are removed from a two dimensional L x L sheet. At hole concentration n per unit area, and remainipg area fraction p, if we increase the hole density, then the area that is still available to be removed is pL2. Therefore the additional area removed by changing the hole density from n to n + dn is pLzAdn. On the other hand, the area remaining is reduced to pL2-(p+dp)L2, so pL2-(p+dp)L2= pL2Adn i.e. dp/p -Adn so that p exp(-An) (1) where we note that p = 1, when n=0 and p = 0 when n = CD. This formula has been used previously for circles16’17 sufficient randomness is present.18 This is because the repeated random but is true for all shapes if placement of an additional object or hole serves as a measure of the remaining area. We are of course always thinking of the thermodynamic limit when the system size is very large. Eq. (1) can be generalized to three dimensions where A is replaced by the volume of each individal 21 hole and n is the hole density per unit volume (this can be visualized as Swiss Cheese). In the two dimensional case we are studying, A=11ab and at the percolation threshold nznc, therefore p0: exp(-nabnc) (2) Eq. (2) allows an immediate determination of the percolation threshold pc for a given a and b once no is known or vice versa. This equation is very convenient to use in practice as it only involves countin ; no area evaluation is involved. For circles we find from our simulations that p0: exp(-na2nc) = 0.33 t 0.02 (3) where a is circle radius and ne is density of circles per unit area at percolation. We have given generous error bars on (3). Our result (3) 19 We also notice that a 19 agrees well with other results for circles. carefu1 finite size scaling study gives better results, but this would need huge amounts of computer time if it was to be done for all aspect ratios. Our purpose here is to look at the general trend of how pc changes with b/a. The percolation threshold for ellipses aligned in one direction but with random centers is the same as that for random circles. The reason is as follows. A conformal transformation can al- ways be performed in one direction on the aligned elliptical system in order to obtain the random circular case. Obviously the change in the area A, due to the conformal transformation, exactly cancels the change 22 in the ellipse density per unit area. Therefore there is no change in p=exp(-An) and the percolation threshold is the same as that for the circular case. This is checked by computer simulations in this work. In 2D, the background ceases to percolate when the inclusions percolate. This is because there is no way aggpgg the infinite cluster. Therefore, there is a single percolation concentration. In higher dimensions, this is obviously not the case and there are two separate percolation concentrations for the inclusions and the backgrounds. In continuum percolation involving identical objects, it is useful to introduce the average excluded area denoted as (aex>.20’21 For given relative orientations of two identical objects the excluded area is defined as the area that if the center of one is outside it, the two ob- jects have no overlap at all. Average means over all allowed relative orientations. The excluded area at percolation is defined as (A >=n (H) ex C ex Although the mean coordination number and the critical area or volume fraction are essentially invariant22 in bond and site lattice percola- tion respectively, (Au) is not quite such a quasi-universal invariant 12 quantity and it has a small range. Our results will be discussed in detail in the next section, but we see from Table 1 that for randomly oriented ellipses, 3.A < (A > < “.5 (5a) 23 For ellipses that can only lie in two directions, we see from Table 2 that, 2.8 < < 11.11 (5h) Taking account of the error bars noted in the table captions, both these sets of results for (Aex) are probably monotonic in the aspect ratio. The excluded area of two identical ellipses can be defined as (aex) = Anabk (6) where k is a geometric factor that is chosen as above so that k = 1 for circles. For randomly centered and oriented ellipses,23 k: 1/2 + s2/812ao (7) where s is the perimeter of an ellipse (this involves an elliptical 11r- tegral which can be evaluated numerically). For randomly centered ellipses that can only lie in the two principal directions the k factor is different from (7) and not available in a closed form for general b/a. For parallel ellipses k = 1 but must be computed (using for ex- ample the contact function described in the next section) for ellipses at right angles. These two results are then averaged. The values of k for both these cases have been calculated and are given in Tables 1 and 2 for various aspect ratios. 24 From Eqs. (2), (4) and (6), we notice that p0: exp(- (Aex>/4k) (8) 16-18 In which reduces to p0: exp(- (Aex>/4) for circles when k = 1. the needle limit where b/a is small, and the ellipses are randomly oriented, using (7) we have s = 4a and k : 2a/(12b) so that from Eq. (8) pc = exp[-3.4/(4k)] = exp(-O.42512b/a) z 1 - 4.2 b/a (9) where we have used (Aex> z 3.4 from Table 1. We note that in this limit, the result (9) is independent of the precise shape of the needles. For example they can be elliptical or rectangular. As b/a be- comes very small, only a few needles are needed to cross the sample and these have essentially no area so that pc + 1 as given by (1). A similar limit is obtained for needles that can only point horizontam or vertically for which k = a/(21b). Using Eq. (8), with (Aex) ~ 2.8 from Table 2, we find that p0 = exp[-2.8/(4k)] = exp(-1.4nb/a) z 1 - 4.4 b/a (10) 25 which we notice is close to the result for randomly oriented ellipses given 1J1 (9). Indeed because the values of are only known numerically; the error bars are sufficiently large that Eqs. (9) and (10) could be identical. (2) Continuum Percolation of Spheriods in 3D Various concepts, just described above, for two dimensional per- colation can be easily generalized to three dimensional case. As we mentioned earlier, in 2D the background ceases to percolate when the inclusions percolate. This is because there is no way around the in- finite cluster and therefore there is a single percolation concentration. In higher dimensions (e.g. three dimensions), this is obviously not the case and there are two separate percolation concentra- tions for the inclusions and the background. For percolation involving identical insulating spheres the percolation thresholds are p0 = .31 for inclusions” to percolate and pa = .968 for background cease to percolate27. Here, of course, pc refers to volume fraction. It is also known that the excluded volume (Vex) in 3D is bounded by 12 0.7 < < 2.8 (So) .. 33. where 2.8 refers to spherical inclusions and 0.7 to very long thin rods. 26 Now we consider the following two cases of inclusion percolation for spheroids of revolution. [see Fig. 1.2 ] Prolate (a > b, b=c) Oblate (a > b, a=c) 2b (b) ,Fig. 1.2 (a) A prolate with a>b and b=c; (b) A oblate with a7b and a=c. 27 We can also make predictions on the dependence of the percolation threshold on b/a for small b/a by using the fact that (Vex) = nc using nc found from the computer simulation and Eqs. (4) and (6) for both randomly oriented and two direction oriented ellipses. The error bar in our computer simulation in determining "c is about 1 51 ; therefore the error bar in Mex) is about t 0.2. Tables 1 and 2 list pc, n k, and (Aex> for various 0’ 31 aspect ratios b/a for the two cases and they show that in both cases decreases very slowly as b/a decreases. In the following discussion, we only consider quantities for the randomly oriented case. The case of only two orientations would give essentially indistinguishable results. Note that although nc and hence pc are virtually indistinguishable for a fixed aspect ratio b/a in the two cases, the quantities k and (Aex) are different as can be seen by comparing tables 1 and 2. In Fig. 1.4 we plot f =11abnc and f = 1 -f = 1 - exp(tabnc). The 1 2 quantity f is the total area in the ellipses for a sample of unit area, 1 not allowing for the overlap effects, whereas f2 is less than f1 be— cause overlap effects are included. Pc 32 0.2 F- .3 0 1 l 1 i 0 0.2 0.4 0.6 0.8 1.0 b/ a Fig. 1.3 Percolation threshold pc for various aspect ratios b/a. Squares are for randomly oriented ellipses while triangles are for ver- tically and horizontally oriented ellipses. Every point is averaged over 25 - 30 samples each containing ~ 2000 ellipses. The solid curve is the interpolation formula (28) for pc. The dashed curve is pl,1niich gives the initial slope, from Eq. (22). 33 Table 1.1 bl‘ pc nc k (“ex 1.0000 0.33 2.8 1.000 11.11 0.9000 . 0.33 2.8 1.002 11.11 0.8000 0.33 2.8 1.009 11.11 0.1000 0.311 2.8 1.021 1 5 0.6000 0.35 2.1 1.050 11.5 0.5000 0.31 2.5 1.093 11.3 0.1000 0.111 . 2.3 1.111 11.2 0.3333 0.1111 2.10 1.2511 11.1 0.2500 0.50 1.16 1.1132 11.0 0.2000 0.511 1.51 1.618 11.0 0.1500 0.62 1.22 1.931 11.0 0.1000 0.10 0.90 2.592 3.1 0.0661 0.18 0.62 3.589 3.1 0.0500 0.83 0.119 g 1.592 3.5 0.0100 0.86 0.110 5.599 3.5 0.0333 0.88 0.311 6.609 3.5 0.0250 0.91 0.26 8.629 3.5 0.0125 0.919 0.133 16.111 3.5 0.0050 0.919 0.0511 111.06 3.5 0.0025 0.990 0.021 81.12 3.11 Table 1.1 Values of pc, no, k, and (Rex) are listed for randomly oriented ellipses for various aspect ratios b/a. The value of k is ob- tained by evaluating the perimeter s of the ellipse from an elliptic integral and using the formula in the text. The value or 11C is obtained from the simulation with the area of the ellipses fixed at I/8 and then (hex) is obtained from formulas (4) and (6). The error bar in (Rex) is t O .2. 34 Table 1.2 b/a pc nc k (Au 1.0000 0.33 2.8 1.000 9.» 0.9000 0.33 _ 2.8 1.002 11.11 08mm , 033 28 um» mu 0.1000 0.39 2.1 1.029 9.3 0.6000' 0.36 2.6 1.049 4.3 0.5000 0.37 2.5 1.091 1.3 0.!000 0.91 2.3 1.162 9.2 0.3333 0.45 2.06 1.237 4.0 0.2500 0.50 1.80 1.391 3.9 0.2000 0.51 1.59 1.598 3.9 0.1500 0.62 1.21 1.812 3.5 0.1000 0.68 0.97 2.342 3.6 0.0667 0.78 0.65 3.137 3.2 0.0500 0.82 0.51 3.933 3.2 0.0100 0.85 0.42 1.129 3.1 0.0333 0.87 0.37 5.525 3.2 0.0250 0.90 0.28 1.111 3.1 0.0125 0.991 0.115 13.18 3.1 0.0050 0.916 0.061 32.58 3.1 0.0025 0.989 ' 0.028 64.41 2.8 Table 1.2 Values of pc, no, 11, and (Ae > are listed for ellipses with x two allowed orientations. The values of k is obtain by evaluating the‘ excluded area (aex) numerically and then using Eq. (6). The values of neare from simulation with the area of the ellipses fixed at I/8 and (hex) are obtained by using Eq. (9). The error bar in is t 0.2. 35 1-2 1 1 l 1 A A A A A 1.0 '_ A —1 A 0.8 - A ~ A .. q: 0 6 - + + A++ 0.4 )- A+ -- + 0.2 Té " 1.0 b/a Fig. 1.4 The quantities f2: 1 - exp(uabnc) and f1: 11abnc are plotted against the aspect ratio b/a for the case of randomly oriented ellipses. Crosses are for f2 and triangles are for f1. 36 It can be seen that for small b/a these two quantities are the same as the overlap area for needles is negligibly small. In the circle limit f1 = 1.09 1 0.02; that is the area in the circles at percolation, before they are thrown down, is greater than unity. Note that if f2 is ex- 1 and the corrections for r body overlap are given by the coefficient of the n: panded in powers of the density nc, the first term is f term. IV Critique of Effective Medium Theories There are extensive discussions in the literature on EMT for 3 and elastic moduli11 dielectric constants of composite materials with circular or spherical inclusions. A strong assumption is always re- quired, in deriving these approximations, that the inclusion concentration is sufficiently low that the overlap of inclusions can be neglected. However these approximations are often used over the whole concentration range where they are of dubious validity. In order to judge how good various EMT results are when applied to completely perme- able objects, we see how close their predictions of pc are compared to our exact (numerical) results. Physical properties, like the conduc- tivity and all elastic moduli, should vanish at pc when holes are punched in the medium. Similarly the resistance and all the elastic compliances should vanish when infinitely hard inclusions are present in the medium. Infinitely hard means superconducting in the electrical case and infinitely rigid or undeformable in the elastic case. All these pc should be the same as it is a geometrical property of the 37 material. However different EMT give very different estimates for pc. These various EMT predictions for pc can be used as a figure of merit, when compared to our exact results, to judge how good the EMT is away from the dilute limit. In what follows we will examine two versions of EMT for each physical property. Depending on whether we treat the in- clusion and background symmetrically or asynmetrically, two versions 3 (i.e. synmetric or asymetric) of EMT can be derived. Thus we have _8_ cases to consider, electrical or elastic, symmetric or asymmetric with inclusions that are either holes (Swiss cheese model) or hard inclusions (reinforced model). Sen, Thorpe and Milton6 have summarised these results for the electrical case. These results can also be obtained from ref. 3. The critical pc for the dilute (Swiss cheese) case are s pc-1/2 (14) p::(82+b2)/(a+b)2 (15) where the superscripts s and a refer to the symmetric and asymmetric cases respectively. The results for the reinforced case are identical to (14) and (15). 7 Similar results have been obtained by Thorpe and Sen for the elas- tic case. For the dilute (Swiss cheese) model, all the elastic moduli vanish at 2 1/2 -1 0: =2 1 1 + [2(a+b)2/(82+ b )1 1 (16) 38 pg: [1 + ab/(a2+b2) 1" (17) For the reinforced model, all the elastic compliances vanish at pg: 1 - 2 { 1 + [2(a+b)2/ (a2+b2)]1/2}-1 (18) (1-p2)-1:2{1+(1-0)(a+b)2/[2ab(1+0)])/(3-0) (19a) (1—p:)“=[(a+h)2/[at(3—o)1 + 1/[1- ab(1+o)/(a+b)2]}/2 (19b) where p: is f0und from Eqs. (19a) and (19b) by eliminating 0, the value of Poisson's ratio at the critical point. If these were exact theories, all the results (14) - (19) would be identical. Note that there is no difference between the symmetric and asymmetric cases in the circle limit for all these results. The above results are shown in Fig. 1.5 as a function of the aspect ratio and we can see that only curve 6 which is the result (19) is reasonable. Indeed all results except for the reinf0rced elastic model fail to get even the circle limit correct. These two approximations for the reinforced elastic model [ Eqs. (18) and (19)] give pc = 1/3 which is the correct result for circles within numerical error as can be seen from Tables 1 and 2. These results show that EMT is inademiate, when strong disorder is present, except in the one special case. In other cases we believe a better procedure is to develop interpolation formulas. Pc 39 1-3 ‘ 1 1 1 1 ‘~s ‘\~ 1!, x 4 ‘8 \.~ 3 0.13--‘i'.I , \~~J -— .\.N ~~\--.~-- n ‘ ‘1‘.==-—-— 0.5 _ 1:1 121 1 ._ .—————— B.— ...——_T—": B --1 0.4 -— n .. _ (5 ___;__ " -35 e; z; ; ,/—/‘/.f’ .- 0.2 7”,” 5 l l l l 00 0.2 0.4 0.5 0.8 1.0 b/a Fig. 1.5 Percolation thresholds pc predicted from various effective medium theories and our computer simulations ownn Fig.1.3. The curves are marked 1 for Eqn. (14), 2 for Eqn. (15), 3 for Eqn. (16), 4 for Eqn. (17), 5 Eqn. (18) and 6 for Eqn. (19). The squares indicate the exact percolation thresholds from computer simulations taken from the results for randomly oriented ellipses in Fig. 1.3. 40 We note that Eqs. (19) could be used as a useful parametric ap- proximation to pc when required. It gives p0 = 1/3 (compared with 0.33 in Tables 1 and 2) in the circle limit and pc 5 1 - 16/3 (b/a) (20) in the needle limit. This should be compared with Eqs. (9) and (10) V Interpolation Formulae (1) Semi-phenomenological Formulation As we mentioned in introduction, the Clausius—Mossotti equation for the conductance of a two phase system is exact when one phase has a very low concentration. All attempts to extend these equations beyond this region are rather uncontrolled and many versions exist in the literature. As we discussed in the previous section, all are unsatis— factory for the electrical case. We therefore develop a simple interpolation formula that gets all the known limits for the dilute (Swiss cheese) model correct. We believe that this should be of con- siderable utilitarian use. Similar fermulas can be written down for all other cases. For a M number of holes in a material with conductance 20 the effective conductance E is given by Z = 20[1 - (i-p)/(1-pl)] (21) 41 where DI = (a2 + b2)/(a + b)2 (22) The quantity pI is where the initial slope for a small number of 3,6 defects would eventually cross the E = 0 axis when extrapolated and is a convenient way to express the initial slope. The relation II‘~ (p- pc)t holds only in the small critical region around pc. Our interpolation formula is designed to link these two limits by assuming the conductivity has the following form 2 = 20(1.0 + co + 8oz)A (23) where c = 1 - p and a, B and l are constants to be determined from the following, E~(p-p)t asp+p (24) c c E = £0[1 - c/(1 - pl) + O(c2)} as c + 0 (25) and to is the conductivity of the sample without any inclusions (c = 0). Of course one would like to include higher order terms in c in Eq. (23). but since we have no more information other than (24) and (25) it is not possible to do better. After some simple algebra we find, 2/20= |1 - c/[t(1-p1)] - 02[t(1-p1) — (1—pc)]/[t(1-pl)(1-pc)2]}t (26) 42 It is rather inconvenient to use the expressions (19) for pc and so we make a simpler approximation (27) for pc that is correct in the two limits b/a : 1, when p0 = 1/3 and b/a small, when pc z 1 - 9/2(b/a) pc : (1 + 4y)/(19 + 4y) (27) where y = b/a + a/b is symmetric in a H b. The result (27) is vir- tually indistinguishable from the computer simulations in Fig. 1.3 and is actually superior to (19) as can be seen by comparing Fig. 1.3 and Fig. 1.5. Of course there is no basis for (27) except that it is cor- rect in the two limits and fits the simulation data for all aspect ratios. By taking t = 1.3,1 and using Eq. (22) for p1: y/(2 + y) and Eq. (27) for pc, we can determine the effective conductance E, (z/xo)‘/t = 1 - c(2+y)/(2t) + c2(19 + 4y)[9(2+y) -(19 + 4y)t]/(324t) (28) which we recommend for use in practice (with t = 1.3) as it reproduces all the known results (i.e. the value of 2 for the pure system, with p = 1; the initial slope for small 1 - p; the value of p = pc where 2 vanishes with critical exponent t) to within numerical accuracy. Note that the term in c2 is always small and positive. This is because pI is always larger than pc for all aspect ratios (1 i pI/pc .<_ 1.5) as can be 43 seen from Fig. 1.5. In Fig. 1.6, 2/20 is plotted against p = 1 - c for various aspect ratios b/a and shown as the solid lines. Also in Fig. 1.6, Z/Eois plotted against 0 but with t = 1.0 for the same aspect ratios and shown as the dashed lines. The two plots are very similar and only differ a little in the critical region. Clearly the EMT described in the previous section would give very different results as the pc are so different. Because of the equivalence of the problems, the interpolation for- mula (28) can be used for, the electrical conductivity of sheets containing holes, the thermal conductivity of sheets containing holes or the dielectric constant of a medium with holes. In all cases, p = 1 - c is the fraction of material remaining after the holes have been punched and y = b/a + a/b where b/a is the ratio of the minor to major axis of the ellipses. If the inclusions are superconducting, rather than insulating (i.e. holes) then the result (28) still holds if we replace 13./II on the left 0 hand side with R/RO where R is the resistance of the sample and lb is the resistance when there are no inclusions (c = 0). These two problems 25,26 map on to one another and are exactly equivalent. Note that pI and pe given in Eqs. (22) and (27) and the critical exponents are the same.1 z/zO 1111 1.0 0.8 - 0.6 - 0.4 - 0.2 - 00.2 Fig. 1.6 Electrical conductance from the interpolation formula (28) for various aspect ratios as indicated. The solid curves are for t - 1.3 and the dashed curves are for t = 1.0. 45 Finally we note that the interpolation formula (28) has two inter- esting limiting forms. Using the limiting forms for pI and pc, we find that for circles, (2/20)‘/t : 1 - 2o/t + 302(4 - 3t)/(4t) (29) and for needles, 1/t 2 (2/20) = 1 - nflLZ/(8t) + n 02Lu(9 - 4t)/(1296t) (30) where we have used Eq. (1) with A = 11ab and put c = 1 - p = 1 - exp(- 11abn) 5: n11ab and the length of the needles is L = 2a. The result (30) is independant of the width b of the needles as would be expected. Here n is the number of inclusions per unit area. A generalization of above interpolation formulae to 3D case is straight forward. One only needs to know pI (can be easily found from one defect problem), pc (from com- puter simulations), and t or s (superconducting diverging exponent near pc). However p1, pc, and l (i.e. t or s) in the interpolation formulae for conductor—insulator (c-i) system and superconductor-conductor (s-c) system are different now because there is no similar duality, which ex- ists in 2D, exists in 30. In fact pc's for c-i and s-c systems are 0.968 and 0.31 respectively. 46 (2) Comparision with Experiments In the previous section we have discussed the formulation of the interpolation formulas. We think these formulas are useful in the predictions of overall electric conductivity for practical applications. The reason is as follows. The actual conductivity of a random composite system, such as the one studied here, in general is a smooth monotonic decreasing function of p. If one can develop a formula which is a smooth nomotonic decreasing function with correct limits in the low and critical p, then the overall behavior predicted by the formula cannot be too far away from the real situation. In our case, as we will see later, the agreement between interpolation formulas and experiments are 2 ’28done on the rather good. There have been two recent experiments conductor-insulator system for which our interpolation formulae have been developed. Therefore these experiments are a direct test of how good these interpolation formulae are in predicting the macroscopic con- ductivity in the whole concentration range. There are no other predictions, to our knowledge, made by any other theory fer overall con- ductivity in the whole concentration range in conductor-insulator systems. In the experiment of ref. 2, which we refer to as experiment I, about 3300 circular holes were drilled on each of two steel and two molybdenum sheets(these materials are used instead of copper or aluminum to avoid deformation of holes during the drill). The size of the sheets are 16cm X 16cm and the radius of the holds is 0.32cm. The thickness of the sheets are 0.13mm, 0.25m and 0.38mm respectively. Effects due to the thickness of samples are negligible after observing no change among 117 the samples with different thickness. Also finite size effect only ex— ists when percolation correlation length is comparable to the system size L(in the experiment this happens when (p-Dc)/pc ~ 21). The macro- scopic conductivity is monitored while holes are drilled. The area fraction is estimated by using the equation p = exp(-nna2) which is explained in section II of this thesis. In the experiment of ref. 28, which we refer to as experiment 11, about 650 ~ 700 circular holes are cut from an metalized myler foil. The size of system is 10cm X 10cm and the radius of hole is 0.4cm. Effects due to sample thickness can be neglected. The area fraction is measured by weighing the cut fragments on an accurate scale. Also the macroscopic conductivity is monitored as holes are cut. The normalized conductivities from inter— polation formula are plotted in Fig. 1.7 , using Eq. (28) in section V for circular holds, against area fraction p. We use t=1.3 in the solid curve and t=1.0 in the dashed curve. The squares and triangles are ex- perimental results from experiment I and experiment 11 respectively. We can see that apart from some fluctuation the general features of the two sets of experimental data are very close to the two lines in a wide range of concentrations. The data from experiment 11 have serious deviations near the percolation threshold. Experiments near the per- colation threshold need to be done with extreme care because of critical and brittle features of the system. Since the percolation threshold is 48 wrong in experiment 11, we attribute the deviation due to some unex- pected failure in the experiment. We have sent our comments to authors in ref. 28. Also we notice that number of holes in experiment 11 is much less than that in experiment 1. While a serious deviation is ob- served in experiment 11 near the percolation threshold the rest of data points are not effected because the experiment is performed step by step. The agreement between experiments and our interpolation formulae is very encouraging and it gives us enough confidence in the interpola- tion formulae. Therefore we believe, with confidence, our interpolation formulae shold be very good in predicting the overall features of con- ductivity in two phase composite systems. Experiments on continuum percolation with randomly distributed elliptical holes are recomended for a fUrther check. Conductance 49 1..0 [ l I t-itC) 0.3.. ——-—-t-L3 a [J Exp 1 J.C. Lobb A Exp 2 E.N. Martinez 0.5 "' .—1 0.4 - - 0.2 - ._ 0 0.2 0 4 0.6 0 8 1 Fig. 1.7 A comparison between interpolation formulae (lines) and Experiments (squares and triangles). .0 50 VI Sumary of Results Our main result has been the numerical determination of the per- colation concentration of randomly centered and oriented ellipses. This has been done by using a contact function to determine if neighboring ellipses overlap, and constructing a connectivity table. 11.13 not necessary to measure any overlap areas to find the areas at percolation; it is sufficient to merely 992p; the number of ellipses. We have also determined the percolation concentration of ellipses when the axes are contstrained to lie in only tw_o_ Cartesian directions. The results are indistinguishable, within our numerical accuracy, from the previous case where all orientations are allowed. We have used these results to critique various effective medium theories that have been developed for the electrical and elastic responses of sheets containing elliptical inclusions. Only one of these approximations is found to give a reasonable percolation threshold while all the others fail to describe the electrical conductivity or elastic properties near the critical point. We have shown that the percolation concentration is described well by the formula pa = (1 + 4y)/(19 + 4y) where y = b/a + a/b. Here pc is the amount of material remaining and b/a is the aspect ratio of the ellipses. We have also developed a simple interpolation formula for the electrical conductance that is correct both for a few inclusions and near percolation. The agreements between interpolation formulas and ex- periments lead us to believe that this kind of formula is superior to effective medium theories and may have useful practical applications. 51 VI I REFERENCES 1. 10. 11 12. 13. 14. B. I. Halperin, S. Feng, and P. N. Sen, Phys. Rev. Lett. 88, 2391,(1985). C. J. Lobb and M. G. Forrester in preprint (1986). R. Landauer in AIP Conference Proceedings No.40, Ed. by J. Garland and D. Tanner (1978). R. I. Cukier, J. Phys. Chem. 82, 246 (1985) S. Kirkpatrick in Ill-Condensed Matter-Les Houches 1978, Ed. by R. Balian, R. Maynard, and G. Toulouse,(North-Holland, Amsterdam, (1979). P. N. Sen, M. F. Thorpe, and G. Milton, unpublished notes, (1985). M. F. Thorpe and P. N. Sen, J. Acoust. Sco. Am., 11, 1674, (1985). R. Hill, J. Mech. Phys. Solids 13, 213 (1965); B.1Budiansky5 J. Mech. Phys. Solids 13, 223 (1965). T. T. Wu, Int. J. Solids Struct. 3, 1 (1966). J. Berryman, J. Acoust. 800. Am. 88, 1820 (1980) and references, therein. J. P. Watt, 0. F. Davies, and R. J. O'Connell, Rev. Geophys. and Space Phys. lfl_(4), 541 (1976). I. Balberg, Phys. Rev. B 31, 4053 (1985). D. Stauffer, Introduction to Percolation Theory, (Taylor and Francis, London, 1985). L. N. Smith and C. J. Lobb, Phys. Rev. B 291 3653 (1979). 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 52 F. Carmona and A El Amarti in preprint (1987). A. S. Skal and B. I. Shklovskii, Fiz. Tekh. Polyprovdn. 1, 1589, (1973) [Sov. Phys. Semicond. 1, 1058 (1973)]. V. K. S. Shante and S. Kirkpatrick, Adv. Phys. 88, 325, (1971). G. E. Pike and C. H. Seager, Phys. Rev. B 1Q, 1421 (1974). C. Mark, Proc. Camb. Phil. Soc. 89, 581 (1954). . We list here the results for pc for circles with references. pc=0.325 (ref. 17); Pa = 0.32 s 0.02 (ref. 18); p3 = 0.33 1+ 0.02, D. H. Fremlin, J. Phys. (Paris) 31, 813 (1976); pc 0.324 t 0.002, E. T. Gawlinski and H. E. Stanley, J. Phys. A 11, L129, (1981); pa = 0.322 1 0.005, T. Vicsek and J. Kertesz, J. Phys. A, 18, L31, (1981). I. Balberg, C. H. Anderson, S. Alexander, and N. Wagner, Phys. Rev. B 88, 3933 (1984). L. Onsager, Ann. N. Y. Acad. Sci. 81, 627 (1949). R. Zallen, The Physics of Amorphous Solids, (John Wiley 8; Sons, New York, 1983). A. S. Roach, Theory of Random ClumpingJ (Spottiswoode and Ballantyne, London, 1968). J. V. Baron, J. Chem. Phys. 88, 4730(1972). J. 8. Keller, J. Math. Phys. §J 548 (1964). K. S. Mendelson, J. Appl. Phys. 88, 917 (1974). W.'I.Enam, A.IL Kerstein and J. J. Rehr, Phys. Rev. Lett.,_5_2_, 1516 (1984). 53 28. J. Sofo, J. Lorenzana and E. N. Martinaez in reprint (1987). Part II Elastic Percolation and Tricritical Point of Stretched Spring Model on Honeycomb Lattice 54 55 I Introduction Percolation of elastic networks has been an interesting subject for the past decade. de Gennes1 was the first to introduce the isotropic force model for the elastic modulus of a gel. The model is defined on a lattice and described by elastic energy _1 _2 V - 2 B (01 ”j) giJ (1) where 8.13 the spring constant and u is the small displacement of site i i. 841::1 when a spring between site i and site j is present with probability p and 81.1 = 0 if a spring is missing with probability 1-p. de Gennes pointed out that the elastic percolation problem of this model is identical to the conductivity percolation problem which has been studied in great detail.2’3 This is reflected by the fact1 that the two problems can be mapped into each other and therefore (i) the percolation thresholds of the elastic and conductivity percolation problems are the same; (ii) the critical exponents describing the vanishing bulk moduli and conductivity near the percolation threshold are identical. Feng and Sen“ proposed the central force model described by elastic energy _ 1 . _ + .“ 2 V - 2 <§J> a {(ui uJ) rij' giJ (2) 56 A where a is the spring constant and rij is the unit vector from site i to site j. Feng and Sen further pointed out that the elastic percolation described by potential (2) is different from that of (1). Again this is reflected by the fact that the percolation thresholds and critical ex- ponents are different for (1) and (2). In fact they belong to different class of the problemsu’s. More specifically, in two dimensions the per- colation thresholds, according to Thorpe's constraint counting method6, are for isotropic force model (3) NIN for central force model. (4) NI: cen Here 2 is the number of nearest neighbors of a lattice. The two models are special cases of the more general Born model7 1 = g (2 {01(81- 5J1-21J12 g,J + e (61-8,)231J } (5) 1j> Feng and Sen also»studied the elastic percolation of Born model" and concluded that even with a small nozero B the elastic percolation cross- overs from the central force model to the isotropic force model. For example, the percolation threshold pcen changes to p0 of the isotropic force model (see Fig. 2.1). 57 1112 ’1) 010L 0 47 .d/ [O 11118L /; ,3/ ’/ ’0 K 006*- lo/ (1.11 = 0.1 “I; \/§/ 004) ,’/ O x / \v - 0 .° /’ )- I 002 ’0; O; ’0’ / ’0’ 0 -°"°’ J 1 oz 04 ‘ 06 06 10 Fig. 2.1 Bulk moduli vs. p for central force model (B = 0, solid line) and Born model (8/0 = 0.1, dash line) for 20 triangular lattice. A strong crossover is shown in latter case. Fig. is from ref. 4. 58 Recently Tang and Thorpe introduced the stretched spring model8 with elastic energy 2 v : X Kuui - L ) (6) <1j> J 0 where K is the spring constant between site i and site j. K 's have 13 1J values of Kij = 1 with probability p and KiJ = 0 with probability 1-p. L0 is the spring's natural length and lij is the distance between site i and site j. Using the relation 1.=|R.-§+G.-Gl, (7) where R1 and R1 are the equilibrium positions of site i and site j, potential (6) can be written as _ 1 _ 2 _ . .n v - 2 (§J>K1J (LiJ L0) + (§J)KiJ(LiJ LO)(uiJ r11) + g { KiJ {(1-L0/L1J) 6112 + LO/L1J(aij.;ij121 + 0(613) . (a) <1j> In (8) L11 = 1R1 - R11 and uij = u1 - uJ. The higher order terms are neglected. Comparing (8) with (1) and (2) one can see that the isotropic fOrce model and the central force model are also special cases of the stretched spring model at L /L1 = 0 and L /L1 = 1 . Therefore J 0 J continuously, the stretched spring model can serve as 0 by changing LO/LiJ 59 a bridging model between the isotropic force model and the central force model. One will also notice that the central force model and stretched spring model described by (2) and (6) are rotationally invariant while the isotropic force model described by (1) is not. Most computer simulations have been done on triangular lattice in the stretched region where O S LO/L1J still in progress. The major conclusion of Tang and Thorpe is that the S 1 while simulations in the compressed region are stretched spring force model can serve as a bridging model between the isotropic force and central force model and the percolation threshold pc changes with L /L. 0 In this paper we study the elastic percolation of the stretched spring model on honeycomb lattice. The motivation of this study is as follows: (i) In honeycomb lattice the number of nearest neighbors (z) is three so the constraint counting method or effective medium theory6’9 will give p , which is meaningless, for the central force model. cen= 3 However, what p = g tells us is that in the central force model the cen honeycomb lattice is unstable against dilution. In other words, if a spring is removed and the network is compressed then it will collapse or have zero elastic bulk modulus. In triangular lattice this is not the case. Therefore by studying the honeycomb lattice we expect to see some new phenomena associated with the instability. (ii) From eq. (8) we notice that there are two quadratic elastic energy terms of u. 1.1 central force components in the stretched spring model. So as LO/L -> 1 which correspond to the isotropic force and the 60 a small isotropic force can help to stablize the network. In the limit, when the strength of'isotropic force goes to zero, we should be able to obtain the correct prediction of the percolation threshold. (iii) A new effective medium theory must be developed to accom- modate the phenomena associated with the instability mentioned before. 'The layout of this paper is as following. In section II we first discuss the elasticity of stretched spring model for the perfect lattice and then use an one defect effective medium theory to estimate the per- colation threshold as a function of LO/Lij' In section III we present computer simulation results on the honeycomb lattice and the explana- tions of the anomalous behavior of the elastic moduli associated with the instability of central force model on the honeycomb lattice. In section IV we use Landau's phase transition theory to establish a tricritical point observed in the computer simulations. We also develop a new effective medium theory which addresses the behavior of the elas- tic constants and predicts the location of the tricritical point. 11 Elasticity of Pure Honeycomb Lattice And Estimate of the Percolation Thresholds We define the stretched spring model on a perfect honeycomb lattice with elastic energy (6) defined in section I. V : g X K. (1. - L 2 (6) (1]) ij ij 0) 61 Here every bond of the lattice is represented by a Hook spring of strength K and natural length L l is simply the distance between ij 0' ij site i and site j. As stated in section I, the elastic energy can be written as (8) for a small external strain. In the rest of discussion we restrict ourselves in the stretched spring region where O S LO/L S 1. The physical meaning of each of the terms in expansion (8) are dis- cussed in ref. 8. Here we just mention them again for reader's convenience. The first term (8a) is the static elastic energy term because L is the equilibrium dis- 13 tance of lattice sites i and j and no 0 is involved. The second term ii in (8) is v2 = <£1110.iJ - L0)(uiJ-rij) (80) which is the elastic energy due to the tension of the springs because it is linear in displacements. The third term 2 A 2} 1 . . v3 = 2 <§J> KiJ {(1-LO/L1J) 0H + LO/Lij(uijorij) (80) 62 is the quadratic elastic energy term of displacements 01 This term JO determines the elastic moduli of the system. By comparing (80) with Born model (5) one can see that a = Kij (1'L0/Lij) (9a) and B = Kij(L0/Lij) (9b) while a + B = Kij . (90) Notice that LO/LiJ = 1 case corresponds to the central force model while LO/Lij = 0 case to the isotropic force model. For an arbitary LOILij , in general, there is a lot of stress in the system due to the stretched and compressed springs. One may also notice that the full elastic potential (6) is rotationally invariant because it only involves the distances between springs. (1) Elasticity of Pure Honeycomb Lattice. Now we first discuss the elastic moduli of model (8) without any dilution (i.e. perfect lattice case). In this case L = L ,the lattice ii Spacing, and K = K. It can be found in the literature7’10 13 strain energy can always be written as that the (10) 63 where 808 are the elements of the stress tensor and COBYT are the second order elastic constants which are directly associated with the elastic moduli. The quantity £08 is called the strain and defined by E :3;- a,B=x,y A where ua and x8 are the components of 0 =(ux + vy) and F. In two dimen- sions -32 -93 Exx ' 3x ny ' 3y -2! -3! 8xy ' 3x eyy ' 31 Therefore an uniform displacement for a Bravais lattice in two dimen- sions can be written as xexx + yeyx (11a) I: II V = x6 + e . 11b xy y yy ( ) Now let us put a small force on the frame which holds the springs in place of the honeycomb lattice. The lattice will respond to this small force by rearranging the lattice sites so that the elastic energy (8) is at a minimun. The response of a Bravais lattice such as the 64 square lattice or triangular lattice is a small uniform displacement ex- pressed in (11a) and (11b). The honeycomb lattice is not a Bravais lattice but can be regarded as a lattice composed of two interpenetrat- ing triangular lattices which are Bravais lattices. Thus the response on the honeycomb lattice will be two identical small uniform displace- ments such as (11a) and (11b) on the two sublattices and an uniform 11’12 In other words the dis- relative shift of the two sublattices placements between any two sites of lattice A and lattice B [see Fig.2.2] can always be written as u : xe + ye + u' (12a) v : xe + ye + v' (12b) where x and y are the x and y components of the distance between two sites. Qantities u' and v' are the relative shifts of the two sublat- tices in x and y directions and can be determined by minimizing the elastic energy (8). '"\ .1 Fig. 2.2 Honeycomb lattice is decomposed into two triangular sublat- tices A and B. Directions of 81 and 62 are shown. 66 Equation (10) can be rewritten as C + S E + S E + S E + S E + __xxxx (62 0 xx xx yy yy xy xy yx yx 2 xx C + C e e + -§l§l(£2 + £2 ) + C e e . xxyy xx yy 2 xy yx xyyx xy yx (13) 111 (13) we have also considered the symmetries of the honeycomb lattice to reduce the number of independent elastic constants13. By minimizing the elastic energy with respect to u' and v' one can find that [see Appendix I] c _ c _ o2 + 508 + 432 _ u - 33 K xxxx ' 11 ’ 2/37(2B + a) I 2 - n 273 02+ 08 - - - ___fl___ K Cxxyy ' C12 ' 2/3 (28 + a) ' 2 — n 273 C _ E _ 132+ 308 _ g<1 - g) K xyxy ' 44 - 2/3 (28 + a) ' 2 - n 273 _ _ aB _ (1-n)(4-g) x nyyx ‘ Cut ’ 273'728'1 a) ' 2 — n 273 and S = S = 0. (14a) (140) (14c) (14d) (14e) 67 In the second column we have used the conventional notation Cxxxx = C11, Cxxyy = C12, nyxy = C44 and nyyx = C44' In the fourth column we have used a = Kn, B = K(1 - n) and n = LO/L. We first see that (14d) is zero when n = L /L = 1. This is the 0 case where the instability, as mentioned in the introduction, rises. It means that even a perfect honeycomb lattice can not resist a shear distortion in the central force limit. We will see later that this in- stability gives an extra zero frequency mode beside the accoustic modes in the dynamic matrix and leads to the wrong prediction for the percolation threshold pcen = g. We can also see that Equation (15) just reflects the fact that honeycomb system, due to the geometric symetry, is isotropic in any direction.8 The bulk modulus is defined by 1 B : 2(C11+ C12) . (163) Using (14a) and (140) we obtain K _ B : 27g . (10b) Equation (16b) just states that bulk modulus is a constant independent of n. The bulk modulus is, however, unstable against dilution in the 68 central force model limit, that is, the bulk modulus will be zero if one spring is removed when n = L /L = 1. It can also be seen that $10,, - 51111) = T. (17) Equation (17) reflects the fact that the strain energy is rotationally invariant.8 (2) Estimate of the Elastic Percolation Thresholds of Honeycomb Lattice We now discuss the theoretical estimate of elastic percolation thresholds of model (8). The method was first developed by Feng, Thorpe and Garboczi9 to calculate the percolation threshold for the central force model. Tang and Thorpe8 generalized it to the stretched spring model. In ref. 9 authors obtained the percolation threshold for the central fOrce model (2) p =— . (18) 11 For two dimensional honeycomb lattice p = 3. Obviously this is un- cen physical because p must always equal or less than one. This unphysical threshold p = g is associated with the instability of the central cen force model on honeycomb lattice. This is clearly seen by the vanishing shear modulus C44 = 0 in previous discussion. The consequence of this vanishing shear modulus will be discussed. In ref. 8 the estimate of the percolation threshold is obtained by removing a sipgle wing and calculating the energy change AE afterward. A straight line is then 69 used to extrapolate AE to a point in the p axis where the energy is zero (see Fig. 2.3). Of course the elastic energy generally does not decrease linearly against dilution but with some curvature. Therefore the approach is just an estimate of the percolation threshoLd. We call this estimate, denoted by p(B), the extrapolated intercept or the ini- tial slope of the elastic energy and p(B) is different from the actual percolation threshold. However, in some cases, for example in the isotropic force model and the central force model, this theory gives very accurate predictions. One may also wish to calculate the similar extrapolated intercepts ( or initial slopes) from tension (T) and elas- tic modulus (B) curves (see Fig. 2.3) denoted by p(T) and p(B) respectively. The extapolated intercept is often used as a guide to roughly locate the percolation thresholds and it is very helpful in com- puter simulations. 70 1.0 0.8 -— £0 / :g: g // //. 0.5 - ’. // //. 0.4 — ,/// // / I./ / ,///, 0.2 "' /'/ / ,jfijpl/ / / 115;». ./ / 0 '4 / 1 58191182) Fig. 2.3 The extrapolated intercepts for elastic energy (E0), tension (T) and bulk modulus (B). It is clear that thses intercepts overes- timate the actual percolation threshold. 71 A derivation of p(E) using the same method in ref. 8 for elastic energy (8) is given below. P(T) and p(B) can also be obtained in the similar manner after differentiating the elastic energy with respect to L. Notice that in the LO/L = 1 limit (8) corresponds to the central force model. During the course of derivation we will point out the crucial step which may lead to the wrong prediction. The previous derivation was on a Bravais lattice and it is now generalized to the honeycomb lat- tice which is composed of two interpenetrating Bravais sublattices. p( E) of the stretched spring model (8) can be written as [see Appendix II] K .. .. iR-R6 .. 1Eo85 p(E) = 285’ X E Tr{—88(A + 0) + 880 e + 888 e I (19) 8 where sum over 6 extends over nearest neighbors and sum over if is restricted to the first Brillouin zone. In (19) the Ra's are the posi- tion vectors of nearest neighbors once a site is chosen and 8's are the corresponding unit vectors in Ra's directions. Matrices A, B, C and D are defined by 4D + D"1 1 _-1 - A ‘ DAB (‘DAB AA BB BA) B : 0‘1 (-0'10 + 0'1 -1 BB AB AA BB BA) (20) c = 11’1 (-0‘10 + D‘1 —1 AA AA AB BA BB) 72 -1 -1 -1 -1 D — DBA (-D D + D AA AB BA BB) where 2 X 2 Matrices DAA’ DAB’ DBA and DBB are the following. 0 _ 3(o + 28)/2 0 0 _ ’31 '32 AA ’ 0 3(a + 2 B)/2 AB “ -a -a 2 3 h I M 0 _ a1 a2 0 _ 3(a + 28)/2 0 BA ‘ ‘ 3 * BB ‘ 0 3(a + 28)/2 a a 2 3 ik /3a ik J33 ik /3a + ik /3a and a1 = 88 1 + (B + g o)e 2 + (B + g a)e 1 2 ik /3a 1k /3a a2 - a]; e 2 (1 - e 1 ) ik1/3a ik2/3a ik 732 + ik /3a a 1 2 a3 :(a + B)e + (B + H)e + (B + u)e . a In above expressions 0 : K(LO/L), B : K(1-iTyTJ and a1 is the complex .. conjugate of a1. k1 and k2 are the k components along the b1 and 8; directions.[see Fig. 2.2] In the central force limit (LO/L = 1), p(B) = pcen should be given by (19) as a special case. By simply substituting B = 0, however, one will find that matices A, B, C and D do not exist because matrices (- -1 —1 -1 -1 DABDAAI DBBDBA) and (-DAADAB+ DBADBB) can not be inverted. One may 73 recall that in the central force limit the shear modulus is zero. It is this vanishing shear modulus which introduces a zero frequency in the dynamic matrices so that the determinants of the two matrices are zero. . . -1 -1 This is a crucial step. If one ignores the fact that (-DABDAf DBBDBA) -1 1 AADAB+ DBADBB) are no longer invertable and symbolically works out the algebra, one may indeed find pee“: g. So now we understand that and (-D the unphysical prediction is associated with the instability of the central force model on the honeycomb lattice. To obtain the correct pm) in the central force limit one shold allow a small 8, which helps to stablize the network, and take the limit 8 + 0. Following this proc- cedure one will obtain the correct extrapolated intercept p(E) = pcen = 1 in the central force limit. Obviously it is hopeless to write down an explicit function of L /L for (18) because it is so complicated. 0 However, the values of p(E) can be calculated by using a computer. In Fig. 2.4(a) P03) is plotted against LO/L. Clearly p(B) goes to one when LO/L = 1. The extrapolted intercepts p(T) and p(B) are also plotted in Fig. 2.4(a) for honeycomb lattice. The same quantities for triangular lattice are plotted in Fig. 2.4(b). By comparing two figures we find (i) on the triangular lattice and the honeycomb lattice both pm) and p(T) approach 32 and 1 respectively as LO/L + 1. However, the slopes for approaching L /L a 1 are different. (ii) while p(B) ap- 0 proaches two thirds as L /L . 1 on the triangular lattice that of the 0 honeycomb lattice stays almost constant. This means the initial slopes for bulk modulus on honeycomb lattice change very little as LO/L changes from zero to one. This is an unique feature of the honeycomb lattice 74 suiich will explain the sudden discontinuities in B in the next section. In Fig. 2.4(a) we also plot the corresponding quantities (in symboles) from one defect computer simulations. There are some small dis- crepancies between the theoretical curves and the simulation data points. They are expected and will be explained in the following paragraph. We can also develop a more sophisticated effective medium theory which is capable of dealing with more defects and we will discuss it in section IV. 75 -P 1.0 (b) Fig. 2.4 (a) the quantities p(B), p(T) and p(B) vs. LO/L (solid lines). The corresponding quantities from one defect computer simula- tions are also plotted; (b) the same quantities for the triangular lattice from ref.9. 76 (3) Parabolicity of Effective Elastic Potential We have discussed above how to obtain the estimate of the percola- tion thresholds (extapolated intercepts) using a one defect effective medium theory. We can also ask what the effective spring constant each spring will have when one spring is removed from the system. Suppose a spring between site i and site j is removed. According to ref. 9 the effective spring constant will be K - ,, K (21) where K = a + B is the original spring constant of the Hooke springs and a = p(E) can be calculated from (19). Notice that K8 is defined in ff the neighborhood where the distance between site i and site j is L (the lattice spacing). We can also consider the case where one spring is removed and the network is allowed to relax so that the distance between site i and site j is Leq a L. Tang and Thorpe8 found a = L - [.03 (22) 1-a L eq In obtaining Leq in (22), however, the authors used the effective spring 11 . l—fii—K for small deformation near Leq which is not a constant Keff = justified. Fig. 2.5 is a schematic drawing of elastic potential when one spring is removed. The whole curve may not be a parabolic, however, 77 for small 81 close to Leq it is indeed parabolic.[see Fig. 2.6] The point indicated by L is the point where a spring is removed and the dis- tance between site i and site j is L. Since the whole curve is not parabolic the effective spring constants deduced from points close to Le q and L may not be the same. We use K: to denote the effective spring ff constant near Leq' In Fig. 2.6 we show some of the computer simulation /L. We can see that b eff are calculated from the potential curves near 81 = O. In Fig. 2.7 we results for the elastic potential for different L0 for small Al the potentials are indeed parabolic. The values of K plot K6 and K: vs. L /L. The solid line is for Ke ff 0 It can be seen that K: ff and the squares are for K ff b . in general is less than Keff‘ eff' ff The difference in Kgff and Keff explains the discrepancies in Fig. 2.4(a) because one uses Ke to calculate the solid line while computer ff simulation are done in the neighborhood where the effective constant is b Keff' We can now conclude that for small 81 close to Leq the elastic potentials are parabolic but with K: is availabe to calculate K: different from Ke No theory ff ff' ff. The second derivative of elastic poten- a tial at L gives Keff = 1-:;E-K which can be obtained from the one a defect effective medium theory. 78 La, L .11 Fig. 2.5 Schematic drawing of the equilibrium positions when a spring is removed and the corresponding elastic potential. Potential “1. 79 -1 v I —1 ‘r 1 I W I I .1 + -+ L0/L-0.5 + + + + + + + + + + + + + + + + + ++ ++ ++ ++ + ++ ++| I |++ . 1 . 1 . 1 . 1 . O -0 6» -0 2 0.2 0 6 i 0 A( U I V I Y I V I 1 + -+ L0/L-.9 + .L + + + + + +1 + + + + + ++ + ++ ++ + ++ +4_ -F+ +++ +++ L l J J L l l I l .0 -4).6 -0.2 0.2 13.6 1 0 A! Fig. 2.6 Computer simulation results of the elastic potentials for L /L : 0.5, 0.9. 0 For small 81 they are parabola. 80 0.8 I I I I l I I I ‘ Kg“- 0 o 6 _ U eff —1 0 . . b . vs. LO/L. Keff is fiom (21) and Keff is deduced . b Fig. 2.7 Ke and K8 ff ff from the computer simulation results of the potential wells. 81 III Computer Simulations And Effective Medium Theory Approach Now we discuss the computer simulation results. Simulations are done in the stretched spring region where 0 S Lo/L S 1 for the model with elastic energy (6) on the honeycomb lattice. The numerical tech- niques to relax the network and obtain the static energy(E0), internal tension (T) and various elastic constants are as follows. Interested readers will find more details from ref. 8. (i) For a given L the honeycomb network is diluted by randomly 0’ removed springs with probability 1-p. (ii) The resulting network is then relaxed using the same methods and criteria as in ref. 8. (iii) After the relaxation is completed the static energy (E0) and tension (T) are calculated. (iv) Various elastic constants are consequently obtained by apply- ing certain designated small strains in two opposite directions in order to eliminate the linear terms in the elastic energy. In our computer simulations we calculate the following major quantities: static elastic energy (E0), internal tension (T) (called pressure in ref.8), bulk modulus (B) and the quantity b = % (CH-C12). For some samples we also calculate some additional quantities such as "s = 5 (C44 + 544) and "r = 31; (C44 - 544) to check the relations (15) and (17). In Fig. 2.8 quantities E T and B are plotted against probabil- O, ity p for L /L = 0, 0.4, 0.7, 0.8, 0.86 and 0.9. The solid lines are 0 the least square fit curves through the data points of the corresponding quantities. 82 '0 1 1 1 1 " 1.0 T , 1 , 4 ln/L-oo / [ML-o: o a - 4 to .. _ + £0 . A 1 / 0'3 A 1 J U U D O / 0.6 - - 0.6 1- //' .. OI. 5 4r i 1 Li. 10 r 1.0 , , “A Lo/L-oa 0 a )- + o.” _ + g. A 1 A 7 U U 8 0 5 — 0 6 _ 0 .1 - 0 ‘ _ 0.2 L- 02 “'- 0 - -l-___. I Refit.— 0 '40 °-5 °5 05 0.8 0.7 . p P ,1, I i 1.0 , , ”4'1““ 1.11/1. -0.9 08 i— + to A I 08 - z to 0 e D I: 05 06 .. 01— 0 a - 02 1. 0.2 )- 0 I I | 05 0.5 0 7 08 0.9 1.0 06 m C m a V Fig. 2.8 Static eneryg E tension T and bulk modulus B are plotted vs. 0’ p for L /L = O, 0.4, 0.7, 0.8, 0.86, 0.9. Discontinuities in B is 0 clearly observed. The solid lines are the least square fits to the corresponding data points. 83 We observe the following features in the graphs (1) Quantities such as E0 and T decrease almost linearly for small O/L . This is 8 similar to what Tang and Thorpe have observed on triangular lattice . p and vanish at the percolation thresholds for all L (ii) The bulk modulus 8 also decreases almost linearly for small p and vanishes at the percolation thresholds for LO/L < 0.8. For Lo/L 2 0.8 there are sudden drops in bulk modulus B at the percolation thresholds and the magnitudes of the discontinuities increase as Lo/L increases. We start to see the discontinuity about L /L 8 0.8. There 0 is some rounding in the discontinuities near the percolation thresholds due to finite size effects and averaging over many different configurations. This is because each different configurations has a distinct discontiniuty at a slightly different percolation threshold and the average Just smears out the distinct discontinuity. In the limiting case when LO/L = 1 we can think of the magnitude of the discontinuity as one.[see discussion below] In Fig. 2.9 we present two sets of typical results that can serve to check the two equivalence relations (15) and (17). Fig. 2.9(a) is for LO/L = 0.6 and Fig. 2.9(b) is for 1.0/L = 0.9. One can see that these two relations are still true in the diluted lattice network where springs are randomly removed. Also the discontinuity in bulk modulus B does not have any effect in the two equivalence relations. 84 Lo/L = 0.6 1 0 I I I I ”1” F B .F _ ‘ g) _ 0.8 [5 b .3 [J Us + O RHUP 0.5 - ‘7 5° 9 'J C) T in _+ v E3 (9 a 0 4 - - E! + V 9' IE 0 2 ‘- a; T 0 | l l 0.5 0.5 0 B 0.9 1 0 p Lo/L a 0 9 1 0 l I I l I r + + B + 0.8 ‘- D b 0 "‘ 0 2"ur + ‘7 £0 -+ ‘V O T 0.6 " A ”B +++ o .— + 4' ‘V 0.4 - C) " 4. C) V 0.2 *- IE 1% E! E! o ‘ ' 0.70 0.75 0 0 95 1 00 Fig. 2.9 Equivalence relations T : 2ur (15) and b = "s (17) are checked for LO/L = 0.6 and 0.9. 85 It is also interesting to see which quantities have discontinuities at percolation thresholds and which do not. From Fig. 2.9(b) and Fig. 2.10 it is clear that quantities such as E T, us, up and b are con- 0’ tinuous while B, C11 and C12 are discontinuous at the percolation thresholds. We will discuss the discontinuity in detail below. From Fig. 2.8 one can determine the percolation thresholds for different LO/L. In Fig. 2.11 the percolation thresholds from the computer simula- tions are ploted against L0 0 not 1 - g = g . This is because our system is not large enough. By /L. One may notice that pc for I. /L = O is studying different system sizes and using an extrapolating method one may obtain more accurate percolation thresholds. Our purpose here is to show that pc changes with LO/L. Fig. 2.11. is called the phase diagram because it indicates two phases. The region above the curve is the region corresponds to the concentration above percolation thresholds and therefore is called rigid phase while the region below the curve is the floppy region and there is a rigid-floppy phase transition cross the boundary.6’1u 86 1.2 I I 1 l I u + 1.0 _ + 011 + .53 C] B + E] 0 C12 E] 0.3 - El 3’ O 0 0 0.5 — "‘ 0.4 L- " dz 0 2 L- Q "‘ a 0 L 1 +4 I J 0.70 0.80 0.90 1.00 F3 Fig. 2.10 The quantities with sudden discontinuity at percolation threshold such as C11, C12 and B are plotted vs. p for LO/L : 0.9. 87 Percolation Threshold vs Lo/L 0 .50 I I I l I I I I I 04m ~ — U] [3 E] E] - E] J a a C10.30~- a a I o - [3 . H020--— a a .. [I] .4 040 — l a _ a. 0 l I l l i I l J l 0 0 2 0 4 0.6 o 8 1 0 Lo/L Fig. 2.11 The actual percolation thresholds for the stretched spring model from computer simulations are plotted against LO/L. 88 (2) Physical Reasons For Discontinuity We have seen above that simulation results indicate that for LO/L 2 0.8 quantities such as C C 2 and B have discontinuities at the per- 11’ 1 colation thresholds. Since the three quantities are related by B = % (CH+ C12) and the deformations associated with CI1 and B are more readily visualized we consider C11 and B here. He only need to discuss the case for B since the arguments applies equally well to C11. Now let us consider model (6) on a honeycomb lattice with LO/L = 1 which is the central force limit. The perfect lattice has bulk modulus 1 B = 573' K. Once the system is diluted , for example one spring is removed, then the bulk modulus is zero. The arguments are the following. Suppose the spring between site i and site J is removed. Since each spring is at their natural length then the positions of all sites will not change at all when the spring is removed. If one tries to measure the bulk modulus by applying small strains at the boundaries one will find that sites i and site 1 can simply pop in with their natural lengths and by doing so to allow other springs to retain their natural lengths.[see Fig. 2.12] This deformation will cost no energy and therefore the bulk modulus is zero. This is verified in computer simulations. In the triangular lattice at the central force limit this argument is not true because site i and site J are connected each by five springs and they can not simply pop in with their natural lengths and by doing so to allow other springs to retain their natural lengths. 89 . / Fig. 2.12 Site 1 and site j pop in with their natural length and by doing so to allow other springs to retain their natural length which results in the zero bulk modulus. 90 Notice that when LO/L : 1, above argument is also not true on the honeycomb lattice because every spring is stretched and removing Just one spring is not enough to relax the remaining springs (which are of order of N) without any change in energy. This is also verified by com- puter simulations. Now we can conclude that in the central force limit ( i. e. LO/L = 1) the honeycomb lattice will have zero bulk modulus if a single spring is removed in the central force limit. One may also recall that shear modulus vanishes even without removing any springs. We will use these instability arguments in the following discussions. Now we proceed to discuss the physical reasons for the discon— tinuities in bulk modulus at the percolation thresholds from the effective medium theory point of view. As we mentioned earlier a more sophisticated effective medium theory which is capable of dealing with more defects can be developed. The procedures of constructing such a theory are as follows. When many spring are removed, we can imagine an effective medium surrounds the remaining and missing springs. He can calculate the energy fluctuation of the remaining and missing spring surrounded by the effective medium. The fluctuation in tension between the effective medium and the remaining and missing springs can also be calculated. In this effective medium lattice all springs will have ef- fective strength Keff and effective natural length L: Obviously the ff' effective spring constant will be weaker. By requiring that the fluc- tuation of elastic energy and tension be zero, a set of self-consistant recursive relations can be derived. The detail of this effective medium theory will be in section IV. In this section we emphasize here that in 91 the proccess of constructing a multi-defect effective medium theory, one will always obtain an effective medium network with an effective spring and an effective spring natural length L: So for a ff’ f/L . 1 , as we argued above in the constant Keff given 1.. /L when the quantity Loef 0 central force limit, the bulk modulus will be zero even if only one ef- fective spring is removed. Since the bulk moduli abouve the percolation thresholds, for L /L 2 0.8 , have finite magnitudes which is roughly 0 controlled by the initial slopes in Fig. 2.11 (a), then the bulk moduli at the percolation thresholds require sudden discontinuities to go to zero. Therefore a more sophisticated effective medium theory should be able to explain the sudden discontinuities and predict the percolation O thresholds when L eff/L - 1. 0 eff medium theory point of view. To see how it is related to real physical The quantity L /L discussed so far is purely from an effective quantities which can be measured, at least in a computer simulation sense, one may do the following analysis. From (8a) and (We), the static energy and internal tension for the perfect honeycomb lattice are 2 2 E20 = W ( 1- LO/L) K (23a) 1 T : T( 1- LO/L)K . (23b) If one uses the physical pictures given by the effective medium theory then one can write down the static energy and internal tension similar 92 to (23a) and (23b) when springs are removed. By replacing K and LO/L 0 with Keff and Leff/L , we have 2 0 2 E0(p) ’ 7‘5 ( 1‘ Leff/L) Keff (2%) T = 1 ( 1- L0 /L) x (zub) (p) 7‘? eff eff where p is the probability that a spring is present. Of course EO(p) and T(p) are also functions of L /L, but for a given L /L they are func- 0 tions of p. Obviously quantities Lgf 0 f./i. and Keff are also functions of . 0 p for a given LO/L. For convenience we denote Leff/L by “eff' What we .are trying to do here is to extract Keff and "eff from computer simula- tion results of E0(p) and T(p) and to see whether they are still meaningful and useful when comparing with other simulation results. From (20a) and (2ub) we have 2 p)“ (25a) Keff : 2T( 003) 1.0 - ”eff = 30(p)/2T(p) . (25b) Now we make least square fits, which are plotted in solid lines in Fig. 2.8 , to discrete data points of E0(p) and T(p) to obtain values of Keff and ”eff against p. In Fig. 2.13 "eff is plotted against p for LO/L = O, 0.2, 0.4, 0.6, 0.7, 0.8, 0.86, 0.9 and 0.95. We notice that maximun 93 value of "eff is one for L /L 2 0.8. As we explained before, the bulk 0 modulus drops to zero when "eff hits one. 911 ,2 e131: 0 A I . 0 . . , Fig. 2.13 ”eff _ Leff/Leff deduced from the computer Simulation iesults and T using (25b). The data sets from bottom to top are for of E0(p) (p) LO/L = 0, 0.2, 0.u, 0.6, 0.7, 0.8, 0.86, 0.9 and 0.95 respectively. .1) 95 One may also recall the equations (111a)-(111d) and (16b) for elasticity of pure honeycomb lattice. By replacing K and n with Keff and "eff one may write down immediately the similar equations for systems with some springs removed. They are as follows u ‘3”err Keff n K _ eff eff C12 - ———2 - new 273" ‘26“ Cuu ‘ 2 - "eff 2) (26°) 2(1 - n ) K _ 1 eff eff K B = 2%: . (268) Equation (26e) indicates that B is proportional to Ke Now we use ff' K and n extracted from B0(p) and T(p) and calculate 611(p) , eff eff C‘l‘HP)’ B(p) and b(p) using (26a-e). The results are plotted vs. p in solid lines in Fig. 2.1“ for L /L = 0, 0.“, 0.7, 0.8, 0.86 and 0.9. 0 Also plotted in Fig. 2.1" are the corresponding quantities from computer simulations (in symbols). [Ill 0 0.0 + m .6 en a n O I I I to}! - o 7 + cu A en 0 I o t '5 I I toll-006 i.2~ + en t. en 0 I o u 0.9 >- L 0.6 - b o.) L- L— L o i L 0.5 0.6 0.7 96 I 1.0/['00 + cu A tn 0 I o u I i 0.7 0.0 p iii/y Fig. 2.1“ C11(p)’ Cflfl(p)’ B(p) and b(p) are plotted vs. p. The symbols are direct computer simulation results and the corresponding solid lines are calculated curves from (26a)-(26e) using only Keff and "eff' I’J G 0.0 1 1 4 i i 1 1 97 In Fig. 2.14 all quantities are normalized to §%§-. From Fig. 2.1“ we observe the following features: (1) The agreements between simulation results of C11, Cu”, b and B and the calculated corresponding quantities using only Ke f and F from E and T(p) are very good. ”err 0(9) ' 1 (ii) Keffs decrease to zero with p for L0 finite at the percolation thresholds for L /L < 0.8 while those stay > O/L - 0.8. Therefore by combining the observations of Fig. 2.13 and 2.111 we can conclude that (i) Keff and "eff are very meaningful and useful quantities to study in the stretched spring model. (ii) As more and more springs are removed, the static energy ED for L /L < 0.8 and and tension T approach zero with vanishing Keff O the bulk moduli decrease to zero continuously. . . _ > (111) E0 and T approach zero with vanishing 1.0 ”Eff for BD/L - 0.8 while K stay finite and the bulk moduli have effs discontinuities. (iv) The two different regions can be distinguished by whether Keff decreases to zero continuously or stays at a finite value at the percolation thresholds. One more interesting quantity to study is the ratio of CH“ to C . From 11 (26a) and (260) we find (1 - n )(4 - n ) Egg ; eff eff (27) 11 u ' 3neff 98 In Fig. 2.15 we plot Gnu/C from computer simulation results, vs. p. 11’ The corresponding solid lines in Fig 2.15 are calculated values of (27) using only K6 and "eff extracted from E0(p) and T(p)‘ Again they ff agree very well. The last point of each curve indicates the ratio of Cull/C11 at percolation threshold. We call this critical ratio of C11/CHQ' We can see that critical value of CHM/C11 changes from one for LO/L = 0 to zero for LO/L = 1. In Fig. 2.16(b) critical values of CM/C11 are plotted vs. L /L. We can see that it is quite close to a 0 straight line. The point where CM/C11 = 0 corresponds to that where discontinuities in B begin to show up. Fig. 2.16(a) is the same as Fig. 12.11. We can locate the point on the phase diagram where quantity C,m /C11 is zero by comparing the two graphs. In phase transition language, it turns out (see next section) that there are two types of phase tran- sitions across the boundary. In the region L /L < 0.8, where Gun/€11 0 are finite, the rigid to flopy phase transitions are of second order, as 8111 studied before ’ , while those for L /L 2 0.8 are first order. 0 Therefore there is a tricritical point about L /L ~ 0.8 in the phase 0 diagram Fig. 2.15(a) to separate the two different types of the phase transitions. Using Landau type phase transition theory we can also study the smoothness of phase boundary at the tricritical point. 99 1.2 1 F 1 I 0.u4 *- (3.2! '- Fig. 2.15 CHM/C11 are plotted vs. p. The symbols are from direct com- puter simulation results of C11 and CM! and the corresponding solid lines are the calculated curves using (27). The data set from bottom to top are for L /L = 0, 0.2, 0.u, 0.6, 0.7, 0.8, 0.86, 0.9 and 0.95. 0 100 0.50 t l v I v I . I . i I) . (l 0 40L- .4 U D C] D L E1 U 121 1 o t ({J H I 0.20 1- | El “ 1- : D 4 (L10 _ ' 4 ' U L ' 0. 0 . i . l . l . l . 0 0.2 0 4 0.6 0 B 10 Lo/L 1.2‘ Y I v I v T I r 1 1 1 1.0 i (b) d I 1 0 a — ' e ' ° ' I C" (I. 6 *- . _i . O I 0.4 ‘- — O 1 0 2 b O l _i 1 o l I j l l l 1 é 96-9-4 O 0 2 0 4 0.6 0.8 1.0 Lo/L Fig. 2.16 (a) is the same as Fig. 2.10. In (b) the critical ratio CNN/C11 is plotted vs. LO/L. A tricritical point about LO/L a 0.8 is located. 101 IV Tricritical Point and Effective Medium Theory We have seen in the previous section that the phase diagram Fig. 2.11 can be devided into two regions in which bulk modulus (or Keff) be- haves differently. For LO/L S 0.8, the bulk modulus (or Keff) decreases with decreasing p and vanishes at the percolation threshold. The phase transition in this region is of the conventional second order rigid e 6’1“. For LO/L > 0.8, on the other hand, the bulk modulus (or Keff) drops discontinuously at the percolation floppy phase transition threshold and the transition in this region is of a first order phase transition. In order to understand this we first study the Landau type phase transition theory and draw an parallel analogy between the phase transition theory and what we have observed in the computer simulations. (1) Landau Phase Transition Theory and the Tricritical Point In a conventional Landau type phase transition theory35’16 the free energy can be written as F : %a M2 + &b M“ + £0 M6 - H M (28) where M is the order parameter which indicates the phase transition, H is the external field, a, b, and c are the phenomenological constants independent of M. In (28) the coefficients of M3 and M5 terms are zero due to the symmetry of the system. In the magnetic model, for example, 102 M is the magnetization and H is the external magnetic field. When con- sidering a second order phase transition, b > O and 0 can be neglected. Thus the free energy has the following form _1 2 1 ‘1 F—éaMirr'bM-HM (29) and gg=aM+bM3-H (30) 2 a—E:a+3bM2 (31) The sufficient condition for the maximum and minimum of the free energy (30) is _35 3 - 3M _aM+bM -H-0 325 Depending on whethere —5 is greater than zero or less than zero the 311 free energy will be at minimum or maximum. For the convenience of the discussion we set a = aO(T-Tc) where T can be regarded as the temperature. The different cases are discussed below. We use 1? to denote the order parameter at the minimum or maximum. (1) When H = 0, T > To (i.e. a > 0) and b > 0, it is easy to see that ii = 0 is the only minimum. Therefore in this phase the order parameter is zero. 103 (ii) When H = O, T < To and b > 0, it can be shown that M = 1 11+ = aO(Tc- T) 1/2 ‘— [ --1;--] are the minima of the system. Thus in this phase as T + Tc ( a critical point) the order parameter M decreases continuously to zero. Combining (i) and (ii) we see that the order parameter decreases to zero and is continuous at the critical point T = To. Therefore the transition is of the second order. In the case when b is equal or less than zero, we must consider the M6 term. So the free energy is FzgaM2+l11bMu+ch6-HM (32) and ggzah+bii3+ci15~ii (33) 320 2 ii ——5 = a + 3b M + 50 M (3“) an For convenience c is set to be positive. When H = 0, using -3% = 0, we expect the free energy reaches minimum or maximum at R = 0i.:.” + and 1 M_ where 2__1__ 2 M? N: - 20 [-b 1_(b — MaO(T - Tc)c) ] (35) Now we discuss the following cases when H 0. (i) If a < 0, it can be shown that i 1 11+ correspond to two min- ima while B = O to a maximum. 10“ (ii) If aO(T—Tc) > O, and b > 0 then 51 = O is the only minimum and therefore in this phase the order parameter is zero. (iii) If aO(T-Tc) > 0 and b < 0 then the minima of the free energy occur at H = 0, 1_M+. It can be shown that F 2 0. So fi = :_M+ may (:M+) not be the lowest minima. This is true in the phase where E = ().is the lowest minimum. Since F( : 0, the location of the other two minima is 0) determined by requiring F (+M ) : 0. Thus we have three co—minima in —+ . . 3F this phase. USing F(1M+) - O and 3M M=IH+- O we obtain Ha (T-T ) 2 _ O c H+ - - -—-—E———— (36) From (iii) we can conclud that in one phase B = O and in the other flaO(T1-Tc)] b 1 der parameter changes discontinuously at the critical point T where T is the b 1 first order phase transition point. The phase boundary for the first 1/2 phase B = H = [- In general T x T , therefore the or- + C 1. The . . . c) 1/2 Jump or the discontinuity in M is AM : [- ] order transition is determined by equating (35) and (36) Ha (T -T ) i.e. 533—1-. 1 (02- 11a0(T1- Tc)c)1/2] = - 0 b‘ ° (37) From (37) we obtain _ 1/2 b - -‘1[caO(T1 - Tc)/3] 105 _ 3b or T —T + 16 (38) In Fig. 2.17 we plotted the phase boundary T vs. b. For b < 0 the phase transition is of the first order while for b > O the transition is of the second order. The two regions meet at b = 0 where the discontinuity AM vanishes. This point is called the tricritical point. We also notice that the first derivative, with respect to b, of the phase bound- ary at the tricritical point is continuous. 106 ; : LseconJ Oder : PM.“ boundary E :71 * 0 1 1 = b < o 5‘0 b > 0 b Fig. 2.17 The phase boundary is plotted vs. b. For b < 0 the transi- tion is first order while b )0 the transition is second oder. b = 0 is the tricritical point. 107 With above phase transition picture in mind, we are ready to draw an analogy between the conventional Landau theory and the observation in the simulations. An analogy between the phase transition theory and the geometric percolation theory was suggested almost two decades ago17. Of course this was restricted to the second order phase transition. The quantity p, the probability that a spring is present, corresponds to the temperature T in the phase transition theory. Thus the percolation threshold corresponds to the critical transition temperature. The quan- tity similar to the order parameter in percolation theory is the probability, denote by P(p)’ that a site chosen at random belongs to the infinite cluster. P‘D) is a much more difficult quantity to study in the elastic percolation problem than in the geometric percolation problems. The bulk modulus, on the other hand, is a much easier quan-- tity to study and has the same behavior as P(p), so we choose the bulk modulus as the order parameter. The analogy is quite obvious now. For LO/L S 0.8, the bulk modulus decreases continuously to zero and is con- tinuous from rigid phase to floppy phase. Therefore the transition is a second order phase transition. For LO/L > 0.8 the bulk modulus drops to zero discontinuously from the rigid phase to floppy phase and the tran- sition is a first order one. There is a tricritical point around L /L ~ 0 0.8 to separate the phase diagram into two regions. Here LO/L is similar to b in the phase transition theory. Our computer simulation results indicate that there is a tricritical point near L /L ~ 0.8, and 0 a more accurate determination of the tricritical point needs more com- puter simulations. This tricritical point can be determined from the 108 phase transition theory once the phenomenological constants a, b and c are known, but unfortunately this is not the case here. Thus we have to determine the tricritical point from an new effective medium theory described below. (2) Effective Medium Theory We start discussing the effective medium theory by first consider- ing the case where only one spring is removed. As shown in Fig. 2.5, when a spring is removed the distance between the sites, where the spring is removed, will change from L to Leq' Let us assume that the resulting spring network will have an effective elastic constant Ke if the distance is restored back to its original length before any spring is removed [also see discussions in appendix 11]. Therefore after one spring is removed the elastic energy is decreased by 2 1 2 1 AB = - 2K(L - LO) - 2 Ke(L - Le ) q (39) where K and LO are the spring constant and natural length of the removed spring. The first term is due to the missing spring and the second term is due to the relaxation afterward. By minimizing AF with respect to L we have 8A8 a L = K(L - L0) + Ke(L - Le ) = 0 (no) 9 From (“0) we obtain K Leq — L + K ( L - L0) (41) e ' a Using Ke : K 1 -*a where a is equal to P(E) as in (19), (01) becomes a a L = L - L03 (“2) eq '_——_—_3. 1 - a a We have discussed in section II that K8 2 K l-:;§- is Just an approxima- a tion, but it is a rather good approximation. Substituting (42) back to a (39) and using Ke = K 1 -.a , the energy change AE can be written as a AB = - 1 K (L - L )2/( 1 - a.) (#3) 2 0 Eq. (“3) indicates that removing one spring decreases more elastic energy than %K(L - L0)2. This is because the network is under tension (L at L0) and the relaxation to Leq releases an additional amount of elastic energy. Now let us add a spring of spring constant K' and natral length L6 in the place of the missing spring. The elastic energy will increase by 2 2 AB = + gK'(L' - L6) + % Ke(L' - L ) (4“) eq 110 Again the first term is simply due to the added spring and the second term is due to the relaxation of the distance from Leq to L' afterward. By minimizing AF with respect to L' we have 3A8 5? = K(L' — L5) + new - Leq) = 0 (“5) So K L + K'L ' . - e as 0 L ' x' +xe (“6) where Ke and Leq are the same as defined before. Substituting (“6) back to (nu) and using (42), the energy change can be written as The net energy change after replacing a single spring of K and Lo with K' and L6 is then - K [ " I 2(1 - a ) K. i i K‘a + K(1 -a ) (L - L a*- L' + L'a*)2 - (L- L )2] (N7) 0 0 0 0 111 If K' = 0, Eq. (47) is reduced to (H3). If K' = K and L6 = LO then AB = 0 as expected. So (“7) is checked. We now construct an effective medium theory using (117). For a given LO/L at p, from the effective medium point of view, we can imagine an effective medium elastic network of N effective springs with spring constant KC1 and natural length Lt? Here the superscript a denots the quantities of the effective medium network. Of course Ka and L: are functions of p for a given LO/L. Let us replace n effective springs with n springs of strength K1 and natural length L Therefore the to- 01‘ tal energy of the resulting effective network, using (117) and replacing * with a, is c x x L - L“ _ _ _NK 0 2 _1_ i e 0 _ c 2 E ' E0 * AB ’ 2 (L ‘ L0) * 2 Z x + x l 0 L01 * L0] 1 1- a 1 Ka(L - Lg)2 - 2 X a (118) i 1 - a a 1- aa a a a where Ke = K a and a = a ( LO/L). The fisrt term, denoted by a EO’ in (118) is the static elastic energy of a perfect effective medium network similar to (2ua). The second and third term, denoted by A8, are the energy fluctuations due to replacing n effective springs with springs of K1 and L . In (118) we neglect the interactions between the Oi replaced springs. The tension is defined by T = 33'; . In 20, T = - ‘ av 511'.- %E Using (‘13) we obtain ‘1 112 a KK L-L BE: 0 0 ie 0 01 T ~ - ——»= T + AT = ~NK (L - L ) - X [ - L + L ] 3L 0 0 i “1* xe 1_ aa 01 0 1_aa 0 a a i 1-a L ii e 1-a (1-a) K a L - L“ K“(L — L“)2 gl (x +1x )2 Eazl g ’ LOi * L812- % ace 1 ("9) i i e a 1-a i (1-a) a Baa a a whereb =--a- andn =LO/ L. In obtaining (119) one should realize an that Ke and aa are also the functions of L. The first term in (119) is the static tension, denoted by To, similar to (24b) and the second term and third term are the fluctuations of tension, denoted by AT, in the diluted effective medium network. By requiring 11 O < AE > (5021) (AT) 11 0 (50b) we obtain two equations describing the effective medium network. In (50), <> implies the ansemble average. The two equations are a a a2 3i—.“:“3 ——3-..+ 32 3: ““23” i i e 1-a i 1-a 113 a a 2 K (L - LO) - i . + l - -' 1 (51b) 802 1_ a“ 01 0 2Ni (1 _ a(1)2 1 -— (———) 2N§ K1+ Ke In our case K1 = K and L01 = LO with probability p, and lg = O with probability 1-p. So (513) and (51b) become a 2 - K (L - L ) ___e __0 _ L . Lg]? . 0 (52a) a 1 _ a - L0 + L0] 0 ' a ’ 1— a 1-a x K L - L“ e 0 - L + La] 0 0 0 L e 1- a (1 - aa)2 Ka(L - Lg)2 a)2 1 (52b) a L - L _ L a 2 1 [ + L ] - a2 1_ as 0 O 2 (1 _ a 5 a Eq (52b) can still be simplified by using (52a) to eliminate the second term in the left hand side of (52b). By using Ke = [(0 -———1 " a 1111 defining 0“ = Ka/K, n“ = Lg/L and n = LO/L (52a) and (52b) can be writ- ten as a a a (1 - naaa - n mac)2 = ( 1 - na)2 (53a) a + Q (1-a ) (I. a a O. a (1 c I 1 - n“ + (00- n)(1- a°)l 1 - n“ 21a“ + Q“(1 - a“)] Here superscript a is used to indicate the effective medium network. Ka and na are the quanties characterising the effective medium network and should be compared with Keff and n from section III. Eq. (53a) and eff (53b) are the self-consistant recursive equations. For a given LO/L, Eq. (53a) and (53b) must be solved self-consistantly and iteratively ii 3a a a “LO/L) . We obtain K and n ii ii starting from p=1, aa = a and ba = b = as functions of p and these values should be compared with the values deduced, using (25a) and (25b), from the computer simulations. The per- colation threshold pc can also be calculated as a function of n using (53a) and (53b). We use nt to denote the value of n at the tricritical a/K = 0 at the percolation threshold thus point. For n < nt, Qa = K (53a) and (53b) are solved iteratively for a given n < nt untill Qa = 0. For n > nt, no : Lg/ L = 1 at the percolation threshold, again (53a) and (53b) are solved iteratively for a given n > nt untill na = 1. In the neighborhood sufficiently close to the tricritical point, (50a) and (50b) can be further simplified. After some lengthy algebra we obtain 115 2 [1 + ba(1 - n )l = 1 (sue) aa + Qa(1-aa) a a a a L1 - 0H: " C /2) = ab (10" Q) a (Sub) [ 1 + b (1 - n)! Zia + 9 (1 - a )i 0 32a“ 0 At the tricritical point 0“ = 0 , n = 1 , a?1) = where c = 301“)2 1, 83(1) = 0.958, and c?1) = 0.16. Thus from (548) we obatin 1 - nt = 0 ba 02 = 0.186 2b + c —b or nt = 0.8111 Using 0“ = 0, "t = 0.811, a?1) = 1, b (1) = 0.958 and (51a) we obtain or 1 - pt = 0.280 Thus the effective medium theory gives the location of the tricritical point at po = 0.28 and nt = 0.81. In Fig. 2.18 we plot the phase bound- ary obtained from the effective medium theory and from computer simulation results (squares). The tricritical ponit is also marked in the figure. The value of "t is very close to our estimates LO/L ~ 0.8 from the computer simulation result. We notice that pc, determined from 116 the computer simulations, at L /L = 0, is 0.38 instead of 1-Zsin(i'-8-) = 0 0.316 which is the exact result3. This is due to the finite size effect. As the system size becomes larger and larger one should be able 3 to extract the exact value by using a scaling rule . It is more impor- tant here to known how pc changes with L /L. The most interesting 0 feature in Fig. 2.18 is that the first derivative of pa, with respect to LO/L, is not continuous at the tricritical point. In Landau phase tran- sition theories, as we have discussed, the first derivative at the tricritical point is continuous. The discrepancy may be due to the negligence of the interactions between the springs in the effective medium theory. An effective medium theory capable of dealing with the interactions is very complicated and technically impossible to apply. After obtaining KCl and "a as functions of p, we can use equations (26) to calculte various quantities such as E0(p)’ T(p), and B(p)’ Notice that Eq. ('18) and (‘19) are the same as (211a) and (2%) when (AB) 0 and = 0, which determine the effective medium network. In Fig. 2.19 we plot na against p for L /L = 0, 0.2, 0.11, 0.6, 0.7, 0.8, 0.86, 0 0.9, 0.95. The corresponding results from the computer simulations using (25b) are also plotted (symbols). The overall agreement is satisfactoy. The agreement are very good for small 1-p, as we can ex- pect, and above the tricritical point. In Fig. 2.20 we compare the quantities E(p), T(p)’ and B(p) obtained from the effective medium theory and the computer simulations. We can see that the agreement, in general, is also quite good although the effective medium theory overes- timates the percolation thershold in some cases. Therefore we conclude P96 117 general, is also quite good although the effective medium theory overes- timates .the percolation thershold in some cases. Therefore we conclude that the new effective medium theory works quite well for the honeycomb spring network under tension although it has some shortcoming. It predicts a tricritical point and the different behavior of K“ and 0° in the first and second phase transition regions. These are the most sucssesful features of the effective medium theory. 0.5 I I f l I I I I I (34:5 C] a [3 _ a I: 0.3~— . - EL i 0.2 *- : .3 I 1 l 0.1 - : ~ ' I 0 . i i i i I . ii i 0 0.2 0.4 0.6 0.8 1.0 Lo/L Fig. 2.18 The phase boundary from the effective medium theory ( solid line) and computer simulations (squares) are plotted against n. The tricritical ponit is marked. 118 1.0 0.8 - 0.6 - 60.4; Fig. 2.19 The quantity na is plotted vs. p for LO/L = 0.0, 0.2, 0.11, (L6, 057,(L8, 0.86, 0.9 and 0.95 (solid lines). The corresponding quantities from the computer simulations are also plotted (symbols). .0 119 Lol’L '0.9 0.8 ~ L. 0'8 + to A T U B 0.6 '- 0.‘ r- 0.2 3- 0 0.5 Fig. 2.20 The quantities E0, T and B calculated [WmMIthe effective medium theory (solid lines) are compared with the corresponding quan- tities from computer simulations (symboles) for LO/L = 0.0, 0.”, 0.7, 0.8, 0.86 and 0.9. 120 V Conclusions We have studied the elastic percolation problems of the stretched spring model on the honeycomb lattice and found two types of rigid s floppy transitions. In the region where L /L S 0.81, the rigid s floppy 0 transition is the conventional second order phase transition studied before8 while in the region where L /L > 0.81 the transition is a first 0 order phase transition. The difference between the two types of transi- tion is that the bulk modulus is continuous in the second order phase transition and discontinuous in the first order phase transition. Landau phase transition theory can be applied to explain the two types of phase transitions and the tricritical point observed in the computer simulations. A new self-consistant effective medium theory has been developed to calculate various elastic constants, phase boundary and the tricritical point. The agreement between the effective medium theory and the computer simulation results is quite good. We have also found that the the elastic percolation threshold of the stretched spring model is 1 rather than 3 in the central force limit on the honeycomb lattice. VI 10. 11. 12. 13. 14. 121 REFFERENCES P. G. de Gennes, J. Phys. (Paris) 31, L-1,(1976). S. Kirkpatrick in Ill-Condensed Matter-Les Houche 1978, Ed. By R. Balian, R. Maynard and G. Thoulus, (North-Holland, Amsterdam, 1979). D. Stauffer, Introduction to Percolation Theory (Taylor and Francis, London, 1985). S. Feng and P. N. Sen, Phys. Rev. Lett. 5g, 216 (1984). A. B. Day, R. R. Tremblay and A. -M. S. Tremblay, Phys. Rev. Lett. 2Q, 2501 (1986). M. F. Thorpe, J. Non-Crys. Sol. 51, 355 (1983). M. Born and K. Huang, Dynamical Theory of Crystal Lattice (Oxford, Univ. Press, New York, 1954). W. Tang and M. F. Thorpe, Phys. Rev. B 31, 5539 (1988); W. Tang, Ph. D. Thesis, Michigan State Univ (unpublished). S. Feng, M. F. Thorpe and E. Garboczi, Phys. Rev. B 31, 276 (1985). T. H. K. Barron and M. L. Klein, Proc. Phys. Soc. 85, 523 (1965). C. Kittel, Introduction to Solid State Physics (J. Wiley 8: Sons, Inc. New York, 4th Edition, 1971). P. N. Keating, Phys. Rev. 145, 637 (1966). A. E. H. Love, A Treatise on the Mathematical Theory of Elasticity, (Dover, New York 1944). H. He and M. F. Thorpe, Phys. Rev. Lett. 54, 2107 (1985). 122 15. L. D. Landau, Phys. 2. SowJetunion ll, 26 (1937); V. L. Ginzburg and L. D. Landau, Zh. Eksp. Teor. Fiz. 29, 1064 (1950); R. B. Griffiths, Phys. Rev. Lett. 24, 715 (1970). 16. Kerson Huang, Statistical Mechanics (John Wiley 8: Sons, New York, 2nd Edition, 1987). 17. R. Zallen, The Physics of Amorphous Solids, (Jhon Wiley & Sons, New York, 1983). Part III The Vibrational Spectra of Random Superconductingy- Normal Networks 123 124 I. Introduction It is well known that the vibrational excitation density of states for homogeneous systems obeys the Debye law g(m)dm ~ md - 1d0) (1a) in the low frequency or long wavelength limit. Here in is the vibra- tional frequency, g(m) is the density of states and d is the dimensionality of the problem. Naturally one can ask such question as what happens to (1a) when disorder or inhomogeneity is introduced. Obviously Eq (1a) will breake down for inhomogeneous systems. However, for certain disordered systems, Eq (1a) simply takes the following form1 ~ g(w)dw ~ wd - 1d0) (1b) where 8 is called the spectral dimensionality which is different than d. A fractal network is a suitable inhomogeneous system for studying quan- tity 8. There have been extensive studies on dilute random percolation networks, which are fractals, in the past five years1-8. If we think of the vibrational model as an array of springs, then the vibrational equa— tions of motion are very similar to those of diffusion and a tight- binding Hamiltonian [see (2a) below]. Because of the similarities among the diffusion equation, equations of spring vibrations and the equations of motion for a system described by a tight binding Hamiltonian (or an 125 elastic Hamiltonian) it is clear that there are two approaches to evaluating 8. One approach is to examine the low frequency dynamical response of the system described by the Hamiltonian _ +_+ H - <1§> ViJ(aiai aiaJ) (2a) where the V1 are randomly distributed on the nearest neighbor bonds ac- J cording to the probability distribution P(V1J) = 05(V - V13) + (1 - p) 8(Vij) (2b) and we restrict the sum to sites on the infinite cluster. The ai, a; are either Bose or Fermi operators (the results are the same when only one particle excitations are considered) and the diagonal term is present in (2a) to preserve translational invariance. This dilute tight binding Hamiltonian describes the physics of several different random systems; dilute magnetic systems, random resistor networks and random scalar elastic networks. By scalar elasticity we mean for example dis- placements perpendicular to the elastic network. Symbol E has the unit of energy and when h:1 we call it frequency for convenience. The low frequency density of states at pc is predicted to behave as 5-2) g(E) ~ E (3) NIL—o 126 where 8 is the spectral dimensionality that can be expressed a31-2 E = ———————- (4) and t is the conductivity exponent. The quantities B and v are the geometric exponents8 which govern the probability of being on the in- finite cluster, and the correlation length, respectively. The fractal dimensionality a = d - B/v. Direct numerical evaluation of the density 3 of states in two dimensions on dilute spin systems has confirmed equa- tion (4) and that a is close to 4/3 which Alexander and Orbach1 have conJectured to be the value of 8 in all dimensions except 1D. An alternative approach is to study the anomalous diffusion of a random walker on the infinite percolation cluster. Scaling arguments"2 ,4 suggest that the mean square displacement (r2) scales with time as (r > ~ t (S) and k = -—-£-B- where t, B and v were defined previously. By relating the probability of return to the origin with the density of states,1 one obtains a = 8k in agreement with Eq. (4). Numerical simulations of random walks on dilute systems are a good indirect way of obtaining the spectral dimensionality 3 using Eq. (5) 127 and the relation 8 = k8. Values of 8 obtained in this ways 3 agree with direct simulations of the low frequency density of states at pc. The randomly dilute problem discussed above can be considered as a random system of strong bonds V3 and weak bonds Vw in the limit VS+V, szo. The other limit, Vsom, VH¢V, is also of considerable interest and is the subject of this article. It would describe a random supercon- ducting network or the elastic properties of a system with rigid grains weakly coupled by soft springs. In the general problem, the density of states is divided into bands: a high frequency band associated with in- ternal modes of the clusters of strong bonds and a low frequency band associated with vibrations of the soft regions. In the limit sz 0 Vs» 1, the low frequency band becomes a delta function at the origin and the spectral dimension describes the low frequency edge of the high frequency band. In the limit Vw » 1, Vs .. on, the high fre- quency band is driven off to infinity and the delta function broadens. There will be a different spectral dimension associated with this low frequency edge. In diffusion language, the superconducting limit is the problem of 9 10-12 termite diffusion which has been the subject of considerable work and some controversy. In principle, it should provide information about the density of states and the spectral dimension 8 for the problem but in practice this is not the case because the random walk is completely dominated by the distances moved on the superconducting clusters. As has been pointed out by Hong et al.,11 ~ tk because even at time zero, (r2> ~ g2-B/v’ the average radius there is no region where 128 of the superconducting clusters. The low frequency excitations we are interested in are associated with the small, time dependent contribution to superimposed on the divergent £278”. Clearly this is impos- sible to obtain numerically and we therefore believe that direct simulations of the low frequency response is the best way to evaluate the spectral dimension. From ones physical intuition and experiences one will also expect that (a) the contributions to the lower frequency part of the vibratitui density of states are associated with the collective vibrational modes of clusters and because such collective modes involve large masses therefore result in lower frequencies; (b) those of high frequency part are due to the internal vibrational modes of clusters or local vibra- tional modes which are associated with small masses and therefore with higher frequencies. While the lower frequency part can be explained using scaling hypothesis, the major features in high frequency part can by explained by Just consider the local vibrational modes. The layout of this paper is as follows. In section II we describe the procedures of obtaining Density Of States (DOS) for Hamiltonian (1) by using Equation Of Motion technique (EOM) and present the computer simulation results for spectral dimensionality 8. In section III we present the scaling cross-over hypothesis to describe the cross-over from phonon region to anomalous region for systems above two dimensions and in sec- sion IV one dimensional case is treated seperately because of unconmon features in one dimension. 129 11. Density of States In this article we obtain the density of states and hence the spectral dimension 8 for the superconducting case. The system is described by the Hamiltonian (2a) but with the interaction strength distribution ) + (1 - p) 5(V - V ) (6) P(vi ) = p 5(Vm' v1 J J 1J where Von = w in the superconducting limit. The wavefunction of the sys- tem at time t can be written as |T(t)> = X ck(t) ak*|0> (7a) k where ID) is the vaccum state, ak+ ck(t) the amplitude of excitation at site k at time t. By substituting the creation operator at site k and |T(t)> into Schr8dinger equation ih éigéfill = H lT(t)> (7b) with Hamiltonian H of (1) then (7b) becomes ih X a+l0> ac (t)/3t = Z V (a+ a - aTa k k k ijk ij i i i j We are interested in how amplitudes ck(t) evolve with time t. By using ) a;io> ck(t) <7c) either commutation or anticommutation rule ] = 5 (7d) the equations of motion for ck(t) in (70) are now ih 5;;— = z; vk,(ck- e2) (7) Here the ck's are the amplitudes of the wavefunction on the k site. Two sites k and 2 connected by a bond V on have the same k2: amplitude so that ck = c2 can be considered as a single site of mass 2. This process is repeated for all superconducting bonds so each supercon- ducting cluster has one degree of freedom and a mass M equal to the number of gitgg in the cluster. We consider a site connected to no sue perconducting bonds to be a cluster of size 1. The equations of motion for the clusters so defined are . k 1h MK 77': = {V nk2(ck - CQ) (8) where k,2 are now cluster indices and n is the total number of normal k2 _ 1 bonds between the k and the 2 cluster. By defining le ck = ck, we have the Hermitian equation of motion C ih-aT-sXVnkQ[-M-:-- ] (9) 131 13 14,15 which we solve numerically using the equation of motion technique. If the original system has N sites,there are only Nf clusters so we have reduced the number of degrees of freedom by a factor of f. For small p f:1-%2+... (10) where z is the number of nearest neighbors and p is the probability that a bond is present. The number of clusters in equation (10) actually is the sum of number of monomers, dimers, trimers and etc. normalized to total number of sites. Here monomers, dimers, trimers and etc. have mass of one, two, three and etc. respectively. For square lattice, as shown in Fig. 3.1 , the probabilities of having monomers, dimers and trimers are (1-p)”, 2p(1-p)6 and 6p2(1-p)8 respectively. In above ex- pressions powers of p represent the presence of bonds and powers of'(1- p) represent the absence of bonds. The coefficients in the front indicate number of different configurations one can possibly have with fixed number of bonds. The probabilities of having clusters of less than five bonds can also be written down by just inspection and those 15 for more than five bonds can be obtained by using graph theory . Therefore by definition the number of clusters at p on square lattice is f : (1-p)” + 2p(1-p)6 + 6p2(1-p)8+ ... (10a) 132 Extending (10a) to a general lattice with z nearest neighbors we can write down immediately the contributions from monomers, dimers and trimers etc. 2‘2 *2(z'1)+ ... (10b) f = (1-p)z + 3 p(I-p)22’2+ 2(2—1)pz(1-p) So f is a function of p and 2. For given 2 in order to find f(pc) we have to add the contributions from higher order clusters which can be found in ref. 15. For small p equation (10) can be easily obtained from (10b). 133 Fig. 33.1 Diagrams of monomer, dimers, trimers. The solid lines repre- sent springs and the dashed lines represent missing springs. There are: one configuration for monomer, two for dimers and six for trimers (omitted here). 134 Note that f decreases to a very small value at pc. In the square net and simple cubic lattice we find it is 0.098, and 0.27 respectively.16 Thus setting V0° strictly equal to infinity reduces the effective size of the system, especially in 2D, although there is some initial time cost in identifying all the clusters. This advantage cannot be acheived by setting Va large and taking a limiting process. Note also that f is finite at pc for d l 2 but it is 22:9 in 1D which affects the scaling arguments for d as we discuss later. The equations of motion were integrated forward in time from 0 + T‘ 13 using standard methods. The first few points were obtained using a fourth order Runge-Kutta algorithm and subsequently a fourth order Adams 17 integration scheme was used. The density of states g(E) is obtained from Fourier transforming the result of the time integration. +0: g(E) = %;,‘ Z M1 c1(t)c:(0)eiEtdt (11) i -0 13 This integral is converted to one over positive times only. Random 13 phases are put on the clusters so that only the on-site terms required in (11) are retained in the limit of a large system. This is a standard way to obtain densities of states in the equation of motion technique. In Fig. 3.2 we show the superconducting clusters of one sample at p = 0.475 on 40x40 square lattice. On each lattice site there is a single mass. The solid lines represent rigid springs with Vs“: while missing lines represent the weak springs with VN+1. Dots are the masses 135 surrounded by weak springs. In Fig. 3.3 , we show results for the square lattice at p0 = $3 The results are from averages over 25 samples of 100 x 100 networks. Because of the reduction in size due to the f factor only ~ 1000 amplitudes had to be monitored. The insert shows a log-log plot at low frequencies. In Fig. 3.4 , we show similar results for the simple cubic lattice8 at p0 = .2491. The results are fnmn averages over 25 samples of 21 x 22 x 23 networks. Because of the reduction due to the f factor only ~ 2800 amplitudes had to be monitored. The insert shows a log-log plot at low frequencies. We ex- amined a single sample with about 16 times as many sites in both 2D and 30 and found no significant changes in the slopes in the log-log plots. There are two ways to obtain the spectral dimensionalities. One is simply measuring the slopes of the log-log plots in the inserts. We use a least square fit program to find the slopes. The slope in Fig. 3.3 is 1.05 1 0.1 which leads to 3 = 4.1 t 0.2 from Eq. (3) in 20. Similar mesurments of insert of Fig. 3.4 gives a slope of 1.9 i 0.15 and hence 8 = 5.8 t 0.3 in 3D. Another method is associated with the crossover scaling hypothesis which will be discussed in next section therefore we postpone the discussion of the method. 136 3333+ ' 33.333333 '7‘] Fig. 3.2 A sample of bond percolation network on square lattice at p::CL475. The solid lines represent rigid springs and the missing lines represent soft springs. The dots are the masses surrounded by soft springs. 137 4. I 1 I 06 L—- oIoI— I T I ——-4 00a - ' 3 a O 4 ~— n -i 1 001 L l 5 02 06 io n 2.6 n 0.2 *- - 5.4 0 l I 0 3 9 12 E Fig. 3.3 The density of states for a random square lattice of supercon- ducting and normal bonds at pc 2 -2-. The insert shows the low energy data on a log-log (base 10) plot. The straight line drawn is a least squares fit to the data points shown from which 8 is obtained via Eqn. 3. The energy is in units with V : 1. The energy values of major peakes are also indicated. g(E) 0 . 0 . 138 1 1 1 1 6 0'0 I V V T 00!)- 1 004- 3 _. .2 S 2 .._ L i i I —-l 10 15 20 25 30 1 - .. 0 l i l 0 4 12 16 E Fig. 3.4 Same as Fig. 3.3 except for a simple cubic lattice at pc : .249. The insert shows the low energy data on a log-log (base 10) plot. The energy is in units where V = 1. The values of peaks are also indicated. 139 As we mensioned earlier in the introduction, the higher frequency (or energy) part of spectra are due to local vibrational modes. Therefore the sharp features at higher frequencies in Fig. 3.3 and Fig. 3.4 can be explained by just considering these local vibrational modes. Since all sites are always connected in superconducting case (i.e. stm and Vus1), at percolation threshold, we can imagine that (1) there is an infinite cluster (largest spaning cluster in computer simulations); (2) there are clusters of different sizes interconnected by weak springs ; (3) there are 'pockets' inside the infinite cluster and large clusters where small masses are surrounded by weak springs and connected back to big clusters. The latter situation is very similar to the case where small islands are surrounded by inland lakes and the lakes are sur- rounded by continents. Here small masses correspond to smll islands, weak springs to inland lakes and big clusters(including the infinite cluster) to continents. The collective vibrations contribute to low frequency part of spectra. Now we consider the vibrational modes of 'pockets' . We assume that big clusters to which the small masses at- tached are very massive and therefore almost stationary. As shown in Fig. 3.5(a) a single mass is attached to big clusters (represented by the shaded area) with four weak springs of strength Vw which we can al- ways assume to be one. The vibrational frequency is simply E = 4 for four neighbors in the 2d square lattice. In Fig. 3.5(b) each site has four springs and there is one spring between the two sites so the secular equation can be written as 140 = 0 (12a) Equation (12a) gives eigenvalue of E = 3 and E = 5. In Fig. 3.5(c) two masses, connected by a massless infinite rigid spring, give a frequency of E = g = 3. Configuration such as one shown in Fig. 3.5(d) gives a secular equation -1 4 - E -1 = 0 (12b) which gives eigenvalues of E = 4 - l2 , 4, 4 + /2. Fig. 3.5(e) gives a frequency of E = g = 2.67. Of course we can go on to exhaust all pos- sible configurations and obtain the corresponding eigenvalues. However the above local vibrational modes are enough to explain the sharp fea- tures in Fig. 3.3 at frequencies E = 2.6, 3, 4, 5, 5.4. The relative heights of the peaks are related to the relative abundances of these lo- cal vibrational modes which can not be estimated by a simple argument. 141 Fig. 3.5 EHagrams showing the vibration modes of 'pockets' in 2D. The solid lines represent rigid springs and the dashed lines repnesent soft springs. The shaded areas represent infinite masses. 142 Similarily the maJor peaks in 3D case (seeiFig. 3.4) can be ex- plained in a similar manner. For example, Fig.3.6(a) gives E = 6 for six neighbors in 3D simple cubic; Fig. 3.6(b) gives 6-E -l -1 6-E :0 (120) which yields E = 5 and E = 7; Fig. 3.6(c) gives E = 1% = 5; Fig. 3.6(d) gives 6 - E —1 0 -1 6 - E -1 : 0 (12d) 0 -1 6 -E Idiich results in E = 6 — J2, 6, 6 + J2; and Fig. 3.8(e) gives E = %fl : 4.67. So peaks at frequencies 5, 6, and 7 can be explained. Again the relative heights of the peaks depend on the relative abundances of dif- ferent modes. We also notice that the relative frequency differenc/es form secular equations in 3D is smaller than those in 2D. Therefore the frequencies in 3D are more close to each other and evenly dis- tributed which results in less spiky spectrum than in the 2D case. 143 cubic lattice. opt for 30 simple Fig. 3.5 exc Fig. 3.6 Same as 11m 111 Scaling Arguments The scaling argument for 8 is given below following the ideas developed for the dilute case.7’8 It is easiest to use the language of phonons, i.e. vibrations on a network of normal and infinitely rigid springs. The results apply equally well to any system described by the tight binding Hamiltonian (1) with the probability distribution (6). Close to pc and for systems of size L much greater than the correlation length i ~(pc-p)-v we expect the system to appear homogeneous with an effective elastic modulus Y ~ (pc-p)'8 where s is the superconducting exponent. Because the system is fully connected, all sites contribute to the inertia in the long wavelength limit, the total mass M~Ld and thus the mass density p is noncritical, i.e. it is independant of E. The low frequency excitations are phonons with a sound velocity 02(g) ~ Y/p ~ gS/” . (13) For wavelengths l~§ the fractal structure of the lattice becomes important. The excitations are no longer phonons and we have a fracton- phonon crossover. The crossover frequency is given by (A) ~£££l ~ E-(Z "' S/V)/2 . (1”) co g Using a Debye type theory we can write the normalised density of states for the phonons as 145 d p(o) ~ C(g)’d (£3) od" /Nf . (15a) To alltnv for the fracton-phonon crossover we introduce a scaling func- tion h(m/mco) into (15a) and obtain d- d 1 ) w h(w/mco) /Nf . (15b) p(m) ~ cm'd {-3,- For d.3 2, f is noncritical and so the critical behaviour of p(m) is -ds dS/V p(w) ~g21) wd-l h(g ) ~ wCOZ-S/V (dd-1 h(lfi )) (16) CO CO where we have used (13) and (14). In the low frequency (phonon) regime x = ”/wcoi<<1 and h(x) + constant, so that p(w) ~ md-1 as expected. For m > ”00’ the fracton regime, we rewrite (16) as - ds/v 2d 2-s/v -1 ~ 2-..). i3 l he wad“ r(-“1 1 mi woo woo woo p(w) ~ w ds/v - 2 - s/v where r(x) = x h(x) is expected to tend to a constant for large x because p(m) should be independent of g or equivalently of “co in this limit. Equation (17) defines the spectral dimension, ~ _ 2d d---——-2_s/v (18) 1146 We have discussed scaling relations for vibrational modes. Exactly parallel arguments can be made in all dimensions for tight binding Hamiltonians where E ~ me ~ 02k2 leading to g(E) ~ E . (19) Comparing Eq. (18) with the result for dilute systems Eq. (11), we notice the following differences. a) The numerator involves d rather than a. The excitations in the superconducting case can explore the whole system whereas in the dilute case, Eq. (11) is obtained for excitations on the infinite cluster. If the finite clusters are included in the dilute case6 the d is also re- placed by d in Eq. (u). b) The conductivity exponent t in Eq. (11) is replaced by the su- perconducting exponent -s in Eq. (18). c) The exponent B does not appear in Eq. (18). In the dilute sys- tem, the probability of being on the infinite cluster scales as (p - pc)8. In the low frequency dynamical response of the infinite cluster, this is the inertia term and appears in (14). In the superconducting system, all sites are connected and so contribute to the inertia which is independent of p. Good estimates of s/v for 2D and 3D are given in Ref. 18. Using s/v = 0.977 :t 0.010 and 0.85 t 0.011 in 2D and 30 respectively, gives 1117 a = 3.91 2 0.014 in 2D and '5 = 5.2 1 0.2 in 3D. This is in good agree- ment with our results in 20 and reasonable agreement with our results in 30 where the numerical results are less accurate because the lattices had smaller linear dimensions. The agreements between the predicted values of d and the computed ones is of course evidences for the existence of the cross-over scaling hypothesis. He now use the scaling property of the system to obtain the best value of 3. By using the equivalence of E ~ m2, equation (17) can be written as 5/2 -1 g(E) : E r(E/Eco) (20a) where r(E/Eco) is the cross-over scaling function and (2-s/v)v 2 co “’co lp - pcI (20b) Eco is obviously different for different p close to pc. For different p, g(E) vs. E curves are expected to be different. However, if we res- cale the E axis by 1'/Eco then all curves g(E) plotted against E/Eco for different p's must fall into a single universal curve if the scaling hypothesis exists. From equation (20a) we can see that if we divide g(E) by Ed/2 ’1 we will get the scaling function r(E/Eco). Since the p's are away from pc we expect some phonon excitations in the density of states curves. In the anomalous region, g(E) ~ BC“2 ‘1, so the scaling 1fl8 function r(E/Eco) has to be a constant in that region in order to have the correct g(E) ~ E d/2 '1 behaviour. 0n the other hand, in the phonon region g(E) Ed/2 '1 so r(E/Eco) has to behave as -ds/v (E/Eco)2(2-S/V) . Therefore we can write constant E > Eco r(E/Eco) ~ -ds/v (200) (E/E )2(2'3/”) E < E co co For convenience we define y = 5/2 - 1 (21a) and x = (2 - s/v)v (21b) so equation (20a) and (20b) become _ y g(E) - r(E/Eco) E (223) x Eco - Eolp - pcl (22b) In Fig. 3.7 g(E)/Ed/2 ' 1_ . - r(E/Eco) is plotted vs. E/Eco using the best values of x and y for p = 0.1175, 0.1180, 0.1185, 0.1190 and 0.1195 where p0 = 0.5 is the percolation threshold for the 2D square lattice. It can be 1‘19 seen that all curves fall into a universal curve extremely well. The procedure of obtaining the best values of x and y is as follows: we first plot g(E)/Ey vs. E/Eco by using rough estimates of x and y, so curves scatter a little bit. Then we fine tune x and y so that all the curves fall into a universal curve. The best values of x and y are 1.37 1 0.03 and 0.95 I 0.03 respectively. Therefore from (21a) we obtain a = 3.90 t 0.06. We see that this value ofd is closer to the predicted value of 3.91 i 0.04. The previous method of obtaining d is simply measuring the slope of log - log plot of low frequency part of spectra at percolation threshold and basically it only uses the informa- tion at percolation thershold. Clearly the latter method of extracting a is better because we have more information and we can fine tune the «mirves to result in the best value of 3. However, this method requires much more computer time. In 3D case we can not get as good results as in 20 because the system sizes are not large enough to produce smooth curves for different p's, but it can be done easily by using a. supercomputer. 150 TABLE 3.1 d 5 (This work) 5 = 2d/(2 - s/v) 2 3.90 i 0.06 3.91 z 0.0“ 3 5.8 1 0.3 5.21 1 0.2 Table 3.1 Comparing the spectral dimensionality 5 obtained directly 111 this work with the scaling relation (18). Values of s/v are taken from Ref 18. 151 * l i 1' 1 p=.475 0.6 ps.480 -‘ p-.485 p=.490 p-.495 CH>+ A O 0.3 0.6 0.9 1.2 1.5 Fig. 3.7 r(E/Eco) is plotted against E/ECO for p = 0.975, 0.H80, 0.U85, 0.M90 and 0.995 where p:0.5 is the percolation threshold. It shows good scaling behavoir. 152 IV. One Dimensional Case The scaling in 10 is rather different because f = (1 - p) ~ E-1 is critical. We start with the Debye form for the integrated density of states (See Eq. (15)). 1(u1) ~§ C(gl’d (ald de[£- ] (23a) 00 Setting d = 1 and keeping only the critical terms we have 1 E. - 2v m In») ~g m 11(5—1 CO ~ if; I Hifi’; 1 (23b) 00 CO ~ Rtfi 1 00 where R(x) = xH(x) is expected to tend to a constant for x large. In this limit we have I(w) ~ md which leads to d = 0 in 10. This result is expected because in the limit p = p0 = 1 we have a perfect superconduct- ing,chain which is completely rigid and thus has only one degree of freedom. This can be considered as a zero dimensional object for which 5 = d = 0. In 1D there are still excitations associated with rigid clusters of length less than g. Because the probability of neighboring clusters 153 having the same mass M is small, we can treat each cluster as an Einstein oscillator of frequency m2 = 51'. The probability of a cluster having mass M is M-i P(M) = (1 - p) p (293) where we have chosen the mass of a single site to be one, so M is an integer. The density of states P(M) lgfll p(w) m 80 E 3» M-1 9) p v p(w) = (1 2 = (1 - p) p(2V/m '1) 5% . (Zub) 0.) For (1-p) small, we can write p Z e-(I-p) and thus 1 - 9 e-2V(1-p)/m?_51 p(w) : p m which to leading order in (1-p) gives p(m) = “V(1-p)/m3 154 or equivalently in tight binding language 2 g(E) : 2V(1 - p)/E . (24a) The exponential can be replaced by unity as there are other terms 0(1- p)2 that we have neglected in deriving Eq. (240). This result is special to 1D. In higher dimensions if we tried to consider the superconducting clusters as Einstein oscillators they would have a frequency 2 niV mi = i— (25) i where ni = Z "ii is the number of surface bonds connecting cluster i to Hi other clusters bonds and M is the mass of the i cluster. However, for 1 n1 superconducting clusters we find numerically that H— ~ constant of order i 1, unlike in 1D where n = 2 always and H can take any value. In 1D, i i the likelihood of having two adjacent clusters with the same mass, and hence the same frequency, is small. Therefore the modes do not hybridize and remain localized with d = 0. However for d 12, the likelihood of adjacent clusters having the same Einstein frequency is high, as all the heavier clusters have essentially the same frequencies. Therefore the Einstein modes from the heavier superconducting clusters hybridize forming extended low frequency modes with 3 > 2. 155 In 1D we cannot do simulations at pc = 1 so we have worked at small (i-p). For frequencies u1< «$0 or equivalently for wavelengths A > g ~(1-p)"1 which is the typical length of a rigid cluster, we expect phonons with a constant density of states -1 p(w) ~ C(p) /f‘ . (26) Because the linear chain Just involves adding springs in series in the static limit, we have 1 ————- (27) /V(1 - p)E g(E) :.%; In Fig. 3 .8 we have plotted p(m) against 1» [or equivalently (fl? g(E) against JE] for p = 0.9, 0.99 and and 0.999. The results were obtained from chains of 20000 clusters (Nf = 20000) using a transfer 19 matrix technique. The results have been rescaled by plotting /E(1-p) g(E) against i/E/(i-p) to show the scaling behavior. In Fig. 3.9 we have plotted JE g(E) against JE for p = 0.999 and the two limiting -1 curves JEg(E) = [2n/V(1 - p)] and ng(E) 3/2 2y(1-p)/E for low and high energies respectively. 1 "g(E) [EH-p)] 156 0.15 O h—‘b O 0 . 05 - 1‘ \‘JR - \ WWW k_ 0 2 4 6 8 [Elm—p)]" Fig. 3.8 The density of states for a Hmukm1linear chain of SUpercon- ducting and normal bonds of p = 0.9, 0.99 and 0.999. The results are rescaled as indicated to show the scaling behavior. The energy is in units where V = 1. JEQE) 157 5 1‘.) Fig. 3.9 forms Eqs. V=1. 0.05 0 The linear chain results for p = 0.999 and the two limiting (2") and (27), are shown. The energy is in units where 158 V Conclusions He comment on our results compared to previous work on the "termite" diffusion problem. Alder et al.10 ~ tk with k = 1+s/[2v+(t-B)] which, using‘d = kd is at variance have predicted that with our results. However it has been pointed out previously that this result is wrong 12 and Hong et al.11 argue that there is no regime where ~tk because (r2) is dominated by diffusion on the superconducting clusters. Random walks on the superconducting clusters are related to eigenstates in the high frequency band which we are not considering. The lxnv frequency states in which we are interested are related to ran- dom walks on the normal clusters and their dependance on the normal- superconducting interface. These walks were studied by Coniglio and Stanley‘?0 who used scaling arguments to find (r2> ~ tk with k : 2/(2- s/v). Using 3 = dk this agrees with our equation (18). It is appropriate to use d = 3k because the number of sites on normal clusters is not critical. Numerically, it would be very difficult to focus on these walks and exclude the random walks on the superconducting cluster. In summary, we have evaluated the density of states of random su- perconducting normal networks at the percolation threshold and found the spectral dimension d which governs the low frequency density of states. Our values of 3 agree very well in 2D and reasonably well in 3D with our predictions for 5 using scaling theory. In this problem, as in the ran- dom resistor network problem, a is also related to anomalous diffusion but in contrast to the random resistor network problem, it cannot be ob- tained from numerical simulations of random walkers. 159 VI References 10. 11. 12. 13. S. Alexander and R. Orbach, J. de Physique Letters 53, L625, (1982). R. Rammal and G. Toulouse, J. de Physique Lettres, 55, L13 (1983). S. J. Lewis and R. Stinchcombe, Phys. Rev. Lett. 5g, 1021 (1989). S. N. Evangelou, Phys. Rev. B 33, 3602 (1986). Y. Gefen, A. Aharony and S. Alexander, Phys. Rev. Lett. 59, 77 (1983). L. Puech and R. Rammal, J. Phys. C 35, L119? (1983). P. Argyrakis and R. Kopelman, Phys. Rev. B 22, 511 (1984). S. Alexander, C. Laermans, R. Orbach and H.M. Rosenberg, Phys. Rev. B 28, 9615 (1983). A. Aharony, S. Alexander, 0. Entin-Hohlman and R. Orbach, Phys. Rev. B ‘31, 2565 (1985). I). Stauffer, Introduction to Percolation Theory (Taylor and Francis, London and Philadelphia, 1985). P.G. de Gennes, J. Phys. (Paris) Colloq. £1 C3-17 (1980). J. Adler, A. Aharony and D. Stauffer, J. Phys. A 18, L129 (1985). D. C. Hong, H. E. Stanley, A. Coniglio and A. Bunde, Phys. Rev B 33, 9569 (1986). F. Leyvraz, J. Adler, A. Aharony, A. Bunde, A. Coniglio, D. C. Hong, H. E. Stanley and D. Stauffer, J. Phys. A 12, 3683 (1986). R. Alben, H. Blume, H. Krakauer and L. Schwartz, Phys. Rev. B 12, 11090 (1975); M. F. Thorpe and R. Alben, .1. Phys. c g, 2555 (1976). 19. 15. 16. 17. 18. 19. 20. 160 M. F. Thorpe, J. Non—Cryst. Sol. 51, 355 (1983). M. F. Sykes, D. S. Gaunt and Maureen Glen, J. Phys. A 19, 287 (1981). The quantity f can be evaluated exactly in 2D for bond percolation on the square net at pc. The value of f = 0.0981 ..... can be ob- tained from H.N.V. Temperley and E. H. Lieb Proc. Roy. Soc. A 332, 251 (1977). This value was given incorrectly in Ref. 19. We have confirmed this value of f by directly counting the clusters in the samples studied in this work. A similar counting of clusters in 3D leads to the value of f = 0.27 for bond percolation on the simple cubic lattice at pc. S. E. Koonin, Computational Physics (Benjamin/Cummings Publishing Company, Inc., Menlo Park, CA, 1986). H. J. Herrmann, B. Derrida and J. Vannimenus, Phys. Rev. B Q, 9080 (1989) and J. G. Zabolitzky, Phys. Rev. B 39, 9077 (1989). M. A. Lemieux, P. Breton, A.-M. S. Tremblay, J. Phys. (Paris) Lett 59 L1 (1985). A. Coniglio and H. E. Stanley, Phys. Rev. Lett. _5_2_, 1068 (1989). APPENDICES 161 162 APPENDIX I Elasticity of Pure Honeycomb Lattice The elastic moduli of model (8) in part II of this thesis is derived below. Without any dilution (i.e. perfect lattice case) L1 = L J , the lattice spacing, and K11 = K. To avoid any confusion we change uiJ to 911 taureserve the notation u for x component of fiij' The elas- tic energy can be written as v = g X K(L - L0)2 + X K(L — L0) (fiiJ-riJ) + (13> <13) 1 2 ‘ 2 + 2 <§J> K{(1 — 1.0/LNiJ + LO/L (fiiJ riJ) 1 + ... (1) As mentioned in the thesis the strain energy can always be written as 1 S * 2 Z CaBYI EQBEY‘I (2) where 308 are the elements of the stress tensor and 00871 are the second order elastic constants which are directly associated with the elastic moduli. The quantity 608 is called the strain and defined by 163 EU 2 = -—9 a B = x y 08 3x ’ ’ B where ”(1 and x8 are the components of samll displacement 61.1 and posi- tion F. In two dimensions E -29. ,3 .22 xx ' 3x yx ' 3y .- -21 E -21 xy ‘ 3x yy ' 3y If x and y are the x and y components of the distance between site i and J then it is conventional [see Kittel] to define the uniform displace- ment fiij in two dimensions as So u = xexx + yeyx (3a) v = xe + ye (3b) 169 (3a) and (3b) are the uniform distortion on a Bravais lattice when a small force is applied. Now let us put a small force on the frame which holds the springs of a honeycomb lattice. As discussed in the thesis the displacements (i.e. the distortions) between any two nearest neighbor sites of the sublattice A and sublattice B [see Fig.2.2] can always be written as u = xe + ye + u' (9a) v : xe + ye + V' (9b) where x and y are the x and y components of the distance between two sites. Qantities u' and v' are the relative shift of the two sublat- tices in the x and y directions and can be determined by minimizing the elastic energy (1). In Fig. 2.2 the three neighbors of site 0 are labled by 1, 2 and 3. If the coordination of site 0 is chosen as (0, 0) then the coor- dinates of 1, 2 and 3 will be 0: (0,0) 1: (J13 a, - ga) 2: ( 2 a, - ga) 3: (0, a) 165 where a is the lattice spacing. So the relative displacements of site 0 to site 1, 2 and 3 will be - B. a . U01 - 2 aexx + ény+ U v =13 a e + as + v' 01 2 xy 2 yy u = - 13 as e + u' 02 2 xx 2 yx £3 a . v02 - - 2 a Exy + 2Eyy + v u03 = - aeyx+ u' - _ 1 vO3 - aeyy+ v (5) If we take the parallelogram in Fig. 2.2 as a unit cell then the elastic energy per unit cell for perfect lattice can be written as 0 - 1 1 V _ V0 + 2NdB [2(J3 ”01 + V01) + é(-/3u02 + V02) - V03) + + N I 8(u 2 + u 2 + u 2 + V 2 + v 2 + V 2) d 2 01 02 03 01 02 03 2 + 31 &(/3 u01 + v01)2 + &(-/3u02 + v02) + yogi} (6) 166 where B = K(1 - LO/L) and a = DO/L. W) is just the static elastic energy independent of 61.1. In (6) we have also used r A A 01’ r02 and r03 which can be easily calculated. Nd in (6) is a constant which is as- sociated with the area of unit cell. Now we minimizing the elastic energy (6) with respect to u' and v' i.e. 3V 55] = 0 (7a) 3V _ and 5;, — 0 . (7b) Substituting all equations in (5):hux>(6) and solving (7a) and (7b) we obtain a(£yx - Exy) “ ' 2(2a+B) ‘83) 1 0(Exx - Eyy) V ‘ ' 2(2a+B) (8") Now substituting (83) and (8b) back into (6) and compare (6) with (2) which can be written as V - S E + S E + S E + S E + xx xx yy yy xy xy yx yx Cxxxx 2 2 Cx x 2 2 --- (E + ) -l-l(€ + E ) + E + C E E + 2 xx yy xxyy xx yy 2 xy yx 167 C e e , (2') xyyx xy yx it is straight forward although a little bit lengthy to show that _ _ _ B _ _ K T - SXX - Syy -/§ - (1 U)/§ c _ c _ a2 + 5aB + 932 _ u ‘ 3" K xxxx ' 11 ' 2/3 (28 + a) ' 2 - n 273 c _ c _ a2+ aB _ __g__ x xxyy - 12 ' 2M; (28 + a) ' 2 - n .273 c _ 5 _ 932+ 3aB _ (1 - x xyxy ’ 99 - 275(28 + a) ' 2 - n '273 - - GB _ (9—nl(1-n) K nyyx ’ Cuu ' 2/§ (23 + a) ' 2 - r1 273‘ In the second column above we have used the conventional notation C xxxx = C11’ cxxyy = 012, nyxy = C99 and nyyx = C99° In the fourth column we have used a = Kn, B = K(1 - n) and n = LO/L. We have also considered the sysmmetries of the honeycomb lattice to reduce the number of inde- pendent elastic constants. 168 APPENDIX I I Percolation Thresholds of Stretched Spring Model on Honeycomb Lattice We give a derivation of the percolation thresholds p for the stretched spring model in part II of this thesis. The method here is the same as that of Tang and Thorpe but it is genralized to the honeycomb lattice which is a no-Bravais lattice. A honeycomb lattice is composed of two interpenetrating triangular sublattices. We start by writing down the forces on the sublattice A and sublattice B. [see Fig. 2.2] Using the elastic energy (8) in part II of this thesis we have AA = - gr?“ = - 3 g (GA - 65) - (1)8 [(EiA - 58).;A51;A5 (1a) #8 = - 3&3 - 3 g (uB - as) - (1% [(6B — 68).},3511338 (1b) Here the subscripts A and B are used to denote the quantities of sublat- tice A and B respectively and the symbole 8 is used for nearest neighbors of A and B. Now substituting the Fourier transformations of displacemnets and forces into (1a) and (1b) and have [for each k] where Ak Bk Ak Ak C1 3’ ll Zl—b WM H» D 11 Zl—s WM "216 w 11 ZI—n WM 169 (2a) (2b) (3a) (3b) writing the results in the matrix form we will - DAA(k)uA - DBA(k)U Ak ' ’ DAB(k)uBk DBB(k)uBk Bk Bk ' (9a) (9b) 170 Here the superscripts x and y denote the x and y components. Subscript k denotes the kth Fourier component. 2 X 2 Matrices DAA’ DAB’ DBA and DBB are the following. D _ 3(a + 2B)/2 o D _ 'a1 ’a2 AA ' 0 3(a + 23)/2 AB ' -a -a a a D - ’ai '32 D _ 3(a + 2B)/2 0 BA ’ * * BB ' 0 3(a + 2B)/2 -a -a 2 3 . . ik1/3a 3 ik2/3a 3 ik1/3a + ik2/3a and a = Be + (B + a)e + (B + a)e 1 9 H a - a“; eik2/3a (1 - eik1/3a) 2 ‘ H ik /3a ik J33 ik /3a + ik /3a a 2 a 1 2 a3 :(a + B)e + (B + a)e + (B + 9)e n where a is the complex conjugate of a and a = KL /L, B = K(1-LO/L). k 0 1 and k2 are the k's components in the 61and 62 directions. [see Fig. 2.2] From (9a) and (9b) we can solve uAk and uBk in the following fOrm D-1F - D F — (—D"B + 0‘10 )u (5a) AB Ak BB Bk ' AB AA BB BA Ak 1F D F = (-D - 1 AA Ak ‘ BA Bk "B + B‘ D D AA AB BA BB)uBk (5b) 171 By defining -1 A - DAB (-D -1 -1 -1 ABDAA+ BB BA) -1 BB - -1 -1 ('DABDAA* DBBDBA) (6) -1 -1 AA (‘DAADAB* DBA BB) 0 n U -1 -1 -1 -1 BA ('DAADAB* DBADBB) ' D n :3 (5a) and (5b) can be written as UAk : AFAk + BFBk (7a) uBk : CFAk + DFBk (7b) Here symbols for matrices A, B, C and D shold not be confused with the subscripts. After removing a spring the two sites move along the direction of the missing;spring (see Fig. 2.9). In order to calculate the effective spring constant of the resulting network in this direction one needs txi pull the two sites back along the same direction. So the net force is A in the direction of unit vector rAB‘ Thus 13A5 f rABSAa (8a) F35- r 'BA5B5 = —r ”A3535 (8b) The Fourier transformation of FAk and FBk are ik-R fiAk = X PAse A6 (93) 6 ik-R er . z is“. 35 . (9b) 5 Substituting (7a) and (7b) into (2a) and using (8a), (8b), (9a) and (9b) we have ’ 1 ik-RA * uA ‘ R E 8 ”AK iE-A -iE-B -iE~A . _ f A A B _ N {X e (Ae - Be ) rAB (103) k5 Similarily 1kofi -1Eofi ~1koR . + f B A B uB - n g e (Ce — De ) rAB (105) 173 $0 + f “ AB ‘ éuAB - rAB.(uB - uA) - fl Eg rAB[-(A + D) + Be + Ce ) ]rAB In the summation RAB = R8 and rAB : 00 Therefore following Tang and Thorpe , iE-AAB -1EoAAB 1 l p = a = 2%; {X Tr {551-(A + D) + Be + Ce (12) k6 174 APPENDIX III Published Paper 175 PHYSICAL REVIEW B VOLUME 37, NUMBER I0 1 APRIL 1988 Spectral dimensionality of random superconducting networks A. R. Day, W. Xia, and M. F. Thorpe Department of Phyrr’cr and Astronomy and Center [irr Fundamental Materials Research. Michigan State University. East Lansing. Michigan 48824 chceived 2 March l987; revised manuscript received 2 July I987) We compute the spectral dimensionality 3 of random superconducting-normal networks by directly examining the low-frequency density of states at the percolation threshold. We find that a? -. 4. HO. 2 and 5.8 l 0.3 in two and three dimensions, respectively, which confirms the scaling rela- tion d=2dll2--s/r'). where .r is the superconducting exponent and v the correlation-length ex- ponent for percolation. We also consider the one-dimensional problem where scaling arguments predict, and our nunrerical simulations confirm. that a .-.-0. A simple argument provides an expres- sion for the de'nsity of states of the localized high-frequency modes in this special case. We com- ment on the connection between our calculations and the "termite" problem of a random walker on a random superconducting-normal network and point out dilllculties in inferring d from simulations of the termite problem. I. INTRODUCTION 1 he question of dynamics on dilrrte random systems at the percolation threshold has been the snbi'ect of consid- erable attention in the past live years." The infinite cluster at percolation is a fractal, so that the excitations are qualitatively diiTcrent to those on a regular Euclidean lattice and are governed by an additional dimension, the spectral dimension :7.” liccause of the similarities be- tween the dillusion equation and the equations of motion for a system described by a tight-binding Hamiltonian (or an elastic Hamiltonian) it is clear that there are two ap- proaches to evaluating (7. ()nc approach is to examine the Ionifrcqrrency dyrrorrri- col response of the system described by the llamiltonian . t l (1.)) where the I". are randomly distributed on the nearest- nrighbor bonds according to the probability distribution rrr;,1~p5(r-'_r,,) (11416113,) . 121 and we restrict the sum to sites on the infinite cluster. The 0,. a,’ are eitirer llose or Fermi operators (the results are the same when only one particle excitations are con- sidered) amt the diagonal term is present in (I) to preserve translational invariance. 'l’lris dilute tight- bindiug llamiltonian describes the physics of several dill'crcnt random systems; dilute magnetic systems, ran- dom resistor networks, and random scalar elastic net- works. liy scalar elasticity we mean for example dis- placements perpendicular to the elastic network. The lmv-frequcncy density of states at p, is predicted to behave as gtI-Ii~ 1;" 4'”. (.1) where J is the spectral dimensionality that can be ex» pressed :rs' I 37. = ___-e_._ (o and t is the conductivity exponent. The quantities I} and v are the geometric exponents“ which govern the proba- bility of being on the infinite cluster. and the correlation length, respectively. The [racial dimensionality is J=d ——/3/v. Direct numerical evaluation of the density of states in two dimensions' on dilute spin systems has -coniirmed Eq. (4) and that d is close to t which Alex- andcr and ()rbach' have conjectured to be the value of d in all dimensions except one dimension t 1 l1). An alternative approach is to study the anomalous dillusion of a random walker on the infinite percolation cluster. Scaling argumentsm" suggest that the mean square displacement (r’) scales with time as (r’)~r* (51 and where t. B, and 1' were defined previously. liy relating the probability of return to the origin with the density of staies,‘ one obtains r7 =r7k in agreement with liq. (4). Numerical simulations of random walks on dilute sys- tems are a good indirect way of obtaining the spectral dimensionality :7 using liq. (S) and the relation (7:417. Valrre old obtained in this way" agree with direct simula- tions‘ of the lmv~lreqncncy density of states at pr. 'l‘hc randomly dilute problem discussed above can be considered as a random system of strong boruls l', and weak bonds 1'", in the limit 1", —- i’, l',,.—-~t). ‘lhc other limit, V, ~ em. l’m-il’. is also of considcrablc interest and is ilrc subject of this article. It would describe a ran- dom superconducting network or the elastic properties of 49.10 (€1l988 'l Ire American Physical Society 3‘? SI‘IECI'RAI, DIMENSIONALI'I‘Y ()l’ RANDOM . . . a system with rigid grains weakly coupled by soft springs. lo the general problem. the density of states is divided into hands: a highlrcqueucy banrl associated with inter- nal modes of the clusters of strong homls and a Iow- frean-ncy baud associated with viluatiorrs of the sott re- gions. In the limit I’wurll, V,—»I', the low~frcquency band bccorrres a 5 function at the origin and the spectral dimension describes the low-frequency edge of the high- lrcanney hand. In the limit I’, ~> l’, l’, -.. on, the high- frequency band is driven off to infinity and the 6 function broadens ‘I here will be a different spectral dimension as- sociated u itlr this low-frequency edge. In diffusion language. the superconducting limit is the problem of ternrite diffusion. which Iras been the subject of considerable work'" ti and sonrc controversy. In "i ll“ iple. it slrorrlrl provide irrforrrratiorr about the density of states amt the spectral dimension of for the problem bill in pram-tire this is not the case because the ramlorrr walk is completely dominated by the distances rrrovcd on the su- pereorrrlrrerrng clusters. As has been pointed out by llong et of. " there is no region where (r’) ~t‘ because even at time 7cm. (r’) ~ ‘6’ n". the average radius of the su- 'rr'rr'ntttlttt‘littg clusters. 'I Ire low-frequency excitations we are interested in are associated with the small, time- dependent contribution to (r’) superimposed on the rliverpent _U "“1 ('learly this is impossible to obtain nu- merically and we therefore believe that direct simulations of the low frequency response is the best way to evaluate the sper tral dimension. II. IIENSI I Y OI’SI‘A I IsS In this paper we obtain the density of states. arnl hence the spectral dimension d for the srrperconrlrrr-tr‘rrg (use. llre system is described by the llanriltorriarr lll brrt with the interaction strength distribrrtiorr I't l'”) ‘pfil I", .- r") I II plbf V~ "'I‘ , I6) uhcre l‘. er} itt the superconducting limit. 'I he equa- tions ot motion for the llarrriltorriarr (ll are I‘f‘ I“ “If " Z'IIH'I "4",, (7, f f u here the r, arc the amplitudes of the “use frrrrctiorr on the l site. luo sites I. arrrl I. connected by a hood I", '2 m. have the some amplitudes so that r, -.c, can he considered as a single srte of ruass 2. 1 his process is repeated for all so- percorrdrn ting borrrls so each snpcrcondrrr ting cluster has one deprce of freedom and a rrrass M equal to the rrrrrrrbcr of sites in the cluster. We consider a site connected to no superconducting bonds to be a cluster of sire l llre equations of motion for the clusters so rletirred are MM, 't --- Z I'nufr-t --~r,l , till I I where AJ are now rlnstcr indices and n“ is the total number of normal bonds between the Ir and the Icluster. lly defining trlf‘ ”we, erg. we have the Iferrnr'tr'rm equa- tion of rrrotiorr 176 49M '* '” “’Z-Z‘Z;I a (9) which we solve numerically using the equation of rrrotiorr leclrnique.‘ if the original system lras N sites. there are only Nf clustersH so we have reduced the rrrrrnbcr of degrees by a factor of I. For smallp fr|_.zf"... . (Ill) 2 where z is the rnrrrrber of nearest neighbors. Note that f decreases to a very small value at m. In the square net and simple cubic lattice we find it is 0.098 and ".27. re- spectively." 'l‘lrus setting I", strictly eqrral to infinity reduces the ell'eetive’sizc of the system, especially in II). although tlrete is sonre initial tirrre cost in identifying all the clusters. 'I his advantage cannot be arhreverl by set- ting f’, large and taking a limiting process. Note also that f is finite at pr for d _> 2 brrt it is zero in II) which affects the sealing argrrmcrrts for d as we discuss later. the equations of motion were integrated forward in tirrre from 0 . 1' using standard methods. H 'I Ire first few points were obtained using a fourth-order Runge-Kutta algorithm and subsequently a fourth-order Adams in- tegratiorr scheme was used. "' ‘l he density of states glffl is obtained from l’ourier transforming the result of the time integration, . i "' o r'r sir-.1: ii}. 23M, f” mm. (file'rft. on This integral is converted to one over positive times only. " Random phases" are put on the elrrsters so that only the ou-site terms required in It It are retained in the limit of a large system. 'I his is a standard way to obtain densities of states in the cqrration rrl rrrotiorr technique. In Fig. l we show results for the square lattice at p, -.- ;. the results are from averages over 25 sarrrples of -_ — l " II I l l on I V e r r r l' ~- tt to / 005 / a", 0 o i ,. " a (”II t r r e r r 0' ~-’ err m r 2 oz 0,, .....-” l_-__--...l _ -.. I ........ l. - l - ,. o z r a a to re l'l lilti. l. ‘Ihe density of states for a random square lattice of superconducting and rmrrrrnl bonds at p, ~ : ‘l he insert shows the low-energy data on a log-log tbase lftl plot, llre straight line drawn is a least squares fit to the data points shown from which d is obtained via liq HI. 'I Ire energy is in nrrits with l' -. l. 177 4‘”? Ion -' ttttt networks. ltecausc of the reduction in size due to the I factor only ~ ttltltl anrplitrrdcs had to he rnorri- torcd. the insett.shows a log log plot at low lrcqucuces “ltlr'lt gises a slope ot' l.tlS t tl.t leading to rl: it till horn liq Lil. ln Hg. 2 we show sirrrilar results for the simple crrhic latticefl at p. - 1124'”. the results are tronr rut-rage over 2‘ samples of It X22.-£2l networks. He- came of the reduction due to the f lactor only ~2lttltl amplitudes had to he nronitored. 'l he insert shows a log- log plot at low frequencies wlrictr gives a slope of I ‘t I H, t5 arnl hence cl .2 5.8 t it. .t. We exaruirred a single sample with about to tirrres as rrrarry sites in truth 2|) and ill amt found no significant changes in the values of il. l he sharp features at higher energies in Figs. l and 2 will he discussed in a subsequent puhticatiou. H lll. SCAHNG ARGllMliN'l’S the scaling argument for (l is given trelow following the ideas developed for the ditrrte case.U It is easiest to use the language of phonons. i.e., vihratious on a network ol rrorrrurl and irrlirritely rigid springs. The resrrlts apply equally well to any system descritred by the tight-binding ltaruittouian ill with the probability distribution (6). (.‘tose to p. amt for systems ot size I. much greater than ttre correlation length g ~lpr ——pl' ' we expect the system to appear homogeneous with an ellective elastic modulus t' - tp‘ — pl ' where s is the superconducting exponent. tlecause the system is fully connected, all sites contribute to the inertia in the long-wavelength limit, the total mass M -- I. ‘ and ttrrrs the mass density p is noncritical, i.e., it is independent of g. 't he low-frequency excitations are phonons with a sound velocity t."t§l~t'/p~§"". tl2l lor wavelengths X~§ the fractal structure of the lattice he. oures irrrportarrl. The excitations are no torrger plur- Irons and we have a fracton-phonon crossover. 'l'he crossover lrequency is given hy (’t l . (it ~ —..§.~§ l2 slit/2. (l3) é .---.-- ...._____-.'.____- -...._.__-_'-..._...___.-._._._._'-.___.. tron 1‘s 9 r'--"—"V‘ “- 03 00]. ' ‘ r‘ 02 Ill .1 r r r l _--_. a .... V Dll l I 3 Qt 0., l.-./- ....... I-____.____- I _ I l0 IS I lti. 2. Same as Fig. l except for a simple crrhic lattice at r, ll 24'”. the insert shows the low energy data on a log log tlr.r'.c till plot. 'l he energy is in units where l' : l. A. R. DAY, W. .‘s‘tA, ANI) M. F. lllORl‘li 37 Using a Dchye-type theory we can write the rrornralired density ol states for the phonons as r! 2!?" ..." '/~/. 'lo allow or the fracton phonon crossover we introduce a scaling function lr tar/cowl into (N) and otrtaiu cl 1. I’llrll‘Clgl d 2'" "l" 'hlltl/ltl", Irlml~ ('tgr " on t/Nf. usr l-‘or cl 2 2, f is noncritical and so the critical hehasior ol ptml is (ll pilots; “Wm." 't. ",I 0 w ’1: It 2r- - slmd ”I “_a_ ll6l (cl ~( (I! where we have rrsed (Ill amt (H). In the low-frequency (ptronorrl regime x =ar/mm << I and Ir txl- oconstant. so that p(ml~m" " as expected. For ar>nr the tractorr regime, we‘rewrite US) as wJI/lh-ll (0 ur p(a”~wllll/‘I~I/Vl-l ...... h _— a)“ (0,", - (ll Ma" 'r . rm (0", wtrere rtxl=x “"""“"lrtxl is expected to tend to a constant for x large because ptmt should be independent of g or equivalently of to", in this tirrrit. Equation tl7l detines the spectral dimension, tl8l We have discussed seating relations for vibrational modes. Exactly parallel arguments can he nrade in all di- rrrensiorrs for tight-binding llmniltoniunr where 15' ~ ml - elk ’ leading to gtlil~ Is" "’1. on Corrrparing liq. (lit) with the result for dihrte systems liq. (4). we notice the following dill'erences. (a) The rrrrrrrerator involves «I rather than (7. ‘l he exci- tations in the superconducting case can explore the whole system whereas in the ditrrte case. liq. t-tl is ohtairred lor excitations on the infinite cluster. It the finite clusters are included in the dilute case.“ then (7 is also replaced by rl in liq. t4). ltrl The conductivity exponent t in liq. Ht is replaced try the srqwrconductirrg exponent - s in liq. ttttt. tcl 'l tre exponent I! does not appear in liq ttttt. tn the dilute system. the prohahility ol' treiug on the irrlirrite clus- ter scales as in p. l". to the low-[requency dyuarrrical 178 37 SPECTRAL DIMENSIONALITY 0F RANDOM . . . TAlll F. I. Comparison of the spectral dimensionality :7 ob- tained directly in this work with the scaling relation. Values of s/v are taken from Ref. t8. M .7---__.fi.—_s ”-.--.. -. .._~ _._..+._._....-«._..._—__ . d a (This workl ‘_ - 2 4.t 1:0.2 18" I004 3 3.8103 5.2 10.2 I 1 response of the infinite cluster, this is the inertia ternr and appears in (4). In the superconducting systenr, all sites are connected and so contribute to the inertia which is irr- deperrdeut ofp. Good estimates of s /v for 20 and 3t) are given in Ref. l8. Using s /r'.-.-fl.977.l.0.flt0 and 0.85.t0.f)4 in 2D and 3f). respectively, gives (7 =3.9t_l 0.04 in 2!) anrl r7:- 5.) t 0.2 in JD. We see from Table I that this is in good agreement with our rcsrrlts in 2t) and reasonable agreement with our results in JD where the rrrrnrerical re- sults are less accurate becarrse the lattices had srrraller linear diruerrsiorrs. l\'. ONE DIMENSION: A SPECIALCASE the scaling in It) is rather different because I r-(l rl~§" is critical. We start with the Debye fornr for the integrated density of states [sec liq. IlSl] 2 If mm or Itr.rr~-§jct§l"’l!’— lm"11[v"’-l. (20) Setting if <- l and keeping only the critical terms we have “""‘-§' """"...Ill 5". l tr)“, ‘ _“l- l" i .9. Ir)", (rlm ~R 1" l. (2|) mm where R txl- .vlllxl is expected to tend to a constant for x large. In this tirrrit we have [lath-tr)" which leads to rl - ft in II). This result is expected because in the lirrrit p -p, - l we have a perfect superconducting clrrrirr which is completely rigid and thus tras only one degree of free- drrrrr. 1 his can be considered as a zero-dimensional oh- jcct for which if t-rl PO. In If) there are still excitations associated with rigid clusters of Icngthless than 5. lteearrse the probability of neighboring clusters having the same mass M is srrrall, we can treat each cluster as an Einstein oscillator of frequen- cy mi ll'/.tl. the probability of a ctrrster having mass I" is PM“ U -/'lp'" '. (22) where we hare chosen the mass of a single site to he one. so .t! is an integer. 'l he density of states 4933 Plft)7’—’p(M) rial" r (r) so )lril“f l l M '36-’2- I t " l’ I’ V =( I --ptp”"""" “4i, . (23) trl For (l—-pI small. we can write [9 ze ” "' and tlrrrs plmlzi:B-e ~1l'fl-pl/ar’i': ' P or“ which to leading order in f l -- pl gives plrrrl=4Vl l —pl/rrr“ . or equivalently in tight-binding language grrn=2m u-pl/EI. on The exponential can be replaced by unity as there are other terms Oil-17)2 that we have neglected in deriving Eq. (Ml. This result is special in It). In higher dimensions if we tried to consider the superconducting clusters as Einstein oscillators they would have a frequency rr- l’ ("31.-A?" , (25) where n, = 2,” m, is the nurrrher of surface bonds corr- rrectirrg cluster i to other clusters bonds and M, is the nrass of the i clrrster. However. for superconducting clus- ters we find numerically that n,/M, ~const of order I. unlike in II) where n,=2 always and M, can take any value. In It) the likelihood of having two adjacent clus- ters with the same mass, and hence the sarrre frequency. is snrall. Therefore the modes do not hybridize and rerrrairr Iocali7ed with (l r-fl. ltowever for rl ,5 2. the likelihood of adjacent clrrsters having the sarrre liinstein frequency is high. as all the heavier ctrrsters have essentially the same frequencies. therefore. the liiustcin rrrodes from the heavier superconducting clusters lryhridi7e forming ex- tended low-frequency modes with (l > 2. to II) we cannot do sirrrrrlatiorrs at p, '— I so we have worked at srrratl ll-pl. l-‘or frequencies or ("’m or equivalently for wavelengths A>§~It~ pl ' which is the typical length of a rigid cluster. we espect phonons with a constant density of states [tlml~('lpl '//. rzm llccarrse the linear chain just irrsnlvcs adding springs in series in the static limit, we have " . l tell-J: --»-> I '-. t27l 2n [hurtful-:1" In lig. .l we have plotted ,rtml against m [or equivalent- ly t'lz'gflz'l against \' l-.' I for p "it“. 0.9". and ".000. 'l he results were obtained from chains of 20000 clrrsters 4934 p=0.900 pnO 990 pr0.999 [K(1-p)]‘f' g(E) two-p)]” FIG. 1. 'I he density of states for a randorrr linear chain of su- perconducting aml normal bonds of p :09. ".99, aml 0.999. The results are rescaled as imlicated to show the scaling behav- ior. The energy is in units where I' = I. le :r2fltlflfll using a transfer matrix tee/grimif The results Ira‘yeébcerLrescalcd by plotting l”l'.’ll-—p)xll'.'l against V E/ll pl to‘show the scaling bilravior. In Fig. 4 we have plotted Vngl-Il against Vl'.‘ for p =0.99‘) aml the two limiting curves v‘Egmzlzm/WITFH ' and Vl-fgflfl: 2m I flow-"3 for low and high energies, respectively. V. CONCLUSIONS We corunrent on our results compared to previous work on the "termite" diffrrsion problcrrr. Alder et ul.'" have predicted that (r!) ~ (I with I.- =-l ls/I2r I-lt -flll which. using cl =rlk is at variance with our results. How- ever. it has been pointed out previously that this result is wrong” and Ilong et al." argue that there is- no regime where (r’) ~ t‘ because (r’) is dominated by diffusion on the superconducting clusters. Random walks on the superconducting clusters are related to eigenstates in the high-frequency band which we are not considering. The low-frequency states in which we are interested are relat- 179 A. R. DAY, W. XIA, AND M. F. TIIORI‘E _31 V’E g(E) FIG. 4. The linear chain results for p :0. 999 are shown as a solid litre mrd the two limiting forms Eqs. (Ni and (27) are shown as dashed lines. The energy is in units where l" a: I. ed to random walks on the rrorrmrl clusters and their dependence on the nornral-supercouducting interface. These walks were studied by Coniglio aml Stanley” who used scaling arguments to lirrd (r1) ~t" with k =2/(2—s/vl. Using d==dk this agrees with our Eq. (l8). It is appropriate to use 3=dk rather than «7 =Jk. because the number of sites on normal clusters is not crit- ical. Numerically, it would be very ditlicult to focus on these walks and exclude the random walks on the super- conducting cluster. ' In summary, we have evaluated the density of states of random superconducting normal networks at the percola- tion threshold and found the Spectral dimension r7 which governs the low-frequency density of states. Our values of (7 agree reasonably well with our predictions for r? us- ing scaling theory. In this problem. as in the random resistor network problem, i? is also related to anomalous diffusion but in contrast to the random resistor network problem. it cannot be obtained easily from numerical simulations of random walkers. ACK NOWI .I'ZI)GI\I I'ZN'I'S We would like to thank C. .l. I.anrbert for useful dis- cussions and I). C. lloug for bringing Refs. II and I2 to our attention. This work was srrpportcd by the National Science Foundation under Grant No. DMRtlll‘lt-rlo and by the Center for Fundamental Materials Research at Michigan State University. 'S. Alcsander aml R. Orbach. J. Phys. fl’arisl I.ctt. 43, [625 ”982). 2R. Rammal aml (i. Toulouse, J. Phys. (l'arisl I.ctt. 44. LI] "981“. "S. I. Iewis and R. Stinchcombe. Phys. Rev. I.ctt. 52. I02l "98“; S N. Iivangelou. Phys. Rev. ll 33. 36le ”986). ‘Y. ('iel‘cu. A. Aharony, and S. Alesarrdcr. l'hys. Rev. I.ctt. 50, WNW”. 3L. I‘uech and R. Ranrrnal. I. Phys. C 35. H I‘ll fl‘ltl.ll; l'. Ar- gyrakis and R. Kopclrnan. Phys. Rev. ll 29. Sll ”“8“. 0S. Alexander. C. Iacrmans, R. Orbach. and II. III. Rosenberg. Phys. Rev. ll 28. 46l5 ”98.“. ’A. Aharony. S. Alexander. 0. Iintin-Wohlman. and R. Orbach. Phys. Rev. I! .‘II. 2565 tl‘)85l. 'l). Stauffer. Introduction to Percolation I'lrmry t 'laylor and Francis. I.ondon, I985). 180 37 SI‘IECI RAI. DIMENSIONALI I'l' OI: RANDOM . . . "I‘. (l clr ( 3mm“, I. I'llg‘s. ll‘atisl Cnllnq 4|. (Ll-I7 ll‘lRlll. "'1 Acllrl. A. Almmny, and I). Slaull'ct. J. Phys. A l8, I.l2'l ”“8“ "I7. ('. llnnp. lI. l3. Slnnlcy. A. Cmnigliu. and A. llumlc. I‘luys. RH. II 33. 4564 ”WIN. ”l7. I('\\I:U. .l. Atlln'. A. Almmny. A. llumlr. A. Coniglin. l). C. llmlg. II. IE. Slanlcy. and I). Slnull'cr. J. Phys. A I9, 368.7 IInRIII. "R. Allvrn, M lllmnc. ll. Kmlmm'r, nml I.. Sclnwmlz. I’loys. Rn. ll IS, 40'!" "075); M. I’. ‘llmlpc «ml R. Alben, J. Phys. C 9. 2‘55 ll'776l. "M. I". llmupr. I. Nun-Cry“. Sulills 57, 355 "98.". "llu- mummy [ can he evalnnlcd exactly in 2|) lm lmml pct- cnlnlinn nu Ilvc squat: ml at 11,. Ill: vnlm: nlf r—IUNRI. . . can lu- nlvmim'cl l'mm II. N. V. 'lrmprtlty and IE. II. Lid», l‘mc R. Sue. anIon. Set. A 321. ISI "977). 1 his value was 4935 givm im‘mtcclly in er’. l4. Wt lmvc cunllmml llll‘ ”mm M f by (lirrclly cmmling Ilac clnslrls in Ilu: sample: mulircl in this wmlc. A similar ("muting uf ('lIKll‘IS in H7 lrmk In Ilse value "I [-4127 Im lmml pcncnlalinn mu IIIC simple culvic lallicc M p,. . "‘5. I3. Kmmin. ('nmpmmimml l'hyu'n lllcnjnmin/Cummiup, Cnlilmnin. l‘786l. "R. llmlcm. (‘. I.:Iml»cvl, A. R. Day. W. Xizl. mul M. I’. llumrc lummlvli‘lmll. "II. J. Ilcummm. II. l7cui¢ln, mul I. meimcnm. I.I|)S. Rev. II 30. vlllllllll‘lll'll; I. (I. anmliltky. ”ml. 30, 4077 ll‘lS-ll. "M. A. I.cmimx. l'. lltclmu. A.-M. S. ‘l tcmlvln). .l. I‘llys. ll'misl I.cll. 46. H lI'MSl. ' 2"A. (‘cmiglin mud II. Ii. Slanlcy. Phys. Rev. I.cll. 52, IUM (I984). "‘fimilnjllullwfllflr