COMPUTER SIMULATIONS OF CRITICAL PHENOMENA INSYSTEMS WITH LONG RANGE INTERACTIONS: A STUDY OF ISINGDIPOLES AND SELF-ORGANIZED CRITICALITY IN EARTHQUAKESByHuang-Jian XuB. Sc. (Space and Earth Science) University of Science and technology of ChinaM. Sc. (Physics) University of British ColumbiaA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESPHYSICSWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIASeptember 1992© Huang-Jian Xu, 1992In presenting this thesis in partial fulfilment of the requirements for an advanced degree atthe University of British Columbia, I agree that the Library shall make it freely availablefor reference and study. I further agree that permission for extensive copying of thisthesis for scholarly purposes may be granted by the head of my department or by hisor her representatives. It is understood that copying or publication of this thesis forfinancial gain shall not be allowed without my written permission.PhysicsThe University of British Columbia2075 Wesbrook PlaceVancouver, CanadaV6T 1W5Date:TAbstractThis thesis discusses scaling and critical behavior of two different models. One modeldescribes Ising dipoles, originates in condensed matter physics and depicts equilibriumcritical phenomena. The other model, taken from the earth sciences, describes faultinginstabilities and the resulting earthqnakes, and involves self-organized criticality— a non-equilibrium phenomenon. Both models are characterized by long range interactions, witha resulting sensitivity to boundary conditions.The ordering properties of Ising dipoles on lattices are studied in a mean field theoryand by Monte Carlo simulations. The mean field theory is manifestly shape independentin zero external field. In the case of dipoles on a diluted lattice the mean field theorypredicts a critical concentration above which the low temperature phase is ferroelectric (oranti-ferroelectric depending on the lattice structure). Extensive Monte Carlo simulationresults are in agreement with those of mean field theory.We propose a finite size scaling form that includes logarithmic corrections for systemsat the critical dimensionality. In the case of dipoles on a body centered tetragonallattice we found that the finite scaling form significantly improved the data collapse overthe scaling form with mean field exponents. With lattice parameters appropriate to theIsing ferromagnetic compound LiHoF4,we obtain a ferromagnetic transition temperature:ii= 1.51K in excellent agreement with experiment. This indicates that the materialLiHoF4 is dominated by the dipole-dipole interaction; since in the simulations we onlyinclude dipole—dipole interactions.For dipoles on the simple cubic lattice, the ordered state is made up of anti-ferromagneticrows. The critical exponents obtained by finite size scaling are 1/7, 8/7 anda 4/7. These results are in good agreement with those of high temperature seriesexpansions.A model of self—organized ruptures in an elastic medium is developed; and appliedto earthquakes. In the model the local ruptures are represented by double couples to beconsistent with elastic theory. The explicit form of this double couple source is derived.The system is driven by slowly increasing the shear stress. The model evolves towards aself-organized critical state in which the earthquake distribution follows the GutenbergRichter law with an exponent in agreement with observational data. By modeling thelocal static fatigue for the rocks, we also obtained Omori’s law for the rate of aftershocks.The effects of annealing are investigated.iuTable of ContentsAbstract jjList of Figures viAcknowledgements ix1 Introduction 11.1 Ising dipoles 21.1.1 Brief historical review and theoretical background 21.1.2 Experimental realization of lattice dipolar system 91.2 Earthquake Model 102 Dipoles on perfect lattices 152.1 A brief review of the finite size scaling theory 162.2 Ferromagnetic phase transition: Monte Carlo simulation and finite sizescaling with logarithmic corrections 182.2.1 Description of the Monte Carlo simulations 192.2.2 Data analysis: Finite size scaling with logarithmic corrections . 27iv2.3 Critical behavior of Ising dipoles on simple cubic lattice: antiferromagneticphase transition 323 Dipoles on diluted lattices 513.1 Model and mean field theory 533.2 Zero temperature Monte Carlo simulation 623.3 Finite temperature Monte Carlo simulation 653.4 Discussion and summary 684 Self—organized ruptures in an elastic medium 774.1 Self—organized criticality 774.2 Earthquake occurrence phenomenology 794.3 Discretization of the elastic equations 824.3.1 A brief review of previous models for discretizing elastic equations 824.3.2 A Plaquet Representation of Ruptures in an Elastic Medium 854.4 Applications to Earthquakes 944.4.1 The effect of weakening and annealing 964.4.2 Simulation of aftershocks 995 Conclusions 113Bibliography 117vList of Figures1.1 Illustration of boundary conditions 131.2 An example of off—center dipolar system 142.1 The structure of LiHoF4 372.2 Specific heat for body centered tetragonal lattice 382.3 Temperature dependence of < M2 >1/2 f/V for body centered tetragonallattice 392.4 Data collapse for the body centered tetragonal lattice—logarithmic corrections 402.5 Data collapse for the BCC lattice—logarithmic corrections 412.6 Temperature dependence of the “Binder number” for body centered tetragonal lattices of different sizes 422.7 Temperature dependence of the “Binder number” for BCC lattices of different sizes 432.8 Data collapse for the body centered tetragonal lattice—mean field 442.9 Data collapse for the BCC lattice—mean field 452.10 Data collapse for body centered tetragonal lattice 462.11 Specific heat of 12 spins on a simple cubic lattice 47vi2.12 Temperature dependence of the magnetization < M2 >1/2 /N for thesimple cubic lattice 482.13 Data collapse for the simple cubic lattice 492.14 Specific heat for simple cubic lattices of different sizes 503.1 Phase diagram of the SK spin glass model (adaped from Binder and Young[12]. ) 693.2 Ground state polarizations for the diluted BCC lattice at various concentrations 703.3 Ground state local field distribution at c=0.1 713.4 Log—log plot of the ground state field distribution for small E at c=0.1 723.5 Data collapse for diluted BCC lattice at concentration c=0.7 733.6 Data collapse for diluted BCC lattice at concentration c=0.5 743.7 Specific heat for BCC lattice at c=1.0 (right), and 0.5 (left) 753.8 Critical temperature for the ferroelectric transition at various concentrations 764.1 The illustration of the earthquake model 1024.2 The shear stress distribution for a fault calculated by lattice Green’s functions 1034.3 Size distribution of earthquakes in the case of instantaneous annealing 1044.4 Size distribution of earthquakes in the case of weaking and annealing 1054.5 Some typical spatial patterns of simulated earthquakes 106vii4.6 The epicenters of 500 quakes in the case of instantaneous annealing . . 1074.7 The accumulated slip during 500 quakes in the case of it instantaneousannealing 1084.8 The epicenters of 500 quakes in the case of weaking and slow annealing 1094.9 The accumulated slip during 500 quakes in the case of weaking and slowannealing 1104.10 Omori’s law 1114.11 Spatial pattern of a main event and after shocks 112viiiAcknowledgementsI would like to thank my supervisor Birger Bergersen for the support and encouragementthat enabled me to finish this thesis. I also wish to thank Zoltán Rácz for numerousdiscussions on the statistical mechanics of dipolar systems.I had the great fortune to meet Kan Chen from whom I learned about self—organizedcriticality. Without his collaboration, the last part of the thesis would not be possible.Thanks also go to my many friends in the Mathematics Department of the Universityof British Columbia. In the long struggle of working on this thesis, I have benefited agreat deal from discussions with Xiaoming Cai, Huaxinog Huang and Qing Ran.I am grateful to Tom Nicol at the UBC Computer Center and Ron Parachoniak inthe Physics Department for valuable assistance. The help of Nathan Weiss in teachingme how to use the graphic program Super—mongo is deeply appreciated.ixChapter 1IntroductionThis thesis discusses the behavior of two seemingly very disparate systems that originatefrom different areas of science, but which have certain common features. Both systemsexhibit critical behavior although in the first system the context is equilibrium phasetransitions, while in the second it is the self—organized criticality of a dynamical system.Both system are characterized by long—range interactions, although for different reasons.Due to the nature of long range interactions, these systems may be sensitive to boundaryconditions. Finally, both systems are studied by computer simulations. In some otherdynamical systems [104], the effective interaction can be produced by the dynamics itself,but we will not address this issue in this thesis.In the first part of the thesis the system consists of electric or magnetic dipoles orspins. We assume that the dipoles cannot rotate freely but can only take on one oftwo possible orientations and are centered at fixed locations. We refer to the spins asIsing dipoles. In this thesis we ignore the effect of the medium. Therefore, terms like“ferroelectric” and “ferromagnetic” mean the same thing except in the former case werefer to the electric dipoles, while the latter we refer to the magnetic dipoles. The1Chapter 1. Introduction 2cooperative phenomena for Ising dipoles in crystals are studied by mean field theory andextensive Monte Carlo simulations are carried out in situations where the dipoles arelocated on crystalline lattices or on disordered lattices.In the second part of the thesis we develop a model for earthquakes in which a steadybuildup of shear stress in an elastic medium causes ruptures that may lead to a chainreaction because of stress redistribution. In this model a self—organized non-equilibriumstate develops with scaling behavior similar to that of dynamical critical phenomena inequilibrium statistical mechanics. This model can reproduce some well known empiricallaws such as the Gutenberg-Richter law for the distribution of earthquakes and Omori’slaw for the rate of aftershocks.1.1 Ising dipoles1.1.1 Brief historical review and theoretical backgroundThe subject of the ordering properties of dipoles has a long history, dating back to theearly studies of Debye [22]. To explain the temperature dependence of the dielectricconstant in certain liquids containing polar molecules, Debye noted that it was necessaryto consider the effect of the molecular electrical dipole moments. He started from theLorentz expression for the average internal local field inside a cavity which is cut outfrom a medium:(Li)Chapter 1. Introduction 3where P is the polarization, F the macroscopic electric field and e is the permittivity ofthe vacuum. From (1.1) Debye proceeded to construct a simple mean field theory whichpredicted a ferroelectric ordering at low temperatures. For the sake of simplicity, letus consider the case where dipoles only have two discrete allowed orientations: paralleland antiparallel to the electric field. We can then consider the dipoles as Ising spins(S = ±1) rather than the free rotors studied by Debye. The generalization to free rotorsis straightforward. For a dipole in a field given by (1.1) it is easy to derive, from theBoltzmann distribution, the following self—consistent equation for the order parameterm=< Sj > at temperature T. Here <> denotes the thermal average.2 Vmm tanh([pE + (1.2)where 3 is 1/kBT, p is the dipole moment, V is the volume of the sample, and N is thenumber of dipoles. Solving the above equation, one finds that in the case of no externalfield, there is an in U solution at T < T0 where T0 is the Curie temperature given by= 3l/Ek11 (1.3)Thus, the above mean field theory predicts a ferroelectric phase transition at T. However,this theory completely ignores the local field fluctuations due to spatial disorder. It willbecome clear later that local field fluctuations tend to destroy the ferromagnetic orderin the case of dilute dipolar systems. This is why, for example, alkali halide crystalscontaining dipole impurities show no sign of ferroelectricity [47] [29].Chapter 1. Introduction 4In the case of isotropic or nematic polar fluids, ferroelectricity has never been observed,although theoretical predictions have been made [60] [76]. In order to explain why dipolarfluids such as water do not exhibit ferroelectricity Onsager [75] argued that only a partof the internal field (1.1), which he called the cavity field, Eca, gives rise to a torque ona dipole inside the cavity. The remaining part, which he called the reaction field, hasno orientation effect. The basic idea is that the reaction field is induced by the dipoleitself and should be subtracted when calculating the torque on a dipole. In the case ofpermanent dipoles randomly distributed in vacuum, Onsager found, in place of (1.1)EcaE+2 (lÀ)where E is the dielectric constant of the medium. When this constant is determinedself—consistently in a generalized Clausius-Mossotti theory there is always a locally stable solution with zero spontaneous polarization. As a consequence, many people believethat a system of randomly distributed dipoles will not exhibit a ferroelectric ordering.However, Pirenne [81], and later in more detail Zernik [111], have shown that the Onsagerreaction field theory does permit a ferroelectric solution. The possibility of ordering canbe seen by noting that in general the dielectric function is given by e = E0 + c9P/DE.In a spontaneously polarized state near saturation OP/OF will be small, and we return to a situation described qualitatively by Debye theory (except the transition to thespontaneously polarized state is now first order).It is possible that the failure to observe ferroelectricity in polar fluids and nematics isChapter 1. Introduction Slargely due to “molecular association” resulting from the fact that the molecules are freeto move. It may be worth noting, that the only molecular dynamics simulation whichhas found ferroelectric behavior in a fluid, was carried at such high densities that thesystem behaved more like a glass than a liquid [102]. In this thesis we only study dipolesin crystalline and amorphous materials, where the positions of the dipoles are frozen.There are many such materials. For example, LiRF4 (where R is a rare earth element)is a magnetic dipolar system, while K+0H with off center Li+ impurities is an electricdipolar system. Many examples of electrical dipolar system can be found in an extensivereview by Vugmeister and Glinchuk [100].A number of authors have studied the properties of the ground state for dipoles onvarious crystalline lattices. The first such calculation was that of Sauer [90] who directlysummed the contribution of 125 dipoles contained in a small cube surrounding eachdipole. The spin configuration with lowest energy was taken to be the ground state. Forvarious lattice structures he found that the ground state could be ferroelectric or antiferroelectric depending on the sample shape. For body-centered cubic lattice the groundstate is anti-ferroelectric for spheres, while it is ferroelectric for a long needle. Thisanalysis is based on a “near neighbor” method. However, since the dipolar interactionsare of long range, this treatment might be inadequate.Luttinger and Tisza [65] performed comprehensive calculations using the more rigorous “normal coordinate” method. For a number of spin structures, their results for theground states are essentially the same as those of Saner. The above calculations were allChapter 1. Introduction 6carried out for the system without an external field. One of their major results was thatthe ground state energies depend on the sample shape. The origin of this dependencelies in the long-range nature of dipolar forces. One can show [98] [90] that the dipolarsum (assuming that all dipoles point in the same direction)U(i) = —(1/N) E(1 — 3cos2O)/r (1.5)is conditionally convergent, i.e., the sum depends on the shape of the external boundary.In the above formula N is the number of dipoles and is the distance between thedipoles i and j; while is the angle between the vector that connects site i and j andthe spin axis. Note that the average of the angular part of the dipole interaction is zero.However, if this angular dependence were not there ( i.e, the interaction is 1/r3), it isclear that the system would not have a thermodynamic limit ( although such systemsare interesting in their own right [83]). Therefore, the peculiar angular dependence ofdipolar forces plays a major part in determining the bulk properties of dipolar materials.A major theoretical advance was made by R. Grifliths [34], who rigorously proved thatfor dipolar lattice systems without an external field, the free energy is shape independentin the thermodynamic limit (sample volume V —* c’YD while keep the number densityN/V fixed). We shall not repeat the detail of the proof here, but simply discuss theconsequences of this result.First of all, let us look at the case of T=0, when the free energy is the same asthe ground state energy. The work by Sauer [90] and Luttinger and Tisza [65] showedChapter 1. Introduction 7that the ground state spin configuration depends on the sample shape, this conflicts withGriffith’s theorem. Where lies the contradiction? And what is missing in the calculationsof Luttinger and Tisza? It turns out that without an external field, the magnetizationcan not be homogeneous inside a sample — there will be domain formation. The ideaof domains was developed by Kittel [51], but the firm theoretical base, is mainly due toGriffith’s proof of shape independence.We next show that the domain structures are a direct consequence of Griffith’s theorem. Let us choose a sample shape, find the minimum free energy of the sample, andcompare the free energy for various shapes. In the thermodynamic limit, the free energyshould not depend on the sample shape. The immediate conclusion is that the net magnetization has to be zero, M=O. Otherwise the free energy would be shape dependent.The reason is simple — the depolarization field is shape dependent (for example, for aspherical sample with net magnetization M, the demagnetization field is 4KM/3 whilefor a thin slab it is 4KM), and the free energy is directly related to the depolarizationfield (the depolarization field will be discussed in more detail in chapter 3). To demandshape independence , we have to set the depolarization field to zero. However, the factthat total magnetization M=O does not exclude the possibility of a ferromagnetic phasetransition. It is possible that the magnetization inside a sample will not be uiliform —if the sample is divided into domains. The size of the domains is macroscopic. Insidea domain, the magnetization is uniform, and in thermal equilibrium the domains arearranged in such a way that there is no demagnetization field from the sample boundary.IChapter 1. Introduction 8The domain structure may be complicated. For example, to the best of the author’sknowledge, even the domain structure for the ground state of a spherical sample has notbeen solved. This thesis is not concerned with the details of the domain structures. Wejust need to emphasize the fact that one cannot ignore the boundary effects when studying the dipolar system; it is the demagnetization field due to the long range interactionsthat causes the domain structure.Therefore, in studying the existence of a ferroelectric phase transition, without specificboundary conditions, we are not allowed to use the order parameter M = Z .s, sinceit is zero at all temperatures. If we assume the domains are of macroscopic size andhence that the domain wall energy per spin goes to zero in the thermodynamic limit [51],we can employ periodic (Ewald) boundary conditions without a depolarization field.An alternative approach is shown in Figure 1. Figure (1.la) shows a sample with freeboundaries. Assuming the lattice structure is such that it favors the ferroelectric phase,a domain structure will develop in such a way that there is no net depolarization field.Figure (1.lb) shows the same sample but now the surface is wrapped by metal and isgrounded. In this case, the surface charge of the boundary is neutralized and there is nodepolarization field. Therefore in figure (1.lb) , the ferroelectric phase is a single domainphase. As we argued, the free energy will be the same in both cases, if the domain wallenergy can be ignored in the thermodynamic limit. Another way to look at it: if wejust concentrate on a single domain inside the sample, the environment in case (a) andcase (b) is the same for this domain, so is the free energy. However, case (b) is muchChapter 1. Introduction 9simpler to deal with, at least theoretically. Only in this case can the ferromagnetic orderbe characterized by an order parameter M . We will discuss this in more detail vhen wedevelop a mean field theory in chapter 3.1.1.2 Experimental realization of lattice dipolar systemExperimentally, the interest in dipolar systems has concentrated on situations where thedipoles are arranged in a diluted lattice, as impurities. The most typical representativesof electrical dipole impurities exhibiting cooperative effects are off—center ions, substitutional impurity ions, whose equilibrium positions are displaced from regular lattice sites,creating an effective dipole moment. The preferred directions of such dipoles is determined by the crystal symmetry. The reorientation of dipole directions can occr due tothermal jumps and tunneling transitions between the various positions of the impurityions. To give a specific example, let us look at the compoundK1_LiTaO3.The localdipoles are attributed to Li ions which occupy K sites. Since the Li ions are smaller thanthe K ions, this causes a misfit which produces a dipole moment. Figure (1.2) shows thethe (1 0 0) plane occupied by oxygen and potassium ions. In this plane Li has four stablepositions. Together with the two displacements perpendicular to the (1 0 0) plane, thereare six possible orientations of the dipole.An extensive review on electric dipolar systems has been given by Vugmeister andGlinchuk [100]. In a highly polarizable dielectric medium, even the form of the interactionbetween dipoles is not well understood. In this thesis, for simplicity, we will not includeChapter 1. Introduction 10the effect of the medium. Therefore, we will oniy compare our theoretical results withmagnetic dipolar compounds where one can ignore the effect of the medium.The compounds LiRF4 (where R is an rare earth element) are magnetic dipolarsystems at low temperatures. For example, the ground state of a Ho+ ions in the crystalfield of the LiHoF4 lattice is an Ising doublet [39] with g 14. The first excited crystal—field level is 9.4K above the ground—state doublet. At low temperatures around 1K onlythe Ising doublet is appreciably populated. Recently there has been increasing interestin this compound and the compound LiHoY1F4which is a disordered Ising dipolarsystem (the Y ion is not magnetic) [32] [84]. We will discuss these in detail in chapter 2for dipoles on perfect lattices and in chapter 3 for randomly distributed dipoles.1.2 Earthquake ModelTn 1987 Bak, Tang and Wiesenfeld [7] developed an important concept to explain thebehavior of driven extended dynamical systems: the theory of Self—Organized Criticality(SOC). They suggested that many systems naturally evolve to a critical state whereavalanches of all sizes occur, and critical dynamics constitutes a general mechanismfor generating spatial and temporal self—similar structures. Immediately following itsdiscovery, the concept of SOC was taken up by many researchers [9] [19] [21] [95] toapply to the dynamics of earthquakes, where self—similar scaling ranging many decadesis common. There are two famous scaling laws that are well established observationally.Chapter 1. Introduction 11One is the Gutenberg— Richter law for the earthquake moment distribution [36]N(M) (1.6)where N(M) is the number density of earthquakes with moment lvi, and lvi is the moment(or energy release) and b is close to 0.5. Another is Omori’s law for the rate of aftershocks(1.7)where P(t) is the number of aftershocks after time I of a major earthquake, and theexponent p is close to unity. The Gutenberg—Richter law (1.6) can be understood as adirect consequence of the critical dynamics of earth’s crust within various SOC models [9][21]. Our model is based partly on the crack-propagation model of earthquakes proposedby Chen, Bak, and Obukhov [21]. In their model, the instability is caused by a localrupture when the local stress exceeds a threshold. In contrast to the stick-slip models [19],[9] in which earthquakes are generated with blocks sliding along a pre-existing fault, thefault system in this model is generated dynamically by earthquakes. In addition, a long—range redistribution of elastic force is introduced following a local rupture. The b valueof the Gutenberg-Richter law obtained in their model is different from the one obtainedin the stick-slip models, but is in agreement with observations [45]. However, there areshortcomings in their model. First, the tensile elastic force is not correctly represented.Second, the model does not exhibit aftershock sequences. The model we consider hereovercomes such shortcomings with a correct discretization of the continuum elastic theoryand with the introduction of a static fatigue law to the local dynamics. In chapter 4 weChapter 1. Introduction 12construct a novel plaquet representation of ruptures in an elastic medium. This schemehas the advantage of being fast and easy to implement in computer simulations. It isquite accurate in cases where we are able to compare with analytical results. With thismethod various empirical laws in earthquakes can be included in our dynamical model.We show that Omori’s law can be accounted for by including a static fatigue law for therocks in this more realistic crack propagation model. The effects of annealing are alsoinvestigated in the model.Chapter 1. Introduction 13a)AVb)Figure 1.1: Illustration of the boundary conditions. (a) Free boundary, the low temperature phase consists of ferroelectric domains and the total polarization is zero. (h) Thesurface charge is neutralized by grounded conducting plates. The depolarization field iszero and the low temperature phase is uniformly polarized.Chapter 1. Introduction 14Figure 1.2: The structure of the (1 0 0) plane of the perovskite KTaO3. The ionicradii and structure are drawn to scale; the K — K distance is 4A. Of the four equivalentoff—center positions shown (plus two out of plane, not shown), only one is occupied by aLi ion.Chapter 2Dipoles on perfect latticesIn this chapter we study, using the Monte Carlo method, the phase transition properties ofIsing dipoles on a number of perfect lattices, namely the body centered tetragonal, bodycentered and simple cubic lattices. Since size effects will be more profound in models withlong range forces than for short range ones [58], it is non—trivial to extrapolate the finitesize Monte Carlo data to corresponding thermodynamic results. In section 2.1 we brieflyreview the general finite size scaling theory. In section 2.2 we modify the usual finitesize scaling form to take logarithmic corrections into account for dipolar ferroelectricphase transitions. We describe in detail Monte Carlo Simulations on body centered andtetragonal lattices and calculate the transition temperature for LiHoF4 in section 2.3.Also in that section, the data collapse with different scaling forms are discussed. Insection 2.4 we discuss the Monte Carlo simulations for dipoles on simple cubic latticewhere the ground state consists of anti—ferroelectric rows.15Chapter 2. Dipoles on perfect lattices 162.1 A brief review of the finite size scaling theorySince the pioneering work of M. E. Fisher [30], the finite size scaling theory has foundmuch application in numerical simulations. Due to the fact that computer simulationscan only deal with a finite number of particles, the finite size theory has become anextremely powerful tool to interpret simulation data. In this section we briefly reviewthe basics of the theory. For a detailed review see Privman [82) and the articles citedtherein.Let the Monte Carlo cell be a cube with side L. At temperatures close to a secondorder phase transition point T, the basic assumption of finite size scaling theory takesthe following form:(2.1)where the reduced temperature I = (T— T)/T, and is the correlation length. FL(I),P (I) are thermodynamic quantities for a system of size L and for the bulk system,respectively. This assumption is plausible since close to T the only important characteristic length is and finite size effects are controlled by comparing the two lengths of theproblem, the linear dimension L and the correlation length . Brezin et al [16] used afinite field theory approach to prove the validity of (2.1) for dimension d < d*, the uppercritical dimensionality. For d > d*, the critical behavior is Landau like (or the criticalexponents are the same as given by a mean field theory). At d* (2.1) breaks down dueto logarithmic corrections. Generally close to T for a second order phase transition weChapter 2. Dipoles on perfect lattices 17havecx (2.2)andP(t) cx t” (2.3)We then have a finite scaling form that is often used in practicePL(t) txf(Lr11) (2.4)A useful method to estimate the transition temperature T and the critical exponents xand ii is to simulate systems with different sizes L and plot PL(t)/F against Lt. Allthe curves for different L should collapse to a single scaling function f when the correctvalues of T and the critical exponents have been employed. We shall use the above datacollapse technique throughout this thesis.However, the finite size scaling theory is not just L scaled by . In one case, supposeone simulates a system of rectangular shape Li . L2, it is then non—trivial what thescaling form is [11]. Another interesting case is when the interactions are sufficientlylong range that mean field theory applies. Then, “dimensionality” and “length” losetheir meaning. It is found that if d > d*, the system is Landau like, and one has toreplace the correlation length by a “thermodynamic length” A cx t_21d [16] [11]. For3-D Ising dipolar system, it is even more complicated. First of all, 3-D is the criticaldimensionality where logarithmic corrections have to he taken into account. Secondlythe system is highly anisotropic, the correlation length in the spin direction is larger thanChapter 2. Dipoles on perfect latticesin the directions perpendicular to it [2]. We will discuss this system in more detail insection 2.2.2.2.2 Ferromagnetic phase transition: Monte Carlo simulation and finite sizescaling with logarithmic correctionsIt is well established in renormalization group theory that for systems above criticaldimension d*, the critical behavior is mean field like. At d = d*, the renormalization grouppredicts logarithmic corrections to the mean field behavior [2]. Actually, even before theinvention of the renormalization group theory, Larkin and Khmel’nitskii were in 1969 [61]able to derive the logarithmic correction for dipolar ferromagnetic phase transitions witha field theoretic approach. A physical realization of systems at critical dimensionality arethe compounds LiRF4 (where R is a rare earth element) [4] [32]. Experimental resultsfor this system are compatible with the theoretical prediction of logarithmic corrections.For the short range Ising model, d* = 4, and computer simulations have been carried outby Mouritsen and Knak Jensen [71] [72] who found logarithmic corrections to the orderparameter in agreement with theory. A finite size scaling analysis of the same model wasdone by Lai and Mon [59]. The effect of high temperature series in the 4—D Ising modelwas studied by Guttmann [37] and McKenzie [68].We choose to study the compound LiHoF4 with Monte Carlo simulations for tworeasons.(i) The material has been studied experimentally and the transition temperature is knownChapter 2. Dipoles on perfect lattices 19to a reasonable accuracy [32]. We can directly compare our simulation results with theexperiments.(ii) It is generally accepted that this compound is a permanent Ising dipole ferromagnet,while e.g. LiTbF4 has an induced dipole moment [84]. It is computationly more feasibleto simulate the permanent dipolar systems.It is also generally believed that in LiHoF4 the dominant interaction is the dipole—dipolecoupling. Knak Jensen and Kjaer [55] carried out Monte Carlo simulations with only thedipolar interactions for the compound and they found T 1.86K while experimentally1.36K . This led them to suggest the disagreement with the experiment is due toshort range antiferromagnetic couplings [13]. However, in their simulation the boundaryconditions were awkward. They used an ellipsoidal sample shape with free boundary. Aswe pointed out in the introduction the ground state will then have a complicated domainstructure [5] and the net magnetization is hard to interpret. One of the purposes of thepresent investigation is to show that a more accurate simulation with correct boundaryconditions, employing finite size scaling with logarithmic corrections, yields a transitiontemperature that is in good agreement with the experiment for the compound LiHoF4.2.2.1 Description of the Monte Carlo simulationsIn this section we describe the Monte Carlo simulations in detail. The basic Monte Carlocell is chosen to be cubic for the reason of simplicity. As we have emphasized, the boundary condition has to be such that there is no depolarization field. There are several waysChapter 2. Dipoles on perfect lattices 20to accomplish this:(i) wrap the cubic cell with a metal which is grounded.(ii) use periodic boundary conditions.(iii) truncate the interaction but compensate the depolarization field by adding an external field in the simulation which cancels the depolarization field; or use a more sophisticated method such as adding the Onsager reaction field term.In both cases (i) and (ii) we have to include all the image dipoles when calculatingthe effective dipole—dipole interactions in the Monte Carlo cell. The spatial arrangementof the image dipoles is different in the two cases. We found that in the former case theeffective D—D interaction does not preserve the translational symmetry and this is to thedisadvantage of the numerical simulations, since we need more internal memory to keeptrack of all the couplings between pairs of dipoles. With the periodic boundary conditions, standard techniques that preserve the translational symmetry exist for calculatingthe effective D—D interactions. Kretschmer and Binder [58] have compared the variousboundary conditions and they found that the periodic boundary condition method, withEwald summation, is clearly the most reliable. Therefore, in this thesis we always useperiodic boundary conditions. We will present the Ewald summation formulas to be usedin the simulation without proof. For detailed derivations see, for example, Adams andMcDonald [1] or Jansoone [43].The Ewald summation method amounts to splitting the slowly convergent sum overall the interactions in the lattice into two quickly convergent sums, one running over theChapter 2. Dipoles on perfect lattices 21reciprocal lattice and the other over the original lattice points. For simplicity, in thefollowing formulas we have let2/4ire0a= 1, where i is dipole moment and a is thelattice constant. The effective coupling between dipole i and j (which is the sum of theinteractions between dipole i and all the image dipoles of 3 due to the periodic boundarycondition) is= a3Z[cv2Gi(aIjnI)z?j ——(4ir/L3)E exp(_k2/4a2)cos(Z. )k/k2 (2.5)where we haven (n1,23)and k 2ir/Ln; andz = — z + n3L(2.6)The functions Q and C2 are defined as:/ — 3H(z) + (2//F)(3 +2z)exp(—zG2(z) — H(z) +(21/exp(—z2 (9 7)where the function H(z) isH(z)= Z’ f exp(—z2)dz (2.8)In (2.5) cv is chosen as ,/TE/L to get good convergence.Chapter 2. Dipoles on perfect lattices 22The Ewald summation technique is quite time consuming if the spatial configurationof dipoles is not fixed, for example, in simulations of polar liquids where the dipoles areallowed to move. For lattice systems, the effective D—D coupling is fixed in the process ofsimulation, we only need to do the Ewald summation once and for all. Therefore, in thisthesis we adopt Ewald technique and avoid any other approximations involving varioustruncation methods.The final Hamiltonian can be written asH= (2.9)One check of the convergence of is to note the fact that for lattices with cubic symmetry, the local field, with all the spins aligned in one direction, is 4ir/3. We mentionthat our numerical value of the local field with n2 = 100 in equation (2.5) is accurate toio—.The lattice structure for the body centered tetragonal lattice is shown in figure (2.1).The lattice parameters of LiHoF4 are c = 10.75A, a = 5.175A [48]. The g—value for theHo3 ion is about 13.3 [39]. We take our unit cell to have c = 2a, and choose latticeconstant and dipole moment so that our ground state energy agrees with Hansen et al[39]. Throughout chapter 2 and chapter 3 the temperature unit is= 4wE0a3kB (2.10)where ci is the dipole moment.Chapter 2. Dipoles on perfect lattices 23Monte Carlo simulations are carried out in order of increasing temperature. We startat low temperatures and the initial spin configurations are set to be ferromagnetic. Aswe gradually increase the temperature, the last spin configurations at the previous temperature is taken as the initial one for the current temperature. For the body centeredtetragonal lattice, systems with N=128. 250, 432, 686, 1024 and 2000 spins are simulated.In the case of the BCC lattice we did not simulate the largest size system. We use a conventional Metropolis algorithm. Since the interaction is long range (every dipole interactswith all the other dipoles in the MC cell), each energy update is very time consuming.Rough estimates show that the computer time is about a factor of approximately N/6of the simulation time for the short range Ising model, if parallel processing is not used.However, we could take advantage of parallel computation on the UBC IBM 3090. Byvectorizing the algorithm, we are able to increase the speed by about a factor 8. At lowtemperatures, when the acceptance rate is low, most computer time is wasted in calculating the trial energy. We can speed up a factor of 10 by simply keeping track of the localfields at each site and updating the local fields. This allows us to avoid recalculating thechange in the energy at each trial step which is very time consuming for a system withlong range interactions. This trick is very useful when we simulate disordered dipolarsystems in the next chapter, where simulations are more computationally intensive.Primary runs with about 2000 Monte Carlo per spin (MC/S) are carried out to locateapproximately the transition region. Then the temperatnres to be simulated are chosento cover the entire phase transition region. The program is written in such a way that itChapter 2. Dipoles on perfect lattices 24is restartable, i.e., the last spin configuration is always saved in case of additional runs.This can save the computer time for equilibration. Between 30,000 and 60,000 MC/S areperformed depending the size and temperature. Close to T we use about twice as manyMC/S as we use far away from fP, since the relaxation time is larger around T0 (criticalslowing down). The first half of the MC data is disregarded for equilibration.The standard Metropolis algorithm is a one spin flip process, therefore, the statesgenerated may be strongly correlated. We must justify that the number of MC/S is largeenough to lead the system into equilibrium and have a reasonable statistical accuracy.We check the equillbration in two different ways. First of all, since the quantity weanalyze with finite size scaling is < Al2 >, where M is the total magnetization, we followthe time dependence of < > and check that a steady state has been reached; makingsure that there is no drift in < M2 >. Secondly, we measure the specific heat in thefollowing two independent ways. (1) by the energy fluctuation formula211— 2kBT (.(2) by numerical differentiation of the energyc = (2.12)where U is the total energy and c is the specific heat. If the system is in equilibrium,i.e., it has a Boltzman distribution for the energy, then these two values of specific heatshould be consistent with each other. After all, (2.11) is derived from the BoltzmanChapter 2. Dipoles on perfect lattices 25distribution. Figure (2.2) plots the specific heat from (2.11) and (2.12) in the case ofLiHoF4.To estimate the statistical accuracy, we explicitly calculate the relaxation time r forWe define the time correlation function R(t) asR(t)= E(< M2(t + i)M2(i) > — < M2(t + i) >< M2(i) >) (2.13)We fit R(t) with exponential decayR(t) = Ae (2.14)and r is defined as the relaxation time. If we have total X MC/S, then the numberof independent measures of the simulation is n = X/r. The statistical error is in theorder of —%-. Except very close to T0, we choose our MC/S to be large enough that thestatistical fluctuation is small compared to the thermal average.Another trick can be used to generate more data points. The method was proposed byFerrenberg and Swendsen [27]. For lattice systems, the probability distribution of somequantity like total magnetization M at one temperature T1 is related to the distributionat T2. To see this, let us start with a general Hamiltonian— H = + hEs = S + hM (2.15)where-S is the energy, h the external field and M the total magnetization. Let K1=and 1(2=then, the probability distribution for S and M at T1 isPK1 (5, Al) = N(5, M)exp(K1(S+ hM)) (2.16)Z(Ii, h)Chapter 2. Dipoles on perfect lattices 26where N(S, M) is the number of configurations at the point (S,M) in the phase space,and Z(K1,Ii) is the canonical partition function given byzq1,ii)= E N(S, M)exp(Ki(S + hM)) (2.17)S,MWe do a Monte Carlo simulation at T1 and the probability distribution PK1 (8, M) canbe calculated from the histogram values of (S,M) of the simulation. Since a relationanalogous to (2.16) also holds at T2, the following formula of the distribution at EL2 canbe derived straightforwardly.— Pjc1(8,M)exp[(K—Ki)(8+hM)]2181(2(8,1 )—Zs,MPK(8,M)exp[(1i—Iil)(S+hM)1 ( . )Once Fk2 (8, M) is calculated, quantities like < M >, < M2 > , the specific heat and theenergy are easily computed from simple averages. The advantage of the above methodis that a single Monte Carlo run at one temperature can study the entire temperatureregion, and it is easy to implement. The disadvantage is that one needs a very high qualityMonte Carlo data at T1, since any small fluctuations could be amplified in the exponentialfunction in (2.18), and errors are hard to estimate. Nevertheless, if T1 and T2 are closeand the system is not too large (i.e. S and lvi small), and the MC data at T1 is of goodquality, (2.18) should give a reasonable result. Again, this can be viewed as another wayof checking if the simulation at T1 has reached equilibrium. We carried out the followingprocedure to produce more data points. Let T be a simulated temperature and AT thedifference between neighboring Ti’s. We generate two more points at T ± AT/3 adjacentto the simulated temperatures, from the simulation data at T. The smoothness of theChapter 2. Dipoles on perfect lattices 27resulting curve of < M2 > against T serves as another check of the adequacy of thesimulated data. Figure (2.3) shows the order parameter \/< M2 >/N for various sizesof the body centered tetragonal lattice.2.2.2 Data analysis: Finite size scaling with logarithmic correctionsIn this section we propose a finite size scaling form that includes logarithmic corrections;and the transition temperature of LiHoF4 is estimated.Larkin and Khmel’nitskii [61] calculated the total magnetization and susceptibility ofbulk uniaxial ferromagnets. For I > 0<M2 >r% xN Nr’Itnt”3 (2.19)while for I < 0< >n.< M >2 NtIlntI213 (2.20)where I = (T— T)/T , N is the number of particles and x is the susceptibility.In order to have a finite size scaling method which agrees with (2.19) and (2.20) inthe N —* cc limit, we conjecture following finite size scaling forms: For I> 0= f+(N1121 n11’73) (2.21)and for I < 0= f(N”2tln(—t)3) (2.22)We see that (2.19) and (2.20) are reproduced from (2.21) and (2.22) in the largeargument limit provided f..(x) —* x as x —* cc and f+fr) —* l/x as x —* cc. WeChapter 2. Dipoles on perfect lattices 28conclude that (2.21) and (2.22) will be valid as N112t —, cc. Note that the dimensionalityand length scale do not enter in the above scaling forms. A similar scaling form formean field like critical behavior, without logarithmic corrections, has been proposedby Botet et al.[15}. In the spirit of Botet et al., we introduce a “coherence volume”VA(t) Since the Ising dipolar system is anisotropic the correlation lengthparallel, and & perpendicular to the magnetization will be different [2], in fact near thetransition temperature cc 1j. withcc t’I1nt”3 (2.23)The scaling form suggested by Botet et al. then gives in our notationFN(i) N=f(1n ) (2.24)V\it)which is consistent with (2.19) if we set F =< M2 > /N2. This result is, however, notconsistent with (2.20). We need to use different arguments in the scaling function on thetwo sides of the transition.A different approach to ours is that of Lai and Mon [59] in the case of 4—dimensionalshort range Ising model. Their scaling form isAN= fQNh/21nTN) (2.25)where x = 1/6. However, for the 4-D Ising model the renormalization group calculationshows that the magnetization and susceptibility have exactly the same forms as in (2.19)and (2.20). Therefore,(2.25) can only be valid when the argument N”2 is small, unlikeChapter 2. Dipoles on perfect lattices 29our treatment which is valid in N”21 —* 00. The scaling form (2.25) is derived from afinite size renormalization gronp calculation for the d=4 Ising model. In the anisotropiclong range dipolar case the corresponding theory has, to the best of the author’s knowledge, not been worked out.In order to determine the transition temperature T, we least square fit the simulatedvalues of y = ln(< M2 > /N312) to a polynomial of eighth order in the logarithm of thearguments of the scaling function f in (2.21, 2.22). The procedure is repeated for anumber of trial values of T0. The temperature with the smallest statistical chi square ischosen as the best estimated critical temperature. In carrying out the above operationwe have assumed that (1) the scaling form (2.21) and (2.22) are valid and correction toscaling is insignificant. (2) The scaling function is a smooth curve so that the eighthorder polynomial is sufficient to define the function in the region of interest. Assumption(1) can be checked by the quality of data collapse which we discuss next; and assumption(2) seems reasonable from the scaling forms known so far.The best fit for the LiHoF4 structure is shown in figure (2.4). Figure (2.5) shows thesimilar fit for the BCC lattice. We note that the scaling function for the body centeredtetragonal and cubic lattices are identical in the asymptotic region while there is somescatter at small arguments (when the temperature is close to fEe). Theoretically, oneexpect these two scaling form to be the same from the argument of, universality. Thefact that there are slight disagreements close to T0 niight be due the quality of the MonteCarlo data, since around T the relaxation time is very large. In the asymptotic regionChapter 2. Dipoles on perfect lattices 30(x —* oc), we also check from figure (2.4) and (2.5) that the slopes of the scaling functionf are ±1 which is consistent with the discussion after (2.22).As a further consistency check, we also calculate the “Binder number” [12] [11].M4>U=l— < (2.26)3 <LW >2From finite size scaling theory, the dimensionless quantity U must be a function of L/,i.e., U = g(L/) where is the correlation function and L is the size. When U is plottedagainst temperature for systems of different sizes, one sees, that at high temperatures,the distribution of M is a Gaussian and U 0. At low temperatures < M >< P/I2 >2and U 2/3 while at T = T0, since diverges, U = g(0) for all sizes. Therefore,the curves for different sizes will all intersect at T and U*. The theoretical values ofU* from a finite size renormalization group calculation for the four dimensional Isingmodel is 0.27 [16] and we expect U* to be the same for the three dimensional Isingdipoles. A plot of U vs. T is given in figure (2.6) for body centered tetragonal lattice.The same plot for BCC lattice is shown in figure (2.7). As can be seen the intersectiontemperature is close to the best fit temperature of finite size scaling (2.21, 2.22). The U*is approximately the theoretical value. Note that the above method of locating T doesnot involve any fitting process.To further check the scaling behavior we also repeated the analysis with the generalscaling form= f(N2t) (2.27)Chapter 2. Dipoles on perfect lattices 31We note that in the N —> c limit < M2 > /N2 oc tj2 fort < 0 , and < M > /N2 cxfor t > 0. We can easily show that the critical exponents fi and are related to x and ythrough:= x 1 (2.28)= 2—x (2.29)The scaling form (2.27) would have been appropriate if we had critical behaviorwithout logarithmic corrections. The best fit with mean field exponents x = 1.5, y = 0.5is shown in figure (2.8) for the body centered tetragonal lattice. The corresponding fit inthe case of the BCC lattice is given in figure(2.9) [105]. We note that the best fit in figure(2.8) is inferior to that of figure(2.4) especially in the low temperature region where thelogarithmic correction is significant. However, it is possible to get a good fit to the datawith three free parameters x, y, T in (2.27) (see figure (2.10) ), although there seems tobe some systematic deviations at low temperatures.We have shown that for systems with Ising dipolar interactions one can get significantimprovement over the mean field result for the finite scaling fit by including logarithmiccorrections to the mean field expressions. It is also possible to fit the data with nonclassical exponents, but this interpretation is not favored by theory. This suggests that thedata collapse method can give misleading results when theoretical support is not available especially when there are many fitting parameters. The same difficulty also occursin the analysis of experimental data where Griffin et al. [32] could not unambiguouslyChapter 2. Dipoles on perfect lattices 32distinguish the logarithmic corrections from weakly nonclassical critical behavior. Thebest fit for the critical temperature is not very sensitive to the inclusion of logarithmic corrections, e.g. in the case of the 13CC lattice our best estimate is T = 0.78T,which should be compared with the value T0 = 0.79T found by Xu et al. [105] withoutlogarithmic corrections to scaling.We conclude that when we use the scaling form with logarithmic corrections (2.21,2.22), with one fitting parameter T we can determine the critical temperature withgood accuracy without going to extremes in computational effort. In the case of thebody centered tetragonal lattice we find for the critical temperature EP = 1, 145T0. Wecompare this with the experimental transition temperature of LiHoF4 by determiningthe value of T0 in (2.10). Our procedure is to equate our ground state energy (which wefound numerically to be -0.7211kBT)with the calculated value -0.9499kBK of Hansenet al. [39] for LiHoF4. This gives T = 1.51K in good agreement with the experimentalresult 1.538K of Griffin et al. [32]. Our results are thus consistent with the interpretationthat the dipolar interactions are the dominant ones in this material.2.3 Critical behavior of Thing dipoles on simple cubic lattice: antiferromagnetic phase transitionIn this section we report Monte Carlo simulations for Ising dipoles on the simple cubiclattice. The finite size scaling method is used to locate the transition temperature andcritical exponents. The critical temperature obtained is in good agreement with previousChapter 2. Dipoles on perfect lattices 33simulations and high temperature series expansions. The susceptibility critical exponentagrees well with the high temperature series expansion result 8/7. The best estimatesfrom simulation for /3 and a are approximately 1/7 and 4/7 respectively.In the previous sections we have shown that for Ising ferromagnetic dipoles at D=3,the critical behavior is that of mean field theory with logarithmic corrections. The aboveconclusion can also be derived from renormalization group theory [3]. Experimentalresults on ferromagnetic materials of the type LiHoF4 are consistent with the theoreticalpredictions [32], and the logarithmic corrections has been studied by computer simulation[106]. However, much less is known about the anti-ferromagnetic case. Let us considerIsing dipoles on a simple cubic (sc) lattice and assume that dipoles can only point in the± z directions. The ground state of the spin configuration is that of anti—ferromagneticrows. That is, the spins line up in z—direction, while in the xy—planes they exhibit a Neelorder ([65],[58],[99]). This model has been used previously to simulate dielectric fluidssince the exclusion of translational degrees of freedom permits studying larger systemwithout excessive computer time [1]. The model has also been applied to a study theeffect of dipole interactions on nematic ordering [85].Theoretically, for a system with long range interactions one would expect mean fieldtheory to be fairly accurate, in fact it gives an exact solution for a uniform long rangemodel [40]. However, previous Monte Carlo simulation results on a simple cubic lattice([58], [99], [85]) show that the mean field results overestimate the transition temperatureby approximately a factor of 2. Aharony [3] used the renormalization group method toChapter 2. Dipoles on perfect lattices 34show that for antiferromagnetic Neel ordered dipolar systems, the critical behavior is inthe same universality class as the corresponding short range Ising systems, therefore, itshould be governed by non classical critical exponents. However, the ordered phase fordipoles on the sc lattice is not conventional Neel antiferromagnetic ordering. For the sclattice Kretschmer and Binder [.58] found that dipoles are more easily aligned in the zdirection while it only costs a small amount of energy to disorder in the xy—planes. Thisunusual behavior is due to the strong anisotropic nature of dipolar interactions. The onlytheoretical work of the model, of which we are aware, is that of Voronstsov-Velyamonovet al. [99]. They used a high temperature series expansion method to calculate thecritical temperature EP 2.58T9 and exponent 7 8/7. Due to the long range nature ofdipolar interactions, the uncertainty involved in the truncation of the expansion is hardto estimate. Therefore, it is worthwhile to investigate the critical behavior in more detailthrough Monte Carlo simulations.Since it is known that for this model there are many low-lying metastable states [58],the initial configuration is chosen to be the ground state and the system is gradually warmed up. We simulated a total of five sizes L=6 (N=216), L=8 (N=5l2), L=lO(N=1000). L=12 (N=1728), L=14 (N2744). For L=6 to L=12 we used 20,000 MonteCarlo Steps /spin (MCS/spin) and the first 10,000 MCS/spin were discarded for equilibration. For the largest size, N=2744, we simulated a total of 4000 MCS/spin, andthe first 1000 MCS/spin were discarded. The quantity that is monitored closely is themean square < M2 >of the total staggered magnetization. The thermal equilibration isChapter 2. Dipoles on perfect lattices 35checked in two ways. Firstly, we monitor < M2 >, to check that after a certain nnmberof MCS/spin this valne tends to stabilize. Secondly, we calculate the specific heat in twodifferent ways: one by fluctuation of the internal energy c =< U2— < U >2> /kT2 andthe other by c = dU/dT. The results are shown in figure (2.11). Except close to T0, thetwo specific heat values agree within statistical error indicating that thermal equilibrium has been reached. Figure (2.12) shows the saturated magnetization as function oftemperature for five different sizes.The finite size scaling form (2.27) is adopted and the same procedures as in section2.2 are followed. Note that since we do not know the critical exponents, we have tosimultaneously fit the three parameters x, y and T. We find T0 2.60T0, /3 1/7 and8/7. The critical exponent /3 1/7 is clearly different from the mean field exponentof 0.5. The critical temperature which we obtained from finite size scaling is in excellentagreement with the result of Kretschmer and Binder [58]. The best fitted data collapseis shown in figure (2.13). Figure (2.14) shows the specific heat data for various sizes.We see that the transition temperature from the specific heat plot is consistent with thefinite size scaling results.We find that for dipoles on the simple cubic lattice, it is much harder to equilibratethe system than in the ferromagnetic case ( for example, dipoles on the body centeredtetragonal lattice). This is why the quality of data collapse in figure (2.13) is worsethan the corresponding ferromagnetic case figure(2.4). The reason is well explored byKretschmer and Binder [58]. It is found that dipoles are more easily aligned in theChapter 2. Dipoles on perfect lattices 36z-direction while it only costs a little energy to disorder in the xy—planes. Therefore,the system behaves in a quasi two dimensional manner. There are many low energymetastable states similar to a “spin glass” phase. This unusual behavior is due to thestrongly anisotropic nature of dipolar interactions. Therefore, to accurately calculate thecritical behavior, much longer simulations are needed.Ghapter 2. Dipoles on perfect lattices 37_VFigure 2± Positions of Ho3 ions in the cell of LiHoF4. a=b=5.l75Aand c=1O.75A.Chapter 2. Dipoles on perfect lattices 382-1.5C105I I I IAA I I I ( I I- l II I1I I IT/TO1.5I I J I I2Figure 2.2: Specific heat of 686 spins on a bodynumerical differentiation of the internal energy(triangles).centered tetragonal lattice calculated by(full curve) and the fluctuation formulaAA0-AChapter 2. Dipoles on perfect lattices 39I I I I I I I1—0.8 - V0.6 -zA0.4 - . N=428t’ N=686Vci N=20000.2 -0 I I I I I I0.5 1 1.5T/T0Figure 2.3: Temperature dependence of < Al2 >1/2 /N for body centered tetragonallattice.Chapter 2. Dipoles on perfect lattices 40I I I I ‘ I7V A 1100a ao o0 6A N=128ri N=250o N=432o N=686-2- N=1024 °oo N=2000I I K K I I K I I-4 -2 0 2 4xFigure 2.4: Data collapse for the body centered tetragonal lattice using the scaling forms(2.21) and (2.22). The best fit is for T = 1.145T0 The abscissa is the log of the argumentof the function f.Chapter 2. Dipoles on perfect lattices 41I I a I I I I I I I I47C’IC4 flOOAC’l 0 A0AV 0-A N=128oA0t N=250 0*o N=432-2 - ° N686o N=102400I I I I i t I I I I I I t i t-4 -2 0 2 4 VxFigure 2.5: Data collapse for the BCC lattice using the scaling forms (2.21) and (2.22).The best fit is for T = 0.78T The abscissa is the log of the argument of the functionf±•Chapter 2. Dipoles on perfect lattices 42I I I I I I I I I-0.6 -0.4-U.AN=250a N=4320.2 - o N=1024\...\I0 I I I I I I I0.8 1 1.2 1.4T/T0Figure 2.6: Temperature dependence of the “Binder number”, U, for body centeredtetragonal lattices of different sizes. The curves are only meant as guides to the eye.Chapter 2. Dipoles on perfect lattices 430.6U0.40.200.6Figure 2.7: Temperature dependence of the “Binder number”, U, for BCC lattices ofdifferent sizes. The curves are only meant as guides to the eye.0.8 1 1.2T/T0Chapter 2. Dipoles on perfect lattices 44I I I I I I I I I4C4zA .AQ00V A0CO—- o 0UN=128o N=250N=432 0%0-2 - N=686o N=10240o N=200041 I I c-4-2 0 2 4In(Nt)Figure 2.8: Data collapse for the body centered tetragonal lattice, using the scaling form(2.27) with mean field exponents (x = 1.5, y = 0.5). The best fit is for T = 1.17T0Chapter 2. Dipoles on perfect lattices 454. a i aa2-rAC40V 0- 0-000N=128-2 - N=250o N=4320o N=686 00o N=10244. I I I I I I I-4 -2 0 2 4In(N2t)Figure 2.9: Data collapse for BCC lattice, using the scaling form (2.27) with mean fieldexponents (x = 1.5, y 0.5). The best fit is for T = 0.79Tchapter 2. Dipoles on perfect lattices 46I I I 4 1 I I I—‘0--ZAA 0000 0A 0AV o-2A N=128ci N=250 0o N=43200oo N=686 0- o N=1024 0o N=2000—6 I i I-4 -2 0 2 4In(N3’t)Figure 2.10: Data collapse for body centered tetragonal lattice, using the scaling form(2.27) with three parameter fit x = 1.72, y = 0.59 and T = 1.126T0 (order parameterexponent 0 = 0.24; susceptibility exponent y = 1.22).CChapter 2. Dipoles on perfect lattices1.50.5047• Figure 2.11: Specific heat of 12 spins on a simple cubic lattice calculated by numericaldifferentiation of the internal energy (full curve) and the fluctuation formula (triangles).1T/TO1 2 3 4Chapter 2. Dipoles on perfect lattices 48I I I I I I I I I• p10.8zA(.4V0.4 N=128i N=686o N=20000.20‘ I I I I I I t1 2 3 4T/T0Figure 2.12: Temperature dependence of the magnetization < Al2 >1/2 /N for the simplecubic lattice.chapter 2. Dipoles on perfect lattices 49I I I I20 00000000 00D000 QO 0000000AUI’ 0Az 0A0AC’10 AAA 0A—-A00 £A02**A N=256 oo-4- N=512 000o N=10000000o N=1728o N=2744—6 I I Io 2In(Nt)Figure 2.13: Data collapse for the simple cubic lattice with the scaling form (2.27). Thebest fitted parameters are x = 1.8, y = 0.7 and T = 2.60T0 (order parameter exponent= 1/7; susceptibility exponent ‘y = 8/7).Cchapter 2. Dipoles on perfect lattices1.510.5050Figure 2.14: Temperature dependence of the specific heat for simple cubic lattices ofdifferent sizes. The curves are only meant as guides to the eye.2T/T01 2 3 4Chapter 3Dipoles on diluted latticesRecently much effort has been put into investigating the cooperative properties of crystals containing off—center ions (for a detailed review see Vugmeister and Clinchuk [100]).These phenomena have attracted much attention because of a general interest in disordered systems, in particular spin glasses [12]. The system with randomly positioned(quenched) dipoles is one of the possible realizations of the spin glass state and couldprovide a test for the theory. When dipoles are randomly distributed in a weakly polarizable medium (such as Li ions in KOHj they appear to exhibit only a spin glassphase but no ferroelectricity [29]. In contrast, similar impurities in a strongly polarizablemedium (such as KTaO) appear to exhibit a ferroelectric phase [101] [52] althoughthis conclusion is controversial [42]. There are even controversies about whether there isa dipole glass state [89] [84].Theoretically it is expected that the presence of polarizable ions in the lattice may playan important role in determining the nature of the low—temperature phase [67]. Takagi[96], for example, showed in a mean field calculation for dipoles on a simple cubic lattice,with polarizable ions in the body centered positions, that the low—temperature phase51Chapter 3. Dipoles on diluted lattices 52could be ferro— or antiferroelectric, depending on the magnitude of the polarizability.For reasons of simplicity we will not address the question of polarizability in this thesis.In order to explain the occurrence of spin glass phases arising from dipole—dipoleinteractions one must consider fluctuations in the spatial distribution of the dipoles.These were neglected in the Debye mean field theory. Klein et al. [53] made an attemptto correct the spatial fluctuations in the case of pure dipole interactions using the methodof “random mean field”. However, the boundary conditions employed by them precludeuniform ferroelectric ordering. The reason is that they use spherical geometry, i.e., theysum up all interactions within a radius R and let R —b on. A uniformly polarized statewould have to overcome the depolarizing field of the sphere and a ferroelectric state wouldhave a complicated domain structure [5]. Similar considerations apply to the Monte Carlosimulations on pure dipoles performed by Medina et al. [69]. Calculations have also beencarried out [54] [49] with interactions which scale with distance as r3 with random sign,i.e. which are not of the pure dipolar form. The question of whether a purely dipolarrandom system can order ferroelectrically thus remains open.In this chapter we will, in section 3.1, construct a mean field theory using availabletechniques from spin glass theory [50]. We will show that this mean field theory ismanifestly shape independent without an external field. As shall be seen, the theorypredicts ordering, for sufficiently high concentrations of dipolar impurities. This orderingmay be either ferro— or antiferroelectric depending on lattice structure. Below the criticalconcentration the theory predicts a low—temperature spin glass phase. In section 3.2 theChapter 3. Dipoles on diluted lattices 53critical concentration is estimated by zero temperature Monte Carlo simulations. Thepredicted ferroelectric ordering at sufficiently high concentration is confirmed by finitetemperature Monte Carlo simulations discussed in section 3.3.3.1 Model and mean field theoryIn building up a mean filed theory for a random dipolar system, we must keep two thingsin mind: (1) we must take the spatial fluctuations into account. (2) the theory shouldbe shape independent in the absence of an external field.We consider a system of N spins that can only be oriented in the +z directions.The coupling between spins is pure dipole—dipole interaction. Assuming that the sampleshape is such that the depolarization field is uniform, the Flamiltonian can be writtenH= —LJd(’%)ssj— EJ0s= (3.1)i i<jwhere the dipole—dipole interaction is=(3cos2Oj) — 1) (3.2)47rwith ,u the dipole moment, Oj is the angle between the vector ij connecting two dipolesites and the z axis, and the spin s can only be +1. If there is a net polarization, thefirst sum in the middle expression in (3.1) will depend on sample shape, because of thedepolarization field. With metallic boundary conditions J0 is such that it cancels thisfield. For periodic boundary conditions J0 = 0. In general, for a homogeneous system ofChapter 3. Dipoles on diluted lattices 54regular shape, J0 will be proportional to the total polarization allowing us to define aneffective dipole—dipole interaction J (ii,To be more concrete, let us consider a regular sample shape such that the free samplehas a uniform depolarizing field, Ed (as in an ellipsoid, sphere, needle or thin slab) andis of the formp D£4=—-——-vEsi (3.3)where V is the volume of the sample; D3 is the shape dependent depolarization factor,for example, D = 4w for a thin slab, D3 = 4ir/3 for a sphere, and D = 0 for a needle.With metallic boundary conditions we have= —Ep (3.4)From the above relation and the Hamiltonian (3.1) we can calculate the effective dipole—dipole interaction2—4 4-4 p SJ(rj, r) = Jd(r1) +---——--i (3.5)‘tirE0 vOne can see it is the competition between the first term, the dipole—dipole interaction,which can be positive or negative depending on the angle Oj, and the second term whichfavors the ferroelectric ordering, which determines the nature of the phase transition.Intuitively, if the second term dominates there will be a ferroelectric phase transition.In the case of random spatial distribution of spins, the strong fluctuations in the coupling J(, ñ) suggests the application of spin glass theory. We note that the restrictionto Ising spins greatly simplifies the calculations, in that the part of the interaction whichChapter 3. Dipoles on diluted lattices 55depends on the spatial separation enters as a scalar factor multiplying the spins. It isthis feature which allows us to carry out the replica trick to be discussed later.In the case of off—center Li+ impurities, however, either of the the six (100) orientations or the eight (111) orientations are possible. The coupling J() then becomesa tensor with spatially dependent components, and the calculations will become moreinvolved. For the similar model with six (100) orientations, the calculation was carriedout by Xu [107]. Since we are dealing with an idealized situation we will not attempt toincorporate such effects here.We complete the definition of the model by specifying the spatial distributi?n of thedipoles. First, we consider the “amorphous” case in which we assume that the continuouspair distribution function g(ñ) is given. We approximate the N—particle distributionfunction by using the Kirkwood superposition approximationgp..r(ri, ‘2, “N) = [J g(i — Ps). (3.6)permutationsThe probability distribution for the strength of the Ising coupling constant Jjj for arandomly selected pair of dipoles is then given byP(J) = fd3rg(flS(J - (1) (3.7)where the normalizationfd3rg(P = V (3.8)has been used.Chapter 3. Dipoles on diluted lattices 56Averaging over spin configurations can be carried out in closed form using the replicamethod if we approximate the coupling distribution (3.7) by a shifted Gaussian— exp(—N(J— (1/N)J)2/24)PG(J)— 2K(J/N)’/ 3.9where we can calculate the moments J1 and J2 from the distribution (3.7):= Jd3rJ(flg(fl (3.10)As we will show later Ji will be independent of sample shape. We have00 1= NJ dJJ2P(J)—(3.11)In what follows we will neglect terms of the order 1/N, such as the last term on theright—hand side of the above relation, since in the thermodynamic limit they go to zero.Once Ji and J2 are known, the phase diagram can be calculated using the replicamethod as described, for example, in Kirkpatrick and Sherrington [50]. The steps of thecalculation will not be repeated here. We simply state the results. There are two orderparameters in the theory. One of them is the polarization in =<< 5 >t>s where first athermodynamic average (denoted by <>) is performed for given spatial arrangements ofspins followed by averaging over the spatial arrangements (denoted by <>). The secondparameter is the Edwards—Anderson [26] order parameter q =<< s . These twoorder parameters are determined from the following integral equations:1 [00 2rn=— I dze2t1 r°°— I d —z2/2j22w-ccChapter 3. Dipoles on diluted lattices 57wheret = tanh(/3[(Jim +J2/ijz)] (3.13)(3.12) and (3.13) is the replica symmetric solution which is not valid at low temperatures[12]. Figure (3.1) shows the phase diagram of Parasi replic symmetry breaking solution [78] [79] [80] which is believed to be exact. The phase diagram for Kirkpatric andSherrington calculation is the same except the phase boundary between the spin glassphase and ferroelectric phase is slightly curved [50]. The phase diagram is determinedby the ratioJ1/J2. When J1/J2 > 1 there is a low—temperature ferroelectric phase withtransition temperature J1 / kp. Physically, this is the regime where the average couplingJ1 is larger than the fluctuations of the couplings, J2. Note that in the limit J2— 0 andwith J1 the Lorentz field, the self—consistent equation for in is of the same form as in thesimple Debye theory of the introduction. Even if J2 0 but J2 < Ji the ferroelectrictransition temperature, which depends only on J1, will be the same as in Debye theory,but the order parameter will in general be different due to the fluctuation term J2.In order to be more concrete, we consider the case where the dipoles occupy discretelattice sites with concentration c, i.e., a site is occupied by a dipole with probabillty c.The coupling constant distribution can be written asP(J) = 6(J- J()) (3.14)where N3 is the number of sites in the lattice (the number of dipoles is then N = N3c).Chapter 3. Dipoles on diluted lattices 58With the above distribution we can calculate J1 as= N.JdJJP(J)= cEJ(i)=c(Jd(ij) + (3.15)where the sum i runs over all the lattice sites. The local field at the site when the sampleis uniformly polarized (all spins point in one direction) is= El+E2+Ed (3.16)whereB1= 3V (3.17)is the Lorentz field due to polarization charges on the surface of an imaginary sphericalcavity and B2 is the field due to the dipoles inside the cavity. For a lattice of cubicsymmetry we have B2 0. Note that £2 0 for the body centered tetragonal lattice( onecan easily calculate it numerically). Substituting (3.16) into (3.15) noting the expressionfor Ed,, given in (3.3), we find that the shape dependent terms cancel. We have— QTI (3.18)EQwhich is shape independent. Similar expression like (3.15) can been derived for J2straightforwardly and we have—____— (3.19)4cirE0aChapter 3. Dipoles on diluted lattices 59where a3 = V/N anda6 37?=- 1)2 (3.20)This value only depends on the lattice structure and can be calculated numerically. Wefind a 3.65, 2.08, 1.93 for SC (simple cubic), FCC and BCC lattices, respectively. Thereexists a critical concentration c* above which J1/J2 > 1, i.e., there is a low—temperatureferroelectric phase. We calculate the critical concentrations c* = 0.76, 0.25, 0.21 for theSC, FCC and BCC lattices.The low—temperature phase for a perfect crystalline SC lattice is antiferroelectric [90][65] with opposite polarization on sublattices < hkl> with h + k even or odd. We can,in the spirit of Takagi [96], recast this situation by reversing the sign of the spins on oneof the sublattices and we find by evaluating the dipole sum numerically that= 0.426cN5v2 (3.21)Since J1 in this case is larger than the value in the ferroelectric phase, the ferroelectrictransition is preempted. The present model thus predicts that the simple cubic lattice willhave an antiferroelectric ordering at low temperatures above an impurity concentrationof C = 0.45. Below this concentration the theory predicts a low—temperature spin glassphase.Next we consider the case where dipoles are located in the centers of hard spheres ofChapter 3. Dipoles on diluted lattices 60radius d. We assume the simple form for g()0 forr<d1 forr>dFollowing the similar steps in deriving (3.18) and assuming that £2 = 0 when averagingover dipole spatial distributions, we have2’ 2J1 = J 3 (3.22)Kand by evaluating the integral in (3.11) analytically32f/522 4wc0d3where f is the filling fraction f = irNd3/(6V). Equation (3.1) is not a bad approximationfor the system of sequentially random parked hard spheres below its jamming limit (f1/3, see Gotoh et al. [33]). We then find that J1/J2 > 1 for f> f* = 1/10. The presenttheory thus predicts a ferroelectric phase for dipoles on a system of randomly parkedhard spheres at high concentrations.We can show that without external field the theory is shape independent. Since thethermodynamics only depends on J1 and J2, we need only to prove that they do notdepend on shape. We have already showed that J1 does not depend on the sample shape.Now let us look at J2 . We have following relation:J2 fdJJ2P(J) f dJJ2fd3rS(J — fdr (3.24)Chapter 3. Dipoles on diluted lattices 61Since the intergation is over a function which decays like 1/r6, it is obviously shapeindependent. Therefore our mean field theory is manifestly shape independent withoutan external field.To summarize, we have constructed a mean field theory which takes the spatial fluctuations of the dipole—dipole interaction into account. Several approximations are involved.First of all we assume that dipoles are randomly distributed; in the diluted lattice casethis means that the concentration c should not be too big. Secondly we approximatethe distribution of the coupling to be a Gaussian, i.e., we neglect the higher moments ofthe distribution function. It turns out that this approximation can be justified in thathigher moments are irrelevant in the thermodynamic limit, as long as they are finite [12].Thirdly in the case of randomly parked hard spheres we assume F2 0. We only dealwith a sample with regular shape so that we can write down a simple Hamiltonian (3.1),although we believe that the same results should hold for any irregular sample shapes.We would like to test the mean field predictions for various lattice structures and randomly parked hard spheres, with Monte Carlo simulations. However, due to the limitedcomputer power and the intensive calculations involved in simulating systems with longrange interactions, we only simulated the case where dipoles randomly occupy a BCClattice. This is discussed in the following sections.Chapter 3. Dipoles on diluted lattices 623.2 Zero temperature Monte Carlo simulationTo estimate the minimum dipole concentration for ferroelectric ordering, we search forand analyze ground state configurations of the BCC system. Finding the ground state fora spin glass is equivalent to a NP—complete problem [70] [31]. This means the computertime needed grows exponentially as the size of the system increases. Therefore, we canonly approximately locate the ground state numerically.The method we employed is a combination of the steepest descent and simulatedannealing methods for finding approximate ground states [50] [77]. The basic MonteCarlo cell is cubic and periodic boundary conditions are used. As described in theprevious chapter, the effective coupling J is calculated by Ewald summation. For eachsample a list of local fields is compiled for ea.ch initial spin configuration, andthe spin with the largest positive local energy is flipped. The list is updated continuouslyuntil no positive local energy can be found. This steepest descent procedure is a zerotemperature quench and leads to a local minimum in the energy. Further improvementcan be achieved by warming the system, using the Metropolis algorithm, to a temperatureT which we took to be four times the mean field transition temperature. The systemis then cooled down gradually to zero, reheated to T/2 and cooled down again. Thisprocedure is repeated six times. The lowest energy spin configuration found during theprocess is saved. The minimum energy obtained through the above annealing process istypically 3 — 10% lower after the original steepest descent.Chapter 3. Dipoles on diluted lattices 63For each sample we start from 50 different random initial spin configurations and alsofrom an ordered state. The lowest energy is chosen to be the approximate ground stateenergy for the sample; and we generated 50 different samples for each concentration. Thenumber of spins in each sample is close to 450. The resulting ground state polarizationis plotted in figure (3.2), the error bars are the sample to sample standard deviation.We estimate the critical concentration for ferroelectric ordering to be 0.3 + 0.1 which issomewhat larger than the mean field result 0.21 found in the previous section.We have also investigated the local field distribution F(E) in the “ground state”.Kirkpatrick and Varma [49] give an argument that F(E = 0) = 0 for random 1/r3interactions. The argument is as follows: suppose P(0) 0. Then the ground stateis thermodynamically unstable at low temperatures. Flipping the orientation of spinsin zero field does not cost any energy. The number of such spins may be infinitesimal;however, flipping these spins alters the field at other spins and the new field distributionis the sum of the previous field and the new extra field due to the flipping of the spinsin zero field. Some of the spins will now be in the wrong direction and we must in turnflip them and find their effect on the rest of the spins, etc. If the number of spins in thewrong direction due to the extra field of a spin is infinitesimal, the change in energy isinfinitesimal and a finite value of the probability distribution at zero field is allowed. Ifthis number is large the process may not converge and the process initiated by flippinga single spin at zero field leads to an extensive reduction in the energy of the glass. Ifthis is the case the stable distribution must have the value zero at zero field. In the caseChapter 3. Dipoles on diluted lattices 64of random i/z’ interaction system, consider a given dipole in zero field. If we flip it, alldipoles at a distance r from it which had fields less than 2z/r3 must also flip. Theirnumber is proportional tof drr2 P(E)dE f drr2(F(0)) (3.25)which is proportional to P(0) inN where N is the total number of dipoles. These willin turn cause flips, etc. The only sensible way out is that P(0) = 0.However, the above argument does not tell us how P(E) approaches zero aroundF 0. Kirkpatrick and Varma [49] suggest that P(E) ‘-.‘ F° , with a = 0.5, for small F.This behavior is compatible with their simulations, but no decisive conclusion of whethera = 0.5 can be drawn from their data. The particular value of a = 0.5 could accountfor the observed T312 behavior of the specific heat of 0H in KC1 at low temperatures[29]. In general the temperature dependence of the specific heat is given by C[49]. However, the origin of the T312 law still remains controversial [62]. We plot in figure(3.3) our results for the local field distribution for c = 0.1. In figure (3.4) we show alog—log plot of the field distribution for small F. Our best fits are a = 0.37 for c = 0.1;a = 0.46 for c = 0.3. In the latter case we may be above the critical concentration forferroelectric ordering. The exponents presented here should be taken with caution, sincethe size in the simulation is not that big (about 450 particles) and we have not studiedthe size effects very carefully. However, these values are significantly different from theChapter 3. Dipoles on diluted lattices 65value cv = 1.0 that have beeu obtained in the infinite range Sherrington—Kirkpatrick (SK)model (see Palmer and Pond [77]).The local field distribution function for a classical freely rotating dipolar system hasalso been investigated analytically by Gulácsi et al. [35]. The basic idea of their methodis proposed by Ma [66]. If we neglect the spin—spin correlation in the ground state, for atotally random system, the local field distribution is Gaussian [66] from the central limittheorem. By including the correlations to first order they found the distribution to beP(E) r.’ exp(—aE2— bE) (3.26)where a, b are constants. This result is not compatible with our numerical simulationdata at small E. The corresponding theory for Ising dipoles, to the best of author’sknowledge, has not been worked out yet.3.3 Finite temperature Monte Carlo simulationWe have performed simulations at non—zero temperature for concentrations above thecritical one for a ferroelectric phase. The calculations were performed starting from arandom spin configuration. Each sample evolved at successively decreasing temperatureand the quantity monitored closely was < P2 >. The initial relaxation time for thisquantity was about 25 Monte Carlo steps per spin (MCS) for the 800 spin systems andwe used 10—20 000 steps per spin after thermalization. Periodic boundary conditions andChapter 3. Dipoles on diluted lattices 66Ewald summation are employed as before. At each concentration the thermodynamicquantities are averaged over 6 to 10 samples depending on the sample size. The longrange nature of the interactions is found, as observed in chapter 2, to lead to importantfinite size effects. We therefore carry out simulations at each concentration for at leastsix different sizes and use finite size scaling theory for < P2 > to locate the transitiontemperature. In this approach we assume, in analogy to (2.27)<p2 >= Nxf(cN9 (3.27)for the mean square polarization, where e = (T— T)/T and the exponents x and y arerelated to the critical exponents, B (order parameter) and 7 (susceptibility) throughX—7y = 1+ 2/3y = 2 (3.28)In mean field theory 7 = 1, $ = 1/2 and x = 3/2, y = 1/2. The fit to T using meanfield exponents is shown in figure (3.5, 3.6). The data points are arrived at by averagingover 6—10 samples. We find that without such averaging it is not possible to collapse thedata. All error bars represent sample to sample fluctuations and do not include thermalfluctuations which we estimate to be smaller. We are not able to significantly improvethe fit by employing non—mean field values for x and y. Due to the large sample—samplefluctuations it is also not possible to distinguish the logarithmic corrections (if thereare such corrections in diluted dipolar systems). From the reasonably good fit of theChapter 3. Dipoles on diluted lattices 67data, our simulation results suggest that at least down to the concentration of 50% , thecritical behavior is mean field like. It has been noted that there might be cross overs;in particular, the critical exponent 7 may depend on the dipole concentration. Thecrossover from “pure critical behavior” [61] to “impure critical behavior” [3] was studiedby Liebmann et al. [63]. Experimentally the concentration dependence of the exponentfor the dipolar—coupled ferromagnet LiTbY1_F4was measured by Folkins and Griffin[28]. They found for the dipolar concentration p close to 1, 1 which is the mean fieldresult. For p below 0.6, significantly deviates from 1. It increases with the decreaseof the concentration p, at p = 0.2,2.0. Our data fail to show the experimentallyobserved behavior. Apart from the fact that the Monte Carlo simulation data have largesample—sample error bars, the following effects might be responsible for the discrepancy.(1) There might be some corrections to scaling to the form (3.27) for disordered dipolarsystem. (2) LiTbYF4 is not a permanent Ising dipole ferromagnet, i.e., the ground stateis not truely doublet. Thus our Hamiltonian may not be applied to this compound.The specific heat obtained from energy fluctuations and averaged over the samples isplotted in figure (3.7) for c = 1 and c = 0.5 respectively and for different sample sizes.We note that in contrast to c = 1 there appear to be no finite size effects on the specificheat for c = 0.5. The concentration dependence of the transition temperature is plottedin figure (3.8) together with the predictions of the mean field theory of section 3.1.Chapter 3. Dipoles on diluted lattices 683.4 Discussion and summaryIn this chapter we have constructed a simple mean field theory for disordered dipolar Isingspins and compared it with the results of Monte Carlo simulations. Both the theory andthe simulations for the diluted BCC lattice show that there is a critical concentrationabove which a ferroelectric phase exists at low temperatures. As could be expected thetransition temperature of the mean field is higher (‘S.’ 30%) than found in the simulationdata. We have also investigated sore properties of the ground state in the spin glassphase and found a local field distribution P(E) E for low fields with a 0.40. Thefact that this exponent is much smaller than the value a = 1 for the SK model suggeststhat the low—temperature diluted dipolar system is less frozen than the spin glass withinfinite range interactions. It is interesting in this context to note that Reich et al. [84]found that for very dilute dipolar magnets (LiHoY1_F4with x = 0.045) the spinswill not freeze at low temperatures while at higher concentrations (x = 0.167) they didobserve the freezing of spins. However, in this material, the situation might be morecomplicated. At very low temperatures quantum effects may come into play [103]. Thusfor very diluted pure Ising dipolar system, the question of whether there is a spin glassphase remains open.Ghapter 3. Dipoles on diluted lattices169Figure 3.1: Phase diagram of the SK spin glass model (adaped from Binder and Young[12].kBTJ2ParamagnetSpin glassFerromagnet12chapter 3. Dipoles on diluted lattices 70I I I j I I •I I I I I I I I I1 00.80.6<,1>0.40.21I I I I I I I I I I I I0 0.2 0.4 0.6 0.8 1CFigure 3.2: Ground state polarizations for the diluted BCC lattice at various concentrations. Error bars are sample to sample fluctuations.Ghapter 3. Dipoles on diluted latticesFigure 3.3: Ground state local field distribution at c=O.1. Here E is the absolute valueof the local field.I • I I—.—710.020.0150.01P(E)0.00500I I I - I I I -0.2 0.4 0.6EChapter 3. Dipoles on diluted lattices72I I I I 5—3 — I I 1 1 11 1 l I -04,-3.50, -,FFFln(P(E)) [//r-4.5I I II-6 -5 -4-3ln(E)Figure 3.4: Log—log plot of the ground state field distribution for small E at c=O.1. Thebroken curve is a least square fit with c = 0.37.Chapter 3. Dipoles on diluted lattices 73—I•12IIIfV 0N=175 4o N=302o N=480-2 - o N=716I I I-2 0 2 41n(eN)Figure 3.5: Data collapse for diluted BCC lattice at concentration c=O.7, with mean fieldexponent x=1.5, y=O.S (see (equation 3.27). The best fit is for T = O.53T0Chapter 3. Dipoles on diluted lattices 744 I2-jzVI V0-a N=343o N=512o N=729-2 o N= 1331r I I I i I-2 0 2 4ln(EN)Figure 3.6: Data collapse for diluted BCC lattice at concentration c=O.5, with mean fieldexponent x=1.5, y=O.5 (see (equation 3.27). The best fit is for T = O.37T0chapter 3. Dipoles on diluted lattices 752N=1024N=432oN=2501.5I:1-8N1331oN729o N=512(.1) 04I..6o .‘0400-a8004 *00*11*0 ‘ I I I I0 0.2 0.4 0.6 0.8TFigure 3.7: Specific heat for BCC lattice at c=1.O (right), and 0.5 (left).Chapter 3. Dipoles on diluted lattices 76I - I -—I I I0—0.4I I0.6CI I I0.8 1Figure 3.8: Critical temperature for the ferroelectric transition at various concentrations.The “experimental points” are from Monte Carlo simulations while the full line is themean field theory of section 3.1.1210.8T 0.6OA0.24Chapter 4Self—organized ruptures in an elastic mediumIn the previous chapters we discussed equilibrium critical behavior of a system with longrange interactions of the dipole—dipole type. In the present chapter we study dynamicalcritical behavior in an elastic medium. The stress redistribution upon a rupture is alsolong range, and in some respects resembles the long range dipolar interaction. In section4.1 we briefly introduce the concept of self—organized criticality (SOC) [7] [8]. Later wewill apply this concept to a model for earthquakes. A very short account of earthquakephenomenology is introduced in section 4.2. We propose a new method to discretize theelastic equations in section 4.3, and in section 4.4 we discuss our simulation results.4.1 Self—organized criticalityMany natural dynamical systems exhibit a power law behavior. Such a power law usuallyindicates that the system is at a”critical” point, and such systems lack a characteristiclength or time scale. This finds an analogy in equilibrium systems, at the critical pointof a second order phase transition, where we observe many power laws. To exploit theanalogy with dynamical systems, Bak, Tang and Wiesenfeld [7] [8] proposed a general77Chapter 4. Self—organized ruptures in an elastic medium 78theory— Soc — to account for the self—similar scaling in dynamical systems. The conceptof SOC can best be illustrated by a pile of sand. Let us think of building up a sand pilefrom a flat surface. Grains of sand are dropped slowly from the top. At the beginningthe sand grains stay close to where they land; soon they will be on top of each other.The local slope becomes steeper and steeper; and then somewhere in the pile grains falloff an edge causing an avalanche. The sand pile will stop growing when the local slopeapproaches a “critical” value. At this “critical point”, if more sand is added avalanchesof all sizes can occur, i.e., one can observe small avalanches as well as large avalanches;and the distribution of the size of the avalanches is a power law. The lack of a naturalscale in avalanche size has an analogy in systems at the critical point of a second orderphase transition. However, there is an important difference, the critical point in thisdynamical system is reached by the system itself (unlike in the equilibrium case whereone has to fine tune the parameters to get to the critical point). In other words, thedynamical instability will lead the system to a critical point.It turns out that this concept can be applied to many dynamical systems and canquite well explain a number of self—similar phenomena that exists in nature. There area few common features in such nonlinear dynamical systems.(i) “Energy” is fed into the dynamical systems more or less uniformly.(ii) The “energy” is dissipated intermittently, and exhibits very large events as well assmall fluctuations(iii) Spatial and temporal structures are generated at all scales.Chapter 4. Self—organized ruptures in an elastic medium 79(iv) The critical state is reached by a self—organization process— making the phenomenonrobust; no fine tunning of the parameters is involved.Bak and Tang [9] were the first to apply SOC to earthquakes in an attempt to toexplain the Gutenberg-Richter law for the moment distribution. Before we move to ourelastic model for earthquakes, let us first give a very brief introduction to earthquakephenomenology.4.2 Earthquake occurrence phenomenologyIt is now well established that most earthquakes are caused by plate tectonic movements.Earthquakes usually occur along a crack (or a ruptured surface) called a fault. Let usconsider a geological region with a fault. The fault is locked together at the beginning(i.e., the two surfaces of the fault stick together). Due to the plate tectonic movement,the stress is slowly building up in the region, and so is the elastic energy. Once the stressexceeds a threshold (the static friction of the fault), the fault slips— causing earthquakes.The energy is then suddenly released and the fault returns to a new stable point. Throughvarious reactions the fault will slowly anneal and it becomes locked again; this allows thestress to build up and an earthquake will occur again most likely in a region which hasbeen weakened by previous events— in seismology this is called a “loading cycle” and themodel we have outlined is the so called “elastic rebound” theory.From this simple point of view one can see a striking similarity between the occurrenceChapter 4. Self—organized ruptures in an elastic medium 80of earthquakes and the sand pile model of the last section. Here the “energy” is fedthrough the plate tectonic motion; and the energy is released through earthquakes. Thedynamical instability in the earth crust causes earthquakes— which in turn release theenergy and returns the earth crust to a new stable point. Therefore, the earth crust isalways operating at the boundary of instability and stability; in another words it is ata critical point. From this point of view it is not surprising that the size distribution ofearthquakes should obey a power law — the Gutenberg—Richter law [36].N(E) (4.1)Here N(E) is the number (density) of earthquakes with energy E. The exponent b isaround 0.5. There is still some debate on the extent to which b is universal and whethera different power law should be used for large and small earthquakes. Apart from theGutenberg Richter power law, earthquakes exhibit many well defined space—time patterns.Earthquakes are not uniformly distributed over the surface of the earth. Most earthquakes occur in relatively narrow belts or zones, and earthquakes occur in the same zoneagain and again. The “loading cycles” usually take 100—500 years.On a much shorter time scale, one can construct the first order temporal moment bycalculating the distribution of the time intervals between all possible pairs of earthquakesin a seismic region, without regard to their spatial location or size. This distribution fitswell over the range 2 hours to 10 years [46] by a function o 1,/f. This is Omori’s law [74]Chapter 4. Self—organized ruptures in an elastic medium 81of aftershocks. It is believed that aftershocks are caused by the stress redistribution ofthe main shock, however, the precise mechanism for the aftershocks is still controversial[73] [93]. We will come to this point later when we discuss the simulation of aftershocks.With regard to the spatial moments, Kagan and Knopoff have shown from severallocal earthquake catalogs that the distribution of epicenters appears to be fractal withdimension D 1.0-1.3 [44]. This is quite universal since the same exponent was foundin all the catalogs they had studied.Other space—time correlations can also be studied, but they are not as well supportedstatistically.From the above observational facts, one has at least 5 essential factors to considerwhen modeling the complex phenomenology of earthquake occurrence [36].(i) Plate tectonics is required to restore the energy dissipated in earthquakes(ii) A model should be able to produce a wide range of earthquake sizes.(iii) A time delay mechanism is required to provide for the time correlation betweenclustering events.(iv) Stress redistribution upon fracture is required to provide the mechanical couplingbetween earthquake events.(v) A realistic model should be able to take the complex geometry of the faults intoaccount.Chapter 4. Self—organized ruptures in an elastic medium 824.3 Discretization of the elastic equationsTo model the earthquakes adequately, one first has to incorporate the tensor characterof the stresses satisfactorily. Also, we have to take disorder in the elastic medium intoaccount. In this section we give a brief review of previous attempts of discretizing theelastic equations. We will then propose a novel plaquet representation of ruptures in anelastic medium and show that this method is very general and has certain advantagesover the previous models , and should have wide range of applications. The method willfinally be applied to modeling earthquakes.4.3.1 A brief review of previous models for discretizing elastic equationsDisorder plays a fundamental role during a fracture process. All materials are heterogeneous, both in yield strength and bulk modulus. Further fractures enhance any disorderwhich may initially have been present, through the nucleation of new cracks, or simplydue to the heterogeneity of the stress field that results from the complex geometricalarrangement of the existing cracks. Even a small amount of disorder present initially canbe be enormously amplified during fracture. This phenomena makes it extremely hardto do practical analytical calculations for the breaking process in real materials. However, computer simulations have added much insight in the understanding of the fractureprocess.Chapter 4. Self—organized ruptures in an elastic medium 83in the last ten or fifteen years, an exciting development occurred in the physics community which merged the techniques of statistical mechanics with old fracture mechanics.Terms like “percolation” and “fractals” are introduced to the shape of cracks, and arerelated to the local distribution of strain. Many approaches have been developed to dealwith local disorder directly (see, for example, reference [41]). The spirit of these studies isto describe very simply the mechanical behavior of microscopic elements, and to considerthe collective behavior of these elements in the presence of disorder. Various models havebeen proposed from this approach. Basically there are three types of methods: Moleculardynamics, finite element method and lattice models (“central—force” [38],”beam” models[86]).In the case of molecular dynamics, once the interaction between atoms has been chosen(usually a Lennard—Jones potential), one should be able to account all the elasticity,plasticity, viscosity, damage and fracture properties of the bulk material. This wouldbe a real first principle calculation. One does not need to put any rheological behaviorinto the model. However, this approach applied to fracture is so computationly intensethat only results on small systems can be obtained. Therefore no firm conclusions canbe drawn from such studies.In the finite element method, the displacement field is discretized on a lattice (usually a triangular lattice), and the total potential energy is minimized over a set of trialfunctions. The basic elements used in the fihite element method have no direct physical meaning. This method is quite popular among engineers due to its simplicity andChapter 4. Self—organized ruptures in an elastic medium 84flexibility in dealing with very different problems. However, the description of a realelastic solid will be accurate only if the displacement field varies slowly over the size ofthe elements used. If there is a crack in a medium, the displacement field around thecrack will be accurate only if the elements are much smaller that the size of the crackitself. Also, in modeling a heterogeneous solid, the size of the heterogeneities should belarger than the size of the elements. This is a serious drawback if one deals with manymicrocracks in a disordered medium.The most intensively studied lattice model is the so called “central force” latticemodel. It is a network of elastic springs which can freely rotate around the sites of thelattice. Therefore, force can only act in the direction of the bond. On a triangular latticeone chooses two spring constants [64] such that the elastic equations for the displacementfield ii (which is defined on the lattice nodes) are satisfied. This model is not independentof lattice structure. A square lattice with nearest neighbor springs, for example, doesnot support any external shear.Another alternative is the “beam” model first proposed by Roux and Guyon [86]. Inthis case the bond—bending effects are taken into account: a bond cannot rotate arounda site without a cost in energy.In the simulation, two primary steps are involved. In the first step, a bond in thenetwork is selected randomly with a probability proportional to its strain dependentbond breaking rate; or by a deterministic procedure (for example, if the strain in a bondexceeds the threshold); and this bond is removed from the system. In the second step, theChapter 4. Self—organized ruptures in an elastic medium 85elastic network is allowed to relax to a new mechanical equilibrium. Here one assumesthat the time scale for breaking a bond is longer than the time scale for relaxation.The simulation proceeds by a sequence of bond breaking and relaxation steps until thematerial has failed or a large crack pattern has developed. Most of the computer time isdevoted to the relaxation process where one needs to invert a very large sparse matrix.The “central—force” model and “beam” model are “physical”, in the sense that thebasic elements (the elastic springs and beams) have real physical meaning. One of theadvantages of these models is that one can naturally introduce local disorder either inthe breaking threshold or the spring constants (or the rigidity of the beams). However,these models have many drawbacks [87]. First of all, the “springs” and “beams” aresomewhat artificial on a mesoscopic scale. This may not be a severe problem, since itis possible that the detailed structure on the microscopic level may not be very relevantfor the bulk behavior of the material. Secondly, the computer time involved is enormousdue to the fact that at each time step one has to solve a huge array of linear equations.This is a major obstacle to progress in this area. In the next section we propose a newmethod of discretizatiou and show how some of the difficulties in the above models canbe overcome.4.3.2 A Plaquet Representation of Ruptures in an Elastic MediumIn this section we will concentrate on a detailed description of our discretization procedure and the representation of local ruptures [109]. The first step is to construct aChapter 4. Self—organized ruptures in an elastic medium 86discrete derivative. For simplicity let us consider a two-dimensional system. Generalization to three dimensions is straightforward. In our work we use the following form of thediscretized derivative (say for a fnnction g(fl):b) + g(+ d) — g(— b) — g(— d)], (4.2)and1 -1 .4 -[g(”+ b) — g(9+ d) — g(— b) + gQ— d)], (4.3)where ê are unit vectors (the lattice spacing is taken to be nnity), d = (e — ê)/2 andb = (e + ê)/2. Note that if g sits on the nodes of a square lattice, then Dg sits on thecenters of the plaquets formed by the nodes of this square lattice, and vice versa. Thisdefinition ensures that the basic unit of our discretization is a plaquet instead of a bond.Once is defined, the discretization of elastic theory is straightforward. First wedefine a displacement vector ‘ on each node (corner of a plaquet). The distortion of theplaquet is then characterized by the strain tensor (defined at the center of each plaquet)1= (DuQ) + Du(’)) (4.4)When all deformations are elastic (no ruptures exist) the stress tensor is related to thestrain tensor through the generalized Hooke’s law:= Ku1j(f)6j + 2p(u( — 6ju1i()). (4.5)where K, ji are bulk and shear moduli respectively and summation over repeated indicesis implied. The discretized elastic equation is given by the force balance condition, which,Chapter 4. Self—organized ruptures in an elastic medium 87in discrete form, is= = 0. (4.6)We further restrict our consideration to the shear mode of rupture (most ruptures inearthquakes are shear mode fractures). In this case only two stress components areneeded in order to determine where a rupture will occur: a1 = a, the shear stressat the ê, direction, and a2 = (a— a)/2 the shear stress in the directions ofor cZ In this thesis we only consider ruptures in the &, and g or ci directions. Theruptures in an arbitrary direction can be represented as a linear combination of these twocases (The shear stress in the direction which makes angle 0 with x-direction is simplya1 cos2O + a2 sin 20).Let us consider the effect of a local rupture in an infinite medium, say a shear fracturein the ê2. direction at plaquet o. After the break, we reduce the shear stress on thefractured surface to a percentage x of the original stress [93]. The stress reduction causesa force imbalance and a subsequent stress redistribution occurs. We denote the stress afterthe rupture as af” and the stress before the rupture as aj’. We have aJfm = a’j’ + a,where a is the additional stress caused by the rupture. At equilibrium afr must alsosatisfy the force-balance equation, and we must haveDa = 0. (4.7)Let u(ñ denote the additional displacement induced by the rupture. a can berelated to u via the generalized Hooke’s law throughout the system except for a5(j)Chapter 4. Self—organized ruptures in an elastic medium 88at the ruptured site. To separate the violation of the Hooke’s law at the ruptured site,we write= u(i) + (4.8)= (4.9)= a,(fl, (4.10)where u7j is the part of stress that is related to u via the generalized Hooke’s law, andrefers to the non—elastic part. Substituting these expression in (4.7) we obtain theequations for gel.+ aD(6n) 0 (4.11)Da(fl + uD1(6) = 0. (4.12)Given u we can obtain a solution for o1, which then can be used to obtain a self-consistently by applying the boundary condition for the ruptured sitegnew(j*)= xu(Po) (4.13)i.e. the stress is reduced to a fraction x of the original stress.The equation for ge can be cast (substituting in a —i./f) as—= 0, (4.14)whereFb=—6+— 64 + (4.15)Chapter 4. Self—organized ruptures in an elastic medium 89andF=— + (4.16)Therefore, can he viewed as generated by the external source F, which is onlynonzero at the corners of the fractured surface, and is shown in Figure (4.1). This forcedistribution satisfy the condition that its net force and net torque are zero, and thusis a double couple (a term used in seismology to describe the earthquake source). Ourdiscretization is in complete agreement with the well-established fact that a shear rupturecan be modeled by a double-couple force redistribution [18).Next, we want to determine oj, given this double-couple force distribution. Thesolution of this problem can be easily obtained in Fourier space. Since the side of theunit cell is unity we define the Fourier transform asg()= (2)2 j dk f2 dkg()e (4.17)The Fourier transform of the discrete derivative is then given by(Dg)(k) = 2i sin(k/2) cos(k/2)g(k) d(k)g(k), (4.18)(D9g)() = 2i sin(k/2) cos(k/2)g(k) d(k)g(). (4.19)Combining equations (4.4) and (4.5) and replacing o- with and u with u’ we have= (K — 2/3)DD1t4()-F + Du9). (4.20)More explicitly, we haveDu = [(K + 4z/3)D -I- zD]u + (K + ii/3)D.Di4, (4.21)chapter 4. Self—organized ruptures in an elastic medium 90and= [(K +41i/3)D + ttD]u, + (K +1u/3)DDu. (4.22)We perform the Fourier transform on equation (4.14) ; this leads to[(K +41t/3)d H- + (K + t/3)ddu) — vfde 0 (4.23)and[(K H- 4/3)d + d]t%() + (K + — = 0. (4.24)These are coupled linear equations can be solved straightforwardly, and the solution isu) = [(K + 4/3)d — (K — 2/3)dd]e0, (4.25)= [(K + 4/3)d — (K — 2/3)dxd]e_i0, (4.26)where () = (K + )(d + d)2. Now we can get the stress tensor u) using thegeneralized Hooke’s law:= (K— 2/3)(du + di4) + 2tdu (4.27)4dd3 -— f 4 28— (d--d)2 ‘= (K — 2/3)(dz4 + dt4) + 2idu (4.29)4d3d •-= (d Hd)2e (4.30)Chapter 4. Self—organized ruptures in an elastic medium 91and= (4.31)4dde0 + fe0, (4.32)=where f = fv’(K + /3)/(K + 4u/3). We only monitor the shear stress o and Theadditional shear stresses are given by= () = + (4.33)= -=- u)). (4.34)Using the fact that o—v”f we obtain4d2d2y —kr (4.35)= f(d+d)2eand9-’ (d2— ) (4.36)(d+d)2 eIn real space we have(4.37)ando) —fG2( ), (4.38)where G1 and C2 are lattice Green’s functions:r dk r dk9 sin2 k sin2 ke; (4.39)G1(r)= J_r27rJ_r27r(1_coskxcosky)2=dk. rj— 2 J_ 2 (1 — cos k cos e (4.40)Chapter 4. Self—organized ruptures in an elastic medium 92If we assume the boundary condition at the ruptured site, that after the rupture theshear stress a1 on the fractured segment is reduced to a percentage x of the originalvalue, then.7 can be determined by a0 + a(Po) = xa0, where a0 is the original shearstress just before the rupture. This leads to—fGi(O) + (1 — x)ao = 0, and a, a can bewritten as,al,2(r) = —(1 — x)aoGi,2(i—0)/G1(0). (4.41)We have checked that the Green functions reduce to correct continuum limit for thedouble-couple stress redistribution. At long distances the Green functions decay as 1/r”,where d is the spatial dimension. At short distance the greatest stress increase is nearthe tip of a crack.For a fracture in the 5 or d direction the double-couple force distribution is shown inplaquet (b) of figure (4.1), and is given below,= f[—+ — &6+ io+ (4.42)We can obtain similarly the additional stress a2 due to a rupture in the direction of bor ci,=—(1— x)aoG2(P—i0)/(l— Gi(0)); (4.43)= -(1- x)ao(690 - G1(P- b))/(1 -G1(0)). (4.44)Again, a0 is the original shear stress a2 at the ruptured site just before the rupture.We can also calculate the stress redistribution when more than one site has ruptured.For simplicity, let us consider the case that all ruptures occur in ê. direction. Let , 1 =Chapter 4. Self—organized ruptures in an elastic medium 931,. , N, (Nb is the total number of rnptures) be the locations of ruptnres. Again we canrepresent the ruptures by double-couple force distribution. We assume the correspondingstrength of the double couples be fi, 1 = 1,. , N6. Because of the crack-crack interactionwe can not determine fi the same way as in the case of a single rupture. Instead, theboundary condition r”() = xaI1”(Jj) corresponds to a set of linear equations whichhave to be solved to obtain fi:(1— + EfmGiVi—m) — = 0, (4.45)for 1 = 1, .. , N6. After fi has been obtained we can determine the additional stresscaused by the ruptures in the entire medium.As a test study we consider the stress redistribution caused by a segment of consecutive ruptures in the ê direction (with x=0 corresponding to complete stress release). Thecorresponding continnum case a segment in the interior of the medium is ruptured andthe shear stress on it reduced to zero — can be solved analytically [57]. In figure(4.2) weplot the additional stress cr generated by the ruptures from the procedure we describeabove. The agreement with analytical results (see, for example [57]) is excellent. Thisagain gives us confidence that our discretization is a very good representation of theelastic theory.Our scheme of discretization has the following advantages:(i) Each plaquet can contain many micro cracks and is characterized by the strength ofthe plaquet. The length scale is set by the size of the plaquets. This should make itChapter 4. Self—organized ruptures in an elastic medium 94possible to study engineering type problems e.g. crack growth in ceramics, as well asgeological scale problems such as earthquakes.(ii) We are dealing directly with the stress tensor distribution, so the simulation processis explicit in its physical interpretation. This makes the method extremely flexible; andthe modification to incorporate different physical situations is straight forward. We willillustrate this point when we apply this method to earthquakes.(iii) In the case of open boundary conditions (the infinite medium), the calculation of thestress redistribution upon a rupture is straight forward as we have shown in this section.In contrast the “central—force” and “beam” models where the calculation of the stressredistribution requires solving a huge array of linear equations which is extremely timeconsuming.4.4 Applications to EarthquakesEarthquakes can be viewed as sequences of ruptures which release the elastic stress builtup gradually by the movement of tectonic plates. The global dynamics of earthquakesis quite rich and complex. There have been many attempts at earthquake modeling,yet there are still many unsolved problems. Part of the reason is the large amount ofphenomenology needed to model; a more important reason is that one does not know howto model the physics well. One of the outstanding problems is to incorporate the tensorcharacter of the stresses satisfactorily into a dynamical model of earthquakes [56]. Thisis the problem we attempt to solve in our model of earthquakes. We also show that a fewChapter 4. Self—organized ruptures in an elastic medium 95empirical laws regarding healing of fractured surfaces [24] (for ensuring a steady state) andstatic fatigue of rocks [93] [6] (responsible for generating aftershocks) can be incorporatednaturally in our model. In addition to reproducing a number of phenomenological lawssuch as the Gutenberg-Richter law [36] (as in a number of previous earthquake models[17] [88] [9] [10] [20] [97] [21]), our model can also dynamically generate a relatively stablefault structure, made up of regions that have been weakened by earthquakes. We willdiscuss how the annealing law affects spatial and temporal structures of faults and variousphenomenological laws of earthquakes in our models. As we will show below this modelis very flexible, we will just concentrate on showing how to implement the scheme andgive some examples.Before going into the technical details of the model, let us first describe an overviewof the model, which is illustrated in figure(4.1). We are focusing on the entire (potential)earthquake region rather than on individual pre-existing faults. For simplicity we onlyemploy a two-dimensional model. The system is divided into L x L plaquets. Due the theplate tectonic movement, the system is driven by uniformly increasing the shear stress ata constant rate. At the beginning the strength (or the yield threshold) of each plaquet israndomly assigned. Since the shear stress is ever increasing, at certain time the stress at aplaquet will exceed its strength, and this site will be ruptured. The stress in this plaquetis reduced. This will cause the force imbalance and the stress in the system has to beredistributed. The redistribution takes place in the speed of sound. At the positions withincreased stress, the stress may now exceed the critical value causing further ruptures;Chapter 4. Self—organized ruptures in an elastic medium 96thus the local instability may cause a chain reaction. Our model does not describe thedetails of the dynamical process of ruptures (we treat the velocity of sound as infinitelylarge), but only the stress distribution before and after the rupture. An earthquake isdefined as a chain reaction in this model.We need to choose a boundary condition as we only can simulate a finite region.In earthquakes it makes sense to use open boundary conditions: we consider a infinitemedium, but neglect the ruptures outside the region on which we concentrate. This hasan enormous advantage computationally, because we can use the Green’s functions forthe infinite medium as described in section 4.3.24.4.1 The effect of weakening and annealingFor simplicity, we only consider in this thesis the case where shear stress in a2 is increaseduniformly by the amount p at each time step, with zero initial stress. After each rupture,the fractured plaquet is weakened: the stress threshold is reduced to a lower value 6(Initially we assign the critical stresses to be random numbers between 61 and 6). In realitythe threshold for the shear failure of a previously fractured surface is better describedby 6= 0? + van, where a, is the normal stress, 6 represents cohesion in the material,and i’ is defined as the internal friction coefficient. As a first step, we only consider thespecial case z-’ = 0 in this thesis.The weakening of fractured surfaces is responsible for generating a fault zone. In orderfor a steady state to be established, the fractured surface is allowed to anneal slowly. TheChapter 4. Self—organized ruptures in an elastic medium 97annealing process is not well understood. In this thesis, we use the empirical finding [24]that the strength of the fractured surface will grow as= O + AlogQ — t) (4.46)where t denotes the time of the last fracture and 01 is the lower cutoff of the threshold(to indicate the fact that after the rupture, the fractured surfaces are most likely stillin contact, parts are not completely ruptured, and the fractured surface will be able tosupport a small applied stress).Occasionally, the stress on one of the plaqnets will exceed its threshold and it willrupture. The stress redistribution from this rupture is computed and we check if thiscauses the stress on other plaquets to exceed their thresholds. This produces a list ofplaquets which will rupture in the next step. In the simulations presented here the stressredistribution from these plaquets are simply calculated individually and summed, ratherthan using (4.45). This process is the same as was used in [21]. It could be mentionedthat, although the time scale for breaking and stress redistribution, are fast, these processes are not truly instantaneous. One can conceive of different selection schemes forthe order in which the plaquets break. Our scheme in carrying out the simulations isonly an approximation to a real dynamical process. However, we believe that as long asthe force balance equations are satisfied, these details should not affect the power lawsobserved here. The overall critical behavior is mainly determined by the symmetry (alsothe range of interactions) of the system.Chapter 4. Self—organized ruptures in an elastic medium 98In the calculations presented here we restrict ourselves to the case where the stressis increased for the component a2. As discussed in [108] a large fraction of the eventstake place near the edges, if instead we uniformly increase the component a3,, in a squareregion with sides parallel to x and y. We only consider the shear mode of rupture andbegin by considering the model without static fatigue.We measure the size of an earthquake by the total number s of sites that have rupturedfollowing the initial instability. If a given site ruptures more than once during an eventevery occurrence is counted. We also define the “seismic moment” which is proportionalto the sum over sites of the stress drop. We have calculated the distribution N(s) andfound that it is a power law N(s) cx 1/si. The exponent r is approximately 1.4 (seeFigure (4.3 and Figure (4.4) ). We have also investigated the distribution of the seismicmoment as defined above and found it to follow the same power law as N(s). We havefound that this power law is quite robust. The exponent does not seem to depend on thevalue of A in (4.46). Similarly, the exponent is unchanged in the case of instantaneousannealing (i.e, the threshold after the rupture is the same as that before the rupture).Also it doesn’t seem to matter if we put 01 to be constant or random after rupture aslong as the initial value of O is random. The exponent for N(s) agrees with that foundin [21] within computational errors.Figure (4.5) shows some typical large and a medium size events. The pattern ofruptured sites associated with a single event does not depend a great deal on the annealinglaw, but is largely determined by the direction of stress enhancement after rupture whichChapter 4. Self—organized ruptures in an elastic medium 99is mostly at the crack tips.While some features turn out to be insensitive to the details of the model the spatialpattern of where the earthquakes occur depends on the annealing law (4.46). We showin Figure (4.6 — 4.9) the location of epicenters and the pattern of cumulative slip after500 events in both the cases of instantaneous and weaking and slow annealing. It is seenthat in the case of instantaneous annealing the quakes will take place over the wholeregion. On the other hand with A 0 and O small stable faults seem to develop. Theearthquakes in the simulation are clearly clustered, which is also one of the characteristicsof real earthquakes [56). In reality the fault systems evolves over very long geologicaltime scale. Some new faults will be generated and some old faults cease to be active.Our model is consistent with this observation.4.4.2 Simulation of aftershocksThe model which we have described so far only has one time scale: the slow geological timescale which determines the accumulation of shear stress (there is no time scale associatedwith the logarithmic annealing process). In a real earthquake sequence, however, thereappears to be a second and a smaller time scale. Although the main shock lasts onlyfor about a minute, the entire earthquake sequence combining foreshocks, mainshock,and aftershocks can last for many weeks. The aftershock sequences are found to obeyOmori’s law: the rate of aftershocks following the main shock decays as 1/Va’ as a functionof time, with the exponent p close to one. The most popular explanation of Omori’s lawChapter 4. Self—organized ruptures in an elastic medium 100is static fatigue ([91], [6]): the rock strength is time-dependent. When the applied stressa is below the instantaneous breaking strength but is above the stress-corrosion limit 00,rocks generally break with probability per unit time proportional to e, where a is aconstant.To illustrate this general physical picture and to verify Omori’s law, we have modifiedour model to incorporate the physics of static fatigue. Consider a specific plaquet at theposition i in our model. If the shear stress a(i) is above the instantaneous breakingthreshold of the plaquet Oj(io), the plaquet will fracture as described above. However,when the stress is below 0 but is still above 0, the plaquet is set to break with probabilityper unit time given by ea(oi).Now, an earthquake sequence can be defined as follows. When all local stresses arebelow the respective instantaneous breaking thresholds as well as the stress corrosionlimit Oo, nothing will happen; this is in the period between earthquake sequences, whichwe call the rest state. Whenever there is a place where the stress is above the localinstantaneous breaking threshold or the stress corrosion limit 0, ruptures in the systemsare expected, and an earthquake sequence begins. The return of the system to the reststate signals the end of this earthquake sequence. Each earthquake sequence can containmany shocks, separated by periods of no activity with some sites remaining to be fracturedby static fatigue. The largest shock in a sequence is defined to be the main shock; theshocks before the main shock are foreshocks and the shocks after it are aftershocks. Tocheck Omori’s law, we calculate the number of aftershocks following the main shock as aChapter 4. Sell—organized ruptures in an elastic medium 101function of time, and average over many earthquake sequences. The result is plotted inFigure (4.10). The linearity of the log-log plot indicates the Omori law N(t) = A/V2 withp 0.90. The spatial pattern associated with a typical series of events is shown in Figure(4.11). Scholz [92] found that p=l can be obtained by assuming randomly distributedstress levels induced by the mainshock; this can be viewed as a mean field result in whichthe spatial correlations are neglected.The origin for aftershocks is controversial. An alternative explanation of aftershocks isthat the main shock induces fluid diffusion which in turns causes changes in the effectivestress and subsequent fractures [73]. In this picture aftershocks migrate away from themainshock rupture boundary following a fluid diffusion front [94]. In contrast, aftershocksin our model, which is induced by static fatigue, occur on the average all over the rupturezone after the mainshock; this is usually observed in real earthquakes [93].chapter 4. Self—organized ruptures in an elastic medium 102fff(a)(b)--1Figure 4.1: Schematic illustration of a region subjected to external shear stress. Thesystem is divided into many plaquets, and is driven by slowly increasing sheer stress.When a local stress is larger than the corresponding threshold stress, the plaquet fractures, causing a long-range redistribution of elastic forces. The force distribution of localdouble-couple sources are also illustrated: (a) Fracture in è. direction; (1)) Fracture inê + ê direction.Chapter 4. Self—organized ruptures in an elastic medium 103>-iFigure 4.2: The shear stress o redistribution calculated form the Green’s functionmethod in the text. The length of the fault at the center is 100 lattice units. The shearstress is reduced by 100 units uniformly along the fault. This distribution is in goodagreement with the analytical results (see, for example, figure 8 in reference [57])xChapter 4. Self—organized ruptures in an elastic mediumloop100N(s)10110S100104Figure 4.3: Size distribution of earthquakes generated in the stationary state of themodel for a 100 x 100 lattice in the case of instantaneous annealing. 10000 earthquakesare included to obtain the statistics. The linear behavior of the log-log plot indicatesa Gutenberg-Richter law with an exponent approximatedly 1.41. The parameters arep = 1.0 x iO, = 0.25, x = 0.5 (see formula (4.41) ), and = 1.0.:‘ ‘1 1 1 tIl1lI-1-[-II-LIH I:IIjj f____j_ 1_f_li1f I till I—Chapter 4. Self—organized ruptures in an elastic medium 105liii •I I 1111111 I I Illilif I I1000\\\“\\\\\100N(s)10 -1 -;I I IIIIII I I 11111,1 I t1 10 100SFigure 4.4: Size distribution of earthquakes generated in the stationary state of themodel for a 100 x 100 lattice in the case of weaking and annealing. 10000 earthquakesare included to obtain the statistics. The linear behavior of the log-log plot indicatesa Gutenberg-Richter law with an exponent approximatedly 1.38. The parameters arep = 1.0 x 10, O = 0.25, x = 0.5 (see formula (4.41) ) and 8 = 1.0. The annealingconstant A=i0 (see formula. (4.46) ).o 20 40 60 80 100J liii l1ILII.7’— ifiIii,ltii!iilIllr10 20 40 60 80 100J 111151 I—IIIIIILIIIIICaaIII IIIIIt Slit S106Chapter 4. Self—organized ruptures in an elastic mediumJ t14I511’(’’l’’ L—••1I liii ui, ,IIi’’’i’ ‘I’’’1511’ L.: \:..Sa1 i I I i t I i I i i I1008060402001008060402001008060402001008060402000 20 40 60 80 100 0 20 40 60 80 100Figure 4.5: Some typical spatial patterns of simulated earthquakes.Chapter 4. Self—organized ruptures in an elastic medium 107100I I I I I I I I I I I I I I I I I‘° 8O 00 0°8 8 000 0 000 0%600o9d’o 0 8• 0 8 0 0 000 000 Q 0• 00 0 00 0 0 00 0 00 Q0 080- 000 0 o°o0 00000 000 800 0 00 bo° 0 0 000 000OQ 00 0 d o 0 000 0 0 0060- 8 0 o 0 0000 0 000000 0 80•00O: 00 0 00 000°b 0 °0Q 0 oQ 0008 0000Qo040-o 0 00 0 0 0• 00 0 0 00ö’ 0 0 00 0 0-0 0 Qo000• 0 0 0 0 o 0OcP 0a0000oQ 0 00 0020_s 00 0 0 0 000 o0d’ 0c° o°oo’0 oo800 00D000 0 0 0 00000 0 Q 0 0 000000OO 0 00oQ°co 0 0 0 00I I r i0 20 40 60 80 100Figure 4.6: The epicenters of 500 quakes in the case of instantaneous annealing. The sizeof the circles is proportional to the size of quakes. The parameters used in the simulationare the same as in figure (4.3).Chapter 4. Self—organized ruptures in an elastic medium 108100 —::.. ..‘ I .:. ‘ 1. .. : ._1 .. :. L: •.I :: ‘..• : ... • ..•••1 S 555.‘ .: : •5 . •..... •..:.I: :., • 55• •5 55 5 5 54 •0 5.4 •• S 55 •4 S I,, j:..; .;..:;i. . :. ;,:.:( ;:. ,:. • •.•, ••• .•,0 20 40 60 80 100Figure 4.7: The accumulated slip during 500 quakes in the case of instantaneous annealing. The scaled length of the line segment represents the amount of slip at each site. Theparameters used in the simulation are the same as in figure (4.3).Chapter 4. Self—organized ruptures in an elastic medium 109Figure 4.8: The epicenters of 500 quakes in the case of weaking and slow annealing.The size of the circles is proportional to the size of quakes. The parameters used in thesimulation are the same as in figure (4.4).I Id, o000 o090 8 00000000 000o0 000 00 0 00 ©o000 o0 000000 000 00001.00 T0 101 I I I I o co’0 0 CO 00d 8o0 00 00- 0 0 00 (Q0 00 do0 900 0 0 0)800 oq,080 000 000 00000(0 0 00000 0 000O0 00o000 060- 00 0-0o000 0000 000 80% 0 0080°‘kb 00:00:00040-a 0 00000o00 0o0 000 0 000 cP 000000 b1:, 0000800 00 o00 000020-0 0 0 o 0 0 0 -00 0 000°c,° 0000•0 0 0000 90o 0000 0 80• 00 0 0000 0) 0 0 0- a o 0 0 0 cB 0 0 00I I I , Io 80 1000 0000o020 40 60chapter 4. Self—organized ruptures in an elastic medium 110I II II 44 I II II 4I.ItI.—. I• I--I ‘• .—. I• 4 I4 4I II. II II —I 4.Figure 4.9: The accumulated slip during 300 quakes in the case of weaking and slowannealing. The scaled length of the line segment represent the amount of slip at eachsite. The parameters used in the simulation are the same as in figure (4.4).10080I I .•1 ‘.••.— .! .‘• I..4. I.I. £•• I.4. I60 -4 I.4I. •4 I .4.4 4 4I.•..••40L :1• . . II I4 • 4. • III• 4 4 4 4 II.,• • 4 4 II I..4. I— ..4 .•• I. I III 4 I I4. II I I4 l1 I I .44. II• 4 4 II .IIII I I4. II I I41 I I.4 • II I I4 II I I• 4 II I• 4I4II II..4• II .4_I I.II •.• •II I •• 4I £• .• .•I •• •20..;-..:.;: .• .I :0 . II0 20 40I 160 80 100Chapter 4. Self—organized ruptures in an elastic medium 11110Figure 4.10: Omori’s law: the log-log plots of the rate of aftershock; as functions of timefollowing the mainshock. x = 0.5, 8 = 1.0, 8o = 0.5, O = 0.25, arid 15.0 were usedin the simulation of a 100 x 100 system.Chapter 4. Self—organized ruptures in an elastic medium 112I I I I I I I I1000c9aSa.as a80 / a a.a•5a..Uaaa00\a.60 -a./&a,20 % (900000I— I I t I 1 1 i 1 .1 —0 20 40 60 80 100Figure 4.11: Spatial pattern of a main event (full square) and after shocks (open circles).The parameters are the same as in figure (4.4).Chapter 5ConclusionsIn this thesis we have studied the behavior of two seemingly very disparate systemsthat originate from different areas of science, but which have certain common features.Namely, their long range interaction nature and the critical behavior governed by it. Bothanalytical work and extensive numerical simulation technique have been used throughoutthe thesis.In the first three chapters, the ordering properties of Ising dipoles are investigated.The boundary conditions are such that there is no net de—polarizing field and both regular lattices and various random arrangements are considered. We construct a mean fieldtheory employing the replica solution for the spin glass, with a Gaussian approximationfor the distribution of dipole—dipole interactions while a Kirkwood approximation is usedfor the spatial distribution of dipoles. This mean field theory is manifestly shape independent in the case of no external field, in accordance with the proof of R.. Griffith [34].The main conclusion from the mean field is that there exits a critical dipole density concentration above which the low temperature phase is ferromagnetic (or antiferromagneticdepending on the underlining lattice structure). Below this critical concentration there113Chapter 5. Conclusions 114is only a spin glass phase.Extensive Monte Carlo simulations have been performed for the diluted BCC lattice.To find the critical concentration, we numerically ( with a combination of simulatedannealing and steepest descent) search for the “ground” state and calculate the netmagnetization for the ground state. The numerical finding is in reasonable agreementwith mean field theory. We have also investigated some properties of the “ground” statein the spin glass phase and found a local field distribution P(E) cx EY for low fields withcv 0.37. This value is much smaller than the value cv 1 for the SK model. Thisindicates that the low temperature diluted dipolar system is less “frozen” than the spinglass with infinite range interactions.Above the critical concentration, the ferromagnetic transition temperature is locatedby finite size scaling. As could be expected, the transition temperature is higher ( 30%)than found in the simulations. But the concentration dependence (which is linear in thedipolar concentration), is in qualitative agreement with the simulation data.In the case of Ising dipoles on pure lattice systems, we investigated the effects oflogarithmic corrections to the critical behavior. In this study we choose the structure ofthe compound LiHoF4— body centered tetragonal lattice. Two major conclusions can bedrawn from the simulation: (1) The finite size scaling form that includes logarithmic corrections significantly improve the data collapse over the mean field critical behavior. (2)With the lattice parameters appropriate to LiHoF4,the transition temperature we foundin simulation is in excellent agreement with the experiment. Thus strongly indicates thatChapter 5. Conclusions 115in this compound dipole—dipole interaction is the dominant one.In the last part of the thesis we have derived in detail a novel discretization schemeto represent ruptures in elastic media. Local ruptures are represented by double couples.The method has the advantage of being fast and easy to implement, and it is quiteaccurate in cases where we have been able to compare with analytical results. The methodcan easily incorporate local disorder in yield strength and can deal with phenomenologicalcorrections associated with static fatigue and annealing. We have illustrated the methodon a two dimensional square lattice with open boundary conditions. The boundaryconditions can easily be modified to the case where the active region has an arbitraryshape, but the use of other “non open” conditions would complicate the calculationby destroying the translational invariance of the Green’s functions. Similar difficultiesarrive if we attempt to deal with a medium with a heterogeneous elastic modulus. Notethat our method of discretizing the elastic equations can be easily generalized to threedimensions; one only needs to redefine the operator D and D in (4.3) and (4.4) andsimilarly define D2, all other calculation procedures are exactly the same as in the twodimensional case. However, the computer simulations will be much more time consumingin three dimensions.When applied to earthquakes, the model will, just as several other models, reproducesome of the size and temporal scaling laws. The Gutenberg--Richter law exponent isfound to be about 0.4, and we find this result is insensitive to the detailed parametersin the model. The Omori’s law for the aftershocks can be successfully reproduced in theChapter 5. Conclusions 116model by introducing static fatigue for the rocks. The medium is weakened locally byrupture leading to the formation of faults. Over a very long time scale a self—organizedstate develops through the gradual annealing of the faults. Over a much shorter timescale aftershock sequences are caused by local static fatigue. With this model there isalso a possibility of studying spatial patterns. Since the model is rather flexible we haveconcentrated on showing how to implement it and giving some examples. We thereforehave by no means been able to explore the full parameter space.Bibliography[1] D J Adams and I R McDonald, Molecular Physics 32 (1976), 931.[2] A Aharony, Phys. Rev. B 8 (1973), 3363.[3] A Aharony, Phys. Rev. B 13 (1976), 2096.[4] G Ahiers, A Korblit and H J Guggenheim, Phys. Rev. Lett. 34 (1975), 1227.[5] A Arrott, Phys. Rev. Lett. 20 (1968), 1029.[6] B. K. Atkinson, J. Geophys. Res., 89 (1984), 4077.[7] P Bak, C Tang and K Wiesenfeld, Phys.Rev. Lett. 59 (1987), 381.[8] P Bak, C Tang and K Wiesenfeld, Phys. Rev. A38 , (1988), 364.[9] P Bak and C Tang, J. Geophys. Res. B94 (1989), 15635.[10] P. Bak and K. Chen, 1990, in Fractals and their application to geology, edited by C.Barton, and P. LaPointe, Geological Society of America, Denver, 1991.[11] K Binder, in Finite size scaling and numerical simulations of statistical systems, edPrivman, , Singapore, World Scientific, 1990.[12] K Binder and A P Young, Rev. Mod. Phys. 58 (1986), 801.[13] P Beauvillain, J-P Renard, I Laursen and P J Walker, Phys. Rev. B 18 (1978),3360.[14] R N Bhatt and A P Young, Phys. Rev. Lett. 58 (1985), 924.[15] R Botet and P Pfeuty, Phys. Rev. Let. 49 (1982). 478.[16] E Brezin arid J Zinn-Justin, Nuclear Physics B257 (1985), 867.[17] R. Burridge and L. Knopoff, Bull. Seismol. Soc. Am.. 57 (1967),341.[18] R. Burridge and L. Knopoff, Bull. Seismol. Soc. Am., 54 (1964), 1875.[19] J M Carlson and J S Langer, Phys. Rev. Lett. 62 (1989), 2632.117Bibliography 118[20] J. M. Carison, J. S. Langer, B.E. Shaw and C. Tang, Phys. Rev. A 44 (1991), 884.[21] K Chen, P Bak and S P Obukhov, Phys. Rev. A 43 (1991), 625.[22] P Debye, Phys. Z. 13 (1912), 97.[23] J R L de Almeida and D J Thouless, J. Phys. A:Math. Gen. 11 (1978), 983.[24] J. H. Dieterich, J.Geophys. Res. 77 (1972), 3690 bf ibid. 3771.[25] For a review of the fracture properties of rocks as well as earthquake mechanics ingeneral see the excellent book by C. H. Scholz, The mechanics of earthquakes andfaulting, Cambridge University Press, Cambridge (1990).[26] S F Edwards and P W Anderson, J. Phys. F: Metal Phys. 5 (1975), 965.[27] A M Ferrenberg and R H Swendsen, Phys. Rev. Lett. 61 (1988), 2635.[28] J J Folkins and J A Griffin, Phys. Rev. B 25 (1982), 405.[29] U T Fiory, Phys. Rev. B 4 (1970), 614.[30] M E Fisher, in Critical Phenomena, Proc. 51st Enrico Fermi Summer School, Varenna,ed M S Green, Academic Press N Y (1972).[31] Y Fu and P W Anderson, J. Phys. A: Math. Gen. 19 (1986), 1605.[32] J A Griffin, M Huster and R. Foiweiler, Phys. Rev. B 22 (1980). 4370.[33] K Gotoh, W S Jodrey and E M Tory, Powder Technol. 21 (1978), 285.[34] R Griffiths, Phys. Rev. 176 (1968), 6.55.[35] M Gulcsi and Zs Guhicsi, Phys. Rev. B 33 (1986), 3483.[36] B. Gutenberg and C. F. Richter, Ann di Geofis. 9 (1956), 1.[37] A J Guttrnann, J. Phys. A: Math Gen. 9 (1978), L103.[38] see for example, A. Hansen, S. Roux and H. J. Herrmann, J. Phys. France 50 (1989),733.[39] P E Hansen, T Johansson and R Nevald, Phys. Rev. B 12 (1975), 5315.[40] P.C. Hemmer and J.L. Lebowitz in Phase Transitions and Critical Phenomena, ed.C. Domb and M.S. Green 5B, (1976), 108.Bibliography 119[41] A collection of reviews: Statistical models for the fracture of disordered media, ed.H. J. Herrmann and S. Roux, North—Holland, Amsterdam, 1990.[42] U T Hochli and M Maglione, J. Phys.: Condens. Matter 1 (1989), 2241.[43] V M Jansoone, Chem. Phys. 3 (1974), 78.[44] Y. Y. Kagan and L. Knopoff, Geophys. J. R. Astron. Soc. 62 (1980), 303.[45] Y. Y. Kagan, preprint.[46] Y.Y. Kagan and L. Knopoff, Geophys. J. Roy. Astron. Soc. 55 (1978), 67.[47] W Kanzig, H R Hart, S Roberts, Phys. Rev. Lett. 13 (1964), 543.[48] C Keller and H Schmutz, J. Inorg. Nucl. Chem. 27 (1965), 900.[49] S Kirkptrick and S M Varma, Solid State Gommun. 25 (1978), 821.[50] S Kirkptrick and D Sherrington, Phys. Rev. B 17 (1978), 4384.[51] C Kittel, Rev. Mod. Phys. 21 (1949), 541.[52] W Kleemann, S Kntz and D Rytz, Europhys. Lett. 4 (1987), 239.[53] M W Klein C Held and £ Zuroff, Phys. Rev. B 13 (1976), 3576.[54] M W Klein, C Held and E Zuroff, Phys. Rev. B 19 (1979), 1492.[55] S J Knak Jensen and K Kjaer, J. Phys.: Gondens. Matter 1 (1989), 2361.[56] L. Knopoff, in Disorder and Fracture, edited by J. C. Charmet, Plenum Press, NewYork, 1990.[57] B.V. Kostrov and S. Das, Bull. Seism. Soc. Am. 72 (1982), 679.[58] R Kretschmer and K Binder, Z. Phys. B 34 (1974), 375.[59] P-Y Lai and K K Mon, Phys. Rev. B 41 (1990), 9257.[60] Lin Lei (Lui Lam), Mol. Gryst. Liq. Gryst. 146 (1987), 71.[61] A I Larkin and D £ Khmel’nitskii, Soviet Physics JETP 29 (1969), 1123.[62] W N Lawless, Phys. Rev. B 23 (1981), 2421.[63] R Liebmann, B Schaub and H C Schuster, Z. Physik. B 37 (1980), 69.Bibliography 120[64] E Louis and F Guinea, Europhys. Lett. 3 (1987), 871.[65] J M Luttinger and L Tisza, Phys. Rev. 70 (1946), 954.[66] S—K Ma, Phys. Rev. B 22 (1980), 4484.[67] G D Mahan, Phys. Rev. 153 (1967), 983.[68] 5 McKenzie, J. Phys. A: Math Gen. 12 (1979), L.53.[69] M Medina, R Nava and M Saint Paul, Solid State Commun. 50 (1984), 51.[70] M Mezard, 0 Parisi and M—A Virasoro, Spin Glass theory and beyond , Singapore,World Scientific, 1987.[71] 0 G Mouritsen and Knak Jensen, Phys. Rev. B 19 (1979), 3663.[72] 0 G Mouritsen and Knak Jensen, J. Phys. A 12 (1979), L339.[73] A. Nur and J. R. Booker, Science 175, (1972), 885.[74] F. Omori, J. Coil. Sci. Imper. Univ. Tokyo, 7 (1894), 111.[75] L Onsager, J. Am. Chem. Soc. 58 (1936), 1486.[76] P. Palify-Muhoray, M.A. Lee and R.G. Petschek, Phys. Rev. Lett. 60 (1988), 2303.[77] R G Palmer and C M Pond, J. Phys. F:Met. Phys. 9 (1979), 1451.[78] G Parisi, em J. Phys. A: Math. Gen. 13 (1980), L115.[79] G Parisi, em J. Phys. A: Math. Gen. 13 (1980), 1101.[80] 0 Parisi, em J. Phys. A: Math. Gen. 13 (1980), 1887.[81] J Pirenne, Hey. Phys. Acta. 22 (1949), 479.[82] V Privman, Finite size scaling and numerical simulation of statistical systems, Singapore, World Scientific, 1990.[83] Z. Rácz, H-J Xu and B. Bergersen, to be published.[84] D H Reich, B Ellman, J Yang, T F Rosenbaum, 0 Aeppli and D P Belanger Phys.Review. B 34 (1986), 4956.[85] 5. Romano, Nuovo Gim. 7D, (1990) 717.[86] 5 Roux and E Guyon, J. Physique Lett. 46, (1985) L999.Bibliography 121[87] S Roux, in “Statistical models for the fracture of disordered media”, North—Hollanded by H J Herrmann and S Roux, 1990.[88] J.B. Rundle and D.D. Jackson, Bull Seism. Soc. Am. 67 (1977) 1363.[89] M Saint—Paul and J. Gilchrist, J. Phys. C. 19 (1986), 2091.[90] J A Sauer, Phys. Rev. 57 (1940), 142.[91] C. H. Scholz, J. Geophys. Res., 77 (1972), 2104.[92] C. H. Scholz, Bull. Seismol. Soc. Am., 58 (1958), 1117.[93] For a review of the fracture properties of rocks as well as earthquake mechanics ingeneral see the excellent book by C. H. Scholz, The mechanics of earthquakes andfaulting, Cambridge University Press, Cambridge, 1990.[94] Recently Bruce Shaw (ITP preprint: NSF-ITP-91-94) has studied a dynamical modelbased on the Burridge and Knopoff model coupled with a diffusion field. His modelalso gives rise to the Omori law.[95] A. Sornette and D. Sornette, Europhys. Lett. , 9, (1989) 197.[96] Y Takagi, Phys. Rev. 139 (1952), 1010.[97] Chao Tang, to appear in Modeling Complex Phenomena, edited by L. Lam and V.Naroditsky.[98] J H van Vleck, J. them. Phys. 5 (1937), 320.[99] P.N. Vorontsov-Velyaminov and l.A. Favorskii, Soy. Phys. Sol. State 15 (1971),1937.[100] B E Vugmeister and M D Glinchuk, Rev. Mod. Phys. 62 (1990), 993.[101] B E Vugmeister and M D Glinchuk, Soy. Phys. Usp. 28 (1985), 589.[102] D. Wei and G.N. Patey, Phys. Rev. Lett. 68 (1992), 2043.[103] W Wu et. al., Phys.Rev. Lett. 67 (1991), 2076.[104] H-J Xii, B. Bergersen and Z. Rcz, submitted for publication.[105] H-J Xu, B. Bergersen, F. Niedermayer and Z. Rácz, J. Phys.: Condens. Matter 3(1991), 4999.[106] H-J Xu, B. Bergersen and Z. Rácz, J. Phys.: Condens. Matter 4 (1992),2035.Bibliography 122[107] H-J Xu, unpublished.[108] H-J. Xu, B. Bergersen and K. Chen, to appear in J. Phys. A. Letter to the Editor.[109] H-J. Xu, B. Bergersen and K. Chen, submitted for publication.[110] A P Young, in Finite size scaling and numerical simulation of statistical systems,ed V Privman , Singpaore, World Scientific, 1990.[111] XV Zernik, Phys. Rev. 139A (1965), 1010.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Computer simulations of critical phenomena in systems...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Computer simulations of critical phenomena in systems with long range interaction: A study of ising dipoles… Xu, Huangjian 1992
pdf
Page Metadata
Item Metadata
Title | Computer simulations of critical phenomena in systems with long range interaction: A study of ising dipoles and self-organized criticality in earthquakes |
Creator |
Xu, Huangjian |
Date Issued | 1992 |
Description | This thesis discusses scaling and critical behavior of two different models. One model describes Ising dipoles, originates in condensed matter physics and depicts equilibrium critical phenomena. The other model, taken from the earth sciences, describes faulting instabilities and the resulting earthqnakes, and involves self-organized criticality — a nonequilibrium phenomenon. Both models are characterized by long range interactions, with a resulting sensitivity to boundary conditions. The ordering properties of Ising dipoles on lattices are studied in a mean field theory and by Monte Carlo simulations. The mean field theory is manifestly shape independent in zero external field. In the case of dipoles on a diluted lattice the mean field theory predicts a critical concentration above which the low temperature phase is ferroelectric (or anti-ferroelectric depending on the lattice structure). Extensive Monte Carlo simulation results are in agreement with those of mean field theory. We propose a finite size scaling form that includes logarithmic corrections for systems at the critical dimensionality. In the case of dipoles on a body centered tetragonal lattice we found that the finite scaling form significantly improved the data collapse over the scaling form with mean field exponents. With lattice parameters appropriate to the Ising ferromagnetic compound LiHoF4,we obtain a ferromagnetic transition temperature Tc= 1.51K in excellent agreement with experiment. This indicates that the material LiHoF4 is dominated by the dipole-dipole interaction; since in the simulations we only include dipole-dipole interactions. For dipoles on the simple cubic lattice, the ordered state is made up of anti-ferromagnetic rows. The critical exponents obtained by finite size scaling are ß~1/7,y ~ 8/7 and a ~4/7. These results are in good agreement with those of high temperature series expansions. A model of self—organized ruptures in an elastic medium is developed; and applied to earthquakes. In the model the local ruptures are represented by double couples to be consistent with elastic theory. The explicit form of this double couple source is derived. The system is driven by slowly increasing the shear stress. The model evolves towards a self-organized critical state in which the earthquake distribution follows the Gutenberg Richter law with an exponent in agreement with observational data. By modeling the local static fatigue for the rocks, we also obtained Omori’s law for the rate of aftershocks. The effects of annealing are investigated. |
Extent | 2316098 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-12-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085635 |
URI | http://hdl.handle.net/2429/3143 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1992-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1992_fall_xu_huang-jian.pdf [ 2.21MB ]
- Metadata
- JSON: 831-1.0085635.json
- JSON-LD: 831-1.0085635-ld.json
- RDF/XML (Pretty): 831-1.0085635-rdf.xml
- RDF/JSON: 831-1.0085635-rdf.json
- Turtle: 831-1.0085635-turtle.txt
- N-Triples: 831-1.0085635-rdf-ntriples.txt
- Original Record: 831-1.0085635-source.json
- Full Text
- 831-1.0085635-fulltext.txt
- Citation
- 831-1.0085635.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0085635/manifest