o . 111111111113 11‘ . .. 4‘. ‘ 'I I 3:. f ‘ . (:{Hl I1 bq’q; ”31111111131111 1.1 1 ;""1: BM 1?“ l 1 . 1 1. ‘ 'Ll-‘N' ’Rufij .1 1. 1.. 11.- ”1‘1; h";- I . v 4 .. .1 ‘i;rL TV.’ 111‘; ' 111111 I.4 KflIM ”1" 1 5 4141‘»..1.1 1-- 111111111.) 91 “1. k9: "[149“ 5%):le n1: ‘3‘” d I _ I '1- >2- “A... 7' w b .r A ’ ~_ 9 ‘ "‘3 . ., . . g . . .‘ -‘éggm_ . - - ‘. I a n " ‘ It- " IU—J ‘ :2: _. . - , _ L“ V . o c __ l u‘ a. ' f"- - f . ' ... ’ ". ' '.....: - 1’. 4 '4. ‘— ww- ‘ o - __4~ I .7 ' . ‘ . “-zF-UUEJ-‘z'l‘, _ < I ‘ h . I . —-”‘ ’77“. r41; J v -7 c ISA-‘5‘: ‘ "J; H. . 1;: .jzuoyfm . . . ' 1: a r- u . _ . ‘ ’I ’ ‘ - . .. ‘41 .42.?2‘ a - Jazz-7'- s: .7‘ .. n ""c.—- o-L'v‘ .— ., c . >- . . “'u. o . '4 ‘ . . c .-. 6'..‘z-l . .- ‘ . . ' . l0. . v - . ‘ . I I - 0' . . Z- a? , .— .WWV‘ .. f.‘ . “—2, ' 2‘2:- '. ‘$‘~_—:- — . -35? .-_-_-:£ 4;? _ of?! - NW.‘_ “ "~ ~ -‘- . "3. . :7. $40- 2% 11" ’ 111...; «.1th "' '3 1 . J’ £- ‘5 5-43? _, ”via: ' . a ‘1‘ . ‘3 :5..,.__:§; 11‘? 1’ ““"“~a“\" 1111711“ ~ :33 ‘j (5:21!!! I? 2.51"? If" 11? 31??“ fifivf 1H ..‘ ‘o 5;"! ‘ 155;.1‘ ; .1! ~"‘""-1"‘; 11.111111 1’; I." c J"? V Eggs 7. YJW‘Izfi1‘Eifltcfi-BF 2‘ ,' Jilin— .jlrlfllff3‘ "‘1 (51:). 5-. '1; ‘ "3"“ :I‘:'- 1 . . «(31‘ L ‘ . ‘ . - “-121 lfilf’”; 1;: . ‘1 1.1:)? 1.111.111} u A E“ 1 .931... "3 £11511 .1 .. _L _._.‘.. -" w ‘I‘ 5*. 3.2“! 12'." ~‘ I" 4...? .» . .311. _ If?! (.3 .... l' , 7 V ‘ 'A l‘):.';9 ‘1‘)"3LIE‘7‘Jg: -5113}. ";;‘ .' ‘ ‘ 11. I :3." ‘. 5' .43: , «$3 . 1w _ .'~ ‘ .F .WJTLE - “31%;; I . 1‘ .3." L- n . ‘ . .. . u .. 1 ‘. ,~ . o ....,.¢“ ' ‘V—-uo -k- , ' .- _ . a c -. 7'HE'»~;u LIBRARY , Michigan State- University Thisistocertifythatthe (fiawflafionenfiflad Elastic Behavior of Noncrystalline Networks [newnufltw Wenqing Tang hasbamammxpmdtowaflkfhflflhmmn of the requirements for Doctor Physics dqgeehl MAM Major profit»; j July 8, 1987 I)ate MS U it an Affmnaxive Action/Equal Opportunity Institution 0- 12771 MSU LIBRARIES w \— RETURNING MATERIALS: P1ace in book drop to remove this checkout from your record. FINES will be charged if book is returned after the date stamped be1ow. _/ L“ ‘/ ELASTIC BEHAVIOR OF NONCRYSTALLINE NETWORKS By Henqing Tang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Physics and Astronomy 1987 ABSTRACT ELASTIC BEHAVIOR OF NONCRYSTALLINE NETWORKS BY Henqing Tang It is shown that the random resistor problem can be mapped on to a related network of Hooke's springs of natural length zero stretched on a frame. The conductance of the network is equivalent to the pressure on the frame. This new viewpoint leads to a useful visualization of conduc- tivity in random networks. The mapping can also be used on tight binding Hamiltonians. Chapter 1 presents the general formulation of the mapping, as well as the results for the conductivity and superconductivity of random networks in two dimension. In chapter 2, the method is applied to the Penrose Tiling network or two dimensional quasicrystal. The advantage of the mapping is explored to establish a sys- tematic way to calculate the diffusion constant. The result is compared with the diffusion constants given by other techniques, i.e. the Equation of Motion technique for the scattering spectrum S(q,E) and density of states and a direct evaluation of the diffusion constant by examining Random Walks. Chapter 3 extends the underlying ideas of the mapping to a deeper level and proposes a "tennis racket" model which bridges naturally the conductivity and elas- ticity percolation problems. A continuous phase diagram is obtained. The results at two extreme cases corresponding to conductivity percolation and pure central force percolation respectively agree with the previous research. The tennis racket model contributes much more enriched information to the macroscopic elasticity of noncrystalline materials. ACKNOWLEDGEMENTS First of all I would like to thank my advisor Professor M.F.Thorpe for all of his careful patient guidance over the last three years. His knowledge and insight helped guide me through this wide assortment of problems to produce some coherent results. I am indebted to Drs. H.H.Repko, J.Cowen, R.Stein, P.Duxbury who, as members of my guidance committee, must tolerate my writing for another 1A0 pages. I would also like to express my gratefulness to Professors S.D.Hahanti, T.A.Kaplan for a lot of useful discussions. I really appreciate the discussion and coorperation with Dr. R.Day in the area of the elastic network under tension. I also appreciate the communication with Professor D.Sherrington, Dr. T.C.Choy in the research about quasicrystals. Dr. E.Garboczi, Dr. H.Thomson, Dr.H.He and Mr. H.Xia, J.Gales, Y.Li, H.Yan have been good friends and the sources of much help. Especially the JPLOT developed by J.Gales has been a great tool in my research and the preparation of this thesis. I am grateful to the special help from J.S.Kovacs and J.King. ii Finally, I should like to thank Liansheng Cao whose un- derstanding, encouragement and assistance have helped bring this project to completion. iii List of Tables TABLE OF CONTENTS 00......OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOVi List of Figures ooooooooooooooooooooooooooooooooooooooooovii Preamble oooooooooooooooooooooooooooooooooooooooooooooooooo1 Chapter 1. Mapping between Random Central Force Networks and Random Resistor Networks 00000-00008 Section Section Section Section Section 1. 1 1 .2 Introduction ............................10 General Formulation ......................7 Illustration of the Mapping and the Centroid Algorithm 0-000-000-000022 Random Networks .........................32 Conclusions ooooooooooooooooooooooooooooofl3 Chapter 2. Low Energy Behavior of Infinitely Large 20 Penrose Tiling oooooooooooooooooooooooooooecu“ section 2.1 Introduction oooooooooooooooooooooooooooouu' Section 2.2 Numerical Simulation on Penrose Tiling 0-46 Section 2.3 Centroid Configuration Analysis -'°°-'°‘°50 Section 2.“ Equation of Motion Approach -°---°--°-¢0-78 iv Section 2.5 Random Halk Approach 000-00-0000-0-00000-88 Section 2.6 Conclusion oooooooooooooooooooooooooooooog1 Chapter 3. Percolation on Stretched Elastic Networks 0-0-09“ Section Section Section Section Section Section Section 3.1 Introduction oooooooooooooooooooooooooooogu 3 .2 Elasticity of Pure Triangular Net under Stress ........................99 Diluted Triangular Net under Stress ---Simulation Technique --°---¢°-o---°°-110 Diluted Triangular Net under Stress ---Simulation Results coo-o-oovoooooco-0113 Symmetry Analysis ......................121 Effective Medium Theory 0-000-0000-00000132 Conclusion and Future Research .........1uu Appendix A Isolated defects in the static limit 000-00-001u5 Appendix B Absorbing rate in continuum media 0-0000-0000-178 Appendix C Formulation of effective media on triangle net 0....00......OOOOOOOOOOOOOOOOOOOOOOOO000......18“ List of References oooooooooooooooooooooooooooooooooooooo186 Table Table Table Table LIST OF TABLES Solution of the Centroid Condition on the Unit Cells ooooooooooooooooo00000000000058 Degree of Incoincidence of the Cells Boundary 00069 Comparison of the Bond Lengths Belong to Different Unit Cells ooooooooooooooooooooooo00071 Summary of the Diffusion Constant on 2D Penrose Tiling oooooooooooooo00000000000092 vi Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure LIST OF FIGURES 1.1 A square net with 702 of bonds present (a) before and (b) after relaxation with the centroid condition.000-00-0000-0-000025 1.2 Square nets with a fraction p of bonds present after relaxation.00-00-00-0-00000-000027 1.3 Conductance of a square net with fraction p of the sites present compared with the effective medium result of Hatson and Leath.oooooooooeooooooooeooooo00000028 Conductance of square net for site and bond percolations.oooooeooooooooooooooooooooooZQ Conductance of triangular net for site & bond percolation.ooooooooooooooooooooooooooozg (a) A three coordinated random network (b) Dual net of (a).oooooooooooooooo000000000033 Relaxed configurations of Figure 1.6.00000-0034 Illustration of the dual problem on random network.eoeooooooooooeooooooooooo000.0037 Conductance o of the network shown in Figure 1.6a, wAere p is the fraction of normal bonds and 1-p insulating bonds. Also shown is resistance R of the dual net shown in Figure 1.6b, ahere p is the fraction of normal bonds and I-p superconducting bonds. The solid lin is the effective medium theory result.-----°------N1 1.10 Same as Fig.1.9 except 0 and R are substituted for 01 and R2.og-oooooo-ooooo-oo-oou1 Deflation and matching rules in the generation of 2-d Penrose tiling with seed (a) Fat rhombus (b) Thin rhombus.---------A7 2.1 2.2 203 A portion of Penrose tiling lattice.----°---°H7 Penrose tiling with fixed circular boundary (a) before and (b) after the relaxation by centroid condition.00-00000-0000048- 2.” Most stationary vertices during vii Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure relaxation connected appropriately to form a super Penrose tiling.----0----------51 2.5 Basic cells (a) unit-1 and (b) unit-2.000000051 2.6 of unit-1 00000000000000000067 Enviromental configuration and unit-2 in Penrose tiling. 2.7 Basic cells (a) unit-3 and (b) unit-H.000000068 2.8 Illustration of the boundary fitting requirements between unit-1 & unit-2 and unit-3 & unit-u.000000000000000000000000000000070 2.9a Illustration of the assignment of weights to bonds in different region of unit.00000000000000000000000000000000000000000076 2.9b Illustration of the similar topology between adjacent generations of units.°°°0--°°°77 2.10 Density of states of Penrose tiling by the Equation of Motion.00000000000000000000000082 2.11 Low energy 8(QtE) V3. E. neutron scattering spectrum O..000.0.0000...0.0.00.0000000000085 2.12 Low energy neutron scattering spectrum S(q’E) vs. Q.000000000000000000000000000000000086 vs. (qa)2 in spectrum S(q,E), 2.13 E Bitgact the slope D.00000000000000000000000087 to 3.1 C C C .purAItrigfigulgg 3.2 Transformations for uniform external strain.000000000000000000000000000000000000000107 and C ' vs. n=L /L for lattigg.0000000090000000000000105 3.3 B. u , u and b vs. n=L0/L for pure triangulapslattice.000000000000000000000000000109 3.“ Relation between the percolation problem ontennis racket model and conductivity percolation and rigidity percolation.000000000000000000000000000000000011n 3.5 C1 and CAN for diluted triangular network.00000000000000000000000000000000000000115 3.6 B, u ur and b for diluted triangular 9 networE000000000000000000000000000000000000000116 3.7 Illustration of reduced "diode effect".-----118 viii Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure .10 Simulation result of CAN - C .11 Simulation result of C .13 B (square) and u .8 Pressure for relaxed diluted triangular net.000000000000000000000000000000000000000000120 .9 Static energy for relaxed diluted triangular net.0000000000000000000000000000000120 P for 0:000000000000126 diluted networks.°°--¢--ooooo +C'=C -C for diluted networks.000000053000001100001g000129 .12 Simulation result for Cauchy's relation C =C' 12 H4 .00000000000000000000000000000000000130 (diamond) vs. p, averaged over 5 sagples, where u is scaled by the factor of 1/(1-n/2I.0--°----0000133 .14 Phase boundary determined by vanishing of elastic constants.000000000000000000000000000013H a .15.Computer simulation gesults for p E’ 000000000000000000000135 p v p o P 1 P ,P . P B ur us b .16 Sketch to explain the single defect configuration before and after relaxation.-¢°-137 a .17 Effective medium theory result of a . The symbol is the initial slope intercepts on p-axis for static energy in the simulation.000000000000000000000000000000000001M0 .18 Effective medium theory results for p E, p P and p B.00000000000000000000000000001’43 ix PREAMBLE Traditionally, solid state physics has meant crystal physics. Yet, many of the materials we deal with do not have the simplicity and regularity found in crystals. In modern technology, noncrystalline solids play a very impor- tant role, for instance, the amorphous semiconductor in solar cells, the ultratransparent optical fibres in telecom- munications, and the ubiquitous everyday uses for various kinds of glasses. In recent years, research in noncrystal— line solids has been one of the most active fields in condensed matter physics. It is well known that the mathematical amenities of the solid state theory (Brillouin zones, Bloch states, group- theoretical selection rules, etc.) are associated with the long range translational periodicity in the crystalline solid state. It is an intellectual challenge to achieve physical insights in noncrystalline materials which have no translational symmetry. While some of the traditional ap- proaches remain useful in noncrystalline materials, this challenge has been met mainly by new approaches. (Zallen 1983; Elliott et al 197“) Percolation theory is one of the new theoretical tech- niques applied to strongly disordered system. (Zallen 1983; Stauffer 1985; Straley 1983) It deals with the effects of varying a parameter such as the concentration or density in a random system. When increasing (or decreasing) the con- centration, a sharp phase transition is found at whikfli long range connectivity of the system suddenly appears (or disappears). The basic concepts in percolation theory can be understood from the following example. Consider the fluid flow in a porous medium which is modeled as a network of interconnected channels. Some of the channels are blocked at random. Those unblocked channels constitute a path through which the fluid can flow. A typical question in percolation is: what fraction of the channels must be blocked in order to prevent the fluid flowing through the whole medium? 0r equivalently, one can ask, what is the threshold concentration of the unblocked channels that guarantees a finite probability of the existence of a perco- lated cluster of unblocked channels which spans the whole system. Since the elements in consideration are the chan- nels (bonds), the problem is called "bond percolation". Its analogy is "site percolation" where instead of channels being blocked, the vertices (sites) of the network are pluged. Whenever a site is pluged, no fluid can flow through the channels adjacent to this site. Both bond and site percolation problems described here are determined by whether or not there exist a geometrically connected path through the system. Therefore, this kind of percolation problem is named connectivity percolation. The concentra- tion pc at which the dramatic change in tn“; long range connectivity occurs is called the percolation threshold. This kind of percolation transition can be applied to several varieties of transport phenomena occuring in amor- phous solids, for instance, the insulator/metal transition in conductor-insulator composite material; the disconnected/connected transition in resistor networks; the normal/superconducting transtion in composite superconductor-metal materials. (Zallen 1983) One of the most prominent applications of percolation theory in condensed matter physics is the glass transition. It is widely believed that the structure of amorphous solids whose bonding is primarily covalent is well described by some form of covalent random network (Zallen 1983). Examples are Se1_ As , Se _xGex and Se x x 1 Each atom (with coordination n) represents one n-functional 1_x_yAsxGey, etc. unit. Two neighboring atoms are linked by a covalent bond. For example, Se (n=2) atom can form covalent bonds with two other atoms, while Ge has coordination number four. The desired fractions of x and 1-x are achieved through the cor- rect compositions and the melt-quench process. ILt is easy to relate this fraction x to the concentration probability in an appropriate percolation model. Then a natural quesm- tion is how do the elastic properties depend on the fraction x? (Thorpe 1983) The connection between the elastic properties of a net- work and percolation was first established in a special type of glass transition: the sol » gel transition. (de Gennes 1976; Stauffer 1976) It is convenient to use the elastic properties to describe this transition since it naturally captures the essential feature of the condensation. In a gelation process, the elastic shear modulus vanishes for the sol (a liquid), begins to grow at the sol . gel transition, and continues to increase in the gel (solid) phase. One can see that this transition in the macroscopic elasticity is intimately linked to the percolation process involving the cross-linking between monomers. As the polymerization reac- tion goes on and more cross-links form, larger and larger molecules appear. Eventually, and abruptly, at a critical stage in the condensation process, an "infinitely extended" molecule appears, a huge molecule whose extent is limited only by the size of the vessel in which the reaction is taking place. Only this gel macromolecule contributes to the elastic shear moduli. It was pointed out by de Gennes that if the interaction between nearest neighbors is modeled by an isotropic elastic force, then the elastic modulus of the gel is exactly equivalent to the conductivity of the resistor network. Therefore the sol + gel transition belongs to the same universality class as connectivity percolation. A new universality class of percolation was found in an elastic network when the interaction between nearest neigh- bors was modeled by pure central forces. (Feng & Sen 198A) In this model, the bond between interacting sites is re- placed by an ordinary Hooke's law spring. Both numerical simulation and analytical results showed that with purely central forces, the threshold pcen at which the elasticity of the network vanishes is significantly higher than the threshold pc of the corresponding isotropic force problem. (Feng & Sen 198A; Feng, Thorpe & Garboczi 1985) This tran- sition is called rigigity_percolation as opposed to connectivity percolation. In this new class of problem, the elasticity of the network vanishes even when the whole net- work remains well connected geometrically. This phenomena is explained below. In the pure central force network, when bonds are removed from a pure lattice, some local floppy regions are created. These geometrically connected struc«- tures are called floppy since they are not able to transmit elastic forces. Eventually these local floppy regions prevent the rigid regions from percolating and the system loses its elasticity. Actually, what we have considered as pure isotropic forces (connectivity percolation) and pure central forces (rigidity percolation) are two extreme cases of the Born model: V = % (0-8) X [(ui-GJ)-rij]2 giJ iJ where the first term or 3:0 case is the pure central force and the second or 8:0 case is the pure isotropic force. A crossover is found when both central and isotropic forces are included. (Feng & Sen 198A; Schwartz et al 1985) However, the Born model has some troubles when one at- tempts to relate the macroscopic properties of amorphous materials to the microscopic atomic structures. First, the Born model is not rotationally invariant when 8 is nonzero. (Keating 1966) Second, the internal stress is neglected, which is not justified. It is well known that all real materials have internal stresses. One expects the stress- induced scalar elastic energy to be large and even dominant in many cases. It is unclear how the percolation transition is modified when the internal stresses of the material are involved. This inspired us to propose and study a new model —-- the tennis racket model, which is defined as the central- force network under tension. In Chapter 1, we establish a mapping between the resistor network and the central force network under very large tension (which can be imagined as the system having springs with natural length zero). The validity of this mapping is tested on several well under- stood systems. One result of this mapping is that it reduces the electrical problem to a geometrical problem and hence provides a useful visualization of the conductivity process. This advantage is explored in Chapter 2 in the cal- culation of the conductivity of the 2 dimensional quasicrystal or Penrose tiling. (Shechtman et al 198“; Socolar and Steinhardt 1986) Another result of the mapping which may be more important, is that connectivity percola- tion and rigidity percolation, which have been widely accepted as belonging to two different universality percola- tion classes, are intimately tied together. In other words, as soon as the internal stress is considered, these two classes of percolation problems occur naturally in a pote- tiongl inyggient central force spring system. Chapter 3 gives a general discussion of the elastic behavior of such a system. We found that connectivity and rigidity percola- tions are actually two special cases of the central force network under extremely large tension or zero tension respectively. As the tension is varied, the percolation threshold evolves smoothly between these two limits of pc and p We believe that the extensive results presented cen' here for the elastic behavior of the noncrystalline network under tension will be useful in understanding many of real systems like glasses (Thorpe 1983, He & Thorpe 1985) and the newly discovered high Tc superconductors. (Khurana 1987) Chapter 1. Mapping between Random Central Force Networks and Random Resistor Networks Section 1.1 Introduction The equilibrium voltage configuration of a random resistor network with an applied external voltage source can be obtained through the minimization of the electric energy stored in the network with respect to the voltage at each node. This leads to the Kirchoff's equations, which can be solved numerically as accurately asldesired. This problem is equivalent to many other physical problems. One example is the spin waves in a Heisenberg ferromag- net at low temperatures. With magnetic atoms on nodes of network, the voltage 111 is analogous to a tilt angle for the magnetic moment i. and the bond conductance maps onto the ferromagnetic exchange integrals. It was shown that the spin stiffness of such systems is the electric conductivity of the corresponding electric network. (Kirkpatrick 1973) Another example is a gel obtained by polymerization of z-functional monomers. Modeled by an isotropic force con- stant, the elastic modulus of a gel is analogous to the electrical conductivity of an electrical network. The poten- tial V represents one component (say Xi) of the elastic i displacement of the i£fl monomer from its rest position on the lattice. (de Gennes 1976) In this chapter we explore yet another mapping. We will map the random electrical network to a central force network. This mapping is slightly more subtle and leads to some new physical insights. Our purpose in this chapter is primarily to develop and explore the mapping. The organization of this chapter is as follows. In section 1.2, we develop the mapping and formula. We also show the relationship to the low energy density of states of’a tight binding Hamiltonian. In section 1.3, we present the results for dilute resistor networks on square and triangular lattices as an illustration of the centroid algorithm based on the mapping. In section 1.“, we use the new algorithm to study the percolation properties of a 2- dimensional random network. A dual random network is constructed and it is shown how the superconducting/normal network on the dual lattice maps on to the dilute resistor network on the original lattice and vise versa. The relationship of this work to rigidity percolation is dis- cussed in section 1.5. 10 Section 1.2 General Formulation 1.2.1 The ngping To illustrate the mapping, let us consider first a resistor network in which a voltage difference V is applied 0 across the two opposite sides. The current flows in the x direction across the voltage difference V0. The conduc- tivity Gxx is defined by E = Gxxvg where E is the electrical energy stored in the resistor network. Let us label the sites by i,j etc. so that a current Ii flows in the ij bond J which has a conductance oij' Sites that are not connected V have 01.1 = 0. Current conservation at each site leads to {111:0 (1.1a) J which using Iij - oiJ(Vi- VJ), where V1 is the voltage at site i, leads to X 01J(Vi- VJ) = 0 (1.2a) J or E a V iJJ v1-1—Z——- (1.3.) o i J J The total energy stored in the network is Z °ij(Vi' VJ)2 = c v2 (1.ua) (i,j) xx 0 r!) u 11 where the angular brackets denote that each bond is only counted once in the sum. The conductivity is given by _ _ 2 2 cxx - °13(V1 VJ) /vO (1.5a) Now let us imagine a central force (Hooke) spring net- work in which each spring has natural length 5359. The whole system is held on a frame to prevent it collapsing to a point. For a large system the shape of the frame is ira- relevant, but it is convenient to imagine it to have a square shape. The equilibrium condition for each site i is that the total force acting on it vanishes. We have in the x direction X KIJ(R1x- RJX) = o (1.2b) J where R1 is the position of the ith atom and Ktj is the spring constant. From this we see that E Kij ij Rix : (1.3b) 2K1, J There are equations similar to (1.2b) and (1.3b) in the y direction. The total energy stored in the spring system is 1 - 2 [Kxx + yy] L (1.ub) where L is the length of the sample and 2 2 K = K R - R /L 1.5b) xx ij ( ix jx) ( and a similar expression for Kyy' We see that the (a) and (b) equations can be made to coincide through the mapping, oij . K1J vi + nix (1.6) Gxx ‘. Kxx There is an important subtlety in this mapping as- sociated with the boundary conditions. In the resistor problem, the net current flow is in a particular direction; x in the case considered here. This leads to the conduc- tivity G In the spring problem, the frame acts equally xx' in the x and y directions so that xxx and Kyy cannot be ob- tained separately. This is not a problem for high symmetry networks, where the tensors G08 and KaB are proportional to the unit tensor. We shall refer to such networks as electrically isotropic. We shall only discuss such networks 13 in this chapter. These include square nets, random square nets, triangular nets, random triangular nets and the random ' networks considered in section 1.“. We are thus led to _ _ 2 2 - (1%) x” (Ai AJ) /L (1.7) The conductance can be calculated from the spring network via eq.(1.7) and the equilibrium condition (1.2b) that may be written in vector form {K(§-A)=o (1.8) J11: j In ZD, the conductance and conductivity of a square sample are the same. For a sample of hypercubic shape in d dimensions (volume Ld), the current flows between parallel hyperplanes (area Ld'1) so that the conductance o is given by d-2 o = G/L (1.9) where 0 comes from eq.(1.7). Note that eq.(1.9) holds for electrically isotropic networks. Of course eq.(1.7) must be modified so that the trace is over all C! dimensions and L2 1“ in eqs. (1.“b), (1.5b) and (1.7) must be replaced by Ld so that Tr’c’ . K. (A - A )2/Ld (1.7a) IJ 1 J and Eq. (1.“b) can be extended to read E = % Tr 'K’ L2 . % Tr ’c‘ L2 = g c L2 = g o Ld (1.“c) The pressure P can be obtained from p = 2E3 . -E%§ = o . (1.10) 3L L Thus the conductivity is obtained from (1.7a) while the egg- ductance is equivalent to the the outward pressure needed to prevent the system collapsing to a point. For 2 dimensional systems, this pressure is independent of the size of the system L. For 1 dimensional systems, the pressure increases with size linearly, whereas for d > 2, the pressure decreases as the size gets larger. The total 32123 required to prevent the system collapsing, ilflili increases with size. 15 It may be convenient in some cases to have the spring network stretched on a circular frame or tennis racket, so .that the boundary condition is isotropic. 1.2.2 Tight Binding Hamiltonian A useful alternative viewpoint is given by the tight binding Hamiltonian H = l I ‘1} [11><1| — 11><31] (1.11) 1.1 th where the li> is a localized 3 state on the 1 site (Alben et a1 1972). The overlap integral K is the same as the 1J spring constant in Eq. (1.2b). It is convenient to take a (large) unit supercell containing n atoms. This cell is repeated periodically. Each atom is designated by a label- ing (2,n) where 2 designates the cell and n is the atom within the cell. This label pair is conveniently denoted by i .. (2,n). We can use Bloch's theorem on the supercell translations R2 to transform to a new basis 16-h. lq,n>:—: {e 2+“ l9.+n> (1.12) 2 1 (No The total number of atoms N = n Nc where Nc is the number of supercells. The Hamiltonian (1.11) is Block diagonal in this basis and within a 5 block the matrix elements are given by 16 ' ' ' N , 11' c 22 (1.13) where atoms 1 and i' are connected by the hopping matrix element K11,. These two atoms may either be in the same unit cell or in adjacent unit cells. The Hamiltonian (1.11) deliberately has related diagonal and off diagonal terms so that the band edge is at E = 0. The eigenvector corresponding to zero energy is denoted by |0> IO) = -l— i 16 = O,n> (1.1“) n Labeling the other n-1 eigenvalues at q = 0 by E we use e! the matrix analogue of hop perturbation theory (Harrison 1980). Given <5.n|H|5.n'> = %; 1Z1,K11'I'15‘(§1' fii') .31: (5411- 11.112} (1.15) then to second order 1 . 2 £9 = 2N 1X11 K11'[Q.(§1- fii)] 6411-11,.) 2 - I Z x , _ (5:0,nlE > /E (1.16) e l1,1' 1‘ Nc/n e I e 17 There is no term linear in q in (1.16) as this vanishes when the sum over all sites is done. If we are interested only in the density of states p(E),then distances are irrelevant. We imagine forming an auxilau lattice in which the second term in (1.16) vanishes. This can be achieved if {x (A-fi)=o (1.17) 31.1 1 J this is the same condition as (1.8) and defines the same set of new atomic positions R1. Thus the sites in the auxilary lattice are coincident with the positions of the sites in the central force network discussed in subsection 1.2.1 In the auxilary lattice sq = g3 .{ xiJ [q-(Ri- 1J112 (1.18) 1»J when i,j go over all sites. The unit supercell is not ex- plicitly appearing in the summations in Eqs. (1.17) and (1.18). Note that (1.17) can be obtained by minimizing the energy in Eq. (1.18) with respect to the R1 if the R1 are initially regarded as arbitrary. From (1.17) we may define principal axes where a :1) q2+D q (1.19) 18 The area of the ellipse containing states up to an energy E is IE/lenyy so that the low energy density of states is L I p(E) = [5;] —:::::: (1.20) nannyy where L2 is the area of the sample. This can be written as L 2 9(8) = [5;] ' (1.21) lbet'b’ where .D’ is the diffusion matrix which from (1.18) and (1.19) has matrix elements 1 Da8 ' 2N izj Kijmia - Rja)(RiB - RjB) (1'22) ’ Note that 0-9 L2 0-9 Tr D =N—TI‘G (1.23) where Tr‘G‘is defined by Eq. (1.7). For electrically isotropic 2d systems we also have Jnet‘n‘ = % Tr D’ L2 0 -0 : a T!‘ G (1.2“) where o is the conductance. Hence p(E) = 5%; (1.25) This density of states is normalized over all energies to the number of sites N, and the low energy part depends only on the conductance 0. Note that this result is very general and applies to any network that behaves as a continuum in the long wavelength limit. This excludes fractals. The geometry of the network enters only through 0. These equations are easily generalized to higher dimensions. For Eqs. (1.11) to (1.18) there are no changes. We have to modify Eq. (1.19) to 2 sq - X Daaqa (1.19a) a where a goes over the d principal«directions. The density' of states given in Eqs. (1.20) and (1.21) becomes 2 2 p(E) = L (L E] (1.21a) uuT(d/2)/Det’p’ where P(d/2) is the Gamma function. Equations (1.22)-(1.23) are unchanged but (1.2“)-(1.25) become 20 «pet'p‘ = % Tr‘p’ L2 0- .0 =‘d—N TPG (1.21‘3) - 1:1 ' N and hence d/2-1 N 1 E p(E) = n10 r(d/2) [3;] . (1.25a) Thus as in 2D, the low energy density of states is deter- mined by the conductance o and the total number of atoms N in the network and p(E) is an extensive quantity. An atom is counted as being in the network if it is coupled to the backbone or conducting path. This includes regions that are connected to the backbone but carry no current. They nevertheless contribute to the inertia and so must be counted. Isolated regions that are not coupled to the back- bone are ignored. This method of using the tight binding Hamiltonian (1.11) to derive the centroid condition (1.17) is less satisfactory than mapping onto central force springs of natural length zero. This is because we had to introduce a "supercell" in order to have a a vector to work with. The final answer did not explicitly contain reference to this 21 supercell and so we regard the result (1.25) as being quite general and true for any conducting network that is homogeneous on distances greater than some correlation length g. The derivation that lead to (1.25) is only valid for q<<§'1 and hence for vanishing by small energies as 5 ~ 0. We believe the result is actually much more general as the relevant length g is probably not the supercell size but rather the distance that characterizes the structural correlations in the random network. A similar quantity g is commonly used in percolation theory. (Kirkpatrick 1973) There is no easy way around this problem in a dynamic ap- proach as given here. In the static approach in 1.2.1 of this section, every- thing is rigorous. In order to get the dynamic result (1.25) one has to assume an Einstein relation between the conductance and the diffusivity. (Kirkpatrick 1973) It is here that the problem is glossed over. This area needs fur- ther study. 22 Section 1.3 Illustrgtion of the Mapping and the Centroid Algorithm 1.3.1 Centroid Algorithm Based on the mapping described in last section, we developed a new algorithm to calculate the electric conduc- tivity of’the random resistor networks, which is called as Centroid Algorithm. This name comes from the fact that the condition E KiJ (A1 - AJ) = o for each 1 (1.8) is just the condition that every site is at the centroid of its nearest neighbors that are present if K takes only two 1i values K (with probability p) and 0 (with probability 1—p). We summarize the centroid algorithm as the following: (i) From the topology of a given resistor network, es- tablish the central force elastic network by replacing a Hooke's spring with natural length gegg at the place of each resistor with the spring constant K1J~ oij’ (ii) Choose an appropriate boundary condition to prevent the elastic network from collapsing. Two boundary conditions usually used are the fixed frame and the periodic boundary condition. 23 (iii) Relax the elastic network until the energy stored in the springs is minimized. Specifically each site in the network is iteratively moved according to: . (n) ‘11 fi) R (n+1) = 1" (1.88) 1 )_ 1cU J For the isotropic case: Kij - K with probability p Kij = 0 with probability 1-p (1.8a) is simplified to: a (n) (n+1)_ j J A1 - 21 (1.8b) where 21 is the co-ordination at site 1. (iv) Finally, the conductivity G is proportional to the mean square nearest neighbor distance in the relaxed network: 6 ~ Z 513 / L2 (1.26) where 61J = I R1 - RJI is a nearest neighbor distance. 2“ 1.3.2 Simulation Results on Square and Triangular Net Before we apply this method to any new structure, we use it on two well studied 2-dimensional systems: square net and triangular net. Let us consider the network shown in Figure 1.1a, where all the nearest neighbor bonds have resistance 061 and a fraction 1-p have been randomly removed. In the percolation problems one really wants to study the system in the thermodynamic limit (Noon). Periodic boundary conditions are better than any other kind as every atom is properly coordinated. He therefore periodically repeat the "supercell" in Figure 1.1a. Bonds are removed randomly in the reference supercell and also in all.<3thers. Rather than put the network on a frame when we go to the central spring model, we hold the supercell repeat vectors cgnstgnt. This is equivalent to an external pressure and more convenient in practice. The network is relaxed itera- tively, site by site, until the energy stored in the springs is minimized. Figure 1.1b shows the relaxed network of“ Figure 1.1a. It is observed that isolated islands do not contribute to the conductivity and so relax to points and the side groups on the backbone that do not carry current also relax to points. Thus the conducting backbone is clearly and simply exposed by this algorithm. Those bonds that carry most current are stretched the most and so make the largest contribution to G. 25 (b) (a) Figure 1.1 A square net with 701 of ham present (a) before and (b) after relaxatim with the centroid midition. 26 In Figure 1.2 we show some typical results for various values of p. Note that many of the long straight connec- tions in Figure 1.2 are actually made up of many bonds and collapsed side groups. Below pc, the whole network is broken and collapsed to many independent points, whikfli cor- responds to the vanishing of the conductivity. The networks for p = 0.70 in Figure 1.2 and Figure 1.1a are examples of different random configurations. In Figure 1.3 we compared our simulation results for the conductance for site percolation on the square net with the result of Waston and Leath (197“). It can be seen that our results lie close to the effective medium theory which agrees well with a direct solution of Kirchoff's laws except very close to pc. The tail that extends slightly below p0 = 0.59 in our simulations is due to finite size effects. As our method is quite general, it works equally well for bond or site depleted networks. The simulation results on a square net and on a triangular net, for both site and bond percolation problems, are shown in Figures 1.“ and 1.5 respectively. The thresholds pc agree well with previous work. Appendix A collects together the known results for the effect of a single defect into a single compact form that can be used in any lattice in the static limit. The result is used to calculate the effect of single defect, either bond or site, on the conductivity of various lat- tices, as well as the corresponding pc. 27 .8393? . 351.. ads 9685 8:8 no a 838.: a 53 Be: 95.3 m; a: 851a means coeua 28 .58: Ba 83m: .8 you? 58.. 938st «5 53 8.8.8 uzommco 83m 93 co o 838.: 55 um: cameo o no 85396.8 w; a“. n. . — a O 1 m6 0 O U 1 v6 0. n O 1... m I m o O a 1 m.o . — . Cofi Conductance Conductance 1.0 0.8 0.6 0.2 1.0 0.0 0.6 0.2 29 A. E] r . r ' .3" Site Percolation EJA Bond Percolation EIEIA " B A Cla B .1 B A A A -1 A B .. G B G A I 0.4 0.5 0.3 1.0 1 I I r I I £ Site Percolation (3'35 Bond Percolation B A - E1 E A EFI A G A _ A A A d3 A ‘ B B B BB .— 63 A9 A _ G A =EIEEI-'--I':-l-l——A#A ' ' 1 0 . 4 0 . 6 0 . B 1 . 0 p Figige 1.5 Gonmctanoe of triangular net for site A bond percolation. 30 Before leaving this section, we consider a general two dimensional network with nearest neighbor conductors of mag- nitude 00. The relaxed network and tight binding model [using Eqs. (1.17) and (1.18)) with the nearest neighbor K iJ . 00, leads to q2 2 Eq - “N 00 .Z 61J . (1.27) 1.J From this the density of states is p(E) = NLz/(noo X 5.2] (1.28) 1J i.J and hence from (1.25), the conductance o is given by )1 513 0 = 00 l‘J-§—— (1.29) “L This is a convenient form. For the square net with N atoms and nearest neighbor distance a, we have L2 = Na2 so that 0 = 00 (1.30) Similarly for the triangular net 0 = J3 00 (1.31) 31 and for the honeycomb lattice 0 = 0 / J3 . (1.32) 32 Section 1.“ Rendom Networks 1.“.1 Conductivity of 3-fold Co-ordinated Random Network and its Dual Network We consider a random network where every site is three fold coordinated as shown in Figure 1.6a. This 800-site network was constructed by successively disordering a honeycomb lattice while maintaining the coordination. (He 1985) The lattice was then relaxed to make lengths as equal as possible. There is a small (~51) variation in the bond lengths. We give all the bond conductor equal magnitude 00. We could have made the resistance proportional to the bond length but this seems unnecessarily complicated. In order to calculate the conductivity, every atom is moved to the centroid of its three neighbors. The relaxed network is shown in Figure 1.7a. The changes are surpr- isingly large. The large polygons have grown and the small polygons have shrunk. There is also much more variability in the bond lengths than in Figure 1.6a. Because the supercell of the sample shown in Figures 1.6a and 1.7a is not a square some additional care must be taken. The supercell shown has repeat distances Lx and Ly’ where Lx/Ly = 2/(3 = 1.16. These are held fixed during the relaxation. Minor changes in the previous formalism can easily be made for such a network. It can be shown that 33 (a) / \\ ‘11K'§/\\\ \ 6. / \ \ / \1 , /// LA / .1 .41 (b) Fim 1.6 (a) A three coordinated randon network (b) Dual net of (a). 311 (b) Fm 1.7 Relaxed configlrations of FIQJPB 1.6. 35 Eq. (1.29) holds for any electrically isotropic network, in- dependent of the shape, if L2 is replaced by the area A. U o = o 14-J— (1.29a) We use the relaxed network of Figure 1.7a, to compute a = 0.555 00. This should be compared to o = 0.577 00 for the honeycomb lattice Eq. (1.32). We also constructed the dual lattice of the random net- work as shown in Figure 1.6b. The dual transformation is simply assigning a dual vertex at the center of each polygon, then connecting those dual vertices that belong to adjacent polygons. An elementary application of Euler's theorem shows that the mean size of a polygon in the network of Figure 1.6a is six. Therefore in the dual lattice, the mean co-ordination is six, although it is not constant but. varies from site to site; the minimum co-ordination is 5 and the maximum is 9 in Figure 1.6b. The dual lattice shown in Figure 1.6b is relaxed using the centroid algorithm to obtain the relaxed lattice of Figure 1.5b. Using Eq. (1.29a) we find that o = 1.81 00. This should be compared to o = 1.73 00 for the triangular net in Eq. (1.31). 1.“.2 Duality Relationship for 2D Resistor Networks 36 It was known that there is a relationship between the conductivity of a 2D lattice with conductances chosen at random and the conductivity of a related distribution of conductances on the dual lattice. (Keller 196“) Figure 1.8 can be used to explain the dual problem, where solid lines denotes a portion of direct random network and dashed lines shows the corresponding dual lattice. Let the potential at the vertex i on direct lattice be V1 and the conductance of the bond between i and its neighbor j be 011. 11 : V1 - Vj’ which correspond to the continuum conditions va V = 0 We can write down two conditions involving V and 7-3 = 0: Z V11 = 0 summed along any single loop (1.33a) Z 011 vij = 0 summed over all bonds involving site i (1.33b) Now we wish to introduce a function Wi similar to Vi , defined on the dual lattice. We label the dual lattice sites by i', j' etc. under the convention that bond (i',j') crosses bond (i,j). Then we define: w = o v. (1.3“) 37 FM 1.8 Illustration of the dial problen on randan network. 38 That means the W difference between regions seperated by a. lattice bond is the current flowing along that bond. It is easy to see that W around any closed loop on the dual lat- tice sums to zero, because it is equivalent to the condition that no net current flows into any vertex of the direct lattice. Actually condition (1.33b) gives: 2 Hi'j' = O (1.353) To get an equation on dual lattice equivalent to (1.33b), we define the conductance of each dual bond n1, as the 1' reciprocal of the conductance of the direct bond it crosses, i.e. o. = 1 (1.36) 111.1. 11 Substituting (1.3“), (1.36) into (1.33a) we have: Z VIJ : l 011J1 "11J1 = 0 (1°35b) Now it is clear that W is the solution to Kirchoff's equa- tion for this dual problem. It was proved (Hendelson 1975) that the conductivity of direct lattice odirect( q,a; p,b) of a concentration q of conductances a and concentration p of conductances b are re- lated to that for dual lattice by: 39 odirect(q’a; pvb) : I odual(p’b ; Qta ) l (1.37) where q=1-p. This relation holds for any planar graph. 1.“.3 Simulation of Insuletor/normal and Superconductor/normal Percolation If we substitute 0 (insulator) N 11 b = 1 (normal) into (1.37), we get: odirect(q’o; p.11 = 1 odu311p.1; q.~) 1“ (1.38) Note that the distribution of insulator in direct lattice is related to that of the superconductor in the dual lattice. This relation tells us that the gonguctiyity in insulator/normal percolation is equal to the resistivity in superconductor/normal percolation on its dual lattice. We proceed four sets of computer simulations: (1) Insulator/normal percolation on random network. (Figure 1.6a) (2) Superconductor/normal percolation on dual network. (Figure 1.6b) (3) Superconductor/normal percolation on random net. (Figure 1.6a) “O (“) Insulator/normal percolation on dual network. (Figure 1.6b) We expect (1) and (2), (3) and (“) would satisfy the duality relation (38). Quantities referring to the direct lattice (Figure 1.6a) are denoted with subscripts 1 and quantities referring to the dual network (Figure 1.6b) are denoted by subscript 2. In the insulator/normal percolation problem (1) and (“), we randomly remove bonds from these networks with prob- ability 1-p. The depleted networks are relaxed using centroid algorithm and the conductances are computed by the method we have discussed in this chapter. The results for 01, against p1, and 02 against p2 are shown in Figures 1.9 and 1.10. These results have been obtained be averaging and 0 over 15 random depletions and o are set equal to one 1 2 for the undepleted systems for convenience. As expected for such dual lattices, the critical points are related by p10 = 1-p2c. Note that for the honeycomb lattice (also coordination 3) the bond percolation concentration pc : 0.65. From our simulations this pc 1 p1‘: as might be expected. In Figures 1.9 and 1.10, the solid lines are the effective medium theory results (for details see Appensix A) 0 1| (392-1)/2 (1.“0) N I 41 01, R; 0 . 1 L l 0 0.2 0.4 0.5 p1, p2 F‘gure 1.9 Conductance c1 of the network Sun in Figure 1.6a, mere is the fraction of normal bonds and 1-p insulating bonds. Also shown is resistance ofthemalnetstminI-ligme1.6b,mere isthefi'acticn of nomal and 1-p2 superconducting bonds. The so id line is the effectiveuadimtheoryresut. 1.0 r 1 T 0.8 - 0- 0.6 '- D: .7 b 0.4 - 0.2 - o l L . ¢ 0 0.2 Pup, Elm 1.10 Sane as Fig.1.9 except 02 and R1 are substituted for 01 and R2. “2 which are there as useful guides to the eye. In the superconductor/normal percolation problems (2) and (3), we randomly change a fraction 1-p bonds from normal to superconducting. The resistance of the normal systems is normalized to 1. The resistivity for these networks, averaged over 15 samples, are also shown in Figures 1.9 and 1.10 and are seen to have the behavior expected of dual networks. In order to facilitate numerical simulations we have used a ratio of conductivities of 1000:1 between super- conducting and normal bonds. The error introduced by this is mainly near pc when there is a small incremental enhance- ment in the finite size tail. It would not be hard in: take a larger conductivity ratio. Note that in the relaxed net- work, the superconducting bonds are fused together to a point so that consequently the normal bonds are elongated as the sample must remain fully connected and is not allowed to collapse. It can be seen that the dilute resistor problem on a network maps onto the superconducting/normal network on the dual network as expected. “3 Section 1.5 Conclusions We have introduced a new mapping for the random resis- tor problem. Its utility is that it reduces the electrical problem to a geometric problem and hence provides a useful visualization of the conduction process. We have applied these results to a few well studied problems. We have also studied depleted random networks and their dual networks for the first time. The results for the depleted three-coordinated random network are very close to those for the depleted honeycomb lattice and the results for the dual networks with mean coordination 6 are very similar to the depleted triangular lattice. This is to be expected since the co-ordination is the single most important parameter characterizing a percolating network. Also, the duality relationship is checked through the simulation of insulator/normal and superconductor/normal percolation. The depletion of central force networks has recently been studied by many authors. These networks are such that every spring has its natural length in the absence of exter- nal stresses. This is often referred to as rigigity percolation. The elasticity vanishes at p > pc because cen connected paths are inelastically ineffective. This is not the case in the work here as all springs in the conducting backbone are stretched by the frame. Hence the conductivity and elasticity are intimately tied together as we have shown. Chapter 2. Low Energy Behavior of Infinitely Large 2D PENROSE Tiling Section 2.1 Introduction Unlike crystalline and glassy structures, which have one and an infinite number of unit cells respectively, the quasicrystal structure can be generated from two or more (but finite) number of unit cells. One of the advantages of the "Centroid Algorithm" developed in Chapter 1 is the geometrical visualization of conductivity problem on random networks. It is noticed that on crystalline lattice the realization of the centroid con- dition is trivial since there is only one unit cell. 0n the other hand, for glessy structure, since there are an in- finite number of units, we can only realize the centroid condition by numerical simulation as shown in Chapter 1. As a phase between these two, a quasicrystal has more than two (but finite) unit cells. The question then arose as to the possibility of solving the centroid condition on Penrose Tiling exactly, i.e. not to depend on the numerical simulation. 44 “5 The so called Penrose Tiling is actually a two dimen- sional quasicrystal with "thin" and "fat" rhombi [ with angles (1/5, “1/5) and (21/5, 31/5) respectively] as two unit cells. The tiling invented by British mathematical physicist Roger Penrose is important here because (1) it gives us a non-trivial structure with purely Bragg diffrac- tion and non-crystalline symmetry; (2) decorations of Penrose tilings are good structural models for the ex- perimental realization. In this Chapter, we report our numerical simulation on Penrose tiling at first in section 2.2. Then in sectixni 2.3 we explore the self-similar properties and establish a sys- tematic method to analyze the centroid configuration from which the conductivity is extracted. Section 2.“ presents the spin stiffness given by the density of states (DOS) and intensity spectrum S(q,E) by using the Equation of Motion technique. Section 2.5 presents the diffusion constant ex- tracted from Random Walk simulation. Section 2.6 gives a comparison of these three techniques and conclusions. “6 Section 2.2 Numerical Simulation on Penrose Tiling There are several methods to generate the Penrose Tiling: (1) the matching and deflation rules, (2) grid projection, (3) direct projections, and (“) generalized dual method. (Levine and Steinhardt 1986) We applied the first method in our computer modeling. The main idea of our algo- rithm is illustrated in Figure 2.1 . We start from one of the two basic unit cells: rhombi with acute angles 36' and 72°, called "fat" and "thin" rhombus respectively. During the process of deflation, each fat rhombus is deflated into three fat and two thin rhombi, and each thin rhombus is deflated into two fat and two thin.‘The arrows are used to explain the matching rules. Solid arrows are associated with seed rhombi which give a unique direction of deflation. Correspondingly, the dashed arrows are added to current rhombi to denote the direction of next generation of deflation. The matching rule embedded in these arrows guarantees that there will be no holes of mismatch appearing during deflation. Figure 2.2 shows a fragment of Penrose Tiling generated by this algorithm. The centroid condition illustrated in Chapter 1 is used to relax the Penrose Tiling through numerical simulation. We took the circular piece as shown in Figure 2.3a, which con- tains 788 interior sites and 100 sites on the boundary. Recall that the bonds between nearest neighbors are imagi- nary central force springs with natural length zero, 47 \ \ \:— / \ / / 2:_ \ \ ) / / (b) (a) / / I Figure 2.1 Deflation and matching rules in the gaieraticn of 2-d Pamose tilirg with seed (a) Fat m (b) Thin Hum. a e .7.‘.§’. ewe-y ”$53.35;. 3:“? 1953:; a": I “u".‘P-fl'a‘ week)" ’0; :‘910 ‘\..7.‘.O.* . O... ’ affinity-awn?" O.‘ I ./.$ I‘z ’ I11" I z 3. “ "A‘aa'lfi'e‘eyingelggliaaaaflge 0 \e. 01.xe c“ 0 Q 1 1" ‘3 e :6 l PM 2.2 A portion of Peruse tiling lattice. 118 . ‘ .1.‘39g:”4'¢;e '01}.'\1‘I’:"‘II’:"‘III‘.'\\I , firmwares To. I'e‘caee“‘u€0’f" (assesses-«:1- e l - e :1 a. elflt'.‘1‘..".f’.§‘.'1:el ‘fl’v‘i’l‘l‘afil‘z’a‘fifi’tézvfifiwfifi ’ “\O’fiéfi‘a‘fi’euefi’fle:afi’fi cultiglelafiglolelello I. as. ‘94. \‘ .,.e;§..‘040.e\ white-:2 straws-av ‘tn‘smsrsigrv I «excess Elm 2.3 Pevose tililg with fixed circular boundary (a) before and (b) after the relaxatim by centroid cmditicn. “9 the whole sample must be put on a frame to prevent it from collapsing to a point. This is actually realized by fixing those boundary sites during the relaxation. Each interior site is moved to the center of its nearest neighbors iteratively. Figure 2.3b is the relaxed Penrose tiling. By using the formula given in Chapter 1, we obtain the conductivity on the Penrose Tiling, o = 0,141! 11 = 0.936 0, (2.1) where a is the unit length in the unrelaxed Penrose Tiling and z is the average coordination number which is four for Penrose tiling. It is known that o is equivalent to the diffusion constant D through Kirkpatrick relation G=PD with the probability of a site belong to the backbone P=1 (Kirkpatrick 1973), noting that the conductance is equal to the conductivity 0 in 2-dimension. An interesting phenomenon is observed that there are some points kept stationary during the relaxation (in range of our precision). The significance of this observation is that if these vertices are really stationary, then we must be able to find out some basic cells and from which to real- ize the controid condition analytically. This possibility is explored in next section. 50 Section 2.3 Centroid Configuration Analysis 2.3.1 Fixed Vertices Hypothesis It is not too surprising to notice that those points which seem stationary numerically during relaxation happen to be those vertices with the local five fold symmetry. In other words, in the unrelaxed Penrose. structure, these points are automatically located at the centroid of their nearest neighbors. What inspired us is that if we connect these points in an appropriate way, we get a pattern with the Penrose Tiling topology as shown in Figure 2.“, where on the relaxed Penrose Tiling, we superimposed this pattern by heavy lines. We know that in the original unrelaxed Penrose tiling, this pattern is exactly the same as the original one but scaled by a factor 13 ,where T = (1+/5)/2 is the golden ratio. However, on the relaxed Penrose Tiling where each site has been moved to the centroid of its nearest neigh— bors, we only can say this superlattice has the Penrose Tiling topology. This is a necessary but not sufficient con- dition for the existence of basic unit cells. The sufficient condition is one of the following three arguments which are equivalent each other: (1) The vertices with local five fold symmetry are sta- tionary points in the centroid relaxation. (2) The superlattice on the unrelaxed Penrose Tiling is identical with that of the relaxed one. 51 Figure 2.“ Most stationary vertices during relaxation connected appropriately to form a alper Pavose tiling. (b) Fm 2.5 Basic cells (a) unit-1 and (b) unit-2. 52 (3) There are only two basic unit cells in the centroid configuration as shown in Figure 2.“. Figure 2.5 presents the topology of these two cells. These cells must match each other preserving the centroid condition on the boundary between them. These equivalent arguments are called as "Fixed Vertices" hypothesis. Our purpose was to find out the exact centroid con- figuration of the infinitely large Penrose Tiling and from which to extract the average bond length (which is related to the conductivity or diffusion constant through the map- ping introduced in Chapter 1). Aiming at this, we took the third argument as our working hypothesis, i.e., we started by studying the two unit cells suggested by numerical simulation. The validity of the hypothesis will be self- evident from the results. 2.3.2 Method of Centroid Configuration Analysis We consider the unit cell of the superlattice with isosceles triangle shapes and internal structures as shown in Figure 2.5. Notice that unit-1 in Figure 2.5a has the same topology as part of unit-2 in Figure 2.5b. We assume the shape of these two units are [72“, 72°, 36°] and [36“, 108', 36°] respectively and the ratio of the sides is 1:1. We take unit-2 as an example to illustrate the procedure of centroid configuration analysis. 53 (1) Degrees of Freedom and Variables To analyze the centroid condition at each vertex, we assign the variables shown in Figure 2.5b. There are four kinds of vertices: F's: The three fixed vertices of the triangle, with zero degree of freedom. Q's: The vertices on the boundaries. They are only al- lowed to move along the boundaries of the triangle so that. with one degree of freedom each. I's: The interior vertices, with two degrees of freedom each. C's: The interior vertices which could be expressed in terms of the F's, 0'3, 1'3 through the centroid condition. Now we introduce variables to describe these degrees of freedom. We use 11's to express the position of 01 relative to a vertex along the boundary and use x1, y1 as the horizontal and verticle coorditions of 11' In summary, we list the functional coordinates of all sites in unit-2 as the following: F1=(0,0) F2=(cos36°, sin36°) F3=(T,0) (2.2a) Qi:(li’o) where i=1,2,o--,5 Qi=[5 + licos36°, (1-11)sin36°] where i=6,7,8 5“ Oi:(licos36°, lisin36°) where i:9,10 (2.2b) I1=(x1, y1) (2.2c) C1=(QZ+Q3+09)/3 C2:(I1+Q3+09)/3 C3:(I1+F2+09)/3 C“=(Q“+Q6+Q7)/3 (2.2d) (ii) Symmetry and Centroid Conditions Let us consider the site interior of the unit cell at first. Remembering that for C's the centroid conditions arebuilt in their definition, we only have constraints from the centroid condition for I's sites, i.e. ) = 0 (2.3) for each i, where j are nearest neighbors of i. It is a little bit more complicated for those sites on the boundary. Let us go back to the global structure shown in Figure 2.“. We can see that all of the three sides of (”thz are the reflecting symmetry planes. One comes from 55 the interior symmetry of the cell, the other two come from the tiling. This observation enable us to be able to write' down the centroid conditions for those sites on the bound- aries by assigning an appropriatelweight to each bond. The constraint at site i which is on the boundary is: A 1k - E ( 01 - AJ) wiJ = o (2.“) where j's are nearest neighbors of i, Yk is the unit direc- tion of the side to which i belong to. Particularly for unit-2, we have: Y=(1,0) .< u ( cos36°, sin36°) ..< 11 (-cos36°, sin36°) (2.5) and w is the weight of bond between i and j: 11 H 1 if j is on the same boundary side as i, i.) w 2 otherwise. (2.6) 1.1 (iii) Equations and Solutions 56 We input these variables and conditions into a computer program by algebra language SCHOONSHIP to derive the set of linear equations for variables 1 and x in terms of the 1 1'31 golden ration T. We express our results in terms of the matrix relation: A X = B (2.7) where for unit-2: -3 T -10 2 1 3 3 3 ‘ 2 :ll 1 £1 2 3 3 6 3 3 i 1 1 1 3 1 3 3 2 1 3 T -_I _, l3 -_2. 3 3 3 3% 3% 1% -T 2836 3 -T 2836 1 21 _ g; 11 T T 1 -5 -_1. _, :1 1 1 13 3 2 2 3 3 s36 s36 €336 1%- 57 and ET: ;g_r_ 7.11 “(1-1) _ _-_2 _7__r_ 1 (0 0 O 3 -1 3 3 (1 I) 3 O 6 3536) where 336 = sin 36° = /(3-1) / 2. By substituting the solution of this equation back to eq.(2.2), we can calculate the coordinates and the bond lengths, which are listed in Table 1(b). 2.3.3 Discussion about the Fixed Vertices Hypothesis We have been concentrated on unit-2 and deliberately did not take unit-1 simultaneously. Now we will consider the inteference between unit-1 and unit-2. At first, let us study the symmetry associated with unit-1. In Figure 2.5, we denote the three sides of unit-1 and unit-2 by a, b, c and A, B, C respectively. Careful ex- amination of Figure 2.5 tells us that the side b of unit-1 is always adjacent with the side 8 of unit-2, while the side c of unit-1 is adjacent to c of unit-1 25 C of unit-2 alternatively. This implies that unit-1 does not have exact enviroment as unit-2 so that it is not suitable to apply the same procedure described as above to solve unit-1, because not all three sides are reflecting symmetry planes any more. However, as we mensioned before, the unit-1 has the same topology as the shaded area of unit-2. Taking this fact into 58 Table 1 Solution of the Centroid Condition on the Unit Cells (a) Epit;1 i Xi Variable 1 0.30833“5“ l, 2 0.61296709 12 3 0.83629833 1. 0.30351959 1. 5 0.626“9066 l, 6 0.312“6202 l, (b) Uplt;g 1 X1 Variable 1 0.16“71666 l, 2 0.3899002? 1, 3 0.69890“36 l, “ 1.0187573“ l. 5 1.25“3825“ l5 6 0.696126“? 1.' 0.3922“76“ l. 0.16665261 1, 9 0.629618“0 1, 10 0.305“01“9 11o 11 0.90719375 X1 12 0.29762““3 Y1 Table 1 (cont'd.) (c) Unit-3 O‘U’I Q 10 11 12 13 1“ 15 16 17 dOOO fl GOOD 0 0. xi .16379“00 .38788255 .6960599“ .01615“60 .25276893 .30525666 .61072082 .83913176 .00978“10 .30786311 .16“03“98 .38771008 .69235305 .3672619“ .50803990 .76697177 56379“5“ 59 variable Table 1 (cont‘d.) (d) Unit—“ 10 11 12 13 1“ 15 16 17 18 19 20 21 22 23 2“ X. l .16331“57 .38681679 .69“““0“5 .01“65665 .25235023 .30988959 .60737915 .7780“903 .00660513 .31239882 .25230592 .01537787 .69“63300 .38“61“50 .15“90521 .982“1525 .62387806 .30280186 .69621370 .2328““15 .836“1“07 .09678353 .O“2171“5 .156“9“22 6O variable X1 Y1 X2 Y2 Y3 Table 1 (cont'd.) i 25 26 27 28 29 30 31 32 x1 1.11907083 0.3692308“ 1.21627382 0.658762“2 1.50“9889“ 0 .36993531 1.90282232 0 .299“8289 61 variable X. Yu X: Y: 62 account, it will not be meaningless to solve unit-1 inde— pendently under the assumed constraints that all three sides are reflecting symmetry planes and see what happen. Table 1(a) lists the results about unit-1. Now let us come back to the global structure (Figure 2.“) to investigate the tiling of the whole space by these centroid unit cells. We have the following observations: (1) Two sides of unit-2 are not only self reflecting but also should match with two sides of unit-1 respectively. (2) Two sides of unit-2, B and C, are not self reflect- ing simultaneously, i.e. if the B or C side is reflecting, the other one must be adjacent with unit-1. (3) Only the side 0 of unit-1 is self reflecting. The other side b is always adjacent to B of unit-2. We express these observations in terms of the Adjacent Probability Matrix AP: a b c A B C al 1 O 0 0 0 0 bl 0 0 0 0 pbB 0 AP: c' 0 0 pcc 0 0 pcc (2.8) A1 0 0 0 1 O O 3| O pBb ° 0 pBB 0 C1 0 ° pCc ° ° pcc where p08 is a number less than one denoting the porbability of corresponding adjacency of cell a and B. 63 Figure 2.6 shows all four possible adjacent situation of unit-1 and unit-2. The double shaded triangles are the host unit-1 and unit-2. The triangles around the host are the possible adjacent cells. Considering in combination the host cell and enviroment we find that the shaded area in each case forms two configurations as shown in Figure 2.73 which are called unit-3 and unit-“ respectively. We apply the same procedure used before on these bigger cells. ,The results are listed in Table 1(c) and Table 1(d) respectively for unit-3 and unit-“. Now we are at a point to check the validity of the "Fixed Vertices" hypothesis. If it is correct, then we would expect: (1) In Figure 2.6, those points marked by "A" in middle of one side should divide the corresponding side into two parts with length ratio 1:1. (2) The two independently solved configurations unit-3 and unit-“ should have sides matchness as shown in Figure 2.8. In Table 2, we showed the degree of incoincidence of the boundaries of the adjacent unit cells which are supposed to matched each other, as shown in Figure 2.8. In Table 3(a), we listed all the bond lengths solved independently from the four unit cells, with the bond se- quence number as shown in Figure 2.9a. The higher order unit cell includes the topology in the lower order one, as unit-“ includes unit-3, etc. We listed the corresponding bonds in an the same row for reader's convinence. The number between is the relative deviation of the two bonds adjacent in the table. In Figure 2.9b, we can see that the larger unit could be chopped into smaller parts with the same topology of the lower order unit cells. we compared these smaller parts in Table 2(b) and (c). Unfortunately, none of our proposed checks has been verified exactly although the deviation is small. We are forced to admit that these points suggested by numerical simulation to be "stationary" are not exactly fixed. 2.3.“ Conductivity of Infinite Large Penrose Tiling Although we showed the "Fixed Vertices" hypothesis not valid exactly, we are still rewarded by this research. We calculated the bond lengths in each unit cell. Through the mapping described in Chapter 1, we know that the average bond length gives the conductivity. To calculate this quantity, each bond in the unit cell is associated with a weight. Depending on where the bond is located in the cell, there are five weights as listed below: (also see Figure 2.9a) P1 = T ( A - inside ) P2 = 1 ( B - inside ) P3 = 1/2 ( B - boundary ) P“ = 1/2 ( A - boundary ) 65 P5 = 1/2 + 1/2 = 12/2 ( A-B boundary ) (2.9) where A is the shaded and B is unshaded area in Figure 2.9a. The formula we used to calculate the average bond length is: < 52> ={9(1) Di /{P(1) (2.10) i i where i is summed over all the bonds in the cell and P(i) has the weight listed in (2.9) depends on where i-th bond is. We finish this section by listing the results about the conductivity as the following,: G (unit-2) = 0.9362 G (unit-3) = 0.9333 0 (unit-“) = 0.9316 (2.11) Result for unit-1 is not included since unit-1 cannot. be futher chopped into two parts which is necessary to tile the whole space. We can see that the results are convergent in going to bigger and bigger unit cells. This can be un- derstood through the fact that when tiling the infinite space by two bigger cells, the error due to the unmatchness between cells, i.e. the unsatisfactory of the centroid con- dition on board between cells is reduced progressively. In 66 principle, to get a result as precise as required, one can always build up bigger and bigger unit cell systematically, and solve the centroid configuration as we have shown in this section. 67 (c) K Figure 2.6 Enviromental configuration of unit-1 and unit-2 in Penrose tilirg. ' 68 Fime 2.7 Basic cells (a) wit-3 and (b) unit-“. 69 Table 2 Degree of Incoincidence of the Cells Boundary (a) Comparison of Unit-1 and Unit-2 a Unit-1 a Unit-2 site [AL (1) a 0.626“907 0.629618“ 0.3128 b 0.3035196 0.305u015 0.1882 c 0.1637017 0.1666526 0.2951 d 0.3870329 0.3922“76 0.5215 e 0.6916655 0.6961265 0.uu61 (b) Comparison of Unit:3 and Unit-“ site Unit—3. Unit-“. IA! (1) A 1.3078631 1.3081““5 0.0281 B 1.0097811 1.01065“9 0.0871 C 0.8391318 0.83998“9 0.0853 0 0.6107208 0.611“288 0.0708 E 0.3052567 0.3056351 0.0378 9 1.2527689 1.2523502 0.0u19 0 1.01615“6 1.01u6567 0.1“98 H 0.6969599 0.69“““05 0.1619 I 0.3878826 0.3868168 0.1066 J 0.16379“0 0.16331“6 0.0179 * Number in the column is the distance from the site to the corresponding vertex in the unit cell. 7O H“."' C. o ,, 11'“ -- .. Table 3 Comparison of Bond Lengths Belgpg to Different Unit Cells 71 (a) Bond No. Unit 1 Unit 2 Unit 3 Unit 1 1 .30833151 3°736 .31985298 '076 .32009166 '203 .32071187 2 .16370167 '620 .16171666 '560 .16379100 '309 .17218996 3 .32297107 '386 .32121691 '38” .32297207 '81” .32560216 1 .21291151 '9‘6 .21513561 '2‘“ .21161007 '052 .21173820 53 .22198558 '302 .22265608 "72 .22303828 '08“ .22322501 6 .20117033 1°367 .20696165 '1‘3 .20672893. .215 .20723138 7 .18730620 '757 .18872109 '300 .18815819 '5"1 .18917661 8 .21611936 '557 .21782309 '089 .21760227 .010 .21710207 9 .22769156 '7"8 .22939720 '385 .22851323. '"32 .22950013 10 .19628166 '620 .19750167 '385 .19671089 '053 .19881272 11 .31216202 "732 .31781701 '555 .31961115 '09' .31993551 12 .21631506 '378 .21552818 '027 .21558575 “170 .21595227 13 .20328655 '7“7 .20180191 '730 .20629913. '060 .20612365 11 .23812550 '272 .23717681 .181 .23632762 "53 .23668902 15 .18361256 '3‘2 .18212236 "62 .18123796 "27 .18117191 16 .23562520 '“20 .23661138 "33 .23692805 17. .30387353 '“55 .30525666 "2" .30563§1§ 18 .16665261 °"°° .17065231 '010 .17066988 19 .23562581 '“30 .23663817 "30 .23991615 20 .19566756 '81; .19725803 “1‘9 .19719211 21 .22711872 '503 .22826175 "25 .22851621 22 .18116311 '5‘7 .18511713 '109 .18561917 23 .19976031 '“62 .20068212 “‘22 .20092813 Table 3(a) (cont'd.) 72 Bond No. Unit 2 Unit 3 Unit 1 21 .21396773 .111 .21191871 "27 .21519078 25 .16103198 "2‘ .16751191 26 .19633892 '998 .20026236 27 .22779152 '82‘ .22592235 28 .18722631 '3“2 .18786608 29 .21576815 '505 .21152781 30 .21268912 '2‘“ .21216976 31 .20381071 '373 .20308063 32 .22091070 '389 .21897710 33 .23792211 '13" .23821238 31 .20136111 .117 .20106871 35 .21580831 “‘68 .21511558 36 . 31017087 '09‘ .30988959 37 .30761695 '323 .30665175 38 .17939189 “03‘ .17911667 39 .32107620 “0 .23769358 11 .32021620 n; .16331u57 13 .19582051 11 .22759118 15 .18721711 u5 .21563965 17 .21382133 18 .23353502 Table 3(a) (cont'd.) Bond No. Unit 1 58. .17516676 50 .20168290 51 .21158133 52 .22092206 53 .20597108 51 .31686367 55 .19313323 56 .23006690 57 .18180567 58 .23763567 59 .20179370 w .mwmfi 74 Table 3 (cont'd.) (b) Unit 1 .308331511 .163701672 .322971073 .21291151“ .221985585 .201170336 .187306207 .216119368 .227691569 .19628166‘0. .31216202“. .21631506‘2. .20328655‘3. .23812550‘". .18361256‘5. .319852981 .161716662 .321216913 .21513561” .222656085 .20696165 .18872109 .21782309 .229397209 Unit 2 6 7 8 1975016710. 3178170111. 2155281812. 2018019113. 23717681‘“. 1821223615. .29961007 .20672803 .18815899 .24760227 Unit 3a .320091661 .163791002 .322972073 1 .223038285 6 one .228513239 196710890 319611151 215585752 206299133 23632762“ 181237965 Unit 3b .3076169537. .1610319825. .322972073 .2126891230. .2209907032. .2038107131 .1872263128 .2157681529. .2277915227 .1963389226 .3101708736. .2158083135. .201361113”. .2379221133. .1793918938. .20723938 . Unit “3 320711871. 172189962. .325602163. 21173820”. 223225015. 6 .189176617. 297702078. .229500139. .198812720. 31993551‘. 215952272. 206123653. 23668902”. 181171915. Unit 9b 32021620“‘ 16331157"2 3210762039 21382133“7 2209220652 2059710853 18721711u5. H6 11 21563965 22759118 19582051“3 316863675“ 2115813351 50 8 20168290 23353502" 17516676u9. .225922352 .30988959 .2151155835 Unit No 3066517537 .1675119125 .325602163 .2121697630 .2189771032 2030806331 1878660828 .2115278129 7 .2002623626 36 .201068713” .2382123833 1791166738 75 Table 3 (cont'd.) (e) Unit 2 Unit 3 Unit “a Unit 1b .319852981 .320091661 .320711871 .32021620‘n .161716662 .163791002 .172189962 .16331157“2 .321216913 .322972073 .325602163 .3210762039 .21513561” .21161007" .21173820" 21382133“7 .222656085 .223038285 .223225015 .2209220652 .206961656 .206728036 .207231386 .2059710853 .188721097 .188158197 .189176617 .18721711“5 .217823098 .217602278 .217702078 .21563965“6 .229397209 .228513239 .229500139 .22759118"“ .1975016710 .1967108910 .1988127210 .19582051"3 .3178170111 .3196111511 .31993551“ .316863675” .2155281812 .2155857512 .2159522712 .2115813351 .2018019113 .2062991313 .2061236513 .2016829050 .23717681‘“ .23632762‘“ .23668902‘“ “ 23353502“8 .1821223615 .1812379615 .1811719115 .17516676”9 .2356252016 .2366113816 .2369280516 .23769358“0 .3038735317 .30525666‘7 .3056351617 .3066517537 .1666526118 .1706523118 .1706698818 .1675119125 .2356258119 .2366381719 .2369161519 .2376356758 .1956675620 .1972580320 .1971921120 .1931332355 .2271187221 .2282617521 .2285162121 .2300669056 .1811631122 .1851171322 .1856191722 .1818056757 .1997603123 .2006821223 .2009281323 .2017937059 .213967732" .211918712” .21519078”4 .2030806360 76 Figur_e 2.9a Illustration of the 355th of weights to bonds in different region of wit. 77 Figure 2.2b Illustration of the similar topology between adjacent gmeraticns of units. 78 Section 2.“ Equation of Motion Approach In this section, we apply the well established Equation of Motion method on Penrose Tiling to calculate the DOS and S(q,E) on purpose to compare with the Centroid method result. 2.“.1 Equation of Motion Method Let us consider a spin system on .a Penrose Tiling lattice. On each lattice point there is a spin and the in- teraction between spins is of the Heisenberg form: H:—z J g ‘g (2.12) where the sum is over pairs of nearest-neighbor sites and J1] is the exchange constants between site i and J. We make the linear spin approximation via the lowest- order Holstein-Primakoff transformation 8+ . (25)1/2 a 1 1 s; . (23)”2 af (2.13) 1 Sz.¢S-a+ a 79 where S is the value of the spin (the same for all atoms) and a1 and a; are Bose destruction and creation operators. From (2.12) and (2.13) we obtain the quadratic, linear spin- wave Hamiltonian: H s J (a; a. - a* a 1 (2.11) =2 (1.1) ‘1 1 1 where pairs occur twice in the summation. In principle, the quadratic Hamiltonian (2.19) can be transformed to a set of normal modes q with energies Eq. What we are interested is to obtain the distribution of the Eq's (density of states) and also the scattering intensity S(q,E). We define a set of quantities gi (t): 15.83) g1q(t) = < a1(t) X a*(0) e (2.15) J J where 83 is the position of site J and < > indicates the ground state expectation value. Note that the sum is over all sites. The giq's obey the following Equation-of-Motion: 11 (giq- sjq) (2.16) with the initial condition 8i (t=0) = eiq' (2.17) The normalized form for the scattering intensity is defined as (Alben et al 1977): S(q,E) - ——l——— 1*“ eiEt/h < (I a (c) e'15’§1)({ a+(0)eia.fij )> dt ‘ ZINh -w i i J J = lim lim -—%fi—— Re 13 X 6-15-81 31 (t) em“h e‘”2 dt i»0* Tow ' i q (2.18) where N is the number of spins and E is the energy. The nor- malization is such that the energy integral of S is unity. The essential point of the Equation-of-Motion procedure is that if we are interested only in broadened spectra, we can do the transform in eq.(2.18) with a nonzero value for the damping constant I and a finite value for the time interval T. Specifically, we proceed the following procedures: (1) Initially assigning N values of giq(t=0) by eq.(2.17) for a given a. (2) Integrating the N simultaneous differential equa- tion (2.16) forward in time to obtain N giq over a time interval from 0 to T. 81 (3) Choose the damping constant A to give an ap- proximately Gaussian broadening function and compute S(q,E) by doing the integral in (2.18). The Density of States can be obtained by using the ex- 16-81 actly same procedures as above Just replacing e by an ia. random phase factor e 1, where the 01's are random angle specified for each site. The average over many sets of 01's gives exactly the density-of-states spectrum. 2.1.2 Calculation of Density of States Our purpose to calculate DOS is to extract the conduc- tance of the resistor network with the corresponding topology. He have shown in Chapter 1 that for a tight- binding Hamiltonian 1 H = - Z J [ li> 13 its density of state at long wavelength limit can be ex- pressed in terms of the conductance p (8»0) = -——!——— (1 25) 110 ’ In Figure 2.10 we show the result of equation of motion calculation of density of state for a 888-site circular Penrose Tiling. We applied the fixed boundary conditdxni. He found out the initial DOS is: 82 ._ 4m 1.. -.. ll LJ 1 [1 l L] l I l 0 cm Naommrxmmvmmuo «dfiV‘V‘OOOODOOOO SOD Fifle 2.10 Da1slty of statw of Penrose tiling by the E‘qation of Votion. 83 0 (8+0) = 0.085;0.002 from which we have 0/0. = 0.9310.02 which is consistent with the result in section 2.3. 2.9.3 Calculation of the Intensity Spectrum S(q,E) It is known that in an ordered crystal, S(q,E) would show sharp peaks 'at energies corresponding to selected ex- citation with the proper momentum. Hhile in the disordered alloy, all modes usually contribute for any q. However, for different q's, different modes contribute most strongly. The spectra are usually best thought of as the density of states with different regions of energy weighted differently as q varies. Because of the total spin rotation symmetry of the Hamiltonian for the1disordered alloy S(q,E) for qu has a rather special behavior, even in the absence of the transla- tional symmetry. It is easy to show that for q=O, S(q,E) is given by a 8 function of unit weight at E=O. For q slightly different from zero, there is always a regime of very small q such that S(q,E) is sharply peaked at an energy S(q) given by 81 S(q) = D q (2.19) where D is the spin stiffness. The width of the peak varies as q”. In Figure 2.11 we show a series of spectrum S(q,E) for different value of a in x direction. Since we used a sample of finite size, we observed multiple discrete peaks. To analysize the contribution of a particular mode at energy E,: we replotted the spectrum S(q,E) vs. q for a value of E at which there is a peak in the S(q,E). The result is shown in Figure 2.12. From which, we find out the position of q cor- responding to the peak through fitting and established the relation between this mode E and the q-value which makes this mode contribute most to the spectra. Figure 2.13 plots this correspondence, i.e. E vs. (qa)2. Its slope gives the spin stiffness constant D=O.9H;0.01. S(q. E) S (a. E) 5111.6) 51‘]. El 85 ' ' ' 1 l (qa)**2=O.1O j, 0.15 ‘ G 20 g, (n 10 '1 0 RNA. d 1 .2 u A .. a 0.: 04 0.: ca 10 0.6 09 E E 40 1 1 1 ‘0 I 1 1 r 0.20 0,25 30- - 30.. NP ~ 20- Slq. E) ( oz 04 0.5 a: 10 c 02 04 us 0.; E E ‘0 1 r I r ‘0 l I r 0.30 0.35 30r- - 30- 3 2°— " U.’ 20!- 5 mu 4 10- A | o 1 J “11.311. 0 1111 ‘ 1.11. c o 2 o 4 0.0 0.9 1.0 c o 2 a 1 0.0 0.: E E 40 r 1 , 1 4° 1 , m 0.40 11 0.45 a 1 20'- - 'j 201- 6 l 10'- -1 20’ 1F 1 I! ) : ' . l' ‘ , 113:1 Mun 1 ,' 1 “C1141 H: C 02 ‘34 05 03 10 7 J2 )- J5 3: E 5 RM 2.11 Low mergy neutron scatterirg spectnm S(q,E) vs. E. 86 .c .n> Amévm .550QO 2.83% 8,58: $.98 33 NTN mama m 5. 35 a... a... Yo me a. -1 1‘ .1— d — ¢ q d C 4 Q Q 4 . _ r . . _ _ m xx .mc. -IIIqIII—IIIIJ1111141( G _ _ d < 0 4 4 o— A: m— ON mw (3 '5) S (3 'D) S m .2... BE .. a... a... To m... o q - q ‘ q q q IIHIQIJIIII a 4 1 4 1 o 1 o L . o T 4 l o o omomd - n b n p h > p b m 5. BS . a... a... To m... o 11.1.111— 114N111... 1.111.111 o 4 o .. 1 o o .1 4 1 mmomdnm . . p . . _ pl _ L O. m. cm mm C— am an O? (3 '915 (3 'D) S 87 .a 22m 65 805.8 B .6.ch .5588 5 was .9 swam rd Mama m xx $3 5.0 m.o n.o «.0 m.o N.o «.o o . _ . . . _ \ 1. \ J \ .. \\ . 1 \ 4 .l k .l \ I \4 1 W\ \ 1 \ 1 \ 1 \ 1 \ _ a a u a u 3 01980) 88 Section 2.5 Random Halk Approach It is well known that the differential equation of the probability distribution in the Random Halk problem is ex- actly the diffusion equation in an appropriate limit. In this section we attempted to simulate the random walk on the Penrose Tiling and try to extract directly the diffusion constant. In random walk problem, a particle undergoes a sequence of displacements F1, 32, ..., Fi’ ..., the magnitude and direction of each displacement being independent of all the proceeding ones. More specifically, in the background of the Penrose tiling, the particle can move to any of its nearest neighbors in one time step. Because we are working on a finite size sample, we have to decide what will happen to the particle at the time it reaches the boundaryu. 'rwo spe- cial boundary conditions are: absorbing and reflecting boundary, as illustrated as follows. Reflecting Boundary Let us imagine the boundary of the region under con- sideration is a perfect reflecting wall. In the isotropic two dimensional case, it is*a reflecting circle at r=r0, which simply has the effect that whenever the partical ar- rives at r0, it has a probability unity of coming back to the interior region when it takes next step. This condition can be expressed mathematically as: 3 P r t _ (2.20) where P(r,thkfi is the probability that the partical will arrive at rSrO after t steps. Absorbing Boundary Now we interpose a perfect absorbing wall on the bound- ary. The presence of this absorbing wall at r=ro means that whenever the particle arrives at r0, it at once becomes in- capable of suffering further displacement. This-boundary condition is written as: P(r=r - O (2.21) .t; r0) - 0 where P(r,t;r0) is defined as the probability of a partical at r at time t. One interesting question which is charac- teristic of the present problem concerns the average rate at which the particle will deposit itself on the absorbing screen. In our computer simulation on the Penrose tiling, we will use the absorbing boundary condition. The statistical quantity we actually observe is the deposit rate, i.e. the reciprocal of the average time (t) for a particle to reach the absorbing boundary. In Appendix B, we derive the rela- tion between the diffusion constant D, the sample radius ro and the deposit rate 1/08». In 2-dimension, with circular sample, and initial distribution of the form 9O ‘ P (r,0) = 8 (r) (2.22) we have r 2 0 = 0 (2 23) (t) ' He will use this relation to extract D from the simulation data. He work on the circular piece of the Penrose tiling of 888 atoms which has been used in the Equation of Motion calculation. Initially, the particle is confined in the small region around the center of the sample. He measure the time (number of steps) t between from the particle starts to move until it is captured by the detectors on the boundary. During this period, at each time step, the particle can jump from its current position to any one of its nearest neigh- is the bors with equal probability 1/z where z i’ i coordination number at site 1. Our result averaged over “0000 runnings gives D = 0.91:0.01. 91 Section 2.6 Conclusion In this chapter we applied three techniques to study the low energy behavior of 2D Penrose tiling. Our results for the diffusion constant D, which is equivalent to the conductivity, of an infinitely large Penrose tiling are sum- marized in Table 1. All three techniques give the consistent results, although with different ranges of error. Centroid Algorithm takes the advantage of the combina- tion of self-similarity property of the Penrose tiling and the geometrical character of the conductivity problem through the mapping. A systematic routine is established to obtain the result which can meet any precision requirement. By solving the centroid configuration of a unit cell with a figltg number of atoms, one can obtain the centroid con- figuration of the infinitely large Penrose tiling. In other words, the conductivity or the diffusion constant of the in- finitely large Penrose tiling can be obtained from unit cell with finite number sites. It should be emphasized that it is the intrinsic self-similar property of the quasicrystal which makes this possible and necessary. In a crystal there is a unique basic cell such that one does not need the series unit-cells, while in glassy structure there is no unit cell in any order. It provides a method to study physical properties, especially those depending strongly on the geometrical structure, by exploring the intrinsic struc- ture of the quasicrystal. 92 Table 1 Summary of the Diffusion Constant on 20 Penrose Tiling CENTROID METHOD numerical simulation 0.936 centroid unit cell 0.9366(2) 0.9333(3) 0.9316(u) EQUATION OF MOTION S(q,E) 0.9“ :_0.02 DOS p(EoO) 0.93 1 0.01 RANDOM WALK 0.91 i 0.01 93 The other two techniques, the Equation of Motion ap- proach and the direct Random Halk approach, have been used widely in dealing with the random system. Because of their numerical simulation characteristics, the results are af- fected more or less by the sample size. Moreover, since for quasicrystals it is unsuitable to use the periodic boundary condition as one usually does for random system, one always encounters a finite sample and its boundary sites are not properly co-ordinated. Chapter 3. Percolation on Stretched Elastic Networks Section 3.1 Introduction In Chapter 2, the visualization of the conductivity problem based on the mapping introduced in Chapter 1 was ex- plored on a new structure --- quasicrystal, which distinctly lies in between conventional crystal and glassy structure. In this chapter, we will explore the implication of the map- ping on a deeper physical level. He will propose and study a generalized "tennis racket" model. It is widely accepted that conductivity percolation and rigidity percolation belong to two different universality classes. (Feng & Sen 1981) The percolation threshold p of a central force network, at which the elasticity cen vanishes, is much higher than pc at which a resistor network loses its conductivity. These two systems also have quite different scaling exponents in the critical region, near p and pc respectively. Mean field theory, as well as cen numerical simulation results strongly support this classification. Hhat is interesting is that our mapping seems to suggest the existence of a bridge between these two 94 95 classes. Our working frame, the mapped tennis racket, com- posed with Hooke's law springs (natural length zero), is certainly a central force network. The problem it solves is identical to the conductivity problem. This can be under- stood by noticing the fact that the percolation in either resistor network or central force network with natural length zero is determined by the geometrical connection. As long as there is a percolating cluster, the current always can flow from one side to the other. Similarly, there is always energy stored in the mapped elastic network through the backbone. This feature is different with the central force network studied by many other authors. The networks they used are such that every spring has its natural length in the absence of external stress. Based on the above consideration, a new model --- the tennis racket model is formed naturally, which seems can fill the gap between conductivity percolation and rigidity percolation. The model is defined as follows. On a lattice background with space constant L, we stick at the place of each bond between nearest neighbors a Hooke's law spring with natural length LO, which need not to be equal to L. The shape of the lattice is kept through a frame or periodic boundary conditions. In the equilibrium state, no net force exerting on each site, but the springs adjacent to this site might be stretched, which could be considered as internal strain. The physical behavior of such a system must be very interesting since it reflects the nature of a lot of real 96 systems. For instance, the percolation threshold p. is one of the interesting quantities. Hhat we already know through constraint counting (Thorpe 1983) and the mean field theory' as well as numerical simulations (Feng, Thorpe & Garboczi 1985) are two special cases: L0 = O, which mapped to the conductivity percolation and L0 = L, which is the pure central force network. The effective medium theory gave the thresholds corresponding to these two cases: I' p (LO:0) 1 'O u and 'O A I" O u l'" V n '0 u cen where z is the coordination number and d is the dimensionality. These are very chose to the numerical simulation results p0 = 0.317 and pcen = 0.65. The question is how p. changes when LO takes arbitrary value ? It was noted that Feng & Sen (1981) observed a cross? over in their numerical simulation by using the Born model. They found that when the contribution of the central force part is dominant, there is a strong crossover from isotropic-force like behavior near p:pc to central-force like behavior near p:p Tris crossover was used as cen' strong evidence to classify the two universality percolation classes. As we will see in the next section, this case is 97 corresponding to the case of, in our tennis racket model, that the natural length Lo tice space L. However, it should be pointed out that Born approaches from below to the lat- model and the tennis racket model are different at least in two respects. First, Born model is actually the contribu- tion to elastic energy in the second order of the displacement, where the isotropic term is not rotationaly invariant. Hhile in the tennis racket model, we are con- sidering the total central force potential in the system. It is rotationally invariant since only the distance r1J in- volved there, of course we have to rotate the frame together with the lattice in the simulation which is done by keeping the periodic boundary conditions. Secondly, all the dis- placements in Born model are relative to the lattice points, while in tennis racket model, the displacement is defined relative to the new equilibrium position. There is also an important conceptual advantage in the tennis racket model because the pressure enters in a natural way into the determination of internal stresses and there- fore into the determination of elastic constants. He believe that a careful research of the proposed ten- nis racket model would be very helpful to understand the macroscopic elastic theory through the microcopic structural modeling. This chapter is organized as the following. Section 3.2 presents the elasticity analysis on a pure triangular 98 network under tension. Section 3.3 describe a diluted ten- nis racket with internal stress and details of the simulation techniques. Numerical simulation results for the overall behavior of the diluted tennis racket are presented in section 3.1. Section 3.5 attempts to understand the results from the symmetry point of view and reduce the num- ber of independent second order elastic constants. Section 3.6 presents an effective medium theory on the stretched network. Section 3.7 gives the conclusion for the current work and suggestion towards further research. 99 Section 3.2 Elasticity of Pure Triangular Net Under Stress He start from an elastic network on a pure lattice where all bonds are the Hooke's law springs with force con- stant K. The natural length of the spring is L0 and the lattice space is L. Except for the case LO=L, the network is stretched (O S LO/L < 1). The elastic potential for this model is: (3.1) where the sum is over all nearest neighbors and rIJ=IR1-RJI is the distance in between which is just L. He introduce 51 to express the small displacement of site 1 relative to its equibrium position. It is noticed that in the pure lattice case, the equilibrium positions are exactly the correspond- ing lattice points, while in diluted case, these two might be far away from each other.) The potential can be expanded ij=ui'uj' To about the equilibrium position in terms of G the second order of u. ij’ we get: 1 2 V=-2-x£ (L-LO) L L . 1 0 2 0 + 5 K X { (1-E- )u (3.2) (i,j) ‘1 A where r is the unit vector between points i and j. 1.1 The first sum in (3.2) is the static energy of the stretched springs and can be written as 1 2 VO - 1 K N z (L-LO) (3.2a) where N is the total number of sites in the system and z is the coordination number. The second sum in (3.2) is the linear order contribu— tion, which is written as: A v1 = K (I.-.LO) X (611- r ) (3.2b) 1.] The value of V1 could be positive, negative or zero, depend- ing on the details of the displacement u or the strain in ij’ the network. He note that the equilibrium condition for the full potential is: Q) < 11 O Q) '10 and for the expansion of the potential in terms of displace- ment, we have: 101 Q) < l O 0) CO both of them lead to the same equation : X K1J(r1J-Lo)rij :0 The last sum in (3.2), i.e. the quadratic term in uij is: 1 0 2 2 L (51) 13 + l K 29 2 (u or )2 (3 2c) 2 L <1,J> ij ij ' Comparing V2 with the Born model potential: v = g (0-8) X 1 (61-0 )-$ 1* born (1’1) j ij 3 z (a -6 )= (3.3) He can see the correspondence between the Born model parameters a, B and the tennis racket parameters LO/L and K: L0 K(1- E_ ) = B (3.33) 102 L0 and K L_ = (a-B) (3.3b) which implies: a = K (3'30) Generally, the strain-energy is defined as: (Barron & Klein, 1965) 1 2 ath $03671 60811 (3.1) where 808 is the stress tensor, 60871 is the second order' elastic constant and the strain is: £03 = a (3.5) For isotropic external strain, we have the expression for $08 and CQBYT as the following: as =1<(1.-1.o){ a 1- (3.6) 08 <1,J> ij and A CaBYt ) X R? r? rY 6 (3.7) where A is the area of the sample, bor distance L and r is the unit vectors. 11 For the triangular net, we have: Rij is the nearest neigh- A:( /3 L2/2 ) N (3.8) and A r = (1,0), (-1,0) 1J (1/2, /3 /2). (1/2, -/3 /2), (-1/2, /3 l2), (-1/2, -/3 l2). (3.9) A From which the moments of the r can be worked out: 13 a - a B - 1 < rij rij > - _2—'8aB (3.10) a 2 B 2 1 1 < ( ”13 ) ( ’11 ) > = ‘8’ * "1"508 Substituting back to (3.6) and (3.7), we obtained the stress tensor and the elastic constants. The results1are listed as the following where we used n to denote L /L: 0 Stress tensor: 101 U) 11 m 11 xx yy /3 K ( 1 - n) (3.11) Sxy = Syx = 0 (3.12) where Sxx and Syy are actually the external pressure. Second order elastic constants: cxxxx . C11 = 13F! 1 1 - n ) (3.13) cxxyy = 012 = 1335 n (3.11) cxyxy = C11 = 1335 ( 1 - 3 n ) (3.15) nyyx = C11 = $135 n (3.16) These results are plotted in Figure 3.1. As expected, these elastic constants include the known results for the pure scalar case and pure central force case. For scalar case: (n = 0) and for central force case: (n = 1) C11:3/3K/u 105 106 and C11 = /3 K / 1 He checked these results through computer simulation, they all agree well. Since neither C11 nor C comes out 12 independently, we design other four cases of initial stress and combine their results to extract C11 and C12. The transformations corresponding to these four cases, in addi- tion to those for computing C and C11 , are sketched in 11 Figure 3.2. The corresponding second order elastic constant are respectively: Pure shear: (Figure 3.2c) -J. . - 112 Pure rotation: (Figure 3.2d) u = g ( c,“ - 0,1 ) = 1335 ( 2 - 2n ) (3.18) 1‘ Bulk compression: (Figure 3.2e) - 13—5 (3.19) -.1. ' B'2(C *012)' 2 11 Compression/expansion: (Figure 3.2f) )-i3—5(2-n) (3.20) 11 12 7 1 107 11-" o ,, I 0.. x'=x(1-8) .. X"X I y'--y 0 y'=y+8x '1 f L...-——-°“""""'_ (a) (b) _____...—-—""‘I Ti,- 1 I . Ix x'-x+cy ‘ x'=x+cy IY ’Y"‘8X __________..4I . y'=y-sx (c) , .(d) 1. 1 v 1 I x=x(1 e) I x'=x(1-c) 'L y y(1-8) I y'=y(1+s) (e) : ‘11) A Figu_rg' 3.2 Transformations for uniform external strain. 108 Figure 3.3 plots these results. From these combinational elastic moduli, we have: C12 : B - b (3.213) and C11 = u - u (3.21b) Note that C and C11 are also given by ”s , ur , B and 11 109 1.0 1V1 1 1 1 c1 c; o ° unpow 0119913 B 6 .4 02 0.6 0.8 - 1.0 0.4 0.2 113 and b vs. rFLU/L for pure trialgular lattice. Figure 3.3 13, up, 110 Section 3.3 Diluted Triangular Net Under Stress ---- Simulation Technique Special care is needed in studying the diluted network under stress. For a pure lattice where no bond is missing, the lattice node serves as the equilibrium position because the forces exerting on a site through adjacent stretched springs are balanced each other. However, as soon as ther diluting process starts, this balance is destroyed and the network deforms around the missing springs. Related sites relaxed coherently to reach their new equilibrium position. It is around these new equilibrium positions, which might be far away from the corresponding original lattice nodes, that the sites in network are vibrating. The elastic constant and the stress tensor should be determined from the struc- ture of this new equilibrium configuration. The key point is that before applying any extra exter- nal stress, we have to obtain the relaxed configuration of' the diluted network. Only internal stress plays the role in this relaxation. A certain amount of static energy is stored in this relaxed network. The details of the numerical simulation are described as below. (1) Generate a diluted triangular lattice by removing bonds randomly with probability (1-p) on a triangular lattice. Note that the resulted network is unrelaxed for the cases L0 = L. 111 (2) Relax the network obtained from (1). The movement. of site i is proportional to the force exerted on it. Specifically, during the procedure of iteratdtni, the posi- tion of atom i is determined by the previous configuration through: x(n+1) = (n) (n) {1 x £1 + a ng (3.22) where F£1 is the force component given by: - ‘ § ng - X Kij (r1J - LO ) rij (3.23) J and a is a number which could be adjusted to control the simulation.’ Principally, the iteration process should stop when the force on any site is zero. Hhile practically, we choose an appropriate small number which is good enough to satisfy the requirement of the final result precision. For example, Fmax = 10'8 was used corresponding to e = 10"2 in our simulation. The energy of the relaxed network is calcu- lated and written as Eo(p,LO/L). (3) Now, depending on whether C11 or C11 or other elas- tic moduli is being computed, the coordinates of the sites are transformed corresponding to a certain external strain. Some uniform transformations are given in Figure 3.2. (1) After exerting the extra strain, one then proceeds to relax the network again. To seperate the linear and 112 second order contribution, we reverse the sign of external strain e and repeate the procedures (3) and (1). The energy of the relaxed configuration is: m u E (p, LO/L, e>0) I") 11 E (p, LO/L, £ 11 (R13 11 Note that the coefficient of the second order term is the (r) bulk modulus. 0n the other hand, 5 is corresponding to the "pure rotation" transformation defined in Figure 3.2d: 125 3‘ II x + m y y' = y - w x (3.39) from which we get the second order elastic constant, / R. ) R ‘3 (3.240) and the corresponding linear order term is always zero. Figure 3.10 presents a check for eq.(3.28) via the numerical simulation on the diluted network. We conclude that the rotational invariance is preserved in the diluted case. Isotropy for the Sound Velocities Now let us consider the eq.(3.29). With the given energy of deformation (3.”), the wave equation is: 2 . o. 9 Z X C a uY (X) (3 ‘41) p u (x) = ----- . a > Y 31 aBYt 3 x8 3 x1 For the plane wave solution: 6 exp ( 1 E . I - i w t ) (3.u2) we get 126 n n n n ‘1" 06 1.0 2.0 1.5 1.0 0.5 \ ' \\ \ \\| v \ h (0 a\ \ \‘ I \ \\V I— ' h o " \ 0 0.0 \ [30 N — d 0 O I l 1 o O = P for diluted (Banks. .10 Sinulatlal remit of C1111 - C1111 Fl 127 p “2 u = § ( X caBYT k8 kt ) uY (3.13) Specifically, for 2-D system, (3.“3) is reduced to: 2 2 O 2 u C11 kx + CH“ ky (C12 + can )kxky u p m = 2 2 v (C12 + Cu“ )kxky Cuu kx + CH ky V (3.HA) In the (1,0) direction, we have two eigenvalues: p m / k = C (3.05) MM Corresponding to the longitudinal and transverse wave respectively. In the (1,1) direction, i.e. kx = ky : k/l2, the two eigenvalues are c + c i ( c + c 1) p “2 / k2 = 11 nu 2 12 nu (3.u6) Since for isotropic system as the triangular net, the sound velocities are independent of directions, by equating (3.“5) and (3.06), we have C11 + C0” + ( C12 + Cuu' ) 2 C11 (3.07) 128 C + C 11 " ( C + C ' ) = 2 Cu“ (3.u8) NH 12 MN both of them leads to equation (3.29). Figure 3.11 presents a numerical check for (3.29) on the diluted system. He con- clude that eq.(3.29) based on the isotropy is valid in the diluted stretched system. Cauchy's Relation Finally, we study the validity of the Cauchy's relation on diluted network. From the numerical simulation results plotted in Figure 3.12, we can see that there is an obvious discrepancy between C12 and CAM. except for LO/L = 1. This result may be not so disappointing since the hypothesis for the Cauchy's relation to be valid is totally destroyed in the diluted tennis racket. Assumptions for Cauchy's rela- tion are (1) central force between pair of atoms, (2) inversional symmetry at each sites. (Love 190“) Garboczi (1985) claimed that the second assumption should be modified since that the numerical simulation suggested the existance of Cauchy's relation for the diluted central network, 1H1ich is the case L /L = 1 in this work. As soon as the network 0 is stretched, i.e. for L /L not equal to 1, we lost the O Cauchy's relation. One explanation is for L /L = 1, the 0 equilibrium configuration still has the lattice background, while for LO/L = 1 and p<1, the equilibrium positions are no longer at the lattice node. 129 1.0 “K I \\ \\ \\\ \\ ‘15 11 e119. .0:\ 0s\ aback CD — "J ' 6h\\6\n\ \\ Q\\ o a any \ "\ \ \ - 3‘03 3‘- ‘3 \ \ \ \ 113111 200 \ \\ v - t '3 x. a 0'0 ‘ 1:10 _ N " 0 l L I o 1.0 0 ID 0 .3 .4 o' 2.0 = C11 — C12 for diluted networks. . 0 Figure 3.11 Sirmlation mlt of CmoC1m 130 . N— ...=0 .1. 0 8539. 903:8 .50 01%! 533250 mp.” .91me 0 0 a... o o» q o m o no 9.. o o o... v o m... o - u u q - a u - ohm 0 an. 1 o 1 . o 1 o ..o a... . 1...... 0 B O T O G J “.0 l O n— n— N.0 o B a o n. 1 o a 1 n o 1 o n. a n... a m a , 1 To . To fl 00.01: 0.01: n p p n “.0 n p — p “.9 0 m o a o. ad a... v... «a o . 0 _ q a . 0 o . O B . 1 n o T o m m B u o o m a 1 u o 1. a a... I 1 a... 1 n... 0' ° 1 1 To 1 2w 0 To 5.0.10 0.0:: p n L - ”.0 b p - n “.0 131 Through above analysis, we conclude that in the diluted tennis racket, eqs. (3.28) and (3“29) are still valid, but not the Cauchy's relation (3.27). This reduces the six quantities E, P, B, ”s’ ”r’ b to four independent ones. 132 Section 3.6 Effective Medium Theory Now we have four independent quantities E (static, energy), P (pressure), B (Bulk modulus) and us(Shear modulus) to describe the elastic behavior of the tennis racket. Figure 3.8 and Figure 3.9 plot P and E vs. p for various n, where n=LO/L. To see the transition effect clearly, both E and P are scaled by their corresponding value on pure lattice. Similarly, Figure 3.13 plots B and “3 vs. p for several values of n, which are averaged over 5 samples. Note that for n=0, B and 118 are exactly equal. Obviously, p., at which the elasticity vanishes depends on n. Figure 3.1“ plots this phase boundary of p. vs. n in the stretched region, i.e. 0 < n < 1, where two extreme cases n = 0 and n = 1 were known before. It can be seen from Figure 3.13 that B is almost a straight line for whole range of p down to p”, while as , as well as E and P which are shown in Figures 3.9 and 3.8, have more or less a tail near p.. This tells us that the phase boundary in Figure 3.1” may be relative to the phase bound- ary determined by the initial slope of B. We introduce a single defect in a network, then let it relax. From the drop of the corresponding elastic quantities due to this' defect , we extract their intercepts on the p-axis. Figure 3.15 presents the computer simulation results of the single a u a n a a defect case for p E’ p P’ p B’ p u , p u and p b' Again, 1" 8 133 . .5113: ..o ..Bmé 9: 3 828m 3 a: 22: 3.95.» m co>o 809.96 .a .n> 329.23 : Em 3.5.63 0 m_.m 9a.... 0.0 0 N.0 . 0 . 1 a... 1 I. 'oo I1 1 0.0 1 I ..o l p . _ p 0.“ . p . . a 0 0.0 0.0 0.0 0.0 «.0 0 0.. 0.0 0.0 0.0 «.0 d _ 0 . . fill-1| .l I ~¢° .l. I I. m I a I .I a I. '.° rl . I. a I a I D I a I I: 0 TI 1 o o I 0 D 1 d “cal-.— coca—.- . . . L 04 . b p . v.0 134 3:398 0:86 ..o 9.629 E 8.3.88 >958 8m: =3” an E 0.0 0.0 v.0 0.0 5.0 135 nu .Qta . nan ...JIQ .le .mna .mta 1.0..“ 31.8. Exam—.350 5083—82 a A 0.0 0.0 v. N0 % _ _ _ a 11 01m ADV Q .1 O W 4 mm D s . .._ + 1 m .1 4 0 . D 4. . 4 .. a a a D Q D Q a a a a D _ m . o 136 C 1‘ d . u,’ pus an pb. I there are two identical pairs, p P and p as the results of the rotational invariance and isotropy which are proved in section 3.5. The differences between the four independent quantities are consistent with what shown in Figures 3.13, 3.8 and 3.9. In the remain part of this section, we will present a theoretical analysis of the single defect case and predict the initial slopes via the effective medium theory. Let us consider a tennis racket with one bond, for ex- ample, in between sites 1 and 2, removed as sketched in Figure 3.16. After relaxation, the distance between 1 and 2 changes from L to L The effective spring constant be- eq‘ a tween site 1 and 2 are K/a , which includes the effects of I both the direct spring and the lattice. We know that 0 < a < 1, since the effect of the lattice background. Then after removing the spring between 1 and 2, the effective spring of the lattice background is K = —-; - K (3.119) I From Thorpe and Garboczi (1985), a can be calculated exactly from the dynamical matrix of the pure lattice, ¢ a = -'i X Tr { X [ 1 - exp (1LE-3) ] (33) . 2"(k) } k a (3.50) 137 .539832 58m Ea b.5000 539500.050 ”.8000 2?? m5 53on B 9390 07m an 138 where g (k) is the Fourier transform of the dynamical matrix: Q (E) = x n X [ 1 - exp (1LE-8) ] 33 a + K (1-n) Z [ 1 - exp(iLE°8) ] l (3.51) 5 where nzLo/L and I is the dxd unit matrix. The formula to calculate a. of the tennis racket model on triangular net are shown in Appendix C. If the potential about the equilibrium position of this effective spring is pure parabolic (checked by simulation for stretched network), we find Leq is given by: = I (3.52) Substitute back to the expression of the energy, we find out the amount of energy that must be added to the one defect case to recover the no defect lattice energy is: I L - L a 1 o 2 1 2 111-:=§1055 208» .80 3.8.0 8 $89.35 83» 3:5 as a. 08% w: ...a .8 3.8.. rows 58.. 9385 a...” Mam... c 0. a 141 a and p P is given by i p P =1-1/A1 (3.58) For the Bulk modulus, we know 2 2 3 E “A 3L where A is the area. We have B = BO [ 1 - (1-p) A2 ] (3.60) where A2 = 1 a [ 1 _ n ( 1 _ n2) 1 . Ga 1-a 1-a 3n 1 * 2 + [ n (1‘n) g ] 1-a 1 2 2 1 32 ' +5 n (1-n) . 32 1 (3.61) 1-a 3 n which gives i p = 1 - 1/A (3.62) 142 I I I All of these results, p E’ p P and p B’ as well as the corresponding simulation results are plotted in Figure 3.18. 143 .msa Em man .msn 8,. 3:5. 2.59: .532. guuootm QTM mar. : m... m... To m.o . _ J _ m.o 144 Section 3.7 Conclusion and Future Research He found that two kinds of percolation problem, conduc- tivity percolation and rigidity percolation can be combined together in the stretched tennis racket model. Our main conclusion is that the threshold for the new model is not discontinuously changing from p to pc, as predicted by cen the simulation on the Born model. In the intermediate s region, i.e. for O1) is certainly needed which is an important part in the complete phase diagram. Secondly, from the ex- perimental point of view, a measurement of the longitudinal and transverse sound velocities in tennis racket must be very interesting. Of course, in a real system, there would be additional contributions to the elasticity which are not due to the network. Finally, determination of the standard exponents for the tennis racket model would be useful but also time consuming. APPENDICES Appendix A. Isolated defects in the static limit 145 Appendix A. Isolated defects in the static limit Submitted to J. Phys. C September 18, 1986 Isolated defects in the static limit H. F. Thorpe and H. Tang Department of Physics and Astronomy Michigan State University East Lansing, MI 48824 Abstract We collect together the known results for the effect of a single defect, into a single compact form that can be used in any lattice in the static limit. The result is used to calculate the effect of single defects, either bonds or sites, on the diffusion constant and conductivity of various lattices. It is shown that the result of Izyumov for the spin wave stiffness can be simply generalized to any lattice. PA Classification: 7155, 6350 146 1. Introduction The effect of a single impurity on some desired response function is well understood (see e.g. Elliott et al 1974). Practical calculations are possible if the range over which the interactions are changed by the impurity, is very short. In general, the detailed symmetry around the defect site must be taken into account in order to complete the calculation. An exception is a mass defect, where a simple formula can be written down that is valid for any lattice. The lattice only comes in when the particular on-site Green function needed is calculated. Another similar situation is provided by a site diagonal impurity in a tight binding model (see e.g. Economu 1979). In this paper, we look at two other cases that can also be carried through without worrying about the particular lattice until the pure lattice Green fUnction has to be evaluated. The first of these is the bond defect in a tight binding model; a result given by Kirkpatrick (1973). A more interesting and complex example is site defect when a.site and all its bonds to its nearest neighbors are altered. 2. General Formalism For convenience we consider a system described by a vibrational potential 0 = X A “B u? u? (1) iJaB iJ 147 where all the atoms are mass points with mass m; i,j are sites and 0,8 are Cartesian components. The vector 31 is the displacement of the 1th site. The Ag? are subject to various symmetry restrictions that we do not need to consider here. He will switch to tight binding language later. The dynamical matrix for this system is defined from the equation of motion. 2 a 0.8 3 m u = {A u (2) 1 JB 133 or 111.123 = 1:32 (23) where u is considered to be a db! dimensional vector. The dimensionality of the system is d and there are N atoms. The dN by dN matrix 1:) is the usual dynamical matrix (see e.g. Born and Huang 1966). If a fOrce £8, also considered as a dN vector, is externally applied to the system, then Ee:-£=22. (3) He will be considering only cases where 2e (and hence _t_‘_ the force produced by the system to oppose 3e in equilibrium) are non zero at the surface. Inverting equation (3) we find that (‘1) II: N "D 1*» where g is the zero frequency limit of the usual Green function 90.12) given by 148 2) = [mwz- 21-1 . (5) 110 (w As we are only concerned with the static limit, the masses are irrelevant and so we will set them all equal to m. If the system described by (1) consists of a perfect crystal with a dynamical matrix 130 plus ardefect, where the defect has a (localized) potential g, then we have the usual Dyson equation §:P+P1:I(;.. (5) where Q = 90 + y. Here P is the perfect crystal Green function and G is the Green function of the system plus defect. Note that go P = D G = -1. We may solve (6) fOrmally “ . (7) "O u A ...o 1 ""0 n< V "'13 The static energy E' of the system produced by the application of the external force 2e is E' = (8) 0 I: N n vol—n I": no I": f —e Nl-a We see from (6) and (7) that C) 1 "U u ""0 "<1 "0 u ""0 u< A -b I ""0 II< V ""0 The change in the energy A5: = E' - E due to the presence of the defect is therefore 149 AB fol—- I": A NC) I "'0 v I": In mo < :(1-gg)g§. (9) Nil—- This form is quite convenient as usually 2 is only nonzero at the surface of the sample. However a better form for our purposes is obtained by writing the equation similar to (A) but for the perfect system L10=E£ (10) where So are the displacements produced by the force g before the defect was introduced and using this we can write AE:-%u°¥(1-E=’\=I) 20. (9a) This result can be rewritten using the 1 matrix 1' = z“ - 22>" (m as A8 = -‘l u T u . (12) 2 -o = -O This is a very convenient form for developing effective medium theories of many defects as well as the single defect problem we are considering here. Equation (12) is simple because the range of 1' is the same as y and so is localized around the defect. In order to calculate AF, it is therefore only 150 necessary to know go around the defect. Remember that the_u) are the displacements produced by the external force fie in the absence of the defect. When the defect is introduced, the force 3 is held constant and the strain changes. The situation is particularly simple in a Bravais lattice if 2e is a uniform external stress (e.g. shear, compression). Then the lattice distortion is homogeneous and the l-lo is very easy to write down. This is the only case we shall be considering here. 3. Isotropic Born Model We apply the work of the preceding section to the isotropic Born model (Born 191”), where the perfect lattice is described by the potential (13) and the sum is restricted to nearest neighbors. This model decouples in the d orthogonal Cartesian directions into d separate tight binding Hamiltonians H = a X (aTa.- afa.) (14) 1’1 1 1 J 1 where the sums in (13) and (14) go over J which are nearest neighbors of’iq which is also summed over. He will work in the language of the isotropic Born model and "translate" into tight binding language later. This is convenient as stresses and homogeneous strains are easier to visualize than the equivalent in the tight binding model. th Let us consider a uniform strain (go) in which the i atom has a displacement. 2': ER. (15) where e is the magnitude of the strain. This homogeneous strain (15) is appropriate for any Bravais lattice. Because a_l_l_ homogeneous strains (i.e. shear, compression, etc.) for the potential (13) are trivially related, we use the simple form (15) where c is a scalar. The methods outlined here become a little more complex for lattices with a basis where the strain is not usually homogeneous. Inserting (15) into (13), we find that the strain energy E is given by E = %a£2Nza2 (16) where the a is the separation of the z nearest neighbors. Note that the origin used for the external strain (15) is irrelevant as only differences in displacements contribute toward the strain energy (16). A. Bond Defect He now consider a very simple defect. A single m is changed from Q ~ ab. He will number the atoms 1,2 as in figure 1. In general the strain (15) will produce displacements 21 = - _u_2 = £12312 where 312 = 31-- 32 is a nearest neighbor vector. Because the defect is localized, no other displacements need be used as basis functions in (9a). It is convenient to choose the center of the bond 1,2 as the origin for the strain. This excludes any uniform motion of the defect bond, which may be included, but does not contribute to (9a). In this basis (Feng, Garboczi and Thorpe 1985), the y matrix may be written 152 11: (a0- a) [ ) (17) - -1 1 .1 and u = ea [ 2 l (18) ‘° -1 '5 u = E-a- [3) w — /2 (19) J. where I3) = [ l2 ] -1 /5 so that the state vector Is) is properly normalized. Note that g |8> = 2(ao- a) |s> (20) so that 2(a - 0) AB = - $22.12 ° (21) 1 - 2(ao- c) where we have also used the fact that Is) is an eigenvector of 1;. This happens because the two states on the atoms 1,2 transform as a symmetric and an antisymmetric vector under the symmetry group of the bond. Because there is a single representation of each kind, there is no mixing so that the antisymmetric vector ls) must be eigenstate of g. The symmetric representation corresponds to the uniform translation and so is irrelevant. 153 This kind of argument becomes much more powerful in the site defect case to be considered later. Therefore 1 - -2-[P11+ 922- P12 - 921] (22) = P11- P12 by symmetry, 1 eprk'fiU) where P ‘- (23) 13 n k 2 2 m(m - wk) and mwz = az(1 - v ) (2n) 5. 1‘. -1 112-i where TE - z E e (25) and g are the z nearest neighbor vectors. The excitations (211) are Just those of the pure isotropic Born model. In general the P must be be 1J evaluated numerically, but in the present case we can make a shortcut that by using the equation of motion fer P11 2 mm P11 - 1 + za(P11- P12) (26) 2 -1 so that as w + 0, (P11- P12) : -(za) and hence (a -a) AB = - 92a2 0 (27) NI-b 1 + 2(ao- a)/za ] 154 Remembering that when the defect is introduced, the stress is held constant, it.is the quantity 1/a that is renormalized as the energy ~ fZ/c. Hence we may write (28) a E' a - E — 1 + - 1 +- e - : using equations (9) and (11). This can be rewritten in terms of the 10 produced by the f, . 2. I u. (1— : 1 - U D U (29) e -o = -o This result is quite general. We find that in the present case 1 2 2 (1°- 0 a E-Noe a 1 + 2(ao- a)/zd F = 1- 4 (30) e %’Nezza2a c(ab- a) ' c + 2(ao- a)/z where No: é-ch is the (small) number of noninteracting defects. It is more usual to invert (30) to give c(ao- a) a =a+ (3021) e 1 - 2(00- a) 155 c(oo- a) : a + a + 2(ao- c)/z to first order in the concentration of defects c. This has been obtained previously by Kirkpatrick (1973) but serves to illustrate the method. A particularly interesting case occurs for a missing bond (00: 0) when 06 = a - ca/(1 - 2/z) (31) which extrapolates to zero at or pl: 1 - cI = (32) This is the intercept in a plot of a against p when the initial slope is extrapolated to where it crosses the horizontal axis. This is discussed in more detail in the next section. Going back to (21), an alternative expression for pI is given by pI = - Za (32a) where Is) is the state vector characterizing the local strain around the defect, but in the absence of the defect, and g is the perfect crystal Green function. This form is useful as a similar but different form holds for a site vacancy. The place where ce actually goes to zero is of course the percolation concentration pc. Kirkpatrick (1973) has shown that the 156 conductance of a network of wires behaves as ae if we identify the tight binding parameter (or in our case the spring constant) a with the conductance of a wire joining adjacent sites. Kirkpatrick has also shown that the diffusion constant D can be obtained from G via G : PSD (33) where P3 is the probability that a site is connected into the infinite cluster. Throughout this paper we will use units such that G and D are normalized to unity at p = 1. This saves having to write 6/60 and D/Do everywhere. It is easy to calculate Ps for small c in the present case. To isolate a site, it must have all 2 bonds cut, therefore P8 = 1 - c2 + 0(022—1] (3n) and so makes no contribution to 0(c). Therefore the initial slope of the diffusion constant D is the same as that of the conductance G and the extrapolated intercept pi is just given by pi = pI (35) In table 1, we show pc,pI and pi for some common Bravais lattices. The quantities G, D, P8 all go to zero at pc. If G ~ (p - pc)t, D ~ (p - pc)t. and P3 ~ (p - pc)B, then t' = t - B. In 2D, t =1.3 and B = 0.111 whereas in 30, t = 1.9 and B = 0.11 from Stauffer (1985). These indices depend only upon dimensionality and are the same for both site and bond percolation. 157 Finally in this section we note that the relation G = PsD is the analogue of E = pv (36) where E is an elastic modulus, p is the mass density and v3 is the sound velocity. It is clear that such a mapping exists because of the equivalence of the tight binding model and the isotropic Born model. u. Site Defect Using the general ideas of the previous section, we can obtain an expression for a general site defect. It is convenient to use a uniaxial, rather than a hydrostatic strain as in the previous section, so that the 1th atom has a displacement a, = eifii- m (37) A A where x, y are any directions; these need not be orthogonal and they can be the same. Also y may be perpendicular to the lattice(e.g. the z direction for a 2D network in the x,y plane). In that case the isotropic Born potential (13) must also contain terms in the z direction. All this arbitrariness is because of the isotropy which leads to a single diffusion constant or elastic stiffness. This formalism can also be used with anisotropic interactions, like central forces, in which case careful track must be kept of directions (Thorpe and Garboczi 1986). The strain energy of the system, using (13) and (37) is 158 E = 11- aezNzaZ/d (38) The vector go is formed from the homogeneous strain (37). All the 2 bonds around the defect site have force constants do, while all the other bonds have the original force constants a. It is convenient to choose the origin at the center of the defect (see figure 2). As in the bond case, this does not affect the result but simplifies the algebra. The nearest neighbors of the central site have displacements 26 = e(§_-x)y (39) where g goes over the z nearest neighbor vectors each of length a. The normalization is given by X u: = :2) (fi-bz = czaZ/d (110) 8 5 if all d directions are equivalent. This is true for all the cases to be considered here and listed in table 2. Thus we are restricted to Bravais lattices in which a second rank tensor is a constant times the unit matrix. The formalism becomes a little more complex in other situations but can still be used. He may write the state vector Is) as _ Z 20' ea/d ls> (‘11) 159 which is a little different from (19). The vector ls) varies depending on the defect and the direction of the strain. For the strain in the triangular net, shown in figure 2. NI-a NI-e _ 1 (‘12) uH-fl '3) :/ nn_. Nl—e Eortunately we will not need to write Is) down in detail. A little thought / will convince the reader that for any 13> derived from (39), Li 13) = (00- a) Is) (113) which differs from (20) by a factor 2, because each atom in the shell connects only to the stationary center in this case. Hence — g_ 2 2 do ' a ] [1 - (ao- c) (nu) where we again use the fact that Is) is an eigenstate of P. The argument in this case is a little more subtle. The state vector Is) transforms as a 160 vector on the shell of the defect. As P is a scalar, and there is only single vector representation on the shell, the result follows. Combining (38) and (44), we see that 2 2c(ao- a) E+AE=71-£Nza2/d [o- 1 - (ao- a) (45) The factor 2 in the numerator occurs because of the counting, and comes about because a bond is removed if either of the two sites at its end are removed. From (28) we can define ae as '2 = 1 - 2c(ao- a) de 1 - (co-a) Inverting to first order in c, we have 2c(aO - a) “e = a * 1 - (ao - a) <31 2 13> (46) (46a) Notice how the factors 2 occur differently in (30a) and (‘16). The problem is now reduced to evaluating . The components of Is) are just proportional to 5110 x as seen for example in (112). If normalized properly to 1, they become / 1%(Bi-x)/a. Therefore 161 . 2 . x 150g . . 'ifi'é (17 ) -d (.8.1 x)e 1 (£1 ’Ue J 5 IKM In 2 Nza k mmk (47) where 3 ) (118) a 3 ‘ [S(kxa1131kya)°°°° and it is convenient to include the nearest neighbor distance a in the definition of the V operator. For the special case of a missing site, co = O and ae extrapolates to zero at p1 = % - g (81:13) (”9) (vvk12 1 1 - :§+§fi£ 1‘Yk ("93) These integrals for pI have to be computed numerically in most cases. They are listed for a number of lattices in table 2. The square net has 1 Yk - 2 [coskxa + coskya] 2 - cos2k a - cos2k a _ l 1 x Y and 80 D1' 2 * ‘BN E 2 - coskxa - coskya u 1 [2 - cost- cosZy]dxdy 1 1 ='2 +'8 J 1 2 - cosx - cosy o 162 = 1 - 1/u as previously found by Watson and Leath (19711). The result for the simple cubic lattice was first obtained by Izyumov (1966) in the slightly more complex situation involving spins. His result agrees with ours. Equation (33) still holds for the site vacancy of course where now P3 = 1 - c + 0(c2) (50) as removing a single site, reduces the number of sites in the infinite cluster. Hence G has a different initial slope from D. Using (33) 3,11 ° 1 11- 1-pi [1-c111- and hence pi 2 — 1/pI (51) These values are also listed in table 2. The results of computer simulations for G and P3 (and hence D) for site percolation on the triangular net are shown in figure 3 (Tang, 1986). The dashed lines corresponding to the intercepts pI and pi from table 2 are shown. The diffusion constant D was obtained in figure 3 using (33). 5. Comments We have shown how a number of known results can be put together. Although we have concentrated on bond and site vacancies, knowledge of pI allows (18 to be found for any defect that has strength (10 in a host a. For a bond, using (30a) and (32a) 163 ca(ao~ a) ae : a + c(1 - p1) + copI (52) while using (46a) and (M9) for a site defect 2ca(ao- 0) ae = a + 2a(1 - p1) + ao(2pI— 1) (53) where in both cases pI can be looked up in tables 1 or 2. For example these can be used in the "superconductingylimit" (see Straley 1983) when 00* a, to give ae = a + ca/pI (52a) in the bond case and oe = a + 2ca/(2pI- 1) (53a) in the site case. It is more usual to write these results in terms of the inverses, which to first order in c are: 1 1 ‘5 = ‘E (1 - c/PI) (52b) e 1 1 1 - 20 and - = - [ -—-] (53b) (16 a 2pI- 1 From (52b) and (53b), we can determine a new quantity p? for the superconducting limit where the initial slope crosses the abscissa. This is 164 clearly simply related to pI for the dilute resistor network on the same lattice. It is more usual to introduce of a fraction (P) of superconducting links so that both problems have a common pc as sketched in figure 11. This can be imagined as follows. Take a resistor network with a fraction p of 2 with 82)) R1. The dilute resistor limit is reached as R2. 0, whilst the superconducting limit is reached when R1’ 0. A similar situation holds for site percolation. Therefore using (52b) and (53b) and remembering to put p ~ 1-p for the bonds having resistence R1 and a fraction 1-p having R superconducting case we find that for bo_nd_s_ p? = pI (51) whilst for _sit_es pi = 91'3 . (55) Again like equations (35) and (51), relating pi to p1, equations (54) and (55) do not depend upon the lattice or the dimensionality. Note that p: for sites is given from (49a) and (55) as s 1 He note that i_f the defects are sufficiently far apart, it is 1/ae (or the resistance) that is linear in c the concentration of defects and contains 92 c2 terms. The c2 terms therefore would represent interactions between defects. On the other hand a‘3 (or the conductivity) is also linear 165 c but contains terms in 02 even if the defects are sufficiently far apart that they do not interact. These come purely from the expansion of the denominator i . e . [1 + c/(1 - pI)]‘1 = 1 - c/(1 - p1) + 0(02). (55b) This is rather obvious but we have not been aware of it before. Finally we comment on the spin wave stiffness. Consider a Heisenberg ferromagnet with spins S coupled by nearest neighbor exchange J and introduce a few isolated site impurities S' that are be coupled to their neighbors with exchange J'. Izyumov (1966) calculated the spin wave stiffness Ds by looking at long wavelengths when the excitation energies Ek are given by Bk: Dskz. This is a dynamic calculation which is technically (miite a bit harder than the static strain calculations done in this paper. Izyumov's result for Ds (again in units where D3: 1, when p=1) is - - - - __22__ D3 - 1 c[o 1 1 + Ap] (56) where 0 = S'/S (57) - £L§L p ‘ JS ' 1 and A is a number that is given as a lattice integral over the cubic Brillouin zone. Using (33) we can identify c.-.1+—2—°P—— (58) 166 and P8 = 1 + c(o - 1) (59) The quantity (58) reduces to our G or de in equation (‘16) if we put 8 = S' (to get to the usual tight binding model) and identify p = J'/J-1 = ao/o - 1 and also use in (H9) (60) He have checked equation (60) for the simple cubic lattice using our expression for pI and the expression given by Izyumov for A and shown that the two lattice integrals are identically the same. He therefbre conclude that the spin wave stiffness (56) of Izyumov can be used fOr agy_of the lattices in table 2, if A = 2pI- 1 is used. This may surprise the reader familiar with the intricacies of Izyumov's paper which seem to be very dependent upon the simple cubic lattice. The work here shows that is not so and the result (56) only depends on the lattice through A. The inertia term (59) reduces to the usual Ps given in (119) if S' = 0, when the central spin plays no role in the dynamics. For a 9133 defect m’ in a host with masses m ml Ps-1+c(m-1] (61) which is proportional to the total mass of the system. As long as the impurity mass m is coupled into the lattice (61) is the correct expression. As soon as do: 0, however m' is irrelevant in the propagation of sound waves and therefore should be set equal to zero, hence recovering (50). This is 167 also the case in the spin wave result (56) where 8' should be set equal to zero (i.e. o = S'/S = 0) if J' = O, as the central spin is decoupled from the rest of the lattice and so should not contribute to the inertia term (59) for the propagation of spin waves. If we return to the original isotropic Born model, we may combine (53) and (61) using (36) to give the sound velocity vs for a small concentration of site defects as from 2 - £22[ 2co(ao- a) ] s ' 2d a * 2a(1 - p1) + 00(2p1 - 1) mv / [1 + c(m'/m - 1)] (62) where the pI are given in table 2. 6. Acknowledgments He would like to thank A. R. Day for useml ooments and the National Science Foundation for support under grant DMR 8317610. 168 Table Captions Table 1. Showing the percolation concentration pc, the intercept pI of the extrapolated initial slope fer the conductance G; the intercept pi of the extrapolated initial slope for the diffusion constant D and the intercept p: of the extrapolated initial slope in the superconducting limit. It will be helpful to glance at figures 3 and u to see what these are. All values are for gang percolation. The values of po are from Zallen (1983). Table 2. Same as table 1 except for site percolation. 169 Table 1 Lattice pc pI pi p? linear chain 1 1 1 1 square net .5 .5 .5 .5 triangular net .347 .333 .333 .333 s.c. .247 .333 .333 .333 b.c.c. .179 .25 .25 .25 f.c.c. .119 .167 .167 .167 d = u hypercube .160 .25 .25 .25 d = 0 hypercube O 0 O O 170 Table 2 Lattice pc pI pi p? linear chain 1 1 1 .5 square net .593 .682 .534 .182 triangular net .5 .6A1 .490 .141 s.c. .311 .605 .347 .105 b.c.c. .2fl5 .579 .273 .079 f.c.c. .198 .561 .217 .061 d = A hypercube .197 .573 .255 .073 d = a hypercube 0 .5 0 0 Figure 1. Figure 2. Figure 3. Figure 4. 171 Figure Captions Showing the two displacements u and _1_1_2 produced by a uniform 1 external stress around a single altered bond. These comprise the components of the vector ls>. Showing the displacements produced by a uniform uniaxial external stress in a triangular lattice around a single altered site. These comprise the components of the vector ls>. The symbols show the quantities D,Psand G for site percolation on the triangular net obtained numerically using methods similar to those described in Kirkpatrick (1973). Each sample contained 1600 atoms and the results were averaged over 20 samples. The quantities pI and pi where the initial slopes cross the abscissa are indicated. A sketch of the conductance G versus p showing the intercept pI of the initial slope. This network has a fraction p of missing units (sites or bonds). Also shown is a sketch of the resistance R versus p superconducting limit showing the intercept of the initial slope p?. This network has a fraction 1-p of superconducting bonds. Both G and R are normalized to the same quantity in the respective pure system. 172 1 .115 Figugg A1. Snowing the two displacements g1 and g2 produced by a mifonn exterml stressamniasimlealteredbad. 'Dmecmprisethecmponentsofthe vector-ls>. FiguLe A2. Showing dedisplacamtspzmicedbyamifommiaxialextemalstrms in a triargllar lattice around a single altered site. These conprise the cmpmmts of the vector ls>. 173 o. _ _ 4.3865 8a 838%.. «5 898 88?. #225 £3 98:: .a Em a 81:35.6 m5. .838 ow .88 83.26 953 3:8. «5 Ba 23m 89 8:580 39% 8mm .33: 30:33:: 5 8588.. 89: 3 33m 89:8 was 5835. 3:38 use .3853 m5 6 8:283 BE 8.. o Be .2. 833:3 65 Si 23% 9: .3 a "a ma 3 me _me ..ve who . _ _ \ 13%. \ E. +B N. v. m. m. o. o o o o a 174 imam? 8:. 93896.. «5 5 $3:qu Bum 9: 8 8335.8 we a 1a a 58 .33 9338883 ..o a; 8:8,... a 3. €032 3.; .aa 83» 325 m5 .3 38.35 «5 mini :5: QBoacoonm a was? m 8.528.. 05 so 3891 m 3 3m 82 Ag .6 8:3 3:3 233... ..o a 838.: a as .538 35. .32» 335 us... no a 3%.: m5 wizofi a 39.3 o mo§o§8 on... ..o :39.» < Q . _. HQ on— Mn— 0 175 References Born H 1914, Ann. Phys. Lpz fig, 605 Born H and Huang K 1966 minamical Theory of Crystal Lattices (Clarendon, Oxford) Elliott R J, Krumhansl J A, Leath P L 197" Rev. Mod. Phys. 56, 465 Economu E N 1979 Green's Functions in Quantum Physics (Springer-Verlag, Berlin) Feng S, Thorpe H F, and Garboczi E J 1985 Phys. Rev. B 31, 276 Izyumov, Yu 1966 Proc. Phys. Soc 81, 505 Kirkpatrick S 1973 Rev. Mod. Phys. QQJ 57H Stauffer D 1985 Introduction to Percolation Theory (Taylor and Francis, London) and references contained therein Straley J P 1983 in Percolation Structures and Processes. Eds. Deutscher G,’ Zallen R, and Adler J Ann. Israel Phys. Soc. V5, ch. 15 Tang H 1986 unpublished Thorpe M F and Garboczi E J 1986 to be published Hatson B P and Leath P L 197“ Phys. Rev. B 2, “893 Zallen R 1983 The Physics of Amorphous Solids (Riley & Sons, New York) 176 Appendix The general result for the effective ae, due to the presence of a defect is A“ A“ (29) cl“ 3“ where no is the strain in the system before the defect is added. Although (29) was derived with a local defect in mind, it is quite generally true for any change 1:] in the dynamical matrix. He give two simple illustrations of the use of (29). In the first we suppose that gl_l_ the bonds in a system described by (13) have their force constants changed from a to do. This is not a localized defect, but (29) is general fOr any 2. In this case ab- a 2 (MB. “‘1’ no. a -1 and hence I = Y [1 - E ( a l 20] a - a :¥[1+ 0a 1-1 a =?\:] (A2) 0 where we have used = -1. Therefore we have D 177 a u l- y) u a '0 C10- “‘0 a Go. a a ‘E- : 1 - u D u = 1 -'E- [ a 1 = '5‘ (83) e "0 =0 0 O O which gives the expected result 0e: do. In the second example, we consider a 10 chain described by (13) in which a single bond is changed from a to co. Equation (30) becomes (with c = 1/N and z = 2) ' N 1 + (a0 - a)/a 2_=, (1 ° ' (M) which can be rewritten as (A5) + 9': .1. a 0 which is the expected result for (N - 1) springs a and one spring do in series. Appendix B. Absorbing rate in continuum media 178 Appendix B. Absorbing rate in continuum media In this Appendix. we study the 29§2£E£fl&_22222é£l. problem in continuum media. We start from the general diffusion equation: 329(F,t) (31) ac2 D v29(3,c) = where D is the diffusion constant. Assuming p(r,t) can be seperated into the form: P(r,t) = R(r) T(t) (32) Then we have the differential equation for R and T respectively: fl = , 1.21) T (133) t and (V2 + k2) R : O (3“) where (83) has the solution: T(t)= e"k D‘ (35) 179 The complete set of eigenfunction of eq.(B1) has the follow- ing form: P(r,t) = X Ak Rk(r) Tk(t) (B6) k where the Tk(r) and Rk(t) are the solutions of (B3) and (BA), and k is indicating the discrete eigenvalue for the corresponding boundary condition. If we have the initial probability distribution given by P(r,t=0), by substituting into eq.(B6), we get: P(r,0) = { AkRk(r) (B7) k and assuming the orthoganality of the eigenfunction: I Rk(r) Rk,(r) g; = 3k skk, (38) Then we have: I P(r,O) Rk(r) 95 = AkBk (39) Substituting (B9) into (B6), we write P(r,t) as: I P(r',0) R (r') dr' P(r,t) = X “R B k " Rk(r) Tk(t) (B10) k k We define the following quantities: 180 I,(t) --- Total number of the particles in the media at any time t: I,(t) = Iv P(r,t) d_r; (311) where V0 is the volume of the media. I,(t) --- Number of particles leaving at time t: 3 I t I,(t) - - a t dt _ 3 P r t) ' Ivo 3 r 9: dt = - I V2 P(r,t) dr (B12) vo -— - I [V P(r,t)] ° d; where s is the surface of v.. --- Average time for a particla reach the boundary: I:t1,(t) (it m (t) = a = I, I,(t) dt (B13) 1.1,(t)dt and substituting (B7) and (B8) into (B12) and (B13) we got: 181 IV P(r,0) Rk(r) pp :2; ° 2 2 I RK(r') d_r (311)) D k Iv. Rk (r) pp v, Actually, the recipracal of is the deposit rate of par- ticle on the absorbing wall. For any given initial distribution P(r,0), boundary condition as well as geometri- cal structure of the system which determines the set of eigenfunction Rk(r), eq. (81“) gives the deposit rate or ab- sorbing rate. We summarize several cases corresponding to the uniform initial distribution as the following: case 1). Initial distribution is a 5-function at center: P(r,0) = 5(_:;) (315) then Ivn Rk(r) pp 2 2 k k IvoRk(r)gp Rk(0) (B16) case 2). Uniform initial distribution in the whole region. P(r,0) 1/vo if r=£ (—-) k k2!v R: (r) g; Actually, case (2) is included in case (3) for v'=v.. All of the formula derived above are suitable to any dimension. Let us consider a circular piece in two dimension , i.e., rSr.. The Bessel function gives the complete set of the eigenfunction: Rk(r) : J,(kr) (321) where J,(x) is the zero-th order of the Spherical Bessel Function and eigenvalues are given by: k=x£0)/r, (322) where J,(x£0)) = o 183 He found out that the diffusion constant is given by: and 2 _LA 1 D'2 case (3) is more complicated. for case (1) for case (2) (323) (824) Appendix C. Formulation of Effective Medium on Triangle Net 184 Appendix C. Formulation of Effective Medium on Triangle Net Q The value of a can be calculated exactly through the formula Thorpe & Garboczi 1985): I N K -1 a = m z Tr { Z [ 1 - exp(iLk°8) ] (88) 0 2 ( E ) } k 8 N 2 (C1) where Q (12) is the Fourier transform of the dynamical matrix: Q (E) = Km n 2 i 1 - exp(1LEo3) 1 33 a + Km (I-n) X I 1 - exp(1L§-3) 1 (c2) 5 where n = LO/L and I is the dxd unit matrix. Through some matrix algebra, we get the result: ( n'1 -1 ) S (C3) : N I Nlm a = —— n - where S is the sum in k space: K 3 = fig 1 i X [ 1 - exp(iLE-3) ] Tr [2’ K 8 ‘ (E) 1 1 (cu) For the triangular lattice, we have and where and XX D YY XY 1 N 185 n" g— n“ -1 ) s D + D X X XX yy D D -D D k xx yy xy yx ( 1-n ) ( 6 - 2 cost - 4 cosx cosy ) + n (3 - 2 cos2x - cosx cosy ) ( 1-n ) ( 6 - 2 cost - fl cosx cosy ) + n ( 3 - 3 cosx cosy ) D : yx - 2 cost where x Y J3 n sinx siny A cosx cosy Nl-a X (C5) (C6) (C7a) (C7b) (C70) (C7d) LIST OF REFERENCES 186 LIST OF REFERENCES R. Alben, S. Goldstein, M. F. Thorpe and D. Heaire, Phys. Stat. Solid (b) 5 , 595 (1972). R. Alben, S.Kirkpatrick, D.Beeman, Phys. Rev. B 15, 346 (1977). Broadbent, S.R., and J.M.Hammersley, Proc. Camb. Phillos. Soc. 53, 629 (1957). T.H.K.Barron and M.L.Klein, Proc. Phys. Soc. 85, 523 (1965). R.J.Elliott, J.A.Krumhansl and P.L.Leath, Rev. Mod. Phys. _lI_,1165(19711). J. N. Essam in Phase Transitions and Critical Phenomena, Ed. by C. Bomb and M. S. Green (Academic Press, New York, 1972) Vol. 2, p.197. S. Feng, H. F. Thorpe, E. Garboczi, Phys. Rev. B 31, 276 (1985). S. Feng and P. N. Sen, Phys. Rev. Lett. 22, 216 (1989). P. G. de Gennes, J. Phys. (Paris) 31, L-1 (1976). N. A. Harrison, Electronic Structure and the Properties of' Solids (H. B. Freeman and Co., San Francisco, 1980). H.He, Ph.D. Thesis, Michigan State University (1985). H. He and N. F. Thorpe, Phys. Rev. Lett. 23, 2107 (1985). C.L.Henley, submitted to Comm. in Cond. Mat. Phys. (1986). K.Huang, Proc. Roy. Soc. A, 223, 178 (1950). J. B. Keller, J. Math. Phys. 5, 598 (1969). A. Khurana, Phys. Today, 3g, No. A, 23 (1987). S. Kirkpatrick, Rev. Mod. Phys. £5, 57“ (1973). P. H. Kognt and J. P. Straley, J. Phys. C 13, 909 (1981). D. Levine and P.J.Steinhardt, Phys. Rev. B 39, 596 (1986). A.E.H. Love, A_ngap1§g_pp_ppg_flathematical Theory of Elasticity, Dover, New York (1999). 187 K. S. Mendelson, J. Appl. Phys. BB, 917 (1975); 3B, 9790 (1975). D.Shechtman, J.Blech, D.Gratias and J.H.Cahn, Phys. Rev. Lett. BB, 1951 (1989). J.E.S.Socolar and P.J.Steinhardt, Phys. Rev. B 39, 617 (1986). D. Stauffer, Introduction to Percolation Theogy (Taylor and Francis, London, 1985). J. P. Straley in Introduction to Percolation Processes, Eds. G. Deutscher, R. Zallen and J. Adler, Isreal Phys. Soc. V5, Ch.15 (1983). J. P. Straley, J. Phys. C B, 783 (1978). J. P. Straley, Phys. Rev. B 15, 5733 (1977). M.F.Thorpe, J. Non-Crys. Sol. _1, 355 (1983). B. P. Waston and P. L. Leath, Phys. Rev. B 9, 9893 (1979). R. Zallen, The Physics of Amorphous Solids (Wiley and Sons, New York, 1983).