COMPUTER SIMULATION STUDY OF FERROELECTRIC LIQUID CRYSTAL PHASES B y Gary Steven Douglas Ayton B. Sc. (Chemistry) U . B. C. 1990 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y in T H E FACULTY OF G R A D U A T E STUDIES (Department of Chemistry) We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C O L U M B I A June 1996 (c) Gary Steven Douglas Ayton, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of The University of British Columbia Vancouver, Canada DE-6 (2/88) A b s t r a c t Molecular dynamics and Monte Carlo simulations are used to study ferroelectric liquid crystals. The effect of molecular shape on the stability and nature of ferroelectric phases is examined. The effect of molecular flatness on the stability of uniaxial nematic liquid crystal phases is studied. Most liquid crystal molecules are elongated and have a degree of flatness. However, most models used in computer simulation of liquid crystals are axially symmetric. A series of molecular dynamics simulations, specially designed to isolate the effects of molecular flatness on liquid crystal phases, are performed. Molecular shape is approximated with a generalization of the Gaussian overlap model [ B. J . Berne and P. Pechukas, J . Chem. Phys. 56, 4213 (1972)]. The form presented here has been extended to include ellipsoidal particles with non-degenerate semi-axes. It is found that small amounts of molecular biaxiality can drive an isotropic to nematic phase transition. Simulations of randomly frozen and dynamically disordered dipolar soft spheres are used to study ferroelectric ordering in spatially amorphous materials. Systems where the dipole moment has 1, 2, and 3 components are considered. It is found that the 1 component (Ising) model has ferroelectric phases. The systems with 2 and 3 dipolar components form disordered phases at low temperatures. A ferroelectric phase diagram is constructed for oblate molecules with point dipoles embedded along the particle symmetry axes. The role of particle shape on the stability of ferroelectric liquid crystals is examined with molecular dynamics simulation. Molecular shape is modeled with the generalized Gaussian overlap. Ferroelectric phases are found in systems with weak (short-ranged) columnar correlations. Systems with long-ranged n columnar order are found to be antiferroelectric. It is found that the stability of fer-roelectric phase is very sensitive to details of molecular shape, and exists only in small regions of the phase diagram. A more robust model for a ferroelectric liquid crystal is developed. Monte Carlo calculations are used to examine ferroelectric order in fluids of disc-shape particles with embedded dipoles. The dipoles are uniformly distributed over a circular "patch" of finite size, placed in the central plane of the particle. It is shown that such systems may undergo spontaneous polarization to form a stable ferroelectric discotic nematic phase. m Table of Contents Abstract ii Table of Contents vi List of Figures vii Table of Symbols xiv Acknowledgement xvii 1 Introduction 1 2 Statistical Mechanics and Computer Simulation of Fluids 7 2.1 Equilibrium Statistical Mechanics 7 2.2 Intermolecular Potentials 10 2.3 Computer Simulation Techniques 12 2.3.1 Periodic Boundary Conditions 13 2.3.2 The Treatment of Long-Ranged Potentials Under Periodic Bound-ary Conditions: The Ewald Sum for Dipolar Systems 13 2.3.3 The Monte Carlo Method and Importance Sampling 16 2.3.4 Molecular Dynamics Method 20 3 A Simple Model Potential for Liquid Crystal Molecules 23 3.1 Generalized Gaussian Overlap Integral 25 i v 3.2 M D Simulation Results: The Effect of Molecular Biaxiality on Liquid Crystal Phases 29 3.2.1 Molecular Biaxiality and the Uniaxial Nematic Phase 34 3.3 Summary 45 4 Ferroelectric Phases of Amorphous Dipolar Systems 46 4.1 Model and Simulation Details 48 4.2 Results for Randomly Frozen Systems 49 4.3 Decoupled M D Simulation Technique 53 4.4 A Mean Field Theory: Isolating the Effect of Long-Ranged Interactions . 63 4.5 Summary 75 5 Ferroelectric Phases of Dipolar Discotic Particles 76 5.1 Solid Lattices of Oblate Particles 77 5.2 Ferroelectric Phases of Simple Dipolar Oblate Particles 83 5.2.1 Simulation Details 83 5.2.2 Results 85 5.2.3 Summary 107 6 A More Robust Model for a Ferroelectric Liquid Crystal 108 6.1 The Cut-Sphere Model 109 6.2 Simulation Details 112 6.3 Summary 124 7 Summary and Conclusions 125 Bibliography 129 v Appendices 133 A The Pair Potential for Point Dipoles 133 B The Ewald Sum for Dipolar Systems 136 C The Reaction Field 142 D Quaternion Parameters 144 D. l Quaternion Parameters and the Orientation of Rigid Bodies 144 E Evaluation of the Generalized Gaussian Overlap Integral 149 E . l Evaluation of the Forces and Torques for the Gaussian Overlap Integral . 152 v i L i s t of F i g u r e s 3.1 A sketch showing the molecular axes (x',y',z') and the principal radii o i , a2 and <T3 which characterize the ellipsoid model 27 3.2 (|r(i*) - r(0)| 2) for 2ax : 2a2 : 2a3 = 2.75 : 3 : 9 (dash line) and 2ax : 2a2 : 2<73 = 1.5 : 3 : 9 (dotted line). These systems are isotropic and uniaxial nematic phases, respectively 32 3.3 Order parameters as a function of particle thickness at fixed reduced tem-perature and density. The squares are the uniaxial order parameter QQQ and the circles are the biaxial order parameter Ql2. The error bars repre-sent one standard deviation 35 3.4 Snapshot of an instantaneous configuration for an isotropic phase with 2ax : 2 a 2 : 2<r3 = 2.875 : 3 : 9 36 3.5 Snapshot of an instantaneous configuration for a uniaxial nematic phase with 2<7i : 2<r2 : 2cr3 = 2.125 : 3 : 9. Here, the director is along the z' axis (up and down on the page) 37 3.6 Snapshots of instantaneous configurations for a biaxial nematic phase with 2(7i : 2<r2 : 2cr3 = 0.5 : 3 : 9. Here, the directors are along x' and z'. T h e snapshots are taken along the axes perpendicular to the z' director. . . . 38 3.7 <jf000(r*) for biaxial 2<7i : 2cr2 : 2cr3 = 1 : 3 : 9 (solid), uniaxial nematic 2ax : 2<r2 : 2cr3 = 2 : 3 : 9 (short dash), and isotropic 2o\ : 2 a 2 : 2<r3 = 3 : 3 : 9 (long dash) systems 40 v i i 3.8 < P2(x'1(0) • x'2(r*)) > for biaxial 2ax : 2a2 : 2a3 = 1 : 3 : 9 (solid), uniaxial nematic 2<7i : 2<r2 : 2cr3 = 2 : 3 : 9 (short dash), and isotropic 2(7! : 2a2 : 2<r3 = 3 : 3 : 9 (long dash) systems 41 3.9 < x'(t*) • x'(0) > for a system with 2ax : 2<r2 : 2<r3 = 2 : 3 : 9 (solid), and a system with 2o\ : 2cr2 : 2cr3 = 1.5 : 3 : 9 (short dash). T h e inset is In < x(t*) •x'(0)> 42 4.1 T h e radial distribution function g(r*) for soft spheres at p* = 0.8 and T* = 10.5 51 4.2 T h e polarization < P > at T^n vs lOO/iV for the randomly frozen Ising (circles), X Y (squares) and X Y Z (triangles). Results are included for 108 ( X Y Z only), 256, 500 and 864 particles 52 4.3 T h e root-mean-square spin length, < S r m s >, for the iV = 256 frozen X Y Z model 54 4.4 < P > vs T * (rotational) for dynamically random X Y Z systems. T h e squares, triangles and circles are for m* = 1, 5 and 10, respectively. T h e error bars represent one estimated standard deviation 58 4.5 < P > vs T* (rotational) for dynamically random X Y systems. T h e squares, triangles and circles are for m* = 1, 5 and 10, respectively. T h e error bars represent one estimated standard deviation 59 4.6 < P > vs T * (rotational) for the Ising model. Results are shown for dynamically random systems with m* = 1 (squares) and m* = 5 (triangles) and for the randomly frozen case (circles). The heat capacities per particle, Cv/N, are plotted vs T* (rotational) in the inset 60 v m 4.7 The mass dependence of the disordered to ferroelectric transition temper-ature Tp (rotational). The squares and triangles are results for the XY and XYZ models, respectively. The values of TF were obtained from plots of the heat capacity, Cv/N, vs T * (rotational) and a typical example is shown in the inset. The error bars represent estimated uncertainties in the peak position 62 4.8 < P > vs T* from the reaction field, as found from a simple mean field theory. The short dash, solid and long dash lines are for the XYZ, XY, and Ising models, respectively. The circles are the frozen TV = 256 Ising results 67 4.9 Th e reaction field contribution to the energy per particle < U^F > /TV for the Ising model, as calculated from simple mean field theory. The inset is the reduced heat capacity per particle Cv/N 69 4.10 Th e reaction field contribution to the energy < U^F > /TV for the XYZ model, as calculated from simple mean field theory. The inset is the re-duced heat capacity per particle Cv/N 70 4.11 Th e short and long-range (reaction field) contributions to the energy per particle < U* > /TV for the dynamically decoupled XYZ model. T h e solid line corresponds to < U* > /TV, the short dash is < U^p > /TV and the long dash is < UgT > /TV. The squares, triangles and circles are for m* = 1, 5 and 10, respectively 72 4.12 Th e short and long-range (reaction field) contributions to the energy per particle < U* > /TV for the dynamically decoupled XY model. The solid line corresponds to < U* > /TV, the short dash is < U^F > /TV and the long dash is < U$T > /TV. The squares, triangles and circles are for m* = 1, 5 and 10, respectively 73 ix 4.13 The short and long-range (reaction field) contributions to the energy per particle < U* > /N for the dynamically decoupled Ising model. T h e solid line corresponds to < U* > /N, the short dash is < URF > /N and the long dash is < UgT > /N. The squares, and triangles are for m* = 1 and m* = 5, respectively 74 5.1 A sketch of a tetragonal I lattice. The circles are "top" views of columns of particles. The dark shaded circles are columns of particles that are shifted 1/2 a particle diameter out of the page compared to the light shaded circles. 79 5.2 A sketch of a hexagonal lattice. The circles represent columns of particles. There is no interdigitation between columns of particles 80 5.3 Dipolar energy per particle U^,D/N for different lattices of oblate dipo-lar particles. Triangles are for an antiferroelectric columnar lattice, and squares are a ferroelectric tetragonal I lattice 82 5.4 T* — 1.9 results for the < P2 > order parameter as a function of reduced density, p*, for oblate particles with different aspect ratios. Squares are for b/a = 1.5, circles are b/a = 1.3 and triangles are b/a = 1.2 86 5.5 T * = 1.9 polarization (< P i >) as a function of reduced density p* for oblate particles with different aspect ratios. Squares are for b/a = 1.5, circles are b/a = 1.3 and triangles are b/a = 1.2 87 5.6 b/a = 1.45 p* = 0.792 < P2 > (circles) and polarization < P i > (triangles) as a function of reduced temperature, T* 88 5.7 g000(r*) (solid), g112(r*) (short dash) and gno(r*) (long dash) for a p* = 0.792, b/a = 1.45, T* = 1.35 system 91 5.8 0||(r[j) (solid) and g±(rl) (dash) for a p* = 0.792, b/a = 1.45, T* = 1.35 system 93 x 5.9 A computer "snapshot" of an b/a = 3, p* = 0.792, T* = 2.5 system. The phase is antiferroelectric columnar 94 5.10 g000(r*) (solid), g112{r*) (short dash) and glw(r*) (long dash) for a / = 0.792, b/a = 1.45, T* = 1.9 system 95 5.11 flf||(rfj) (solid) and g±(rl) (dash) for a p* = 0.792, b/a = 1.45, T* = 1.9 system 96 5.12 g000(r*) (solid), # 1 1 2 (r*) , (short dash) and <?n o(r*) (long dash) for & p* = 0.792, 6/a = 1.45, T* = 2.5 system. 97 5.13 <7N(rf|) (solid) and g±(r*±) (dash) for a/)* = 0.792, b/a = 1.45, T* = 2.5 system 98 5.14 A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 1.35 system. The phase is antiferroelectric columnar 100 5.15 A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 1.9 system. The phase is ferroelectric 101 5.16 A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 2.5 system. The phase is isotropic 102 5.17 The reduced energy per particle, < U* > /N as a function of reduced temperature, T*. The inset is the excess contribution to the reduced heat capacity per particle, < Cv > /N 103 5.18 T * = 1.9 , p* = 0.792, < P 2 >, (circles) and polarization < P x > (triangles) as a function of aspect ratio, b/a = width/height 105 5.19 A rough phase diagram for dipolar, p* — 3, oblate ellipsoids. The re-duced density is p* = 0.792. The width of the ferroelectric phase is only estimated. The phases are I (isotropic), F E L(ferroelectric liquid crystal), A F E C L (antiferroelectric columnar fluid), and S (solid) 106 6.1 The coordinate system used to evaluate the interaction potential. The relevant vectors defined in the text are shown. The inset is a cross sectional view of a cut-sphere with point dipoles dispersed over a circular patch of diameter D' I l l 6.2 The mean square displacement, ((R(r) - R ( 0 ) ) 2 ) , for T* = 0.02 (short dash), T* = 0.03 (solid), and T * = 0.05 (long dash). R ( r ) is the position vector of a particle and r represents the number of Monte Carlo moves per particle 114 6.3 The temperature dependence of the polarization < P > and reduced con-stant volume heat capacity per particle Cv/N (inset) for a system with D/L = 10, D'/L = 9 and p* = 0.525. The error bars represent one estimated standard deviation 115 6.4 "Snapshot" of a configuration for a system with D/L = 10 and p* = 0.525. The system with D'/L = 9 is shown in the ferroelectric (T* = 0.025) phase. 116 6.5 "Snapshot" of a configuration for a system with D/L = 10 and p* = 0.525. The system with D'/L = 9 is shown in the unpolarized (T* = 0.050) phase. 117 6.6 "Snapshot" of a configuration for a system with D/L = 10 and p* = 0.525. The system with D'/L = 5 is shown in the unpolarized (T* = 0.0275) phase. Note that this configuration consists of loose polarized columns with nearly equal numbers up and down. O n average, there is no net polarization in this system 118 6.7 gs(r*) f ° r T* = 0.025. No long-range smectic order is detected in this system. 120 xn 6.8 The variation of the reduced electrostatic potential energy, L3ues(12)/p?, as one particle is slid over another with their dipoles directed in the same direction and their surfaces in contact. The inset shows the same curves on a magnified scale. Curves are shown for single point dipoles (D'/L ~ 0) (long dash) and for D'/L = 4 (short dash), 5 (dotted) and 9 (solid). . . . 122 6.9 The contributions to the reduced energy per particle, < U* > /N. T h e squares are the total energy < U* > /N, the triangles are the reaction field contribution < U^F > /N, and the circles are the short-range structural contribution (everything else) < UsT > /N. The reduced density is p* = 0.525 123 A . l The geometry that is used for the point dipole pair potential 135 D . l Clockwise rotation of the vector p to p' about the unit vector n 146 x i n Table of Symbols Symbol Description b/a The aspect ratio of an oblate ellipsoid. b Shape tensor. B Rotation matrix. Cv The constant volume heat capacity. D' The diameter of the charge distribution embedded in the cut-sphere. d The director unit vector. g(r) The radial distribution function. <7||(r||), <7|_(r|_) Parallel and perpendicular pair correlation functions. &B The Boltzmann constant. K Kinetic energy. L Box length. m,m' The mass and reference mass of the particle. M The instantaneous moment of the sample, n A n integer lattice vector . < P > Average polarization. < Pi > The average first rank order parameter. < P2 > The average second rank order parameter. QNVT The partition function for the canonical ensemble. Ql0 Uniaxial nematic order parameter. Q 2 2 Biaxial nematic order parameter. . xiv qo, qi, ?2, 93 Quaternion parameters. R The reaction field. 5 Gaussian overlap "hardness" parameter. The root-mean-square spin length. r The absolute temperature. T Orientation-shape tensor. u(12) A pair potential. [/ Instantaneous configurational energy. <U> IN The average energy per particle. < URF > IN The average reaction field contribution to the energy per par-ticle. < UST > IN The average structural contribution to the energy per particle. V Volume. ZNVT The configurational integral. a The convergence parameter for the Ewald sum . a\ Moment of inertia parameter. 1 / f c B T V Gradient operator. £ The reference unit of energy. e The dielectric constant . A* The dipole moment vector. Number density. PNVT Probability density function for the canonical ensemble. <Kr) A n electrostatic potential . U) Angular velocity in the body fixed frame. XV a The reference unit of distance. 0 i 5 0"2> o"3 Ellipsoid principle radii. x v i Acknowledgement Almost five years ago I met Dr. Gren Patey. I remember vividly a few things of our first meeting. I remember seeing running shoes, a baseball mit, and piles of paper. Gren started telling me about his field of research. I think I caught the first few words; I heard the word "liquid" and "computer", but other than that I really had no idea what he did. In all honesty, my choice to work for Gren Patey was pure gut instinct. A n d it was the best choice that I have ever made. I hope some of Gren's patience, insight, and ability to see through the little details and get to the important physical picture has rubbed off on me. He has taught me so much more than the chemistry of the fluid state, he has shown me the art of science. I consider myself lucky. I would like to thank Dr. Michel Gingras for the opportunity to work with him, and for his infinitely long and valuable comments. Although I am quite sure that Michel and I agreed on possibly nothing, I found our collaboration to be challenging, and very rewarding. Our lengthy discussions motivated some of work in this Thesis. Again, I consider myself lucky. U N I X ? F O R T R A N ? I never realized how computer illiterate I was. Dr. Dan Berard saved my life. In my formative years as a theoretical chemist I met my nemesis, the computer. If it were not for Dan, I would have gone mad. I don't know how many times I said "hey Dan, how do I..." Dan gave me his S M O N G O plot macros, his c-shells, his L A T E X macros, and his thesis macros. In fact, every plot I have done was made possible because of Dan. However, I still hate computers. xvn Final ly , John Shelley not only put finger-prints, but great computer graphics on the screen. A l l of the computer snap-shots were made possible with his M O V I E program. A picture is worth a 1000 words. A n d as I have found out in the course of writ ing my Thesis, one has a picture A N D a 1000 words. I have had the pleasure of working a truly great group of people. There were many, many Friday nights spent at the G r a d Center. I won't forget them. Natalie, Claudio, Richard, M a r k , John, L i a m , and Frank, it's been fun under the pa lm tree. xvin Chapter 1 Introduct ion What is a phase of matter? Sometimes, the phases of matter are so familiar that we forget how to define them. We are accustomed to the solid, the liquid and the gas phase, so when we are given a rock, we call it a solid. When we are given a bucket of water, we say it is filled with a liquid". Likewise, when we are given a balloon, we say it is filled with a gas. But how do we come to these conclusions? What set of criteria do we use to identify a phase of matter? For example, we call a candle a solid and vegetable oil a fluid. But what is grease? Or the goo in a soap dish? It would be valuable to construct a set of simple guides to help us define the phases of matter. If some substance can be seen (even with a microscope), looks uniform through-out, and does not appear to be changing into something else, then it is probably a phase of matter. We can define a liquid as a phase of matter that will take the shape of its container. Conversely, a solid will never take the shape of its container. We can also define a phase of matter by describing the motion of the particles. In the liquid phase, the particles are close together, but are constantly moving, much like people in a crowded room. This explains the ability of the liquid phase to take the shape of its container. Of course, there are different degrees of motion. In oil, the particles are moving fairly rapidly, but in wax, they are moving more slowly. In the solid phase, particle motion is so limited that the molecules only shake back and forth. A phase where the particles are moving V E R Y slowly might be confused with a solid. However, if it takes thousands of years for a particle to travel on average a few Angstroms, for "human" purposes the 1 Chapter 1. Introduction 2 phase is solid. We can further describe a phase of matter by the spatial arrangement of the particles. T h e spatial arrangement, or structure, is determined by the degree of positional and orientational order of the phase. The degree of order indicates whether the positional or orientational structure persists over long (macroscopic) ranges, or short (microscopic) ranges. For example, the crystal has long-ranged positional order, that is, if one particle is "tagged", and the space around it is studied, there will be other particles found at regular intervals. The fluid phase has short-ranged positional and orientational order. At small distances, there might be regions of space where the probability of finding another particle in a certain orientation is greater than in other regions. At farther distances, this structure decays, and particles can be found anywhere, with any orientation. W i t h the idea of long and short-ranged positional and orientational order, we can describe the liquid crystal phase of matter. In this phase there is only short-ranged positional order, but there is long-ranged orientational order. The particles are liquid-like in their translational motion, but crystal-like in their rotational motion. In other words, the particles are free to move about, but they point in roughly the same direction. This phase has been known for over one hundred years and was first identified by an Austrian botanist, Friedrich Reinitzer (although this is still debated) and Otto Lehmann, a German physicist, when studying the melting behavior of cholesterol derivatives in plants [1]. For one sample, there appeared to be two melting temperatures; one at 145.5C and one at 178.5C. A cloudy blue liquid was observed in between these two temperatures. Above 178.5C the sample was a clear fluid, and below 145.5C it was solid. The cloudy blue fluid was similar to crystals in that it interacted with polarized light, but it also had fluid-like flow properties. Around 1916, Born [2] suggested that the dipole moments of the molecules were the dominant interactions in this phase. In 1960 Maier and Saupe [3] developed a mean field theory that used dispersive interactions, that is Chapter 1. Introduction 3 instantaneous dipole moments induced in one neutral molecule due to the presence of another, as the driving force for liquid crystal phases. There are a variety of liquid crystal phases, where each one is characterized with a specific name. The simplest liquid crystal phase is the nematic phase. Here the molecules have translational freedom, that is, the translational motion is fluid-like, but on average the particles have the same orientation. The average direction can be described by a vector, the director d. In the smectic phase, there is long-ranged positional order where the particles arrange themselves in layers, much like a stack of pop cans. If the particles arrange themselves in columns, the phase is called columnar. There are more complicated types of liquid crystals, but they will not be discussed here. Most liquid crystals have a similar structure and are neutral, large organic particles that are significantly elongated and have a degree of flatness [4]. They usually have a "backbone" of aromatic linkages, with functional groups at the ends [4]. Recently, liquid crystal phases of disc-shaped molecules [5] have also been observed. In the discotic nematic phase, the director is associated with the symmetry axes of the molecules. From the observation that most liquid crystal molecules have similar structures, it seems reasonable to consider that molecular shape is important to the stability of the liquid crystal phase. Computer simulations performed by Frenkel et al. [6] showed that hard ellipsoids of revolution can have liquid crystal phases if the elongation of the particles is sufficient. Systems of prolate ellipsoids with aspect ratio a/b > 3 were found to spontaneously order and form a stable nematic phase. Also, oblate ellipsoids exhibited nematic phases for aspect ratios a/b < 1/3. The results of these simple hard models suggest that the anisotropic "shape" of liquid crystal molecules plays an important role in stabilizing the liquid crystal over the isotropic phase. However, this model does not explain the role of attractive dispersion forces, or the the role of long flexible alkane or alkoxy end chains [4]. Most isotropic fluids freeze directly to the solid. However, liquid crystal molecules form Chapter 1. Introduction 4 an orientationally ordered fluid before they freeze. The stability of the liquid crystal over the solid phase depends not only on the "shape" of the molecule, but on details of the molecular structure. The unique properties of liquid crystals, namely their solid-like orientational order, but fluid-like flow properties, generates very interesting behavior in the presence of an electric field. For liquid crystal molecules without a permanent dipole moment, each particle "feels" the electric field through its induced dipole moment. Since the particle is elongated, the induced dipole moment could be larger in some directions than in others. If the largest induced dipole moment is parallel to the long axis, the particle will tend to align with the electric field. However, if the largest induced dipole moment is along some other axis, the particle will align at some angle to the electric field. Since in a liquid crystal all the particles are oriented in about the same direction, A L L of the particles will similarly align. Therefore, we have a means of easily controlling the orientation of Avogadro's number of molecules. This type of control would be very difficult in the solid phase since the particles cannot rotate, and not easy in the isotropic phase since the aligned state is not stable. Once the field is turned off, the ordered fluid would return to isotropic phase. In fact, it is this response of liquid crystals to electric fields that has been employed in all of the liquid crystal displays ( L C D ) used in watches, computer screens, and scientific gauges. B y controlling the orientation of many particles, we can "switch" their orientations, much like opening or closing the blinds on a window. T h e response of a liquid crystal to an electric field depends, to a large extent, on the strength of the induced dipole moment. If on average, the induced dipole moment is small, the response time will be slow and a strong electric field will be needed. Conversely, if the induced dipole moment is large, the response time (that is the average time it takes for a particle to align with the field) will be shorter, and a weaker field will suffice. Chapter 1. Introduction 5 A ferroelectric liquid crystal would consist of molecules with permanent dipole mo-ments where the dipole moments are all aligned. In this case, the liquid crystal phase would be polarized, the response time to electric fields would be less [1], and the strength of the field could be reduced. However, to date, a "true" ferroelectric liquid crystal has not been observed. The closest candidate is the chiral smectic C * phase, where the par-ticles arrange themselves in layers, but only each layer is polarized. Meyer [7] suggested that the "in plane" ferroelectric nature was a result of symmetry requirements, while Tredgold [8] challenged this theory, stating that the ferroelectric state was a result of short-range molecular repulsions. I would like to pose an interesting theoretical question: a true ferroelectric liquid crys-tal phase has not been observed experimentally. Why? There are two possible reasons. Firstly, such systems do not exist. Secondly, none has been found. At least theoretically, both of these statements have been challenged. Recently, computer simulations of dipo-lar hard and soft spheres have been found to exhibit a stable ferroelectric nematic phase [9,10]. The orientational order comes from purely dipole-dipole interactions. From a the-oretical perspective, this is an exciting result, however, from an experimental viewpoint a real ferroelectric nematic phase is still far off. In reality it would be very difficult to construct a near spherical particle with a large dipole moment. Clearly, these preliminary results demonstrate that a ferroelectric nematic phase can exist. However, much more work, not only in understanding the existence of this phase, but in developing new and more realistic models is necessary before any detailed laboratory experiments should be performed. This Thesis will address several questions concerning ferroelectric liquid crystals. What factors favor or disfavor ferroelectric liquid crystal phase formation? Can mod-els other than dipolar spheres be found that have ferroelectric liquid crystal phases? Is Chapter 1. Introduction 6 there a possibility of constructing a ferroelectric liquid crystal in a laboratory? C o m -puter simulations are used extensively in this study. These computer "experiments" will address these questions and offer some insight, explanations and suggestions. This Thesis is divided into 7 Chapters, and a series of Appendices containing details of some relevant computational tools and derivations. In Chapter 2, the statistical me-chanical tools needed to study liquid crystals with computer simulation will be discussed. Chapter 3 examines how the flatness associated with most liquid crystal molecules affects the stability of nematic phases. In Chapter 4, the ferroelectric ordering and disordering interactions that exist in dipolar sphere systems are studied in detail. Systems of amor-phous dipolar spheres will be used to study the effect of short-range spatial structure on ferroelectric phases. Chapter 5 examines the effect of molecular shape on the stability of ferroelectric systems. The ferroelectric phase behavior of oblate particles with dipole moments along the symmetry axes is examined. In Chapter 6, a more realistic model for a ferroelectric liquid crystal is developed and studied. Finally, Chapter 7 is a summary and review of the major conclusions. C h a p t e r 2 S t a t i s t i c a l M e c h a n i c s a n d C o m p u t e r S i m u l a t i o n of F l u i d s Computer simulation of fluids is one application of equilibrium statistical mechan-ics. In this Chapter, we will discuss the necessary tools to study equilibrium fluids by computer simulation. 2.1 E q u i l i b r i u m S t a t i s t i c a l M e c h a n i c s In a classical fluid, the Hamiltonian for a system of N particles can be written as HN(^p) = KN(p) + UN(q), (2.1.1) where AOv(p) is the total kinetic energy and t/;v(q) is the total potential energy. A set of generalized coordinates [11] is given by q = (qi , q 2 , ...q./v) and the corresponding set of conjugate momenta is p = ( p i , p 2 , •••PN)- The classical approximation is valid for particle translations as long as the thermal de Broglie wavelength is less than the mean nearest neighbor separation [12]. For particle rotations, this approximation is valid only at high temperatures, where the ratio T/Qr is large. Here, T is the absolute temperature and Qr is the characteristic temperature of rotation [13]. Of course molecular vibrations and electronic motion cannot be correctly treated classically. A simple solution to this problem is to treat the molecule as a rigid body [14]. The time evolution of is determined by Hamilton's equations [11,12,16] , as given by <9H/v q, = -g j j - , (2.1.2a) 7 Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 8 P i = (2.1.2b) where, i refers to the ith particle. For a molecular system, each particle will require six coordinates and six momenta. Three coordinates describe the particle's position, three specify orientation, three are translational momenta and three are angular momenta. Of course, if the particle is a symmetric top, these coordinates will be reduced to ten, and if it is spherical, they are reduced to six. In the general case, iJ/v(q, p) is a function of 12iV variables. These variables are simply coordinates in a multidimensional phase space [13]. T h e system will evolve in time, sampling points in phase space as determined by Hamilton's equations of motion. A n y mechanical thermodynamic observable, W , can then be written as a time average, However, as an alternative to calculating mechanical observables as time averages, Gibbs suggested removing the time dependence and instead, calculating ensemble aver-ages. A n ensemble is a collection of systems, identical except for their values of (q, p), that are "prepared" such that they are in equilibrium and have the same thermodynamic state parameters. Each member of the ensemble has a corresponding phase point, and the ensemble forms a "cloud" of points in phase space. The thermodynamic state pa-rameters determine the accessible phase points. For example, if the number of particles N, the volume V and the energy E are constant (NVE), then the phase points are on a hypersurface of constant E. If the number of particles, volume V and temperature, T are constant (NVT), then all phase points are accessible, but some are more probable than others. We will only consider the case of constant (NVT), known as the canon-ical ensemble. The probability of finding a member of the ensemble with phase space coordinates between q and q + dq, as well as p and p + dp is given by PNVT(^ p)dqdp, (2.1.3) Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 9 where dc\ — dqidc[2---d(\N, and dp = dpidp2-..dpN- That is, if the number of mem-bers of the ensemble is S, the number in the small volume element dqdp is given by SPNVT{<\-> p)dqdp. The function PNVT{^ p) is the probability density for the canonical ensemble. Mechanical observables can then be calculated from <W>= J W(q, p)PNVT(q, p)dqdp, (2.1.4) where PNVT(<1, p) can be written as [12] ( n ^\ 1 ' 1 e x P ( - ^ A f ( q , p ) ) / f c B T ) PNVT(q,p) = ^ • (2.1.5) Here, QNVT is the partition function and has the form, QNVT = -^j^ J dqdpexv(-HN(q,p))/kBT), (2.1.6) where kB is the Boltzmann constant. The factor x is the number of coordinates required to completely describe the position and orientation of a particle (6 for an asymmetric top, 5 for a symmetric top, and 3 for a spherical body.) The factor of takes account of the indistinguishability of the particles, and the factor of h~xN is from quantum cor-respondence [13]. The expression for PNVT can be simplified by recalling E q . (2.1.1) and factoring out and averaging over terms that contribute to the kinetic energy. This gives exp(-UN(r,n))/kBT) PNVT = ~ , (Z.l.t) ZJNVT where the set of generalized coordinates q = (qi, q 2 , ...qyv) has been split into a set that describe particle positions, r = ( r i , r 2 , • ••YN) and a set that describe particle orientations $7 = ( f i i , f i 2 , . . . O N ) . The term ZNVT is the configurational integral, and is defined as Z N V T = J J drdnexp(-UN(r,n))/kBT). (2.1.8) Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 10 2.2 I n t e r m o l e c u l a r P o t e n t i a l s In real molecular systems, the configurational energy U arises from all molecular inter-actions. In general, for TV particles it can be written as [14] = £«(»)+ £5>(*i) + E E E «&•*) + ..., (2-2.i) i i j>i i j>i k>j>i where the sums include distinct pairs, triplets and higher terms. The first term arises from any external fields. The second term is from pairwise interactions and usually dominates the total energy. The third and higher terms come from many-body interactions. Three-body interactions (eg. the triple-dipole dispersion potential) can be significant in dense systems [12,15] and, therefore, a simple truncation of the expansion at the pairwise term could be a poor approximation. However, if only the most robust features of a molecule, are needed, for example the approximate "shape", then a pairwise additive potential is sufficient. A convenient form for a pair potential is u(12) = u r e p (12) + 1 ^ ( 1 2 ) , (2.2.2) where u r e p (12) is a strongly repulsive term that defines the "shape" of the particle, and itir(12) represents longer ranged electrostatic as well as attractive dispersion interactions. The short-range repulsion originates from the overlap of the charge clouds of the two molecules, and is formally a quantum mechanical problem [16,17]. The exact functional form is complicated, however, if only the approximate "shape" of the molecule is desired, a variety of simple repulsive functions can be used. The simplest is a "hard" potential [18] and can be written as tthc(12) = 0, no particle overlap, Uhc(12) = oo, particle overlap, (2.2.3) Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 11 where "he" stands for "hard core". This potential has been used to model hard spheres, ellipsoids [19-21], spherocylinders [22], and "cut" spheres [23,24]. Here, a cut sphere is simply a sphere where the top and bottom, equidistant from the center of the sphere are cut out. The next class of simple model repulsive potentials are continuous functional forms, written as U c t(12) = e / ( - ) , (2.2.4) where, a is the unit of length and is roughly the distance of closest approach for the two particles. The intermolecular distance is r = |r| , and e is the fundamental unit energy. For spherical particles, a simple example is the "soft" sphere potential, uss(12) = £ ( - ) 1 2 . (2.2.5) r The choice of 12 for the exponent is arbitrary. In this Thesis, three repulsive potentials corresponding to different particle shapes will be used. These are soft spheres, "cut" spheres, and the Gaussian overlap. The generalized Gaussian overlap potential was developed as part of this Thesis as a conceptually simple way of modeling ellipsoidal molecules whose length, breadth and width are different. Shapes similar to "surf boards" and "used bars of soap" can easily be modeled. The it]r(12) contribution is usually a sum of terms, representing dispersion and multi-polar interactions. Dispersion interactions are a result of the instantaneous fluctuations in the charge clouds of molecules and are attractive. T h e best known model that incor-porates dispersion interactions is the Lennard-Jones (LJ) potential written as u L J (12) = As ((^)12 - i^f) . (2.2.6) W i t h appropriate choices for e and a, computer simulations with the LJ(12-6) potential have been used to model liquid argon . Electrostatic interactions due to the charge Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 12 distributions of a molecule can be included as a multipole expansion [16]. For neutral molecules, there can be dipole-dipole, dipole-quadrupole, quadrupole-quadrupole and higher interactions. Using a multipole expansion is usually only successful for simple charge distributions, for example, when a distribution can be modeled with a simple point dipole. The interaction energy of two point dipoles, p,\ and /* 2, separated by r is given by the well known expression, uDD{\2) = / x 2 - V ^ ( r ) p,x • p2 3(A*I • r ) ( > 2 • r) (2.2.7) where <f>(r) is the potential at // 2, due to p,\. A simple derivation of the pair potential between two point dipoles is given in Appendix A . 2.3 C o m p u t e r S i m u l a t i o n T e c h n i q u e s There are two basic techniques used to simulate classical fluids. These are the Monte Carlo ( M C ) and Molecular Dynamics (MD) methods discussed below. Very briefly, in Monte Carlo simulations, configurations are generated according to PNVT-, so that aver-ages can be calculated from i M <W>=uEW^ (2-3-1) i=l where M is the number of configurations selected from the distribution. W i t h molecular dynamics simulations, the particles move according to the equations of motion, so that the trajectory of the phase point correctly samples the distribution function. For example, under (NVT) conditions, the phase point will sample phase space according to PNVT-Averages can then be calculated by i M <W>=^EWt, (2.3.2) 1V1 t=i where Wf. is the instantaneous value at time t. Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 13 2.3.1 Periodic Boundary Conditions In a computer simulation, evaluating the potential energy U is the most time consuming calculation. For a system of N particles interacting with a pairwise additive potential, the calculation of the U involves N(N — l)/2 operations. The size of the simulation is limited by the power of the computer. Generally, the two limiting factors are the number of particles and the "complexity" of the potential. The number of particles is usually less than 1000, which is much less than a bulk system. A finite system of around 1000 molecules would have significant surface effects and would not adequately represent a bulk N ~ 10 2 3 system. To avoid this problem in simulations, the technique of periodic boundary conditions is used. That is, a small system of particles, N < 1000 , is repeated periodically in space. The small system is contained in a "box", or "cell", and the cell is repeated in a lattice. As one particle leaves one side of the cell, it enters on the other side. In this way, the "surface" is removed. Under periodic boundary conditions, a particle in the central simulation cell will interact with the other particles within the cell, with the images of the other particles, and with the images of itself. If the potential between the particles is short-ranged, then the interaction with the central particle with the long-ranged periodic images is negligible. However, if the potential is not short-ranged as in the dipole-dipole case, the interaction of the central particle with the long-ranged images can be significant and must be considered. This problem will be discussed in the next section. 2.3.2 The Treatment of Long-Ranged Potentials Under Periodic Boundary Conditions: The Ewald Sum for Dipolar Systems The Ewald summation method is a means of handling long-ranged potentials under periodic boundary conditions. W i t h this method, we surround the central simulation Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 14 cell with an "infinite" sphere of periodic images. The images form a three-dimensional periodic lattice. Also, beyond the "infinite" sphere is a dielectric continuum characterized by the dielectric constant of the system*. Even though the sphere is infinite, we must specify the surrounding material. Hence, we should think of the sphere as being a V E R Y large sample embedded within another larger continuum. A potential (j>(r) = A/rn is considered to be long-ranged if n < 3, [25]. W i t h long-ranged potentials, the energy of the N particles in the central simulation cell U is written as 1 O O U = oT, 2 n N N E E " ( l r + n^l t=i i=i (2.3.3) where L is the length of the cell, and r = r3; — r;. The vector n designates an integer simple cubic lattice point with n = {nxL, nyL, nzL) and n r , ny, nz integers. The prime on the first sum means that the term corresponding to i = j is omitted when n = 0. T h e pair potential u(|r + nL\) gives the interaction between a particle i in the central cell, and a particle j that is either in the central cell (n = 0) or is a periodic image located at r + n L (n ^ 0). The lattice sum is conditionally convergent [25], meaning that the value of U will depend on the order that the periodic images are summed. The particular order of the summation will determine the "shape" of the periodic array. For example, two different values of U will be found depending on whether the summation order was done by adding up cubic shells or spherical shells. Cubic shells are simply cubes of thickness L surrounding the central simulation cell, while spherical shells are contructed by building "jagged" spheres of periodic images. As the radius of the sphere becomes larger, the shell becomes "smoother". For the case n —> oo, a "natural" order of summation is to add up roughly spherical shells about the central cell. In this case the energy USH (SH = spherical shell) of the central cell corresponds to the central cell surrounded by Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 15 an infinite sphere composed of periodic images. Deleeuw et al. [25] showed that conditional convergence, or the order of summation can be removed by introducing a "convergence parameter" f(s). T h e convergence pa-rameter applies to a particular order of summation, hence the choice of the convergence parameter defines the shape of the array. Here, the convergence parameter was chosen to correspond to the summation of spherical shells. The energy is now written as U(s), and USH can be found from 1 CO ' EI>(|r + n I | ) / ( 3 ) t = l j = l (2.3.4) 2-USH = lim*7(s). (2.3.5) s—>0 T h e convergence parameter makes E q . (2.3.4) absolutely convergent and thus removes the order of summation. This makes evaluating USH much easier, since the summations in E q . (2.3.4) can be taken in any order. We will restrict ourselves to the problem of point dipoles under periodic boundary conditions. This situation has been studied extensively [25,26]. We consider a system of N point dipoles, given by {fix, J A 2 , /*3, ••• MAT}, in a cubic box of length L under periodic boundary conditions. T h e total energy of the dipoles in the central simulation cell, the interactions of the dipoles in the cell with the other other periodic images, as well as the interaction of the surrounding continuum of dielectric constant e , is given by U = ^EIXBDC(12; e'), (2-3.6) where v^^(\% e) is an effective pair potential and has the form u PBC, D D \ 12; e) = -(/ix • V)(/x 2 • V ) * ( r ) 4TT 4 ^ 2 a 3 2(e' - 1) 4TT + ^ • ^ - 3 ^ - ( 2 7 T l ) 3 ^ / i l - / X 2 ' ( 2 - 3 J ) Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 16 where a is a convergence parameter. The term ^ ( r ) has the form (2.3.8) where the first term and second terms are known as the real space and Fourier space sums respectively. Details of the derivation leading up to this result can be found in Appendix B and C. In practice, the the value of a is chosen so that only the first term n = (0,0,0) in the real space sum is used. Convergence of the Ewald sum for different values for a has been studied [27]. For the simulations in this Thesis, the parameters employed were a = 6.4 and n 2 = 27. Other parameter sets were tested to ensure that the results were not affected by the particular choice of a and n 2 . The term that depends on the dielectric constant of the surroundings, e', is due to the polarization of the surrounding dielectric. It is, in fact, the reaction field contribution. The reaction field will be discussed in detail in later chapters, as well as in Appendix C. For now, the reaction field contribution to the pair potential is simply a contribution to the energy from the polarization of the surroundings beyond the spherical array. In the case of e' = oo, this term exactly cancels the term j^Vi • p>2-2.3.3 T h e M o n t e C a r l o M e t h o d a n d I m p o r t a n c e S a m p l i n g The Monte Carlo Method is simply a means of estimating multidimensional integrals. In the context of equilibrium statistical mechanics, it is a method of calculating ensemble averages. For example, the mechanical observable < W > is defined by where /3 = 1/fcfiT. Estimating < W > with the Monte Carlo method could in principle be done by randomly selecting a number of configurations ( M configurations in total), <W> = JJWe-Pudrd£l (2.3.9) (2.3.10) Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 17 and for each one, calculating W ( r , fi) and weighting it by e $ u ( v & ) . The estimate of < W > is then 1 TM w-e-pu< < W > ^ = M (2-3.11) However, calculating ensemble averages with this crude method of Monte Carlo integra-tion is impractical, since a good estimate of ZNVT would require generating a V E R Y large set of configurations. Almost all of the randomly chosen configurations would be highly improbable, and would have very small contributions to the integral. Instead, we will use importance sampling, where configurations are no longer chosen at random, but are concentrated around important contributions to the integral in E q . (2.3.9). If con-figurations are selected from PNVT, then ensemble averages can be calculated from E q . (2.3.1). Now the problem is how to sample configurations from the distribution, PNVT-The method used to sample PNVT was originally developed by Metropolis et al. [28] and involves constructing a Markov chain of configurations. The set of configurations that are generated by the Markov chain will be selected with a probability given by the limiting distribution, PNVT- Briefly, a Markov chain is simply a sequence of configurations (or states) { c i , c 2 , C 3 . . . . C M } such that the probability of finding the system in a state n at some point c; in the sequence depends only on the state m that the system occupied at C j _ i - Th e idea is to continue generating new configurations, where each new configuration depends only on the last one in the chain. Ensemble averages can then be collected by E q . (2.3.1), using the configurations in the chain. Obviously, the trick is to generate new configurations so that they are chosen according to PNVT-Generating new configurations in the Markov chain is done by choosing an appropriate transition probability, irmn, defined as the probability of going from the state m to n in one "step". Then the probability that the system is in a state n, given by pn = PNVT{^n,^n)-, Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 18 is = Pn, (2.3.12) m where the sum over m is a sum over all states. The configuration associated with the nth state is written as ( r n , fin). This equation can be conveniently written in the matrix notation, PNVTT = PNVT- (2.3.13) Here, the limiting distribution is written as a vector with each element k corresponding to one configuration, pk = pNVT(rk,^k)- TT is a stochastic matrix since its rows add to 1, that is Y,m Knm = 1 • There are a variety of methods for finding 7r, however we will only use the solution first suggested by Metropolis et al, known as the asymmetrical solution [14,28]. There are two steps involved in calculating the value of irmn. First we impose the constraint of "microscopic reversibility", given by Pm^mn = Pn^nm- (2.3.14) If both sides of this equation are summed over all m, it is simple to show that this equation satisfies E q . (2.3.13). We then apply the asymmetrical solution, where the transition matrix elements are found from T^mn = OCmn pn ^ pm 771 ^ TI TTmn = « m n ( ^ ) Pn < Pm m^U (2.3.15) ""mm = 1 Z /n^ra ^mn-Here, a is a symmetric stochastic matrix. It is possible to show using the symmetric property of a that this n satisfies E q . (2.3.14) and hence E q . (2.3.13). We now have a means of generating new configurations in the Markov chain. In practice, using the Monte Carlo ( M C ) simulation technique for fluids involves randomly choosing one of N molecules and attempting a trial move. The move is either Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 19 accepted, and the new configuration is the next in the Markov chain, or it is rejected, and the old configuration is the next configuration in the chain. A move can be a combined or separate translation and rotation, or a more elaborate move such as rotating many particles at once. For translations, a particle i is chosen (usually sequentially) and is displaced from r™ to r™, where the new position can be anywhere within a cube of dimensions 25rmax centered at i*™. This corresponds to a stochastic matrix with elements Qt-mn — if MR if r " is within the cube, and a m n = 0 if r™ is outside the cube. MR is simply all the (finite) positions available within the cube. The change in the energy of the system is calculated by considering the energy of the ith particle before and after the move, given by TV N i=i j=i If SU < 0, then the new configuration is accepted. If 6U > 0, then the move is accepted with a probability pn/pm, in this case given by Pn _ Z N V T exp(-/3Un)dr Pm ZNVTexp(-pUm)dr = e x p ( - / W ) . (2.3.17) T h e volume elements dr have been explicitly included, however they are identical in the m and n states, and therefore cancel. To accept the move with this probability, a random number n r a n d o m is generated uniformly on the interval (0 ,1) . If n r a n d o m < exp(—/?<?£/), then the move is accepted, otherwise the move is rejected. For rotations, the general case is when the orientation of a molecule can be described by three Euler angles, <f>, $, and ip [11 ,14] . However, when the required ratio of the probabilities is written, it is clear that the "orientational' volume elements for the particle that was rotated do not cancel. For non-linear molecules, this volume element is dSl7? = am 0?d6? dip? d<p?. (2.3.18) Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 20 T h e ratios of the probabilities can then be written as pn ZNyT exp(-/3Um + 8U) sin Q\ pm ZNVTexp(-f3Um) sin Of1 = e x P ( - W ^ f • (2-3-19) In this form, random orientational displacements have been done by small changes in the Euler angles themselves, 8(f>, 89, and 8ip. However, if random displacements are in cos Of rather than Of, this amounts to taking displacements in a new variable £, such that £f = cos Of, and d£f = — sin0™dcos Of. The ratio of the probabilities is now given by E q . (2.3.17), and the translations and rotations can be combined. Furthermore, the calculation ZNVT is avoided, since we only deal with ratios of probabilities. 2.3.4 Molecu la r Dynamics M e t h o d Molecular Dynamics solves the equations of motion for a system of N molecules inter-acting through a potential U as in E q . (2.1.1). If the trajectories of the particles follow Newtonian equations of motion, it is easy to show [29] that the total energy H, is con-served. That is, the Hamiltonian is a constant of the motion. If the particles interact with a continuous potential, the force on the ith particle is F,- and can be found from Fj = —V£/j, where Ui is the ith particle's contribution to the energy. For pairwise additive potentials, Ui can be found from Ui = E u{ij). (2.3.20) j=l,i^j The total force F j , acting on the ith particle is then Fj = YljLi,i^j f (u) where the force on the ith particle due to j is f(ij) = —Vu(ij) . For non-spherical particles, the torque on particle i due to j , Tij can also be calculated from the pair potential u(ij). T h e method that is used is described in Appendix E . Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 21 T h e idea with M D is to obtain new positions, velocities, orientations and angular velocities at some time t + 6t, given the initial values at some time t. T h e trajectories of the molecules are solved on a step by step basis, where the time interval, or timestep, is given by St. T h e new positions, velocities, orientations and angular velocities as well as higher derivatives are calculated by finite difference methods [14,30,31]. Here we use the Gear predictor-corrector algorithm [14]. M D Simulations that constrain N and V to be constant and use Newton's equations of motion will sample the microcanonical (NVE) ensemble. If we wish to sample the canonical (NVT) ensemble, the temperature T rather than the total energy E must be constant. The equations of motion for a classical particle in the canonical ensemble can easily be formulated by constructing a set of constrained trajectories that will ensure that the temperature T is constant. Obviously Newtonian equations of motion, which imply that the energy is a constant, do not apply. Instead, constant kinetic temperature dynamics is needed [27]. Following Evans [32], we will use the procedure of Gauss's principle of least constraint. Evans [33] showed that the equations of motion generated under the constraint of constant kinetic temperature will have trajectories that sample PNVT- The constraint of constant kinetic temperature can be written as y E | r l | 2 - ^ B T = 0, (2.3.21) which assumes the equipartition of energy. From this, an equation constraining the accelerations is N m £ r , - • £ • = (). (2.3.22) j=i This equation is simply the time derivative of E q . (2.3.21). Following Evans [32], the constrained accelerations r , c , in terms of the unconstrained Newtonian accelerations, r'i" = F , / m can be written as f,-c = r - 7 - A - . (2.3.23) m Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 22 The constraint on the accelerations is then N m V r , - - ( f , - u - A ^ ) = 0, (2.3.24) implying A = t ? = l F : (2.3.25) Here, F,- is the total force acting on particle i. The constrained acceleration on the particle can be written as *• = M F - ( § f e ! * - ) ( 2 - 3 - 2 6 ) For particle rotations, the equation of constraint is if;i:a;,-u;t--{iV^r = 0, (2.3.27) where a>4- is the angular velocity if the ith particle in the body fixed axis frame, and I is the diagonal moment of inertia tensor. The number of rotational degrees of freedom / is 2 for a linear particle and 3 for a non-linear particle. The constraint on the angular acceleration is «? = <J?_Au;V, (2.3.28) N X ) I : w , - ( w i - A u ; i ) = 0, (2.3.29) t=i A = tit*1.""*'. (2.3.30) E,-=i I : « t W i Note that in the non-spherical case, the orientation of the particles is not introduced. Only the angular velocity, and the angular acceleration is needed. C h a p t e r 3 A S i m p l e M o d e l P o t e n t i a l for L i q u i d C r y s t a l M o l e c u l e s Liquid crystal molecules are generally elongated, and have a degree of flatness due to the presence of rigid aromatic rings [4]. The aromatic rings are linked so that they form a flat and rigid "backbone" for the molecule. In the nematic phase, the particles align, and the direction associated with the ordering is the director d. The orientation of the particles about the director is completely random, meaning that the "flat" faces of the molecules are not ordered. In model calculations, uniaxial nematic ordering is usually presumed to arise from particle elongation. Computer simulations of liquid crystals have generally used axially symmetric hard models [18], such as spherocylinders [22], cut spheres [23] and ellipsoids [19]. Continuous potentials such as the Gaussian overlap model, suggested by Berne and Pechukas [34], approximates molecular shape by a three dimensional Gaussian. This method gives continuous interaction potentials which are very convenient for molecular dynamics ( M D ) simulations of prolate and oblate particles. Also, extensions of this method are used to construct the so-called Gay-Berne potentials [35] which have been employed in recent simulations of liquid crystals [36-39]. Again, this approach has been limited to systems of axially symmetric particles. A n important exception is the Monte Carlo study of Allen [40] which showed that for some particle shapes (in particular ellipsoids with semi axes a : 6 : c w 1 : 3 : 10), a biaxial nematic phase exists and lies between uniaxial nematic and uniaxial nematic discotic phases. Allen's work focussed upon locating biaxial transitions in model systems 23 Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 24 already possessing uniaxial nematic order. The use of axially symmetric particles as a simple model for liquid crystals neglects the effects of particle flatness. It would be interesting to know if, and to what extent, introducing particle flatness can create uniaxial nematic order in systems where the particle elongation alone is insufficient to do so. In order to investigate this problem, we have generalized the original Gaussian overlap model. The new form derived here applies to ellipsoids with non-degenerate semi-axes. Molecular shape is approximated by a three dimensional Gaussian of general ellipsoidal symmetry. It is then assumed that the repulsive pair interaction of two molecules is proportional to the overlap volume of the Gaussians. This yields an interaction potential which depends upon the orientation and separation of the two Gaussians. The potential decays as a scaled Gaussian where it's maximum height can be varied while holding it's width constant; in the limit that the height of the overlap Gaussian becomes infinite the ellipsoids become "hard". Of course, as with the original Berne-Pechukas model, the present interaction is physically unrealistic, however, repulsive interactions are very short-ranged and we would not expect the exact nature of the decay to be of crucial importance. Obviously, this is the same assumption underlying the widespread use of hard-core models. The repulsive form presented here was specifically chosen to isolate the effects of particle shape on liquid crystal phases. The strength of the Gaussian overlap model lies in its ability to represent molecular shape and in it's computational convenience in M D simulations. The calculation of the forces and torques necessary in M D simulations is greatly simplified by the Gaussian form. A simple related expression which might be used to model short-range attractive interactions is also suggested. The effect of molecular biaxiality on liquid crystal phases is studied with a series of M D simulations that are carefully designed to determine the influence of molecular flatness. We begin with a system of axially symmetric particles which have no liquid Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 25 crystal phase whatsoever and then slowly flatten the particles keeping all other variables (i.e., the reduced density, reduced temperature and the remaining two aspect ratios) fixed. As the particles become flatter we observe a transition from an isotropic to a uniaxial nematic phase. This transition results entirely from increasing particle flatness and clearly demonstrates that molecular biaxiality can be an important factor in the formation of uniaxial phases. As the particles are made flatter still, a uniaxial to biaxial transition is also observed. The generalized Gaussian overlap model will be used in later chapters to study the effect of molecular shape on ferroelectric liquid crystal phases. 3.1 G e n e r a l i z e d G a u s s i a n O v e r l a p Integral In this section the Generalized Gaussian Overlap Integral is briefly described. A three dimensional Gaussian function centered at the point p with orientation (6,(f>,ip) can be defined as j(r) = 7 e s / 2 - s T : ( r -P) ( r -P) . (3.1.1) The parameter s determines the maximum of the Gaussian when r = p. When r lies on the surface of the ellipsoid defined by 2 T : ( r - p ) ( r - p ) = l , (3.1.2) then j ( r ) = 7 . (3.1-3) The value of 7 will be chosen after the overlap integral is calculated. The orientation and shape of the Gaussian in the space fixed frame are given by the second rank tensor T . It is defined by • T = B ' b B , (3.1.4a) Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 26 where b is the tensor / \ (3.1.4b) l/2o-l 0 0 0 l / 2 a 2 0 V o o i / 2 a 3 2 ; The parameters <Ti, a2, and <73 are the principle radii when the ellipsoid is cast in standard form (see Fig. 3.1). B is a rotation matrix (its transpose is denoted B*) which expressed in quaternions has the form [14] ( ql + q\ - q\ - ql 2(9192 - qoqz) \ 2 ( ^143 + q0q2) where 2(<?i<?2 + 9o9s) 2(5^3 - q0q2) ^ ql - q\ + ql-ql 2(9293 + 9091) 2(9293 - 9o9i) 9o - 9? - 9! + 9l ) 90 = cos ^0cos ^(<f) + ip) , 1 1 qx = sin -0 cos -((f) - ip) , 1 1 q2 = sin -6 sin-(<f> - xp) , 93 = cos ^ s i n ^(<f> + ip) . (3.1.4c) (3.1.4d) (3.1.4e) (3.1.4f) (3.1.4g) T h e choice of quaternions to express orientations simplifies the computations. T h e calculation of trigonometric functions is avoided in the rotation matrix, and the equations of motion have no discontinuities [14]. A review of quaternions when they are used to describe the orientation of rigid bodies is given in Appendix D. We next consider the overlap integral of two Gaussians, i and j , where i is centered at the origin and j at p. The overlap integral can be written as /co /•co rco r _ , m / w , . -CO J — OO J — OO Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 27 Figure 3.1: A sketch showing the molecular axes (x',y',z') and the principal radii <7i, <T2 and <r3 which characterize the ellipsoid model. Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 28 and some details of the integration are given in Appendix E . T h e result obtained can be expressed in the compact form /(P, q , q,) = ^ ( ^ 3 / V [ i + ( B B / 4 A - C ) : p p ] ? ( 3 L 6 ) V / I T ' + T ' I 5 where the scalar A and the elements of the vector B and of the tensor C are given explicitly in Appendix E . In order to associate E q . (3.1.6) with a repulsive pair potential, 7 2 must have the units of energy per unit volume. If we consider identical ellipsoids and choose 7 = 2V2i ( - f/ 4 (2 |b | ) 1 / 4 , (3.1.7) then the interaction energy, u(p,q 8 -,qj), is given by u(p, q,-, q,-) = 4 e / | b , " + b j , | e * [ i + ( B B / 4 A - C ) C T ] - ( 3 L 8 ) y/\Ti + Tj\ For this particular choice of 7 the energy is 4ee s when the ellipsoids are aligned and the intermolecular distance is zero. This is true regardless of the size or aspect ratios of the ellipsoids. The factor of four will scale the well depth when the attractive term is added, as discussed below. Equation (3.1.8) is a purely repulsive potential which can be used to model the hard core or "shape" of a molecule. However, it is possible to add simple attractive terms without loss of generality or the introduction of computational difficulties. For example, we may choose the attractive term « - t ( p , q , - , q , - ) = _ 4 £ V / l b ' + b ^ e S / 2 [ i + ( B B / 4 A - C ) : p p ] ; ( 3 L 9 ) ^/|T,- + T i | such that the total pair potential is given by u ( p , q,., q.) = 47|b' + b j l ( V [ i + ( B B / 4 A - C ) : P p ] _ e S / 2 [ i + ( B B / 4 ^ - C ) : p p ] \ _ { 3 A l Q ) 'lTi + T,-| v 1 Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 29 This attractive term is simply a negative Gaussian with a maximum of e3'2. T h e factor of four scales the well depth so that the minimum for two spherical Gaussians is —e, independent of the scale factor, and comparable with the well depth of a Lennard-Jones potential. Of course, this attractive potential has a Gaussian decay and will not reproduce the correct long-range behavior of the dispersion interactions. This attractive term is only given as a demonstration of how simple attractive interactions might be incorporated into the pair potential. For this study, we use only the repulsive form. The repulsive potential given in E q . (3.1.8) is continuous and it is possible to obtain analytical expressions for the forces and torques. This involves taking the gradient of the pair potential with respect to the intermolecular vector p and with respect to the eight quaternion parameters which define the orientation of the two ellipsoids. The procedure and explicit results are given in Appendix E . 3.2 M D Simulation Results: The Effect of Molecular Biaxiality on Liquid Crystal Phases If we take a as the reference unit of length and e as the.reference unit of energy, then a sys-tem of Gaussian ellipsoids can be characterized by the reduced density p* = 8N0-10-2&3/V (N is the number of particles and V is the volume) and the reduced temperature T * = ksT/s, (ksT is the Boltzmann constant times the absolute temperature). A l l distances are reduced by er, such that r* = r/cr. It is also useful to introduce the reduced energy U* = U/e, the reduced pressure P* = Pa3/e and the reduced t i m e i * = t^Je/o2m, where m is the mass of a particle. Suitable moments of inertia about the molecular axis (x', y', z') (see Fig. 3.1) ), Ia'a', can be calculated by assuming a mass density of the form £ e s / 2 - s b : r r ? where ( is a constant Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 30 and has units of mass per unit volume. One obtains L ai(a\ + erf) , (3.2.1a) y'y' ai(a2 + al) (3.2.1b) (3.2.1c) where Oi = C 2 \ / 2 V / 2 C r 1 C r 2 C r 3 ( 7 r / s ) 3 / 2 . (3.2.2) In M D simulations it is advantageous to have similar relaxation rates for translations and reorientations. That is, the decay of velocity-velocity and angular-velocity-angular-velocity autocorrelation functions should be similar [14]. This can be accomplished by adjusting (. Also, as usual, it is convenient to define reduced moments of inertia of the form I*,a, = Iaiai jma2. Of course, the equilibrium properties of the classical fluid will not depend on the choice of the moments of inertia. In M D simulations, the existence of liquid crystal phases is established by monitoring both dynamic and static observables. The mean square displacement, (|r(t) — r(0)| 2 ), was used as a diagnostic to determine if a particular system was liquid or solid. For all systems studied, the mean square displacement increased with time indicating fluid behavior (Fig. 3.2). Uniaxial and biaxial ordering was detected by calculating the second rank order parameters, Ql0 and Ql2, defined as [40,41] where the Euler angles describing molecular orientation are now defined with respect to the laboratory fixed reference frame (X, Y, Z) selected as described below. QQ0 is the usual uniaxial order parameter and (except for finite size effects) will be zero in the (3.2.3a) Q l 2 2 = (-(1 + cos 2 0) cos 2(f) cos 2ip — cos 9 sin 2(f) sin 2ip) (3.2.3b) Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 31 isotropic phase. Q22 is a biaxial order parameter which will be zero in isotropic and uniaxial phases but non-zero in a biaxial phase. Following Allen [40], we calculate the order parameters by first forming the ordering matrices QXX, QYY, and Q 2 Z with elements defined by Qlp = ^ E ^ V ^ - M . = 1 ,2,3, (3.2.4) (and similar expressions for QYY0 and Qzap) where X;, y,- and z; are unit vectors along the three symmetry axes of particle i. We next determine the largest eigenvalue associated with each of the three matrices. The molecular axis corresponding to the largest of these three is identified as the principal molecular z axis and the corresponding eigenvector becomes the laboratory fixed Z axis. The second largest and smallest eigenvalues are used to identify the molecular y and x axes, respectively. The associated eigenvectors are then orthogonalized together with Z in order construct the Y and X axis of the laboratory fixed frame. In terms of these axes, the order parameters have the cartesian forms Q200 = (Z • C T • Z) , (3.2.5a) Q222 = -(x • QXX • x + Y • qyy • Y - x • qyy • x - Y • qxx • Y) . (3.2.5b) 3 These averages can be easily evaluated in the simulation. It is important to note that the variable molecular x, y and z axes referred to in this paragraph are used only in the definition of the order parameters. They should not be confused with the fixed molecular coordinate system (x',y',z') shown in Fig. 3.1 . The effects of molecular flatness on liquid crystal phase transitions have been studied by carrying out M D simulations of selected repulsive Gaussian ellipsoids. The length 2<7i was varied from 0.5 to 3 while the remaining dimensions were held fixed. The a2 and o"3 radii were chosen such that the ellipsoid of revolution obtained when <7i = a2 Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 32 30 A C\2 J=L 10 v 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 . ' ' ' / -/ ,1 / / — / / -- / / -/ / // // // J / s / — / / ~ y \ i 1 1 1 . . . jr l l 1 l l 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 2 4 6 8 10 12 t * Figure 3.2: ( |r(f) - r(0) | 2 ) for 2oY : 2a2 : 2a3 = 2.75 : 3 : 9 (dash line) and 2cri : 2<72 : 2a3 = 1.5 : 3 : 9 (dotted line). These systems are isotropic and uniaxial nematic phases, respectively. Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 33 has 2<7i : 2cr2 : 2<r3 = 3 : 3 : 9 . These dimensions were selected because prolate hard ellipsoids of revolution with an aspect ratio of 1 : 3 are believed to have no stable liquid crystal phase [21]. This result was also confirmed with the present model. Thus, at fixed temperature and reduced density, any ordering found by varying <7i can be attributed to molecular flatness. A l l results reported are for the reduced parameters T* = 1.35 and p* = 0.8. T h e scale parameter s was set at 20 which gives a very "hard" potential. Some runs were also carried out with s values of 10 and 100 and very similar phase behavior was observed. We note that for hard ellipsoids with principal radii cr l 5 o2 and a 3 , the value of p* used would correspond to p/pcp = 0.8/^/2 — 0.5657, where pcp is the close packing density. Most calculations were carried out with 216 particles in a non-cubic simulation cell starting from a body-centered orthorhombic lattice. For systems with 2&i : 2cr2 : 2<r3 ranging from 1 : 3 : 9 to 2 : 3 : 9 , the initial particle array was 12 x 6 x 3. For systems with 2<7i > 2, the initial particle array was 6 x 12 x 3. Systems with 2ax < 1 were done with 360 particles to ensure that the simulation cell did not become too narrow in one direction. Possible system size and/or cell shape dependence was examined by carrying out simulations of several systems using 216, 288 and 360 particles. T h e results obtained were found to be very similar and there was no significant dependence upon system size. In all simulations, the reduced timestep St* = 0.0012 was used and the moment of inertia parameter, (, was adjusted to keep a\ = 0.5 (see E q . (3.2.2)). Typically, systems were equilibrated for 100,000 timesteps and this was followed by "production runs" of 100,000 timesteps. Standard deviations were estimated by dividing the production runs into ten equal blocks and assuming that the block averages were independent measurements. In Fig. 3.3, the order parameters QQ0 and Q\2 are plotted as a function of the particle thickness 2o\. Again, it should be noted that for the points shown in this plot, the lengths 2(j2 and 2cr3 are held fixed (at 3 and 9, respectively) as are the reduced temperature and Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 34 density. Thus Fig. 3.3 displays only the influence of changing width; as 2ci is varied from 3 to 0.5 we go from axially symmetric ellipsoids of revolution to flat, surfboard-like particles. Snapshots of various systems are shown in Figures 3.4, 3.5 and 3.6, and give a clear physical picture of the observed liquid crystal phases. We see from Fig. 3.3 that, as expected, when the aspect ratios are close to 3 : 3 : 9 the liquids are isotropic (see Fig. 3.4 ) and apart from small finite size effects both order parameters are zero. However, as 2<Ti is reduced the system quickly orders to form a uniaxial nematic phase (see Fig. 3.5) with the z' molecular axes (see Fig. 3.1) aligned along the director. This clearly shows that molecular biaxiality can have a very strong effect even upon the formation of a uniaxial phase. It is evident that as the particles become flatter the tendency to form a liquid crystal phase is significantly increased. At 2<Ti ~ 1.2, Q\2 increases sharply indicating the formation of a biaxial phase. Also, at 2<7i ~ 1.1, Qoo increases rather sharply as well indicating that the principal uniaxial director is now associated with the x' rather than the z' axes of the molecules. W i t h respect to Fig. 3.3, it should be noted that Q\\2 does not approach 1 as 2<7i is decreased. This is because as the particles become thinner the ordering of the z' molecular axes is diminished and, eventually, for 2<7i sufficiently small, we would expect to obtain a uniaxial discotic phase. 3.2.1 Molecular Biaxiality and the Uniaxial Nematic Phase Clearly, molecular biaxiality has a strong effect on the stability of the uniaxial nematic phase. However, the question remains, why does making a particle flatter in one direction result in ordering along another? In this section, we will study the uniaxial nematic phase in more detail. It would be interesting to know to what extent particle flatness effects the structure of the fluid. In Fig. 3.7, the radial distribution function g(r*) is plotted for biaxial, Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 35 1 0.8 0.6 -0.4 -0.2 h 0 1 1 ^1 1 II 1 1 1 1 1 1 1 1 1 i i i | i i i i j_ - T * S M * * at QfJ0> — 1 T Q 2 2 " X -r i 1 — * 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 T 0.5 1 1.5 2 2.5 3 2(7! Figure 3.3: Order parameters as a function of particle thickness at fixed reduced tem-perature and density. The squares are the uniaxial order parameter QQ0 and the circles are the biaxial order parameter Q\2. The error bars represent one standard deviation. Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 36 Figure 3.4: Snapshot of an instantaneous configuration for an isotropic phase with 2ox : 2o2 : 2a3 = 2.875 : 3 : 9. Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 37 Figure 3.5: Snapshot of an instantaneous configuration for a uniaxial nematic phase with 2ai : 2a2 : 2a3 = 2.125 : 3 : 9. Here, the director is along the z' axis (up and down on the page). Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 38 Figure 3.6: Snapshots of instantaneous configurations for a biaxial nematic phase wi th 2<7i : 2<r2 : 2<r3 = 0.5 : 3 : 9. Here, the directors are along x' and z'. The snapshots are taken along the axes perpendicular to the z' director. Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 39 (2<Ti : 2<72 : 2<73 = 1 : 3 : 9 ) uniaxial, (2<7i : 2<72 : 2<r3 = 2 : 3 : 9 ) and isotropic (2ci : 2 a 2 : 2<73 = 3 : 3 : 9 ) phases. In the biaxial nematic there are weak shoulders at r* ~ 1.7 and 3.7, indicating that there is a enhanced probability of finding particles with their flat faces "pressed" together. This.is evident from the snapshot of the biaxial nematic in Fig. 3.6. However, in the uniaxial nematic, there is a peak in g(r*) at r* m 3.2, indicating that the orientations of the particles about the long z axes are completely averaged, and few particles are "pressed" with their flat faces together. Short-ranged orientational correlations can be detected by calculating [18] <P2(x'1(0)-x'2(r))>, (3.2.6) where P2(x) is the second Legendre polynomial and we are considering an average over pairs of particles (ie. particle "1" and particle "2"). This function gives an indication of the average value of P2 as a function of the intermolecular distance. Figure 3.8 shows < P2(x[(0) • x'2(r*)) > for biaxial, 2<7i : 2cr2 : 2<T3 = 1 : 3 : 9 , uniaxial nematic 2<7i : 2<r2 : 2<T3 = 2 : 3 : 9 , and isotropic 2<7i : 2<r2 : 2cr3 = 3 : 3 : 9 systems. We see that in the biaxial phase, this function decays to a large non-zero constant, indicating that there is long-ranged biaxial correlations. In the uniaxial nematic, there is a small peak at r* ~ 2cri, but any biaxial correlations weak and short-ranged. Particle reorientation about the long axis can be examined by the time correlation function < x'(t*) • x'(0) >. If the decay is nearly exponential at long times (Debye-like) [12], then reorientation occurs through many small uncorrelated steps. This type of decay also indicates that the particles are not "spinning" about the z axis, nor are they "wagging" back and forth. Free rotor, or "spinning" reorientation would result in oscillations in < x'(t*) • x'(0) >. Likewise, "wagging", or back and forth motion would appear as < x'(t*) • x'(0) > decaying to a non-zero constant. Fig. 3.9 shows < x (t*) • x (0) > and ln(< x (t*) • x (0) >) for two uniaxial nematic Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 40 Figure 3.7: <7°00(r*) for biaxial 2ox : 2cr2 • 2cr3 = 1 : 3 : 9 (solid), uniaxial nematic 2<7i : 2<r2 : 2o3 = 2 : 3 : 9 (short dash), and isotropic 2o\ : 2<72 : 2u3 = 3 : 3 : 9 (long dash) systems. Figure 3.8: < P2(x'1(0) • x'2(r*)) > for biaxial 2<7i : 2a2 : 2a3 = 1 : 3 : 9 (solid), uniaxial nematic 2<7i : 2cr2 : 2a3 = 2 : 3 : 9 (short dash), and isotropic 2o~\ : 2o~2 : 2<r3 = 3 : 3 : 9 (long dash) systems. Figure 3.9: < x'(t*) • x'(0) > for a system with 2<TI : 2a2 : 2a3 = 2 : 3 : 9 (solid), and a system with 2a\ : 2a2 : 2a3 = 1.5 : 3 : 9 (short dash). The inset is In < x (t*) • x (0) >. Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 43 systems. One is a system near the isotropic to uniaxial nematic transition, corresponding to 2<7i : 2o"2 : 2cr3 — 2 : 3 : 9 , and the other is near the uniaxial to biaxial nematic transition with 2o~i : 2o~2 : 2er3 — 1.5 : 3 : 9. If the reorientational correlation is given by the simple exponential < x'(t*) • x'(0) >= exp[—f/r*] , where r * is the relaxation time, then for the 2at : 2cr2 • 2cr3 = 2 : 3 : 9 (solid line) and 2ox : 2a2 : 2cr3 = 1.5 : 3 : 9 (dash line) systems, r * « 3.5 and r * 8, respectively. Although particle reorientation in both of these systems is reasonably Debye-like at long times, we see that the correlation time increases as the particle is made flatter. Obviously, at the biaxial-uniaxial nematic transition, the correlation time would diverge as this motion freezes out. We can summarize the previous results. First, from the local fluid structure, there is evidence that the orientation of the particles about the long z axis is completely random. This is consistent with the biaxial order parameter calculations (Q\2 ~ 0). There are "biaxial" orientational correlations when the intermolecular distance is near 2<rl5 but this structure decays quickly. The reorientational motion of the particle about the z axis is complicated, and occurs through many small, uncorrelated reorientations (Debye-like). T h e role of particle flatness on the stability of the uniaxial nematic is not obvious from these observations, and is in fact quite subtle. However, we can offer a simple explanation. Consider the argument used to explain isotropic to nematic phase transitions in hard axially symmetric liquid crystal models [19,23,42] as well as in the hard ellipse fluid [43]. In systems that are purely repulsive, the isotropic to nematic phase transition is entropy driven. Under constant volume con-ditions, the transition will happen at a specific density. At this density, the isotropic to nematic transition occurs because the nematic phase has greater entropy than the isotropic phase. Axially symmetric particles will have an isotropic to nematic phase transition if their shape is sufficiently oblate or prolate. At the phase transition, some Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 44 "orientational" freedom is lost when the particles align. However, there must be a corre-sponding gain in translational freedom, so that the entropy of the nematic is greater than the isotropic. B y electing to "line up" the particles can more efficiently avoid collisions and explore the available translational space. This arguement can also be used to explain the isotropic to uniaxial nematic transition driven by increasing particle flatness for a fixed elongation. For the aspect ratios chosen in these studies, the axially symmetric system had no observed nematic phase. This is because if the particles were to align, they would be too "fat" to easily slide past each other, and hence would gain no translational freedom. However, this is not the case in the uniaxial nematic when the particles have a degree of flatness. B y making the particles slightly flatter in one direction, there is a greater probability depending on their relative orientation, that they can "slip" past each other without colliding. Molecular biaxiality allows the particles more translational freedom, and occurs despite the fact that there is no orientational ordering associated with the x axis. As an example, consider the situation where two particles in the nematic phase are nearly aligned (with respect to their z axes) and are end-to-end. T h e orientation of the particles about their corresponding z axes is completely random. One particle is moving towards the other. If their relative orientation about their corresponding z axes just happens to be such that they do not collide, then they will be able to slip past each other. This event would not happen if the particles were axially symmetric. T h e "chance" that the particles will not collide increases with increasing molecular biaxiality, that is, there are more possible orientations about the z axes where they can avoid collisions. A transition will be observed when the molecular biaxiality is great enough that sufficient numbers of aligned particles can slip past each other without colliding. The results of these simulations indicate that particles with aspect ratios around 2<7i : 2a2 : 2a3 = 2 : 3 : 9 have sufficient molecular biaxiality to form a stable nematic phase. Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 45 3.3 Summary In this Chapter, the Gaussian overlap model has been extended to include ellipsoidal particles with non-degenerate semi-axes. This was a valuable extension to the original Gaussian overlap model [34], since the original model was restricted to particles that were axially symmetric. We have studied the effect of particle flatness on uniaxial liquid crystal phases. These are the first M D simulations of biaxial ellipsoidal particles, and they clearly demonstrate the role of particle flatness on uniaxial liquid crystal phases. Detailed constant tempera-ture M D simulations clearly show that even small degrees of molecular flatness (that is, flattening a prolate ellipsoid of revolution into a surfboard shape) can result in a uniaxial nematic phase where the director is associated with a different molecular axis. That is, flatness in one direction can stabilize order in another. This result is particularly relevant to real liquid crystals, which generally tend to be somwhat flat. In the systems selected in this study, molecular flatness is crucial in order to stabilize the liquid crystal phase over the isotropic phase. For example, at a reduced density where particles characterized by the aspect ratios 3 : 3 : 9 are isotropic, those with ratios 2 : 3 : 9 form a strongly ordered uniaxial nematic phase. As the particles become flatter, the expected biaxial nematic phase is obtained. For the present model, particles with the aspect ratios 1 : 3 : 9 formed a fully developed biaxial phase. Chapter 4 Ferroelectric Phases of Amorphous Dipolar Systems Recently, computer simulations of fluids of strongly interacting dipolar spheres [9,10] were found to possess a stable ferroelectric nematic liquid crystal phase. The phase was nematic in that it exhibited long ranged orientational order, but only short-ranged spatial structure. T h e sample was polarized, but was fluid-like in terms of particle translations. The local spatial correlations were similar to those found in the ferroelectric tetragonal I lattice and is the structure believed to be the ground state for dense dipolar sphere systems [44]. The existence of the ferroelectric nematic phase was believed to be due to these specific short-ranged tetragonal I correlations [9,10]. The translational mobility of the fluid allowed specific short-ranged correlations to build up, and the system sponta-neously polarized. In this system, short-ranged spatial correlations are believed to play a large role in the stability of the ferroelectric phase. In a recent paper, Zhang and W i d o m [45] proposed a mean field theory that predicted ferroelectric phases in dipolar systems that lacked any specific spatial correlations. They considered amorphous solids of dipolar hard spheres where the particles were free to rotate, but were frozen in random sites. Here, random refers to the spatial structure of the material (i.e., the radial distribution function g(r) = 1 for r > a, where a is the diameter of the sphere). The prediction of ferroelectric phases in dipolar systems that lack any specific spatial correlations suggests that short-ranged structure may not be necessary for ferroelectric phase formation. 46 Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 47 However, experimental systems of near spherical Fe30\ particles in a frozen non-magnetic solvent [46-48] do not exhibit ordered phases. In these systems the dipolar interaction is magnetostatic. The particles are translationally frozen in random locations, but above a certain temperature they can rotate freely. As the system is cooled, the particles do not orientationally order, but rotationally freeze in random directions. Th e random spatial structure is believed to be responsible for the disordered low temperature phase [48]. The random environment about a particle creates a strong random field, and the direction of the field is different for different particles. At low temperatures, the moments of the particles will align with the field, resulting in an orientationally disordered phase. The experimental results seem to contradict Zhang and Widom's theory, although it is possible that the experimental dipole density was too low, or impurities in the sample were present. The role of short-ranged spatial correlations on ferroelectric phase formation is still not well understood. The question of whether or not dense spatially disordered dipolar materials can have ferroelectric phases is still open. In this Chapter, we will identify and study the fundamental forces that both promote and destroy ferroelectric order in dipolar systems. T h e effect of spatial disorder on the phase behavior of dense dipolar systems will be studied with M D and M C simulations. In the first section, some details of the model and the simulations will be discussed. The second section concerns ferroelectric ordering in amorphous frozen dipolar systems. The third section introduces a new simulation technique devised to study orientational order in spatially random media. The results of the "dynamically disordered" simulations offer significant insight into the observed behavior of both the dipolar fluid [9,10] and the dipolar amorphous solid [45-48]. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 48 4.1 M o d e l a n d S i m u l a t i o n D e t a i l s . We consider dipolar soft spheres defined by the pair potential u(12) = ug B(12) +uDD(l2), (4.1.1) where u s s(12) is the soft sphere potential given in E q . (2.2.5) and uDD(12) = - 3 ( > i • r ) ( / i 2 • r ) / r 5 + m • fi2/r3 (4.1.2) is the dipole-dipole interaction as previously given by E q . (2.2.7). Under periodic bound-ary conditions, the dipole-dipole pair potential is replaced by the effective interaction given by E q . (2.3.7). For the dipolar systems that will be studied, the dielectric con-stants are sufficiently large that the reaction field contribution is essentially independent of the exact value [10,25], therefore the dielectric constant of the surroundings was set to e' = oo. We will study systems where the dipole vector has one, two and three components. Systems with one and two components are usually referred to as the Ising and X Y models, respectively, and for notational simplicity we shall refer to the three component dipole as the X Y Z model. A l l of these systems have three spatial degrees of freedom, that is the position of a particle is specified by a three dimensional vector. However, the rotational degrees of freedom of a dipole is varied. The Ising system consists of dipoles in a 3-dimensional space that can either point in the -\-z or —z direction. The X Y model has dipoles in a 3-dimensional space that are constrained to rotate in the X Y plane. Although the components of the dipole vector are varied, the dipolar interaction is still given by E q . (2.2.7). The existence of a ferroelectric phase can be detected by calculating the average polarization per particle < P > defined as < P > = 1 / / V ( | f > r d | ) , (4.1.3) Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 49 where d is a unit vector in the direction of the total instantaneous moment, M = J2i=i A*;> and N is the number of particles in the system. In ferroelectric systems, < P > is non-zero , while in disordered phases < P > is near zero. There is a possibility that the system could "orientationally freeze" at low tempera-tures. Orientationally frozen, or "spin glass" phases can be detected by calculating the root-mean-square spin length < Srms > [49] given by < srms >= iv-1E(< •»-•>•< m< >)]1/2> (4-L4) i where < m , >= r - 1 &(>'), (4.1.5) T ' =0 and r is the number of M D timesteps, or M C sweeps (N moves). Briefly, systems with both < P > and < Srms > non-zero will be ferroelectric. Systems with < P 0 but < Srms > non-zero will not be polarized, but the particles are orientationally frozen. When both < P > and < Srms > are near zero, the system is rotationally "fluid-like". Dipolar soft-sphere systems can be characterized by the reduced density, p* = Nas/V, the reduced temperature, T* = kBT/e, and the reduced dipole moment, p* = (p2/ecr3)1/2. For this study, the reduced dipole moment and reduced density were constant at p* = 4 and p* = 0.8. This density is well within the range where Zhang and W i d o m predict a ferroelectric phase. For example, the lowest density of the ferroelectric phase is predicted to be p* — 0.31 for the Ising system and p* = 0.55 for the X Y Z model [45]. For the present case with p.* = 4 and p* — 0.8, Zhang and W i d o m predict a ferroelectric phase for the Ising case if T* < 41.6 and for the X Y Z model if T * < 6.4. 4.2 R e s u l t s for R a n d o m l y F r o z e n S y s t e m s We have performed simulations where the amorphous spatial structure was obtained from a "random" configuration of soft spheres. The configuration was generated from a M D Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 50 simulation of a soft sphere fluid at T* = 10.5 and p* = 0.8. At the end of the simulation the final configuration was used as the frozen system. Unfortunately, it is impossible to have a truly "random" and structureless spatial configuration at this density, that is, it is impossible to have g(r) = 1 for r > a [45]. The radial distribution function, g(r) for p* = 0.8 and T* = 10.5 is shown in Fig. 4.1. Clearly, the local structure of the soft sphere fluid is weak and very short-ranged. The configuration of soft spheres determines the positions for the point dipoles. Ini-tially, the dipoles were placed with random orientations at the centers of the soft spheres. Constant temperature M D simulations were done for the X Y Z and X Y models, and M C simulations were used for the Ising model, where one M C sweep (N attempted moves) consisted of randomly flipping the orientations of the dipoles. Fig. 4.2 shows the average polarization versus 1/N for the Ising, X Y , and X Y Z models. Systems sizes ranged from N = 108 to N = 864 particles. The polarization values that are plotted were obtained at the lowest temperature where equilibrium could be achieved, denoted as T ^ n . The values of T ^ n for the Ising, X Y and X Y Z models are 10.0, 4.0 and 3.5 respectively. Below these temperatures, simulations that were started from perfectly aligned and random states did not converge to the same result. Lower values of T ^ i n could possibly be obtained with longer simulations. However, the temperatures that were attained were within the range that Zhang and W i d o m predicted a ferroelectric phase. The Ising model at T ^ i n = 10.0 is almost completely polarized and shows little or no system size dependence. In Fig. 4.6 a detailed polarization versus reduced temperature plot is shown for the p* = 0.8, N = 256 frozen Ising system. We see that ferroelectric order develops spontaneously in this system-at T* K, 25. This transition temperature is much lower than the value predicted by Zhang and W i d o m , where they predict a transition at T* = 41.6. Returning to Fig. 4.2, we see that the X Y model at T^n = 4.0 and the X Y Z model Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 51 1.5 tad 1 -0.5 0 i—i—i—r i—i—r i—i—r i—r j i i L j i L J I I L J L 0 1 Figure 4.1: The radial distribution function g(r*) for soft spheres at p* = 0.8 and T* = 10.5. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 52 1.2 fl • r - t . s "X P-i V 0.8 0.6 0.4 -0.2 -0 1 I I I I I I I I I I I I I I I I r L- i I i i I i i i i i i i i i i i i i i 0 0.2 0.4 0.6 100/N 0.8 1 Figure 4.2: The polarization < P > at T^n vs 100/N for the randomly frozen Ising (circles), X Y (squares) and X Y Z (triangles). Results are included for 108 ( X Y Z only), 256, 500 and 864 particles. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 53 at T ^ i n = 3.5 show significant polarization for the N — 108 and TV = 256 systems. However, the polarization decreases monotonously with increasing system size. Indeed, it appears to approach zero in the thermodynamic limit. The observed polarization of the X Y model is strongly dependent on system size, decreasing from ~ 0.49 for iV = 256 to ~ 0.10 for TV = 864. The polarization for X Y Z model shows a similar, although not as pronounced system size dependence. A n indication of the behavior of the (N = 256) X Y Z model at T ^ n = 3.5 can be seen in Fig. 4.3, where the root-mean-square spin length < Srms > is plotted versus the reduced temperature T*. The non-zero value of < Srms > at T^n suggests that the system is orientationally frozen at this reduced temperature. Furthermore, < Srms > at ^ m i n s n o w s little or no system size dependence (eg. for N = 256, < S r m s >~ 0.37 and for N = 864, < Srms 0.35). To summarize, we find no evidence of ferroelectric states in the thermodynamic limit for the randomly frozen X Y and X Y Z models. This clearly disagrees with the calcula-tions of Zhang and W i d o m , where they predict that the X Y Z model will have a stable ferroelectric phase at T* = 3.5. The Ising model has a ferroelectric phase, however spon-taneous polarization was observed at T* w 25, which is much lower than the transition temperature predicted by Zhang and Widom. 4.3 Decoupled MD Simulation Technique In the simulations described above, a configuration was selected from an M D soft sphere simulation at p* = 0.8 and T* = 10.5, and was then used as a "typical" randomly frozen system. Point dipoles were embedded at the centers of the soft spheres and rotational M D or M C simulations were performed. Ideally, many randomly frozen configurations should be used to obtain accurate values of the average polarization < P >. However, Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 54 Figure 4.3: The root-mean-square spin length, < Srms >, for the N = 256 frozen X Y Z model. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 55 this would result in very long simulations. As an alternative approach to studying orientational order in random media, we have developed a decoupled M D simulation technique. Here, we consider a dynamic random substrate instead of a frozen system. The dipoles are no longer embedded in randomly frozen soft spheres. Instead they are embedded in a soft sphere underlying substrate that moves, but also lacks any specific spatial correlations. The structure of the substrate is not affected by the dipolar interactions. The "equilibrium" state of the dipoles, however, will depend on the underlying motion of the substrate. This model is similar in spirit to those used in recent studies of non-equilibrium phase transitions in magnetic systems subject to Levy flights [50]. It is fundamentally different from dipolar fluid simulations in that the spatial structure is not affected by dipole-dipole interactions, but it is also not an amorphous solid simulation since the particles do move. W i t h this technique, we not only avoid the sampling problem and the system size dependence of the frozen systems, but we can also gain insight into the role of short-ranged spatial correlations on phase behavior. B y controlling the rate that the substrate moves relative to dipolar reorientations, we can effectively "turn on" or "turn off" the specific details of the random spatial structure. The implementation is straightforward. The force between two particles is and the torque between two particles comes from the dipole-dipole interaction. Th e spatial structure of the soft sphere fluid is then only determined by the soft sphere part of the pair potential. The translational and rotational temperatures are decoupled so that the underlying substrate is not affected by changes in the rotational temperature. Finally, the rate that the substrate moves relative to the dipolar reorientations can be varied by changing the mass of the particle. Obviously, adjusting the mass of the particle f (12) V r u s s ( 1 2 ) , Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 56 will have no effect on the equilibrium structure of the soft sphere fluid. For large masses, the substrate changes slowly relative to dipole reorientations. A n extrapolation to the infinite mass case would give results of a randomly frozen system. For light masses, the substrate moves rapidly relative to dipole reorientations. At some finite mass, the substrate motion is so rapid that the dipoles cannot react, and a mean field like limit is reached. Increasing the mass effectively "turns on" the specific structural details of the environment surrounding a dipole. The moving substrate is a means of simulating dipolar systems in a dynamically random media that lacks any specific spatial correlations. T h e idea with the decoupled technique is to study the behavior of intermediate mass systems, then extrapolate the results to the infinite mass limit. In decoupled M D simulations, it is convenient to introduce the reduced mass m* = m/m', where m' is a reference mass defined such that the reduced simulation timestep 8t* = (e/m'a2)1/28t = 0.00125. (4.3.2) In all the decoupled simulation results that follow, the underlying substrate is a soft sphere fluid at p* = 0.8 and T* = 10.5. In Fig. 4.4, we have plotted < P > versus T* (rotational) for the X Y Z model. Systems with particle masses m* = 1, 5 and 10 are included. Clearly, spontaneous polarization develops for all systems and the temperature at which < P > begins to rise decreases with increasing mass. For m* = 5, the transition is T * « 0.55, and for m* = 10, the transition is T* « 0.25. Larger masses were investigated up to m* = 20, but they were all disordered above T* = 0.1. Slow convergence below T* = 0.1 made calculations for very large masses difficult. Systems with N = 108, 256, and 500 particles had similar results, and there was no system size dependence on the polarization. The results of Fig. 4.4 indicate that for any finite mass there will be some temperature where the system will spontaneously polarize, but as the mass becomes large, the transition temperature Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 57 will approach zero. The dynamically decoupled X Y model behaves much like the X Y Z system. In Fig. 4.5, we see that the three different masses, m* = 1, 5, 10 spontaneously polarize with tran-sition temperatures at T* « 6, 2, and 1.8. Again, the transition temperature decreases with increasing mass. The Ising model is more suited to M C simulations where the dipoles can be randomly flipped in opposite directions. Therefore, a suitable M C scheme was devised that allowed the substrate to move independently of the Ising dipoles. This involved combining a soft sphere M D simulation with a M C Ising dipole simulation. The substrate evolved as in the X Y and X Y Z decoupled cases, however after each M D timestep, one M C cycle (N attempted Ising flips) was performed. The behavior of the dynamically decoupled Ising model is different from both the X Y and X Y Z models. Fig. 4.6 shows polarization versus reduced temperature results for m* = 1, m* = 5, and the randomly frozen system (m* — oo). We see that the ordering behavior of the Ising model is essentially independent of mass and that the results for the dynamically disordered systems lie very close to those for the randomly frozen case. Reduced heat capacities per particle Cv = Cv/NkB, obtained by numerically differentiating the average dipolar energy with respect to the rotational temperature are also shown in Fig. 4.6. The randomly frozen and m* = 5 results are very similar and indicate a phase transition at T* « 25. The dependence of the ferroelectric transition temperature, denoted as TF, on particle mass m* for the X Y and X Y Z models is shown in Fig. 4.7. The transition temperatures were estimated from heat capacities obtained from numerical differentiation of the dipolar energy per particle. Results for the Ising system are not plotted because the transition temperature is essentially independent of the mass. As the mass increases the transition temperature drops for both the X Y and X Y Z models. As noted earlier, for large masses Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 58 Figure 4.4: < P > vs T * (rotational) for dynamically random X Y Z systems. The squares, triangles and circles are for m* = 1, 5 and 10, respectively. The error bars represent one estimated standard deviation. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 59 Figure 4.5: < P > vs T * (rotational) for dynamically random X Y systems. The squares, triangles and circles are for m* = 1, 5 and 10, respectively. The error bars represent one estimated standard deviation. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 60 1 -0.8 -A 0.6 \-V 0.2 -0 -"i—i—i—i—I—i—i—i—i—I—i—i—i—i—I—r 0.8 I I I I I I I I I | ^ ^ | I I I I I M * 0.6 E-0.4 h u 0 4 t-0.2 0 1111 m 11111111 0 10 20 30 40 r i i i i i i i i i i i o 10 20 J I I I L 30 Figure 4.6: < P > vs T* (rotational) for the Ising model. Results are shown for dy-namically random systems with m* = 1 (squares) and m* = 5 (triangles) and for the randomly frozen case (circles). The heat capacities per particle, Cv/N, are plotted vs T* (rotational) in the inset. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 61 and low rotational temperatures convergence becomes prohibitively slow, but it seems reasonable to assume that the graph would simply continue with the transition temper-ature approaching zero in the infinite mass limit. This is consistent with the fact that we did not find a ferroelectric phase for randomly frozen systems at finite temperatures. The results for the randomly frozen and dynamically decoupled systems can be inter-preted as follows. It is useful to divide the local field E i o c a i experienced by a particle in an infinite medium into two parts such that, E l o c a i = R + E , (4.3.3) where R is a reaction field contribution [25,51] and E is everything else. The reaction field favors ferroelectric order. It arises from the long-range nature of the dipolar interactions and is independent of the spatial correlations. The remaining contribution E is sensitive to positional correlations and may or may not favor ferroelectric order. Thus if R dominates a ferroelectric phase is to be expected, but if E is important the existence or nonexistence of a ferroelectric phase will depend on the details of the spatial correlations. This simple picture allows us to rationalize the various systems considered. For fully coupled dipolar fluids [9,10] the short-range spatial correlations can adjust and become tetragonal I like. Consequently, the remaining contribution to the local field E adjusts in order to favor (or at least allow) a ferroelectric state. O n the other hand, for randomly frozen systems the spatial correlations cannot adjust and at equilibrium E dominates over R . The random spatial correlations, and hence E , inhibits the formation of a ferroelectric phase except for the Ising case. Apparently, the discrete nature of the Ising model makes it much less susceptible to the development of disordering fields than are the X Y and X Y Z systems. The reaction field dominates in the Ising system and gives rise to a ferroelectric state. The behavior of the dynamically disordered systems can also be understood. If the Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 62 10 E-i 8 -6 -4 -0 i i i i i i i i i i i i i i i i i i i i i O 1 0 . 5 P I I f | I I I Hj I I I I | I I I LJ 2 . 5 E-2 I" / ? 1 . 5 =- > - i i i I i i i i I i i i i I i i i -1 1 . 5 2 2 . 5 3 . . . - - " i i 0 0.2 0.4 0.6 0.8 (1/rn) 1 Figure 4.7: The mass dependence of the disordered to ferroelectric transition temperature TF (rotational). The squares and triangles are results for the X Y and X Y Z models, respectively. The values of T-p were obtained from plots of the heat capacity, Cv/N, vs T* (rotational) and a typical example is shown in the inset. The error bars represent estimated uncertainties in the peak position. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 63 translational diffusion of the substrate is sufficiently fast compared to the dipolar reorien-tational time, then a ferroelectric phase is possible. The rapid movement of the particles l imits the extent that spatially dependent frustrating correlations can bui ld up, allowing R to prevail over E . As the mass of a substrate particle is decreased, the translational motion becomes faster and the above condition is met at higher and higher rotational temperatures. Thus the observed transition temperatures increase wi th decreasing mass. As the substrate particles become sufficiently light, the transition temperature is deter-mined only by the reaction field and hence becomes independent of mass. In fact, m* = 1 gives essentially this l imit ing behavior and reducing the mass further has l i t t le effect on the transition temperature. 4.4 A Mean Field Theory: Isolating the Effect of Long-Ranged Interactions In this section we present a simple mean field theory where a dipole in a spherical cavity interacts only with the reaction field. We wi l l begin with a brief summary of the details of the reaction field. A more thorough discussion can be found in Appendix 3. From this simple theory, we can calculate the reaction field contribution to the total energy per particle. This theoretical result wi l l be compared to the reaction field contribution for the same model (Ising, X Y and X Y Z ) from the simulations. The results from the simulations gives quantitative information on the role of long-ranged and short-ranged interactions on the stability of ferroelectric phases. The reaction field R at a point dipole in a sphere of radius r c , is due to the polarization of the dielectric beyond the sphere. We wi l l restrict ourselves to the case that the dielectric constant of the surroundings, e = oo. If there are N dipoles in the sphere, and the total moment of the sphere is M = f>, (4.4.1) t=i Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 64 then the reaction field is given by [14,52] R = ^ f v l , (4.4.2) '"c where the direction of reaction field is along the total moment of the sample. A point dipole // , within the sphere wi l l interact with the reaction field, and its interaction energy is given by [14] ut = -^pt-K. (4.4.3) Clearly, the low energy state is found when the dipole aligns with the reaction field. The factor of a half takes into account the work of polarizing the surrounding continuum [51,52]. The total moment M can be written as M = Md, (4.4.4) where M = |M|, and d is a unit vector in the direction of M. If the instantaneous polarization P is given by TV E M«: • D (4.4.5) then it is easy to show that the total instantaneous moment is given by M = NPd. (4.4.6) If the N dipoles are contained within a sphere of radius rc, then N = (4.4.7) where p is the number density. The total moment can be written as M = ^ J ^ P d . (4.4.8) The reaction field is then R = ^ P d , (4.4.9) Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 65 and the i particle's instantaneous energy is u% = -^-fii-d. (4.4.10) For the N dipoles within the sphere, the total energy is U — J2iLi uii then 2ir N U = --Pp2P\Z(ir<l- (4-4.11) 6 i = i \i p,i • d is written as cos (9t, then the average value of P for the X Y Z model in the canonical ensemble can be written as < P > _ K C O S 0 1 e x P ( ~ ^ U l ) s i n 9 l d d l (4 4 12) JQ exp(—Bu\) sin 9id9i W i t h the mean field approximation, we replace the instantaneous polarization P with the average polarization < P > resulting in = ^ c o s - ^ e x p f l ^ / V < P > cosf l i )s inMfl i ( A 4 ^ < > _ J ^ e x p ^ / V < P > c o s 0 i ) s i n M 0 i ' This can easily be evaluated, to obtain < P >= ^ ^I'-'V'. (4-4.14) a^e37 — e x) where 2 7 T / V < P > x = 3kBT For the X Y model, the average polarization < P > is (4.4.15) n ffcOsMxptf p/3p2 <P> COS 91)dei %exp(^p8p2 < P > cos 8 ^ ' which cannot be easily evaluated analytically, and therefore, is calculated numerically. For the Ising case, with the orientation of the dipole parallel to the z axis, the energy for the ith particle can be written as 27T u i + = -—pp2 <P> (4.4.17a) Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 66 or 2TT = + ^-pp2 <P>, (4.4.17b) depending on whether the dipole is directed along ( —) or against (+) R. For a system of TV particles, the configurational integral is N £ e x p ( - / 3 [ / ) = n(exp(-/3«,-+) + exp(-/3«,-_)), (4.4.18) j = i t=i where N */ = £ > , - . (4.4.19) t=i The number of possible states for an Ising system of iV particles is easily shown to be 2N. The average polarization can be written as < f > 4 ^ f f . ^ ; ) ^ » , ( 4 , . 2 0 ) N E , = i exp(-/9t/) which is written in the form exp(-/3u+)+exp(-/3u_)' 1 J 2ir = ^fPV2 <P>, (4.4.21b) Figure 4.8 shows the polarization per particle versus the reduced temperature for the X Y Z , X Y and Ising models. The reduced density and dipole moment are p* = 0.8 and p* = 4 respectively. A l l three models are predicted to have an ordered phase, with the Ising model having a transition at T* « 24, the X Y model at T* « 12, and the X Y Z model at T* « 8. The driving force for this transition is simply the reaction field due to the polarization of the surrounding dielectric. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 67 Figure 4.8: < P > vs T* from the reaction field, as found from a simple mean field theory. The short dash, solid and long dash lines are for the X Y Z , X Y , and Ising models, respectively. Th e circles are the frozen TV = 256 Ising results. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 68 The average energy per particle, due to its interaction with the reaction field can be found from E q . (4.4.11) and is given by -2 T T <URF>/N = -pp1 <P>' (4.4.22) The heat capacity per particle Cv/N can be found from d<URF> IN d<P> Cv/N dT 47T PP2 <P> dT (4.4.23) For the X Y Z model we have, d< P > dT <P> T 1 - X :i-x) 3kBT 2TTPH2 (4.4.24) w here X x2(ex + e~x)2 - ((ex - e x2(ex — e~x)2 -x\2 (4.4.25) For the Ising case, a similar expression is found but with X given by X (ex - ^-x\2 (4.4.26) (ex + e~x)2 The heat capacity for X Y model would have to be found from numerical differentiation of the energy versus temperature. Figures 4.9 and 4.10 show the reaction field contribution to the reduced energy per particle < URF > /N as a function of reduced temperature for the Ising and X Y Z models, respectively. The X Y model would have behavior similar to the X Y Z model. For both the X Y Z and Ising case, we see that the reduced heat capacity per particle Cv = Cv/N kg is continuous across the transition, indicating that the transition is not first order. Of course, in the frozen systems this transition was not observed for the X Y Z or X Y models. However, for the Ising model, the transition temperature as predicted from the mean field Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 69 Figure 4.9: The reaction field contribution to the energy per particle < URF > /N for the Ising model, as calculated from simple mean field theory. The inset is the reduced heat capacity per particle Cy/N. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 70 Figure 4.10: The reaction field contribution to the energy < URF > /N for the X Y Z model, as calculated from simple mean field theory. The inset is the reduced heat capacity per particle Cv /N. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 71 theory is very close to the transition temperatures observed from the simulations (both dynamically decoupled and frozen). For a particular system ( X Y , X Y Z , or Ising), the reaction field contribution to the reduced energy per particle gives an indication as to whether the long-range (reaction field) or the short-range structural interactions are dominant. The reaction field con-tribution to the energy can easily be calculated in computer simulations by using E q . (4.4.22). As in the previous section, where we "conveniently" separated the local field F J ' o c a / into E/ o c o ( = R + E , we can separate the total energy per particle into a reaction field contribution < URF > and the rest, < UST >•> such that < U >=< U R F > + < UST > • (4.4.27) To a large extent, < UST > is dependent on the local spatial structure surrounding a particle, hence ST refers to "structural". In Figs. 4.11, 4.12, and 4.13, the total energy per particle < U* > /N, the reaction field contribution, < URF > /N and < U*ST > jN are plotted versus the T*. In each plot, we consider the mass dependence for a particular model ( X Y Z , X Y , and Ising respectively). From these figures, we note that for all three models, < URF > jN dominates in ferroelectric phases, and < UST > /N dominates in disordered phases. In fact, in ferroelectric phases, < UST > /N rises (becomes less negative), indicating that structurally dependent interactions do not favor ferroelectric order. For the Ising model, (Fig. 4.13) we see that the contributions to the energy are independent of the mass. Even for large masses (i.e., approaching the amorphous solid), the largest contribution to the energy in the ferroelectric phase comes from the reaction field. However, we note that for the X Y Z model (Fig. 4.11) and the X Y model (Fig. 4.12), the contributions to the energy mass dependent. As the mass is increased, the structural contribution dominates at lower and lower temperatures. As the specific details Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 72 Figure 4.11: T h e short and long-range (reaction field) contributions to the energy per particle < U* > /N for the dynamically decoupled X Y Z model. The solid line corre-sponds to < U* > /JV, the short dash is < URF > /N and the long dash is < U*ST > jN. T h e squares, triangles and circles are for m* = 1, 5 and 10, respectively. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 73 Figure 4.12: The short and long-range (reaction field) contributions to the energy per particle < U* > /N for the dynamically decoupled X Y model. The solid line corresponds to < U* > /N, the short dash is < URF > /N and the long dash is < U*ST > /N. The squares, triangles and circles are for m* = 1, 5 and 10, respectively. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 74 Figure 4.13: The short and long-range (reaction field) contributions to the energy per particle < U* > /N for the dynamically decoupled Ising model. The solid line corre-sponds to < U* > /TV, the short dash is < URF > /N and the long dash is < Ugr > /N. The squares, and triangles are for m* = 1 and m* = 5, respectively. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 75 of the short-ranged structure are "turned on" (by increasing the mass), the structural contribution to the total energy increases, and the ferroelectric phase disorders. 4.5 S u m m a r y From the randomly frozen systems, we find no evidence of a ferroelectric phase for the X Y Z or X Y models at p* = 0.8 and p* = 4. However, we do observe a ferroelectric phase for the Ising model, with a transition temperature at T* « 25. We note that at this reduced density and dipole moment, Zhang and W i d o m [45] predict ferroelectric phases for the X Y Z model, and a transition temperature for the Ising model that is much higher. From the dynamically decoupled simulations, we find for the X Y Z and X Y models, that as the details of the random structure are built up (by increasing the mass), the systems disorder. In infinite mass limit, both the X Y and X Y Z systems are predicted to be disordered at all temperatures. T h e behavior of the Ising model was significantly different. Ferroelectric phases were found regardless of the mass, and the ferroelectric transition was driven by the reaction field. Structurally dependent interactions that favor unpolarized states had little effect on ferroelectric phase behavior. In all mass dependent systems ( X Y , X Y Z or Ising) that were ferroelectric, the reac-tion field was responsible for stabilizing the ferroelectric phase. In all cases, the random structural interactions did not favor ferroelectric order. Ferroelectric phases were ob-served when the reaction field contribution to the total energy was just large enough to dominate the structural interactions. Chapter 5 Ferroelectric Phases of Dipo la r Discot ic Part icles Liquid crystal phases of disc-shaped or discotic molecules have been observed ex-perimentally [53,54] and with computer simulation [23,42]. The phase diagram of these systems includes isotropic, nematic and columnar phases. In the nematic phase, the director d is associated with the symmetry axis of the particle. In the columnar phase, there is long-ranged spatial order where the particles arrange themselves in columns. Systems of discotic particles with dipoles embedded along the symmetry axis have been investigated with computer simulation [56,57]. The shape of the particle was a cut-sphere, as described in Chapter 2, with length L and diameter D such that L/D = 0.1. Two columnar phases were observed [56,57]. Systems with dipole moments given by p / ^J~(D3 kBT) < 0.125 exhibited unpolarized columns. However, systems with u/^J~(D3kBT) ~ 0.2 had completely polarized columns, but the total polarization of the system was zero. In this "antiferroelectric columnar" phase, equal numbers of columns are polarized in opposite directions. The columns are arranged with hexagonal symmetry where there is no well defined arrangement of polarization "up" and polarization "down". No ferroelectric phases were observed. It is interesting to compare the phase behavior of systems of dipolar cut-spheres (L/D = 0.1) and dipolar spheres [9,10]. At approximately the same dipole moment and reduced density, systems of dipolar cut-spheres are antiferroelectric columnar while dipolar spheres are ferroelectric nematic. Spherical, or near spherical dipolar particles form ferroelectric nematics, whereas, disc-like or oblate dipolar particles with L/D m 0.1 76 Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 77 do not. Evidently, particle shape has a significant role in determining the stability of ferroelectric phases. In fact, the only model up to now that has a ferroelectric phase is the dipolar sphere. For example, fluids of dipolar prolate ellipsoids [55] and spherocylinders [57,58] also do not exhibit ferroelectric behavior. The role of particle shape on the stability and nature of ferroelectric phases is not well understood. As a first step in understanding how particle shape affects ferroelectric behavior, we will investigate the break down of ferroelectric order in systems of oblate particles with point dipoles embedded along the symmetry axes. We will examine how changing particle shape (i.e., changing the aspect ratio) affects phase behavior. The first section involves some preliminary calculations that examines different lattices of dipolar particles. These calculations help explain the stability of the ferroelectric tetragonal I structure for systems of dipolar spheres, and the antiferroelectric columnar structure observed in fluids oiL/D = 0.1 dipolar cut-spheres. In the next section, a phase diagram for dipolar oblate particles is determined and discussed. W i t h intermediate aspect ratios b/a m 1.5 we find interesting ferroelectric phase behavior. 5.1 Solid Lattices of Oblate Particles As a preliminary investigation, we have constructed different lattices of dipolar particles. The lattices were chosen to coincide with the observed local structure in the ferroelectric nematic and the antiferroelectric columnar phases. The local structure in the dipolar sphere ferroelectric nematic is reminiscent of the tetragonal I crystal. This is the predicted ground state for dense dipolar sphere systems [44]. The antiferroelectric columnar phase has columns arranged with hexagonal symmetry. The arrangement of polarization "up" and polarization "down" is completely random so that the total polarization of the lattice is zero. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 78 Both of these lattices can be made by arranging columns of dipolar spheres, where the distance between spheres in a column is the diameter a of the sphere. T h e tetragonal I lattice is constructed by repeating a square unit cell with dimensions (x,y) = (n/,m/), where / is the length of the cell. The location of the cell in the lattice is given by the integers n — {0,1,2...} and m = {0,1,2...} . One column is at (x,y) = (n/,m/) and the other column is at the center of the cell with coordinates given by (x,y) = (l[n + 1/2], l[m + 1/2]). The column in the center is shifted by cr/2 in the z direction. The arrangement is shown in Fig. 5.1. T h e antiferroelectric columnar lattice is constructed by repeating a rectangular cell with dimensions (x,y) = (nl,m\/Zl) with columns at (x,y) = (nl,ml) and (x,y) = (l[n + 1/2],/[m + A/3/2]). Antiferroelectric columns were arranged by randomly orient-ing successive columns either "up" or "down". In this lattice there is no arrangement where every polarization-up column can be completely surrounded by polarization-down columns. Also, unlike the tetragonal I lattice, there is no vertical shifting, or interdigita-tion, of columns. A sketch of the lattice is shown in Fig. 5.2. In these calculations, lattices were "prepared" at the same reduced density, p* = pa3/\/2 = I/A/2, and dipole moment p* = (p2/ea3)1/2 = 3. The value of p* was chosen to coincide with the simulations as described in the next section. T h e results of the lattice calculations are independent of the reduced dipole strength p*, so that the energies reported are simply scaled by (p*)2. Expansion of the different lattices from close packing to this density was done only in directions perpendicular to the dipole. In this way, columns of head-to-tail dipolar particles were not broken. T h e anisotropic expansion is consistent with the dipolar soft sphere simulations [9,10], where strong head-to-tail correlations were found. Again, dipole-dipole interactions were calculated by the Ewald summation method, E q . (2.3.7), where the dielectric constant of the surroundings was set at e' = oo. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 79 Figure 5.1: A sketch of a tetragonal I lattice. The circles are "top" views of columns of particles. The dark shaded circles are columns of particles that are shifted 1/2 a particle diameter out of the page compared to the light shaded circles. Figure 5.2: A sketch of a hexagonal lattice. The circles represent columns of particles. There is no interdigitation between columns of particles. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 81 Figure 5.3 shows how the energy of the tetragonal I and antiferroelectric columnar lattices change as the particle's aspect ratio is increased. Keeping in mind that both the reduced density and dipolar strength are constant, we see that making the particle more oblate has a great effect on the relative stability of the different lattices. For reference, a dipole with LI* = 3 in an infinitely long column of dipoles will have an energy of Upp/N « —21.6 [44]. For the spherical case, the tetragonal I lattice has the lowest energy per particle with U£,D/N « —24.0. The interdigitation between columns results in column-column interactions that are short-ranged and attractive. A particle in this lattice will have a stronger interaction than it would have in an isolated infinite column. For spherical particles, the energy of the antiferroelectric columnar lattice is UpD/N s=s —22.4. The interaction between columns is attractive, but not as great as in the tetragonal I lattice. Again referring to Fig. 5.3, we see that as the aspect ratio of the particles is increased, the energy per particle of the ferroelectric tetragonal I and antiferroelectric columnar lattices rises and converges to the energy of a dipole in an infinite column. At an aspect ratio of b/a = 1.4, the energy difference for the two lattices is less than 0.2, and by b/a = 1.6 the energy per particle has converged to the energy of a dipolar particle in an infinite column of similar particles, Uf)D/N w —21.6. Three comments can be made based on these preliminary calculations. First, for near spherical particles, the ferroelectric tetragonal I lattice is the most stable, consistent with the local structure observed in the ferroelectric nematic liquid crystal. Secondly, at large aspect ratios, b/a > 1.6, there is little or no difference between the energy of the ferroelectric tetragonal I and the antiferroelectric columnar lattice. In fact, the energy per particle has converged to that of a dipole in an infinite chain of dipoles, and the short-ranged attractive column-column interactions have decayed to near zero. The antiferroelectric behavior observed in fluids of dipolar cut-spheres [56,57] can be Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 82 -22 Q * Q 23 -_ 1 1 1 1 1 1 1 .i...L.i..LJ 1 1 1 — t - i i — # — * * " i i i 1 i i i I n , i i i I I I ~ 1 1.2 1.4 1.6 b / a 1.8 Figure 5.3: Dipolar energy per particle UE>D/N for different lattices of oblate dipolar particles. Triangles are for an antiferroelectric columnar lattice, and squares are a ferro-electric tetragonal I lattice. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 83 explained simply in terms of entropy. The column-column interactions in the fluid are so weak that entropy ensures that equal numbers will be polarized in different directions. Finally, for intermediate aspect ratios 1.2 < b/a < 1.5, the suggested ground state is a ferroelectric tetragonal I lattice. However, from these calculations, the behavior of the fluid systems is not obvious. Fluid systems with these intermediate aspect ratios could have characteristics of both lattices. 5.2 Ferroelectr ic Phases of Simple Dipo la r Oblate Part icles In this section we present a partial phase diagram for dipolar oblate particles. In partic-ular, we focus on locating ferroelectric phases for dipolar oblate particles with different aspect ratios. A T * versus b/a phase diagram is constructed from the results of a series of M D simulations. We will examine systems of dipolar oblate particles with aspect ratios ranging from 1.2 < b/a < 3 at various reduced temperatures and densities. T h e observed phase behavior is explained by examining the local structure of the fluid. 5.2.1 S imula t ion Detai ls Dipolar oblate ellipsoids can be easily characterized by choosing the fundamental unit of length a to be the height of the particle a. We then can define the reduced dipolar strength p* = (p2/ea3)1/2, and reduced density p* = pab2. The pair potential for this system can be written as u{\2) = u r e p (12) + u D C ( 1 2 ) , (5.2.1) where u r e p (12) is the repulsive generalized Gaussian overlap potential given by E q . (3.1.8), and tt£>£)(12) is the dipole-dipole pair potential, given by E q . (2.2.7). Of course, in simulations UDD(12) is replaced by the effective pair potential, u^p(12; e'), given by E q . (2.3.7) and the parameters for the Ewald sum are as given in Chapter 2. The reduced Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 84 timestep St* = Sterna2je = 0.0012, hardness parameter s = 20, and moment of inertia parameter a\ = 0.5 (Eq. (3.2.2)) were as given in Chapter 3. Orientational order in ferroelectric liquid crystal phases can be detected by calculating the < Pi > and < P2 > order parameters [14,59]. In models where the dipole is along the symmetry axis of the molecule, the instantaneous value of the P2 order parameter is calculated from the largest eigenvalue of the ordering matrix, Q, with elements QAP = 1/N £ 1/2(3^4 - M , (5-2-2) w here pla is the a component of the unit vector given by fa = (p) 1^t4-. The eigenvector corresponding to the largest eigenvalue is the instantaneous director d . Again, if the dipole moment is along the symmetry axis of the particle, we can define the instantaneous first-rank order parameter P i as [10,59] N P i = 1/N (5.2.3) When the dipole is along the symmetry axis of the particle, P i is identical to the polar-ization P (Eq. (4.1.3)) in ferroelectric phases. For ferroelectric nematics, both < P i > and < P 2 > will be non-zero. For antiferroelectric columnar phases, < P i > will be near-zero, while < P 2 > will be non-zero . The non-zero < P 2 > results from the formation of columns, and the near-zero < P i > is due to columns polarized in opposite directions. For isotropic phases, both order parameters will be near-zero. A l l simulations were initiated with oblate particles on a face centered cubic ( F C C ) lat-tice. The particles were then equilibrated with no dipolar interactions for 20000 timesteps at the desired value of p* and T*. A configuration was then chosen from the equilibrated isotropic system and point dipoles with p* = 3 were embedded along the symmetry axes of the ellipsoids. The reduced density and temperature were held constant. A n initial equilibration run of 100,000 timesteps was then followed by a production run of 200,000 Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 85 timesteps. Most simulations were carried out with N = 256 particles, however, some runs with TV = 500 particles were done to ensure that the results were independent of system size. Calculated values of < P i > and < P 2 > did not exhibit N dependence. We note that systems of oblate ellipsoids with aspect ratios b/a < 2 and no dipole moment have no liquid crystal phases whatsoever [14,19]. A n y liquid crystal phases of dipolar oblate particles with aspect ratios in this range are due the addition of dipolar interactions. 5.2.2 R e s u l t s We begin discussing the results with Fig. 5.4, where the < P 2 > order parameter versus the reduced density p* is shown for systems with different aspect ratios. Here, the reduced temperature and dipole moment are T* = 1.9 and p* = 3, respectively. The results include systems with b/a = 1.2, 1.3 and 1.5. A l l three systems have < P 2 >^ 0 below p* « 0.7 and it begins to increase at p* w 0.8. However, the < P 2 > for the b/a = 1.5 system rises sharply to < P 2 >~ 0.7. This sharp increase in < P 2 > for systems with b/a = 1.5 can be attributed to strong short-ranged columnar "correlations". Strong head-to-tail dipole-dipole interactions, along with the oblate shape of the particles favors the formation of columns. Clearly, increasing the particle's aspect ratio favors < P 2 > ordering. The < P i > order parameter versus p* for the same T * , p* and aspect ratios is shown in Fig. 5.5. Each system has an isotropic to ferroelectric phase transition at p* « 0.8. The b/a = 1.5 system shows a sharper transition, with < P i > rising to around 0.85. This result is interesting in that is suggests that increasing a particle's aspect ratio (making it more oblate) favors ferroelectric ordering. However, if the same density scan is done with a lower reduced temperature, T* = 1.35, the order parameter decreases with increasing aspect ratio. For example, at when T* = 1.35, p* = 0.792, b/a = 1.5, then < Pi 0.15, compared to ~ 0.75 for T* = 1.9. This effect is clearly demonstrated in Fig. 5.6 where Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 86 1 I I I I I I I I i i i r i i i i 0.8 -A DTO.6 v I I -0.4 -0.2 -0 I I I I I:::::::* I I I I I 1 I I I I I I 0.5 0.6 0.7 * P 0.8 0.9 Figure 5.4: T* = 1.9 results for the < P2 > order parameter as a function of reduced density, p*, for oblate particles with different aspect ratios. Squares are for b/a = 1.5, circles are b/a = 1.3 and triangles are b/a = 1.2. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 87 1 I I I I I I I I I I I I I I I 0.8 -A pTO.6 v 0.4 -0.2 0 I I I I I I I I I ' l l I I I I 0.5 0.6 0.7 * P 0.8 0.9 Figure 5.5: T* = 1.9 polarization (< Pi >) as a function of reduced density p* for oblate particles with different aspect ratios. Squares are for b/a = 1.5, circles are b/a = 1.3 and triangles are b/a = 1.2. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 88 Figure 5.6: b/a = 1.45 p* = 0.792 < P2 > (circles) and polarization < Pi > (triangles) as a function of reduced temperature, T*. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 89 a temperature scan is performed on a system with b/a = 1.45, and p* = 0.792. There is a sharp rise in the polarization around T* = 1.9. Below and above this temperature, < Pi > is near-zero. Above T* = 1.9, the system is simply isotropic, as indicated by < Pi >~ 0. However, for reduced temperatures below T* = 1.9, < P2 >~ 0.8. This suggests that strong columnar correlations are present, and that the loss of polarization might be due to antiferroelectric columns similar to those observed in the dipolar cut-sphere systems [9,56]. T h e values of b/a, and p* that were chosen for this plot are not unique. Systems with aspect ratios between 1.4 and 1.5 and reduced densities ~ 0.8 all have similar behavior. In systems where the dipole lies along the symmetry axis of the particle, we can describe the local structure of the fluid by the pair correlation functions [9,10,56] g000(r) = < g ( F j ~ r ' ~ r ) > , (5.2.4a) 1 1 2 , s _ 3 < E,-,j,,?j 8{TJ - rt - r)[3Q,- • v^jfa • r,-,-) - fa • fa) > g^{r) = < ^ ^ ^ - r ^ ' ^ > , (5.2.4c) g000(r) is the usual radial distribution function. T h e average interaction energy between a pair of dipoles separated by a distance r is related to g112(r). T h e function <?110(r) gives an indication of parallel or antiparallel dipolar correlations as a function of the intermolecular distance. Also, we can define parallel and perpendicular correlation functions, denoted by g\\(r\\) and <7_L(rj_), respectively. If the instantaneous director is given by d , then r|| and rj_ are vectors such that r-j| = (r • d ) d and rj_ = r — r||. The correlation functions are defined as = < " ( ; i i - n i + fr|l)> (5.2.5, " l l v < n 0 (r N , r f l + Sru) > K ' Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 90 ( \ < n(rurL + 8rL) > < «o(r±,r ± + 8r±) > In the parallel case, (n) is the average number of particles in a disc-shape shell perpen-dicular to d at r|| of width Sr\\ and radius [(L/2)2 — r 2 | ] 1 / / 2 , where L is the length of the simulation cell. In the perpendicular case, (n) is the average number of particles in a cylindrical shell about d of radius rj_, width Sr± and height 2[(L/2)2 — r 2 ] 1 / 2 . In both cases (n0) is the average number of particles one would expect in the relevant shell calculated using the angle-averaged pair distribution function. Positional order parallel or perpendicular to the director will appear in these correlation functions as peaks, and in an isotropic fluid they will both be 1. Briefly, columnar ordering is detected by a large peak in g± at rj_ « 0, and peaks in g\\ at multiples of r-|| ~ a. In computer simulations, spatial correlation functions are cut off at L/2, where L is the box length. Therefore, long-ranged spatial order is indicated by correlations that do not decay beyond L/2. In Fig. 5.7, g 0 0 0 , g110 and g112 are shown for a system with p* = 0.792, b/a = 1.45, T * = 1.35. The sharp peaks in all three functions at r* « 2.1 and 3.2 indicates column-like structures where the dipoles are aligned head-to-tail. The parallel and perpendicular correlation functions for the same system are shown in Fig. 5.8. The peaks in g\\(r^) at rj| w 1, 2, 3 and the strong peak in g±(r]_) at r]_ = 0 indicates columnar ordering. Peaks in g± at r*x ~ 1.5 and 3 indicate that the columns are arranged with hexagonal symmetry. The peaks in <7||(rj|) at r*| « 1, 2, 3 arise from head-to-tail particles only, indi-cating that there is no significant interdigitation of the columns. If there was significant interdigitation of neighbouring columns (as in the Tetragonal I lattice), there would be peaks in 5||(^ |j) at rj| « 1, 1.5, 2, 2.5 and 3. The peaks in between the integer values are due to particles in adjacent columns that are staggered by half a particle diameter. For the system at T* — 1.35, these results indicate that the phase is antiferroelectric columnar, where the columns are arranged with hexagonal symmetry. For comparison, Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 91 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 0 1 2 3 4 Figure 5.7: g000(r*) (solid), <?112(r*) (short dash) and glw(r*) (long dash) for a p* = 0.792, b/a = 1.45, T* = 1.35 system. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 92 a snapshot of a typical antiferroelectric columnar phase is shown in Fig. 5.9. Here, the columns are well defined, and there is strong hexagonal symmetry of the columns. There is no average polarization. The pair correlation functions <7000(r*), g110(r*) and g112(r*) for the ferroelectric phase at T* = 1.9, p* = 0.792 and b/a = 1.45 are shown in Fig. 5.10. The functions have structure similar to the T* = 1.35 case, except that the peaks are now broader, suggesting that the columns are not as well defined and columnar correlations are short-ranged. One important exception is that giw(r*) decays to a non-zero constant, indicating that there is long-ranged ferroelectric order. The functions and <7x(rI)) a r e shown in Fig. 5.11. It can be seen that the small peaks in g\\ at multiples of the particle height decay with increasing r^. This suggests that the columnar correlations are short-ranged. Although the T* = 1.9 system is ferroelectric, its local structure is radically different from the soft sphere case. Here, the system is ferroelectric with short-ranged columnar order and short-ranged hexagonal symmetry. The only difference between this and the T* = 1.35 system, other than the fact that this system is ferroelectric, is that the columnar structure is "weaker" and shorter-ranged. For the isotropic T* = 2.5 system, g000(r*), gno(r*), and g112(r*) are shown in Fig. 5.12 and <7||(r-|j) and fl'i(r^) are presented in Fig. 5.13. Weak head-to-tail dipole-dipole correlations in g000(r*), g110(r*), and g112(r*) are found at r* 1, but the structure is very short-ranged. As expected in isotropic systems, ^(r-jj) and g±(r]_) are uniform. The structures of the isotropic, ferroelectric, and antiferroelectric columnar phases for the b/a = 1.45 system can be seen in configuration snapshots of these systems. Fig. 5.14 shows a typical antiferroelectric columnar, (b/a — 1.45, p* = 0.792, T* = 1.35) configuration. The system is fluid, and the columns bend and "snake". O n average there is no net polarization. A snapshot of a ferroelectric (b/a = 1.45, p* = 0.792, T * = 1.9) phase is shown in Fig. 5.15. This system has short-ranged columnar correlations, and Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 93 Figure 5.8: ^(rfj) (solid) and g±(r*_) (dash) for a p* = 0.792, b/a = 1.45, T * = 1.35 system. Figure 5.9: A computer "snapshot" of an b/a = 3, p* = 0.792, T* = 2.5 system. The phase is antiferroelectric columnar. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 95 Figure 5.10: # 0 0 0 (r*) (solid), g112{r*) (short dash) and glw(r*) (long dash) for a p* = 0.792, b/a = 1.45, T* = 1.9 system. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 96 Figure 5.11: g\\(r*\) (solid) and g±(r*±) (dash) for a p* = system. 0.792, b/a = 1.45, T* = 1.9 Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 97 I I I I I I I I I I I I I I I l l l l I l I 0 1 2 3 4 Figure 5.12: g000(r*) (solid), g112(r*), (short dash) and g110(r*) (long dash) for a p* = 0.792, b/a = 1.45, T* = 2.5 system. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 98 5-H J_ OO oo 1.5 1 0.5 i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r j i i i I i i i i I i i i i I i L 0 1 Figure 5.13: flrN(r[j) (solid) and g±(r*±) (dash) for a p * = 0.792, b/a = 1.45, T* = 2.5 system. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 99 the < Pi >K, 0.8. A snapshot of the isotropic (b/a = 1.45, p* = 0.792, T* = 2.5) phase is shown in Figure 5.16. In this phase there is no average polarization, and < P2 >~ 0. This contrasts with the antiferroelectric columnar phase, where < P2 > is large. Further evidence of a ferroelectric phase transition for the system with b/a — 1.45 p* = 0.792 is given in Fig. 5.17, where the average energy per particle < U* > /N is shown as a function of the reduced temperature T*. The inset is the heat capacity per particle, < Cv > /N, as obtained from numerically differentiating < U* > /N. There is an isotropic to ferroelectric liquid crystal phase transition at T* ~ 1.9, and at a slightly lower temperature T* & 1.6, there is ferroelectric to antiferroelectric columnar phase transition. The first isotropic to ferroelectric phase transition involves a significant change in local structure and energy, and is most likely first order. The ferroelectric to antiferroelectric columnar transition is less strong and is possibly second order. In the second transition, the columns become longer-ranged and the system is depolarized. Figure 5.18 shows how particle shape affects ferroelectric order at T* — 1.9 and p* = 0.792. As the particle is made more oblate, both order parameters rise. At b/a « 1.45, just as < P2 > is rising, the ferroelectric phase is observed. As the aspect ratio continues to increase, < Pi > drops. This is consistent with the greater columnar ordering tendency of the large oblate particles ( b/a « 2). At b/a ~ 1.45, the columnar correlations are strong enough for < P2 > ordering, but not strong enough for long-range columnar ordering. It is the lack of columnar order that allows the system to polarize. We can describe this ferroelectric phase as a thermal ferroelectric. Here, the low temperature state is an antiferroelectric columnar phase, but at some finite temperature near where the columnar phase "melts", the thermal ferroelectric phase exists. This is a fundamentally different situation than the soft sphere case because now the low temperature state is unpolarized. The breakdown of the antiferroelectric columnar phase results in a ferroelectric phase with short-ranged columnar structure. In this case, the Figure 5.14: A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 1.35 system. The phase is antiferroelectric columnar. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 101 Figure 5.15: A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 1.9 system. The phase is ferroelectric. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 102 Figure 5.16: A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 2.5 system. The phase is isotropic. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 103 A * V 6 rr i i i i i i i i i i i i i r -1 i i i i i i i i i i i i i i 30 - 9 8—ao L<c>/N j -10 -10 r 0 -I I I I I I I •mi 1 1 1 1 1 1 i -12 - 1.5 , 2 T 2.5 / I -14 -16 J I I I I I I I L J I L 1 1.5 2.5 Figure 5.17: T h e reduced energy per particle, < U* > jN as a function of reduced temperature, T*. The inset is the excess contribution to the reduced heat capacity per particle, < Cv > /N. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 104 local spatial structure of the ferroelectric is not reminiscent of a ferroelectric ground state. A convenient way to summarize the previous results is to draw a simple "phase dia-gram". In Fig. 5.19, a T* versus b/a phase diagram is sketched for p* = 3 and p* = 0.792. The reduced density was chosen to emphasize the ferroelectric liquid crystal regions. For lower densities, the isotropic phase would exist at lower temperatures, while at higher densities, the solid phase would occupy a larger region of the phase diagram. The phase boundary lines are approximate, and were estimated by examining the < P i > and < P2 > order parameters for selected systems. Beginning with the near spherical regime, b/a < 1.2, the ferroelectric ( F E ) phase is bounded at low T* by the ferroelectric solid and at high T* by the isotropic fluid. In this case the low temperature state of the system is ferroelectric [44] and the nematic liquid crystal has local spatial structure similar to the solid. At some point above b/a « 1.2, the low temperature phase is an antiferro-electric columnar solid. A n antiferroelectric columnar liquid crystal phase is observed upon melting of the solid. This phase is characterized by antiferroelectric columns of head-to-tail dipoles, where on average the columns are oriented so that the system is not polarized. The unpolarized phase is a result of weakly interacting columns of particles. As the aspect ratio is increased, the columns become more defined as the interactions between adjacent columns become weaker. For some aspect ratios, a ferroelectric liquid crystal phase exists between the isotropic and antiferroelectric columnar phases. There is evidence that the local correlations are not reminiscent of a ferroelectric Tetragonal I solid, but more closely resemble those in the antiferroelectric columnar phase. However, there is no long-range columnar order and the structure of the columns is "loose". As the aspect ratio of the particle is increased, the temperature range of the ferroelectric phase decreases, and the transition temperature from the antiferroelectric columnar to the ferroelectric phase increases. For example, systems with b/a = 2 have the antiferroelectric columnar to ferroelectric transition near Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 105 1 0.8 -0.6 -A P?0A v 0.2 h 0 1 1 1 1 1 1 1 1 i | i i i | i i i | i i i | i * * 1 — I I ' • -/ P*\ : i » i * , % i * : <p 2 > 4 _ i 1 i i i 1 i i i 1 i i i 1 i i i 1 i i i 1 i 1 1.2 1.4 1.6 1.8 2 b /a Figure 5.18: T * = 1.9 , / = 0.792, < P 2 >, (circles) and polarization < Pi > (triangles) as a function of aspect ratio, b/a = width/height. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 106 3 Figure 5.19: A rough phase diagram for dipolar, p* = 3, oblate ellipsoids. The reduced density is p* = 0.792. The width of the ferroelectric phase is only estimated. The phases are I (isotropic), F E L(ferroelectric liquid crystal), A F E C L (antiferroelectric columnar fluid), and S (solid). Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 107 T* ~ 2.3, and the reduced temperature range of the ferroelectric is about 0.2. Systems with b/a = 3 showed no ferroelectric phases at any temperature. The "collapse" of the ferroelectric phase with increasing aspect ratio is due to stronger columnar correlations. Particles with larger aspect ratios form columnar phases more easily, and eventually there is some aspect ratio where only an isotropic to antiferroelectric columnar phase transition occurs. There are no ferroelectric phases above this aspect ratio. 5.2.3 Summary Clearly, particle shape can have a profound effect on the existence of ferroelectric liquid crystal phases. Small changes in particle shape can either support or destroy ferroelectric fluid phases. There is evidence that columnar phases will rarely, if ever be ferroelectric. However, there is evidence that ferroelectric phases can be found in systems other than spherical dipolar particles. Systems with low temperature antiferroelectric columnar phases can have a ferroelectric phase at higher temperatures. This thermal ferroelectric phase is characterized by short-ranged columnar correlations and exists in small regions of the phase diagram for specific p, T and b/a. Chapter 6 A More Robust Model for a Ferroelectric Liquid Crystal From the last Chapter, there was evidence that a ferroelectric liquid crystal phase might be found in systems of dipolar oblate particles that have low temperature antifer-roelectric columnar phases. Usually, the antiferroelectric columnar phase melts directly to the isotropic fluid, however for some combinations of particle shapes, dipole moments and densities, the antiferroelectric columnar phase melts to a ferroelectric liquid crystal. The short-ranged columnar correlations in the ferroelectric phase are not strong enough to form antiferroelectric columns. However, the polarized phase is very sensitive to the details of molecular shape, and only occupies a small region of the phase diagram. It is clear that if one wishes to build more diverse models which might eventually serve as a guide to the discovery or "construction" of new ferroelectric liquid crystals, then a more robust approach must be found. There are a number of criteria that must be met in order for a system to have a ferroelectric phase. In terms of particle features, there are three important components that should be incorporated into the model. First, it must have a dipole. Secondly, the short-range interactions (arising from both the particle shape and the charge distribution) must not strongly inhibit ferroelectric ordering. In fact, if possible the short-range interactions should favor ferroelectric order. Thirdly, interactions that favor columnar ordering should be minimized. In this Chapter we describe one possible approach which readily yields a ferroelectric discotic nematic phase. Monte Carlo calculations are used to examine orientational order in fluids of disc-shape particles with embedded dipoles. However, instead of one point 108 Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 109 dipole along the symmetry axis, now many dipoles are distributed over a circular "patch" of finite size. The "patch" of embedded dipoles acts to breakdown columnar ordering by dispersing the strong head-to-tail point dipole interaction. It is shown that such systems undergoes spontaneous polarization to form a stable ferroelectric discotic nematic phase. This ferroelectric phase exists in a much larger region of the phase diagram, and is not overly sensitive to the details of molecular shape. 6.1 The Cut-Sphere Model For convenience we have chosen to work with cut-spheres although other models such as oblate ellipsoids would likely serve as well. The phase behavior of plain cut-spheres has been studied in detail [24,23]. For our purposes the main observation is that for sufficiently large aspect ratios and moderate densities cut-spheres form a discotic nematic phase with no long-range spatial order. This is true of cut-spheres treated with all rotational degrees of freedom [23] and restricted so that the symmetry axes of all particles are parallel [24]. At higher densities this model possesses columnar and other spatially ordered phases. The model consists of cut-spheres of length, L, and diameter, D , with a uniform distribution of point dipoles embedded in the central plane (i.e., L/2 from a surface). A sketch is shown in Fig. 6.1. The dipoles are directed along the symmetry axis and cover a circular "patch" of diameter D'. If the dipole density in the plane is pes, then the total dipole moment of a particle p is given by p = pesA, where A = TTD'2/4 is the area of the dipolar patch. For this model the pair potential has the form (6.1.1) where u c s(12) is the hard cut-sphere (cs) interaction and u e s(12) is the total electrostatic Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 110 (es) potential. The total electrostatic pair interaction can be obtained by integration over the dipolar distributions of a pair of particles. Explicitly, one has 3(12) = pi J J[-3(AI • s)(/x2 • s ) + A i • fa}/s3dSlds2 , (6.1.2) where the integrations are over the areas A\ and A2 of the dipolar patches and the various vectors are shown in Fig. 6.1. In particular, A i a n d A 2 a r e u m 1 : vectors directed along the dipoles of particles 1 and 2, Si and s 2 are position vectors of points within the dipolar patches, s = s 2 — S i , s = s/s and s = | s | . At long-range, u e s(12) clearly must obey the limiting expression «es(12) ~ - 3 ( / i i • r)(fi2 • r ) / r 5 + /*x • p,2/r3 , (6.1.3) which is the interaction between two point dipoles. Unfortunately, it is not possible to obtain an analytical expression for the short-ranged electrostatic interactions. Therefore, for computational purposes the short-range part of the interaction potential must be calculated numerically and tabulated. In general, a complete description of the potential would require a four dimensional table which is computationally problematic. However, this can be significantly simplified by considering a model which is Ising in nature. That is, the symmetry axes of all cut-spheres must be parallel to each other and the dipoles may point only up or down. For this case the pair potential is a function of r and of the angle between p,\ (or p,2) and r , hence only a two dimensional table is required. Therefore, for computational purposes the total electrostatic contribution to the pair potential is given by •ues(12) = utabie(12) + u ^ c ( 1 2 ; e ' ) , (6.1.4) where «tabie(12) = u e s ( 1 2 ) - f - 3 ( / i i - r ) ( / i 2 - r ) / r 5 + A*i-A*2/r3l . (6.1.5) Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 111 Figure 6.1: The coordinate system used to evaluate the interaction potential. The rel-evant vectors defined in the text are shown. The inset is a cross sectional view of a cut-sphere with point dipoles dispersed over a circular patch of diameter D'. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 112 T h e long-range nature of the electrostatic interaction is correctly taken into account by UDE>C(12; e' = oo), which is the effective dipole-dipole interaction under periodic boundary conditions with a conducting surrounding continuum ( E q . (2.3.6)). The short-range electrostatic interactions due to the dipolar patch are contained in u t a bie(12). 6.2 Simulation Details A set of parameters was chosen to minimize the influence of the Ising approximation. In particular, only highly anisotropic particles with D/L = 10 are considered. T h e phase behavior of plain cut-spheres with this aspect ratio has been studied both at the fully rotational [23] and restricted to be parallel levels [24]. In both cases it is found that a nematic phase is stable at p* = pLD2 = 0.525. Furthermore, for the fully rotational system the second-rank order parameter < P2 > is ~ 0.95 at this density. Hence, the particles are strongly ordered and one would expect the parallel model to give a reasonable representation of the nematic fluid. Of course, for the parallel case < P2 > = 1 by definition and the description "nematic" refers to the absence of any long-range spatial order. In the present calculations, we have deliberately avoided higher densities where parallel and fully rotational cut-spheres may give different spatially ordered states [23,24]. The phase behavior of the model was explored by carrying out Monte Carlo ( M C ) cal-culations in the canonical ensemble. Calculations were performed at p* = 0.525 varying the reduced temperature T* = ksTL3/p2. In most calculations 288 particles were used, but several tests with 540 particles gave very similar results. Typical M C calculations consisted of an equilibration run of 100,000 moves per particle followed by a production run of equal length over which averages were obtained. Care was taken to verify that all results reported were independent of the initial starting configuration; calculations begun with fully polarized and randomly oriented dipoles gave the same final result. Also, to Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 113 confirm that all systems simulated were fluid, we calculated the mean square displace-ment, ((R(r) — R(0)) 2), where R(r) is the position vector of a particle and r represents the number of Monte Carlo moves per particle. This is just the M C "analog" of the M D mean square displacement, and it has similar behavior. However, in this case, it is not related to the diffusion coefficient. From Fig. 6.2, it is clear that in all cases the mean square displacement steadily increased with increasing r indicating fluid behavior. This contrasts with the situation for solids where ((R(r) — R(0)) 2) becomes constant as r increases. The presence of ferroelectric order was detected by calculating the average polariza-tion per particle < P > defined in E q . (4.1.3). Of course < P > will be near zero in an unpolarized system and will approach 1 if the fluid becomes ferroelectric. Also, here the polarization < P > is identical to the < P i > order parameter since the particle dipole vector is along the symmetry axis. Calculations were carried out for dipolar patches ranging in size from D'/L = 4 to 9. For values of D'/L ^ 7, spontaneous polarization occurred and a ferroelectric phase was formed. For smaller values of D' ferroelectric transitions were not observed. For D'/L = 9, the temperature dependence of < P > and of the reduced constant volume heat capacity Cv = Cv/Nks are shown in Fig. 6.3. The heat capacity was found by numerically differentiating the reduced energy per particle < U* > /N. It can be seen that this system has a clear ferroelectric transition at T* « 0.03. For physical reference, note that if T = 298K, and T * = 0.03, then for L = 3A, p = 6.08 Debyes and p e s = 0.011 D e b y e s A - 2 . Likewise, for L = 10A, the dipole moment is 37.0 Debyes and pes = 0.0058 D e b y e s A - 2 . "Snapshots" of fluid configurations for D'/L = 9 and D'/L = 5 are shown in Figures 6.4, 6.5 and 6.6. From Figure 6.4 and 6.5, we see that the D'/L = 9 system has no discernible long-range spatial order in either the ferroelectric or unpolarized state. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 114 Figure 6.2: Th e mean square displacement, ((R(r) - R(0)) 2 ) , for T* = 0.02 (short dash), T* = 0.03 (solid), and T* = 0.05 (long dash). R ( r ) is the position vector of a particle and r represents the number of Monte Carlo moves per particle. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 115 1 A P-i V 0.8 -0.6 -0.4 -0.2 -0 1 1 1 * 1 1 1 1 1 1 1 i \> 2.5 F-2 F-I 1.5 i i i j ii 111111 ii I i i i 111111 i i i I i i i I i i i I i i i I i i i"l -\ 0 2 4 6 8 102T* ^ i i i i i i i I i i i i i i i 0 2 4 6 8 102T* Figure 6.3: T h e temperature dependence of the polarization < P > and reduced constant volume heat capacity per particle Cv/N (inset) for a system with D/L — 10, D'/L = 9 and p* = 0.525. The error bars represent one estimated standard deviation. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 116 Figure 6.4: "Snapshot" of a configuration for a system with D/L = 10 and p* The system with D'/L = 9 is shown in the ferroelectric (T* = 0.025) phase. = 0.525. Figure 6.5: "Snapshot" of a configuration for a system with D/L = 10 and p* The system wi th D'/L = 9 is shown in the unpolarized (T* = 0.050) phase. = 0.525. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 118 Figure 6.6: "Snapshot" of a configuration for a system with D/L = 10 and p* = 0.525. The system wi th D'/L = 5 is shown in the unpolarized (T* = 0.0275) phase. Note that this configuration consists of loose polarized columns with nearly equal numbers up and down. On average, there is no net polarization in this system. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 119 To ensure that the system did not possess long-ranged smectic order, we calculated an appropriate "smectic" correlation function gs(rz) [9,58]. This function identifies long-ranged smectic order in planes perpendicular to the z axis. Since the orientation of the particles are constrained to be along this axis, smectic order in other directions is not possible. gs(rz) gives the probability of finding a particle j with an rz component of the intermolecular vector given that there is a particle i at the origin and the particle j lies within a cylinder defined by the shortest box length (either the box length in the x direction or the y direction). If the height of the cylinder is given by Srz, and its volume is K y i i n d e r = 7rL2c>r2/4, then gs(rz) can be written as , . ( r , + l / M r . ) = <"('y+fr')>. (6.2.1) P ^cylinder Here, n(rz, rz + Srz) is the instantaneous number of particles found in the range rz,rz+6rz, and the average is over the number of particles JV, and the number of M C cycles. Figure 6.7 shows gs{r*z) i " o r the ferroelectric T* = 0.025, p* = .525 system. The function oscillates with a period slightly greater than L, and decays to 1 at long-range indicating that only short-range smectic correlations exists, and that the ferroelectric phase is nematic. In contrast, it is apparent from Fig. 6.6 that at lower temperatures the D'/L = 5 system is developing strong columnar correlations. At this value of D' the columns are "loose" and interpenetrate to some extent. Nevertheless, the loose columns are polarized and on average equal numbers point in opposite directions such that no net ferroelectric order is achieved. As D' is decreased keeping u fixed, the columns simply become better ordered and more rigid as the single point dipole limit is approached. It is also worth noting that the homogeneous ferroelectric state shown in Fig. 6.4 was obtained with conducting boundary conditions. If instead the sample is surrounded by a vacuum, then a domain structure is obtained exactly as was found in previous simulations of ferroelectric fluids of dipolar spheres [10]. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 120 Figure 6.7: gs(r*z) for T* = 0.025. No long-range smectic order is detected in this system. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 121 In order to understand the physical origin of the phase behavior described above, it is useful to more closely examine how the behavior of the electrostatic part of the pair interaction varies with D'. In Fig. 6.8 the dimensionless quantity L3ues(12) / p2 is plotted as a function of the component of r perpendicular to the particle symmetry axes, rj_, as one particle is slid over another with their surfaces in contact and their dipoles in the same direction. Note that when r^/L = 0 one particle completely covers the other and that surface contact is broken when r±/L = 10. Curves are shown for single point dipoles (D'/L ~ 0) and for D'/L = 4, 5 and 9. It is evident from Fig. 6.8 that spreading the dipoles over a circular patch has a dramatic effect on the pair potential. For single point dipoles, the strong head-to-tail interaction results in a very sharp minimum at r±/L = 0. It is the sharpness of this potential which leads to columnar order and, consequently, to antiferroelectric behavior in these systems. As the dipoles are spread over larger and larger patches, the potential looses its sharp focus and one has instead a weaker more slowly varying attraction. Thus, as D' is increased the columnar structure is weakened and is eventually replaced by the ferroelectric nematic phase. The ferroelectric transition for the D' = 9 system can be explained by examining the long and short-range contributions to the energy per particle. As in Chapter 4, we can conveniently separate the total energy < U > into a reaction field contribution < URF > and the rest < UST > given by < u > = < URF > + < UST > • (6.2.2) The long-ranged reaction field contribution < URF >, depends on the total dipole moment of a particle's charge distribution. The short-range structural contribution < UST > depends on the shape and size of the electrostatic "patch". In Fig. 6.9 the total reduced energy per particle < U* > /iV, the reaction field Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 122 Figure 6.8: T h e variation of the reduced electrostatic potential energy, I / 3 u e s ( 1 2 ) / / u 2 , as one particle is slid over another with their dipoles directed in the same direction and their surfaces in contact. The inset shows the same curves on a magnified scale. Curves are shown for single point dipoles (D'/L ~ 0) (long dash) and for D'/L = 4 (short dash), 5 (dotted) and 9 (solid). Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 123 0 A * V -2 -4 i i i i i i i i i r i i i i i i i i • i i i I I - I I 4 102T* 6 Figure 6.9: T h e contributions to the reduced energy per particle, < U* > jN. T h e squares are the total energy < U* > /N, the triangles are the reaction field contribution < URF > /N, and the circles are the short-range structural contribution (everything else) < U*ST > IN. The reduced density is p* = 0.525. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 124 contribution < URF > /N, and the structural contribution < UgT > /N are shown as a function of T* for the p* = 0.525, D' = 9 system. Recall that the ferroelectric transition occurred at T* w 0.03. The total reduced energy per particle < U* > /N becomes more negative as the reduced temperature is decreased. The interesting observation is that both the reaction field and the structural contributions to the energy per particle become lower when the phase becomes ferroelectric. In this case, both the long-range and the short-range interactions favor the ferroelectric state. The distribution of charge not only "dampens" the interactions that favor columnar order, but it also favors ferroelectric over disordered phases. This situation is different from the ferroelectric dipolar sphere case where the short-range structural interactions do not favor the polarized phase. 6.3 S u m m a r y We have described a model discotic nematic liquid crystal which forms a stable ferroelec-tric phase. The essential feature of the model is that the dipole density is spread over a significant portion of the central plane of the disc-shaped particle. Spreading the dipole density has the effect of almost cancelling out short-range interactions that favor colum-nar ordering. This avoids the antiferroelectric columnar phase previously observed [56] for cut-spheres with single point dipoles. The ferroelectric phase formed by this model is not overly sensitive to the details of molecular shape, and exists over a large temperature range. Chapter 7 Summary and Conclusions This Thesis has been mainly concerned with ferroelectric liquid crystals. Computer simulations were employed to study the behavior of dense systems characterized by strong dipolar interactions. Various simple liquid crystal models were used to examine the role of particle shape on ferroelectric phases. A n extension of the original Berne and Pechukas [34]Gaussian overlap model was developed. The generalized Gaussian overlap model was designed to incorporate molecular flatness typical of most liquid crystal molecules. Molecular shape was approximated by a three dimensional Gaussian ellipsoid with non-degenerate semi-axes. The role of molecular flatness in uniaxial nematic liquid crystal phases was examined by a series of molecular dynamics simulations. Spatially amorphous systems of dipolar soft spheres were used to study the effects of short-ranged structure on ferroelectric ordering. Monte Carlo and M D simulations of both spatially frozen and dynamically random dipolar media were performed. The results offer an answer to the long standing question of whether or not spatially random dipolar systems could have ferroelectric phases. The effect of particle shape on the stability of ferroelectric liquid crystal phases was studied with M D simulations of oblate particles with dipoles embedded along the symmetry axes. A phase diagram was constructed, locating ferroelectric re-gions. Finally, a more realistic and robust ferroelectric liquid crystal model was designed, where specific features that favored ferroelectric phase formation were incorporated. M C simulations were used to study the phase behavior of this model. T h e generalized Gaussian overlap potential was designed for use in M D simulations of 125 Chapter 7. Summary and Conclusions 126 simple liquid crystal models. As with the original Berne-Pechukas model [34], the shape of a molecule is approximated by a three dimensional Gaussian. The repulsive energy of two particles is associated with the mathematical overlap of the two Gaussians. Expressions for the pair energy as well as the forces and torques necessary for M D simulations were presented. The form derived here is capable of modeling ellipsoidal shapes that are not axially symmetric. This was a necessary extension to the original axially symmetric Berne-Pechukas [34], and subsequent Gay-Berne [35-39] Gaussian overlap models. In these earlier studies, the liquid crystal molecule was modeled by an ellipsoid of revolution, and the flatness characteristic of real liquid crystals was not incorporated. M D simulations designed to isolate the effects of molecular biaxiality on the stabil-ity and formation of uniaxial liquid crystal phases were performed. Beginning with an isotropic phase of axially symmetric particles, we slowly flattened the particles, holding the reduced density and temperature constant. It was found that systems with small amounts of molecular biaxiality can form uniaxial nematic phases. For example, at a reduced density where particles characterized by the aspect ratios 3 : 3 : 9 are isotropic, those with ratios 2 : 3 : 9 form a strongly ordered uniaxial nematic phase. The role of molecular biaxiality on the stability of the uniaxial nematic is found to be very subtle. In our purely repulsive model, the isotropic to uniaxial nematic phase transition is entirely entropy driven. When the particles align in the liquid crystal phase, the increase in trans-lational freedom is greater than the loss of orientational freedom. In the uniaxial nematic phase, molecular biaxiality acts to increase the translational freedom of the particle. The possible existence of ferroelectric phases in spatially disordered dipolar materi-als was examined with M C and M D simulations of randomly frozen and dynamically disordered dipolar soft spheres. This study addresses the basic questions concerning the existence and nature of ferroelectric order in positionally disordered dipolar materi-als [46-48]. Systems with one (Ising), two ( X Y ) and three ( X Y Z ) dipolar components Chapter 7. Summary and Conclusions 127 were considered. Frozen spatially random systems of Ising dipoles exhibited a transition to a polarized state at finite temperature. However, randomly frozen systems with two and three dipolar components showed polarization that decreased with increasing system size. The results strongly indicate that the systems are unpolarized in the thermodynamic limit. From dynamically decoupled simulations, it was found that specific short-ranged spa-tial correlations act to oppose ferroelectric ordering. Furthermore, ferroelectric order, when it is found, is driven by long-range correlations. For the X Y and X Y Z models, the position dependent short-range correlations dominate over the long-range correlations, and a ferroelectric phase was not observed. The Ising model has a ferroelectric transi-tion which can be accurately predicted by considering only the long-range reaction field interactions. Constant temperature M D simulations of oblate ellipsoids with point dipoles em-bedded along the symmetry axes were performed. A phase diagram was constructed, emphasizing the location of ferroelectric liquid crystal phases. It was found that particle shape has a large effect on the stability of ferroelectric phases. Ferroelectric liquid crystals were found in systems where columnar correlations were short-ranged. The stability of the ferroelectric phase depended on details of particle shape and state parameters. Sys-tems with longer-ranged columnar correlations formed antiferroelectric columnar phases. It was found that in columnar phases the interaction energy between columns decreased rapidly as the aspect ratios of the particles were increased. Antiferroelectric columnar phases were observed in systems where the interaction energy between columns was neg-ligible, and the unpolarized state is simply an entropic effect. Entropy ensures that equal number of columns are polarized in opposite directions. A more "realistic" and robust model that forms a ferroelectric liquid crystal was designed. Systems of large cut-spheres with embedded "electrostatic patches" along the Chapter 7. Summary and Conclusions 128 central planes of the particles were considered. The electrostatic patch consisted of an even distribution of point dipoles such that at long-range the patch behaved as a single point dipole. It was found that this charge distribution not only broke down column formation, but resulted in short-range interactions that favored ferroelectric order. The model spontaneously polarized at finite temperature and a stable ferroelectric discotic nematic was observed. In this model, it was found that both the long-ranged and short-ranged interactions favored the ferroelectric phase. This result contrasts the dipolar soft sphere case, where it was found that only the long-ranged interactions favored ferroelectric phases. As the size of the "patch" was decreased, antiferroelectric columnar phases were observed. B i b l i o g r a p h y [1] P. J . Collings, Liquid Crystals (Princeton University Press,New Jersey,1990). [2] M . Born, Sitz. Phys. Math. 25, 614 (1916); Annalen der Physik. 55, 221 (1918). [3] W . Maier and A . Saupe, Z. Naturforsch., 13a, 564 (1958); 14a, 882 (1959), 15a, 287 (1960). [4] G.R. Luckhurst and G . W . Gray, The Molecular Physics of Liquid Crystals (Aca-demic, London, 1979), chap. 1. [5] C. Destrade, P. Fouchet, H . Gasparoux, N . H . T i n h , A . M . Levelut and J . Malthete, Molec. Crystals liq. Crystals, 106, 121 (1984). [6] D. Frenkel, Mol. Phys. 60,1 (1987). [7] R. B. Meyer , L. Strzelecki and P. Keller, J . Physique Lett., 36, 69 (1975). [8] R. H . Tredgold, J . Phys. 59, 29 (1975). [9] J . J . Weis, D. Levesque, and G . J . Zarragoicoechea, Phys. Rev. Lett. 69, 913 (1992); J . J . Weis, D. Levesque, Phys. Rev. E , 48,3728 (1993). [10] D. Wei and G . N . Patey, Phys. Rev. Lett. 68, 2043 (1992); Phys. Rev. A . 46, 7783 (1992). [11] Goldstein, H . Classical Mechanics 2nd edition (Addison-Wesley,1980). [12] J . P. Hansen, I. R. McDonald, Theory of Simple Liquids (Academic Press, London, 1986). [13] D. A . McQuarrie, Statistical Mechanics (Harper & Row, New York, 1976). [14] M . P. Allen and D. J . Tildesely, Computer Simulation of Liquids (Clarendon, Ox-ford, 1989). [15] P. W . Atkins, Physical Chemistry, 3 r d ed. (W. H . Freeman, New York, 1979). [16] C. G . Gray and K. E . Gubbins, Theory of Molecular Fluids Vol. 1 (Clarendon, Oxford, 1984). 129 Bibliography 130 [17] G.R. Luckhurst and G . W . Gray, The Molecular Physics of Liquid Crystals (Aca-demic, London, 1979), chap. 2. [18] M . R Allen, G . T . Evans, D. Frenkel, and B . M . Mulder, Adv. Chem. Phys. 86, 1 (1993). [19] D. Frenkel, B. M . Mulder, and J . P. McTague, Phys. Rev. Lett. 52, 287 (1984); [20] J . W . Perram, M . S. Werteim, J . L. Lebowitz and G . 0 . Williams, Chem. Phys. Lett. 105 no. 3, 277 (1984); J . W . Perram and M . S. Wertheim, J . Comp. Phys. 58 409 (1985). [21] G . J . Zarragoicoechea, D. Levesque, and J . J . Weis, Mol. Phys. 75, 998 (1992). [22] D. Frenkel, H . N . W . Lekkerkerker, and A . Stroobants, Nature 332, 822 (1988). [23] J . A . C . Veerman and D. Frenkel, Phys. Rev. A 45, 5632 (1992). [24] H . Azzouz, J . M . Caillol, D. Levesque, and J . J . Weis, J . Chem. Phys. 93, 3520 (1990). [25] S.W. de Leeuw, J . W . Perram, and E . R . Smith, A n n . Rev. Phys. Chem. 37, 245 (1986); Proc. R. Soc. A373, 27 (1980); A388 (1983). [26] D . J . Adams and I. R. McDonald, M o l . Phys. 32 931 (1976). [27] P . G . Kusalik, J . Chem. Phys. 93, 3520 (1990). [28] N . Metropolis, A . W . Rosenbluth, M . N . Rosenbluth, A . H . Teller, and E . Teller, J . Chem. Phys. 21, 1087 (1953). [29] H . F . Davis, A . D . Snider Introduction to Vector Analysis 4th edition (Allyn and Bacon) 1979. [30] C. W . Gear Numercial initial value problems in ordinary differential equations. Prentice-Hall, Englewood Cliffs, N . J . (1971) [31] L. Verlet, Phys. Rev. 159, 98, (1967). [32] D. Evans, Phys. Rev. A 28, 2, 1016 (1983). [33] D. Evans, Phys. Lett. 98A, 433 (1983). [34] B. J . Berne and P. Pechukas, J . Chem. Phys. 56, 4213 (1972). [35] J . G . Gay and B . J . Berne, J . Chem. Phys. 74, 3316 (1981). [36] D . J . Adams, G.R. Luckhurst, and R . W . Phippen, Mol. Phys. 61, 1575 (1987). Bibliography 131 [37] E . de Miguel and L. F . Rull , Mol. Phys. 74, 405 (1991). [38] E . de Miguel, L . F . Rull , M . K . Chalam, and K . E . Gubbins, M o l . Phys. 74, 405 (1991). [39] E . de Miguel, L . F . Rull, M . K . Chalam, and K . E . Gubbins, and F . van Swol, M o l . Phys. 72, 593 (1991). [40] M . P. Allen, Liquid Crystals 8, 499 (1990). [41] J . P. Straley, Phys. Rev. A . 10, 1881 (1974). [42] D. Frenkel, Mol. Phys. 60, 1 (1987). [43] J . Vieillard-Baron, J . Chem. Phys. 56, 10, 4729 (1972). [44] R. Tao and J . M . Sun, Phys. Rev. Lett. 67, 398 (1991). [45] H . Zhang and M . W i d o m , J . Magn. Magn. Mater. 122, 119 (1993); Phys. Rev. B 51, 8951 (1995). [46] U . T . H6chli and M . Maglione, J . Phys. Condens. Matter 1, 2241 (1989). [47] T . Jonsson, J . Mattsson, C. Djurberg, F . A . Khan, P. Nordblad, and P. Svedlindh, Phys. Rev. Lett. 75, 4, 4138 (1995). [48] B . E . Vugmeister and M . D . Glinchuk, Rev. M o d . Phys. 62, 4 (1990). [49] J . R. Thompson, Hong Guo, D. H . Ryan, M . J . Zuckerman, and Martin Grant, Phys. Rev. B 45, 3129 (1192). [50] B. Bergersen and Z. Racz, Phys. Rev. Lett. 67, 3047 (1991); H - J X u , B. Bergersen, and Z. Racz, Phys. Rev. E 47, 1520 (1993) [51] D. Wei, G . N . Patey, and A . Perera, Phys. Rev. E 47, 506 (1993). [52] C. J . F . Bottcher, Theory of Electric Polarization 2ND ed. (Elsevier Scientific Pub-lishing Company, 1973). [53] S. Chandrasekar, B. K. Sadashiva and K . A . Sureshi, Pramana, 9, 471 (1977). [54] J . C. Dubois, A n n . Phys. 3, 131 (1978). [55] G . J . Zarragoicoechea, D . Levesque, and J . J . Weis, Mol. Phys. 75, 989 (1992). [56] G . J . Zarragoicoechea, D. Levesque, and J . J . Weis Mol. Phys. 78, 1475 (1993). Bibliography 132 [57] J . J . Weis, D . Levesque, and G . J . Zarragoicoechea, Phys. Rev. Lett. 69, 913 (1992). [58] D . Levesque, J . J . Weis, and G . J . Zarragoicoechea, Phys. Rev. E 47, 496 (1993). [59] F. Biscarini, 'C. Zannoni, C. Chiccoli, and P. Pasini, Mol. Phys. 73, 439 (1991). [60] D . J . Evans, M o l . Phys. 34, 317 (1977). [61] M . P. Allen, Mol. Phys. 52, 717 (1984). A p p e n d i x A T h e P a i r P o t e n t i a l for P o i n t D i p o l e s The pair potential for point dipoles can be easily derived. A point dipole is obtained in the limit that two charges q and —q, separated by a distance 1 are brought together such that the distance / —> 0, the charge q —> co and the product gl = fi. To begin the derivation we must first define the gradient operator, W ) / = J i m l ( ^ ( i ' ) - W ) ) , ( A . l ) where <f> is some scalar function evaluated at points P and Q, and the distance between these points is 81. Consider the two charges q and —q separated by a distance /. The negative charge, for convenience, is at the origin. The situation is shown in Fig. A . l . The potential at some point r + , due to q is <^+(r+) = q/r+, r+ — | r + | . Likewise, the potential at r + due to the negative charge at the origin is <^-(r_) = — g/r_, where r + + 1 = r_ and r_ = |r_|. The potential at r, where r = r_, due to these point charges is < K r ) = M r + ) + ^_(r_). (A.2) " If the potential due to the positive charge is evaluated at a point given by r' = r_ + 1, then the potential at r' is ^ + ( r + ) = - M O , ( A - 3 ) # r ) = - ( M O " M r - ) ) . ( A - 4 ) As / becomes small, the definition of the gradient operator can be employed to give # r ) = - l - V # r ) , (A.5) 133 Appendix A. The Pair Potential for Point Dipoles 134 which is the potential at r due a point dipole at the origin. If we now consider another set of charges, — q at the point r, and the other, q at r ' , we can calculate the energy of these two charges due to the point dipole field at the origin. Explicitly, uDD(l2) = q(-<f>(r) + <^(r')), (A.6) where r' = r + l . (A.7) Again we can apply the gradient operator to obtain the well known result, uDD(12) = <?1- V ^ r ) , (A.8) and if the limit is taken, / —> 0, q —> oo and q\ = p, we can write uDD(l2) = A * - V ^ ( r ) M i ' A * 2 3(/xi • r)(f* 2 • r) (A.9) endix A. The Pair Potential for Point Dipoles Figure A . l : The geometry that is used for the point dipole pair potential. A p p e n d i x B T h e E w a l d S u m for D i p o l a r S y s t e m s In this Appendix a brief derivation of the Ewald summation for point dipoles is presented. This particular derivation is similar to the derivation given for point charges in [25]. We consider a central simulation cell containing N point dipoles. T h e cell is surrounded by an "infinite" sphere of periodic images. We must be careful with the word "infinite", because we must specify the nature of the material beyond the sphere. Specifically, we must specify its dielectric constant e . For this derivation we set e = 1, corresponding to surrounding the sample with a vaccuum. Furthermore, we must specify the order that we sum the periodic images. We will choose to sum approximately spherical shells about the central cell. Each spherical shell is composed of cubic cells arranged to form a "jagged" sphere. As the diameter of the shell gets larger, the sphere gets "smoother". In this case, the energy of the central cell, USHI SH = spherical shell, for a system of N dipoles is given by [25] -USHL* = lim iEEE'(^ • V ) ( M i • V ) r ^ - r , ( B . l ) s ^ u z ,-=i j=i n l r T n l where r = rj — r 4 and n are reduced by the boxlength, L. The reduced lattice distance between a dipole i in the central cell, and a dipole j either in the central cell (n = 0) or in one of the periodic images (n ^ 0) is given by |r + n|. The sum over n ' means that the term n = 0 is not included when i = j. T h e gradient operator V , is rigorously given by V = V r + n , (B.2) 136 Appendix B. The Ewald Sum for Dipolar Systems 137 where |r + n| is the reduced lattice distance. The dipole strength, LI is not reduced by the box length. The sum in E q . (B.l) is conditionally convergent, meaning that it depends on the order of summation. Deleeuw et al. [25] showed that in the spherical summation case, the conditional convergence, that is the order of summation, can be removed by a convergence parameter, f(s). The energy is now written as U(s), and USH c a n be found from USH = HmU(s). (B.3) s—>-0 For dipolar systems, the explicit form of U(s) is given by -U(s)L» = l im i £ £ J > , • V ) ( ^ . • V ) j — ^ - (B.4) '=1 i=i n' 1 1 where e~sn is the convergence parameter for spherical shell summation. We can continue with the derivation. This particular derivation is similar in spirit to the derivation for point dipoles in [25]. Following Deleeuw, we consider the substitution 1 1 f°° x-1/2e-\T+nfxdx. (B.5) |r + n| ypK Jo T h e energy can then be written as s^O Z^TT i = 1 j = 1 n , JO -i N N -oo + H ^ E E E ( * • V)(^' • V ) e _ S n 2 L *-1/2e-lr+n{2*dx, (B.6) s^O Zy/TT i = 1 j = 1 n , Ja2 where the integral has been split by introducing the parameter a2. T h e sum over n on the first term can be extended to include n = 0 as r —> 0 if this additional contribution is subtracted off. -U(s)L* = l im -L= £ £ J > , - • V ) f o • V ) e - 2 / X-We-W*dx ' 0 t = l j = l n ' a Appendix B. The Ewald Sum for Dipolar Systems 138 u Z V ^ j=i j=i n J0 V I — 1 J — 1 T h e last term in E q . (B.7) can be easily evaluated, giving l im {pt • V)(p3 • V ) f * x-V'e-W'dx = -44 ,^ (B.8) T-*0,]-H Jo 3 so that the energy expression is 1 N N -U{s)L* = l im — = ^ £ E ( t t • V ) ( ^ . • V ) c - » 2 / x - ^ e - l ' + n l ^ x + h ^ E E E ( * • v ) ( ^ , • v ) f 2 x - 1 / ^ - 2 - ^ ! 2 - ^ s u ^ V 7 1 ^ i=ij=i n 7 0 1 ^ , V « 3 + Now consider the following identity, e -|r+n| 2 y = ^ 3 / 2 e-ir2n2/y+2*in-r^ (B.10) n V n and the manipulation —sn2 — |r + n| 2 x = —sn2 — r2x — 2n • rx — n2x = —(s + x)n — 2xn • r — x 2 x r sxr 5 + x = —(5 + x ) n 2 — 2xn • r — S + X s + X . . , 2xn • r x 2 r 2 , s x r 2 = - ( s + x ) ( n 2 H H ^ ) V n s + x ( 5 + a;) 2 ; s + x = - ( s + x ) ' X n H r 5 + x 2 s x r 2 s + x ( B . l l ) Appendix B. The Ewald Sum for Dipolar Systems 139 We can now make the substitution to give the fourier space term, 1 N N /-co -U(s)L3 = l i m ^ E E E ( M r V ) ( M r V r n 7 2 ^ V | r + n | 2 ^ V i=l j=l n ' 01 + - f — 2 7 r m - r ^ - r 3 1 w 4 u 2 a 3 + ^ E ^ ' (B.12) iv/vr . = 1 3 Here we have separated the n = 0 term from the second sum. Consider the ri = 0 term, and the substitution u = x/(s + x), = i / s + « 2 u - i / 2 ( i _ S O T 2 + L(sur2y + ...)du s Jo 2 = I feu1'2 ^ - r 2 f - u 3 / 2 S + a 2 + 0 ( 3 ) ) 3 = -^ + ^ r 1 / 2 - 2 ^ + ^r3/2 + o(s). (B.i3) 5 or 3 a We can now write the energy expression as N N 1 J V J V /-co - t / ( , ) I 3 = l i m - i = E E E ( M r V ) ( M r V ) e - 7 Q 2 ^ V ^ d s V 8 = i j = i n ' Q V i — l 7 — 1 + ^ S ^ ' (B'14) When the gradient operator acts on 7(1 + ^ ) _ 1 / ' 2 ) the result is zero, since it is independent of r. The limit as s —> 0 is also zero. For the term ^-(1 + ^ ) ~ 3 / 2 , the gradient operators Appendix B. The Ewald Sum for Dipolar Systems 140 gives The limit can then be taken, resulting in N N r . 2 ( M , - V ) ( , V V ) ^ ( l + -ij)- 3 / 2 = |(1 + i)-"V< • f>,- (B.15) V i = l 7 = 1 1 * N 4TT 4 E E f ^ - ^ • (B-16) » = i j = i This gives, TV N s+xr dx 1 J V J V /•CO - c 7 ( , ) Z 3 = l i m - ^ E E E ( M r V ) ( M r V ) e - 7 2 V i=i j=i n ' 0 1 1 " " 4* Z i=l j = l 0 1 + — f 4 ^ V , (B.17) The limit s —> 0 can now be taken. We note two results relating to the first, and second integrals respectively, erfc(aa) i r . - v 2 ^ ^ , (B.i8) f x-2e~v2n2/xdx = - r ^ e - v 2 n 2 / a \ (B.19) The expression for the total energy of the central cell can be written as - w 3 = i E E E h ' V K . r V ) " ' * 1 ' 1 " 1 1 » = l j = l n ' lr + n l 1 N N A _£ ^ 47T 1 ^ , 4 ^ 3 ^ Q ' (B.20) 3 ' Appendix B. The Ewald Sum for Dipolar Systems 141 where the pair potential u(12) is 1 v - w w . . ^ e r f c ( a | £ + n|) u (12) = 7 E - ( ^ - V ) ( M 2 - V ) -n i p2irin-r/Le—K2n2/a2 + r 7 E - ( ^ i - v ) ( / , 2 - v ) r ' / V f ~ 2 • / 2 ^ n^O n 4vr 4u2a3 ,^ n . + u , a 7 - — . (B.21) 3 £ 3 ^ 1 V2 3 N L 3 \ ) In this expression, r, p and V are not reduced by the box length L. n is a lattice vector with length reduced by L. If the gradient terms are factored out, the pair potential can be written as, u(12) = - ( A * 1 - V ) ( ^ 2 - V ) * ( r ) 47r 4 u 2 a 3 / T , . + 3 L 5 " ' ^ - S b ' ( B ' 2 2 ) 1 e r f c f a l - + n h 1 ^in-r/L - K 2 n 2 / a 2 *w =jE r;-K, " + — £ — • <B-23> ^ n ' IL + n l ^ n#0 n In the limit that L —> oo, *P(r) behaves like 1/r, and we recover the pair potential for point dipoles, E q . (2.2.7). The term ^JA^I ' A*2 i s important and arises because we are considering a macroscopic sphere of dipolar material in a vaccuum. We might conclude that this term is due to the charges at the boundary of the sphere and the vaccuum. However, a more instructive identification of this term arises when we surround the spherical sample with a material with dielectric constant e . This situation will be discussed in detail in the next Appendix. A p p e n d i x C T h e R e a c t i o n F i e l d In this Appendix we consider the situation where a large, macroscopic sphere of dipolar material is placed within a continuum with dielectric constant e . The radius of the sphere is a, the instantaneous moment of the sphere is M s , and the moment is located at the center of the sphere. The dipolar material within the sphere will polarize the dielectric, and in turn M s will interact with the polarized continuum. Specifically, there will be a reaction field R at the center of the sphere due to the polarization of the continuum. For a spherical sample, the reaction field will have the same direction as M s , and is given by [52] _ 2 ( e ' - l ) M , R - ( 2 7 T T ) ^ - ( a i ) If the sphere of dipolar material was constructed from s cubic cells of length L, such that the cells formed a periodic lattice (Ewald periodic boundary conditions), then M s = s M , where M is the moment of one of the periodic cells. Also, the volume of the sphere can be written as V =fa3= sL3, (C.2) and the reaction field R can be written in terms of the moment of the central cell However, there is an energy cost in polarizing the continuum. That is, work must be done in order for the reaction field to form. We must take this work into consideration 142 Appendix C. The Reaction Field 143 when we calculate the reaction field's contribution to the energy of the central cell. We can accomplish this by slowly "turning on" M by the parameter A, A = 0 —> 1, M ( A ) = A M . For this value of A, the reaction field is R(A) = A R . The work required to increase M ( A ) to M ( A + dX) against the reaction field R(A) is dU(d\) = — A M • RdX, (CA) and the total work for the charging process A = 0 —• 1 is the reaction field contribution, URF, given by URF = — f1 A M • RdA Jo = -1M - R . (C.5) Here R F refers to the reaction field. We see that although the polarization of the dielectric results in a lower energy, there is a "cost" of polarization (by a factor of 1/2). If E q . (C.3) is employed, then the reaction field contribution to the energy of the central cell is U R F = - I M . ^ ^ - ^ - M . (C.6) R t 2 (2e + 1) 3L3 y ' The work done in forming R is £/RF, and it is also the free energy change for the process of "turning on" M and hence R . If the moment of the cell is M = A*i? URF can be written as a sum of pair potentials, given by "»<12> = -WTT)m*-»- (C'7) The total energy of the central cell will be the sum of its configurational energy, and the reaction field contribution, given by ^ c ( 1 2 ; e ' ) = u(\2) + u(\2)RF (C.8) where u(12) is given by E q . (2.3.7). A p p e n d i x D Q u a t e r n i o n P a r a m e t e r s In this Appendix, the method used to describe the orientation of a rigid body is discussed. A non-spherical body can be a point dipole within a sphere, a prolate or oblate ellipsoid, a spherocyfinder, or a "surfboard" shape. That is, a non-spherical body requires the its orientation as well as its position to be specified. For point dipoles, or volumes of revolution (prolate/oblate ellipsoids, spherocylinders), two parameters are required to specify its orientation. These are commonly the polar and azimuthal angles. For more general shapes, three parameters are required, and are usually the Euler angles, (M,V0 [14,11]. D . l Q u a t e r n i o n P a r a m e t e r s a n d the O r i e n t a t i o n of R i g i d B o d i e s T h e orientation of any three dimensional body can be described by "inserting" a cartesian axis frame within the body and specifying the orientation of the body axis frame relative to some laboratory axis frame. This can be done by an appropriate transformation. If the cartesian axis vectors in the body fixed frame are given by the columns of the second rank tensor A , then the transformation from body to space axis systems can be written as T = B * A B , (D.la) 144 Appendix D. Quaternion Parameters 145 where B is a rotation matrix (its transpose is denoted B ' ) , and in terms of Euler angles has the well known form [14], cos (j> cos ip — sin <p cos 9 sin ip sin <p cos ip + cos <f> cos 9 sin ip sin 9 sin ip — cos <p sin ip — sin (/> cos 9 cos i/> — sin <p sin ^ + cos <f> cos 5 cos ip sin 0 cos ?/> ^ sin <p sin 0 — cos <p sin 0 cos 9 j (D.lb) Likewise any vector p in the body fixed axis frame can be expressed in the space fixed frame by p' = B*p, (D.lc) where p' has been rotated by (9,cp,ip). Using Euler angles to specify particle orientations is not computationally efficient. Evaluating trigonometic functions is computationally expensive, and the equations of motion when expressed with Euler angles have singularities [14]. These problems can be avoided if the rotation matrix in written in a different form. A n y finite rotation of the vector p can be accomplished by a rotation about some suitable axis of rotation, given by a unit vector, n , as shown in Fig. D . l . This can be written as [11] p' = p cos $ + n ( n • p)[l — cos $] + p x n sin $. (D.ld) Here $ is the angle in the plane normal to the axis of rotation between p' and p. If we define the scalar qo = cos($/2), (D.le) and the vector q = qxi -f q2j + q3k such that q = nsin($/2), (D.lf) the rotation can be written as p' = p(?o - 9i - 4 - ql) + 2 c i ( q • P) + 2P * <w0. (D.ig) Appendix D. Quaternion Parameters 146 Figure D . l : Clockwise rotation of the vector p to p ' about the unit vector n. Appendix D. Quaternion Parameters 147 The elements of the rotation matrix can written in terms of q0 and q by comparing Eqs. D . l c and D . l d and writing the components of p' in terms of p, go and q, to obtain where ^ 9o + ?i -<&-<& % i ? 2 + 9o9s) 2 ( ^ 3 - q0q2) ^ 2(qiq2 - q0q3) ql - q\ + q\ - q\ 2(q2q3 + q0qt) y 2(^1^3 + 9092) 2(q2q3-qoqi) g 2 - q\ - q\ + q\ j qo = cos ^6>cos ^(<f) + xp) (D.lh) fD.li] qi = sin ^6 cos ^ - xp) , (D-lj) q2 = sin i(9 sin ^(<f> - xp) , (D.Ik) g 3 = cos^6>sini(<^ + xp) . (D. l l ) The terms {go, 9 i , 92,93} a r e known as quaternion parameters [14,29,60] and are con-strained such that 9o2 + 9 1 2 + 92 2 + 932 = l - (D.2) Briefly, a quaternion, Q , can be written as the combination of a scalar and a vector, [61] Q = (g 0 ,q), (D.3a) q = gii + g2j + g3k- (D.3b) The product of two quaternions can be defined as Q R = (g 0 r 0 - q • r, g 0 r + r 0 q + q x r) (D.3c) Appendix D. Quaternion Parameters 148 The square of the magnitude of a quaternion is | Q | 2 , where | Q | 2 = Q Q , (D.3d) and Q is the conjugate of Q , given by Q = ( ? o , -q ) . (D.3e) A p p e n d i x E Evalua t ion of the General ized Gaussian Overlap Integral T h e overlap integral of two three dimensional Gaussians can be written as /O O f O O /* O O i I I e-sW,:VV+TyAv-v){r-p)}drzdrydrx ^ ( E 1 ) - O O J — OO J — OO and evaluated by application of standard integrals. However, it is useful to give some intermediate results which are required in the expressions for the torques obtained below. Integration with respect to rz yields f°° e-{arl+br,+c)drz = ^ ^ / A a - c ? (E2a) J-co yja where a- = Tlzz + TUz , (E.2b) h- = 2(Tlxzrx + T3xz(rx - px) + Tlyzry + TJyz(ry - py) - TJzzpz) , (E.2c) -s = T3ZZP2ZZ ~ 2[T3xz(rx - px) + T3yz(ry - py)]pz . (E.2d) Collecting terms in descending powers of ry yields j e~{a r ; + & ry+c = _ 'aa t l e - ( a V S + 6 S W ) d r y = * 6 * / 4 a ' - c ' . ( E . 3 a ) \fa J-oo v aa' where a Syy ^zz S y z S S; zz S S S (E.3b) h ' ~ ^ + ^ , (E.3c) 149 Appendix E. Evaluation of the Generalized Gaussian Overlap Integral 150 b\ 2 S. 9 9 uyz uxz SZz (E.3d) b\ Hi 5 = 2 { ^ ( T J x z P x + riw,ft, + TJzzpz) - (T3xyPx + T^p, + T3zyPz)) , (E.3e) c s = Tjyyp2y-2Tjxy(rx-px)py-(r,-IzrE + r ^ r - , . - ^ ) - TJyzPy - TJ22pz)' <~>zz 2(TJxz(r- - P*) - TiyzPy)Pz + TjzzPl > (E.3f) and the tensor S is defined by S = T i + Tj . Finally, terms are collected in descending powers of rx to obtain / e - ( a r'+b r ,+c )d = * ^ / 4 a - c / i I / i II \/aa. J-co \/aa, a, (E.3g) (E.4a) where <s _ Q c _ C2 ' 15 ^>yy^>zz ^>yz (E.4b) = B p , = C : p p . (E.4c) (E.4d) The components of the vector B are Q Q _ Q2 ^>yyJzz l^yz Tj Tj Tj J xx J xy Jxz ^yx yy yz SZy SZ (E.5a) Appendix E. Evaluation of the Generalized Gaussian Overlap Integral 151 B„ - 2 OyylJzz >->yz Bz q O _ C2 T Ixy T Jyy Tn •> zy q q °yy Syz q uzx q °zy Szz T-x3xz Ti 3yz T-3 zz q uyx q uyy q °yz SZx S z y S z z (E.5b) (E.5c) T h e elements of the symmetric tensor C are Jyy C — T-^xy — J - j xy c. ~ M3xz Cyz — Tjyz = T-3 xx (SyzTjxz SzzTjxy)2 rp 2 1 3 x Z SZz(SyySzz — Syz) >->zz - T ]yy (SyzTjyz — SzzTjyy)2 rp 2 3yz SZz(SyySZZ — SyZ) Jzz = T 3 zz {SyzTjzz — SzzTjzy)2 Tj" S z z ( S y y S z z — Syz) Jzz (SyZTjxz — SzzTjxy)(SyzTjyz — S z z Tj ) Ti Ti Jxz Jyz Szz(SyySzz — Syz) Szz (Sy*Tjxz — SzzTjxy)(SyzTjzz — S z z T>zy) T- T 3 xz 3 zz S z z ( S y y S z z — Syz) S z z {SyZTjyz — SzZTjyy)(SyzTj z z — S z z Tlzy) T- T 3yz 3zz SZZ(SyySZZ — Syz) Szz (E.6a) (E.6b) (E.6c) (E.6d) (E.6e) (E.6f) Appendix E. Evaluation of the Generalized Gaussian Overlap Integral 152 E . l E v a l u a t i o n of the F o r c e s a n d T o r q u e s for the G a u s s i a n O v e r l a p I n t e g r a l Here we give expressions for the forces and torques needed in M D simulations. If a continuous potential is employed, the forces and torques can be found by taking the appropriate gradients of the potential. For the present model the force on particle i due to j is easily obtained using fn = - V p o - u f a ) , (E. la) and fji = -fn- (E-lb) The torque on particle i (at the origin) due to j (at p ) , T\-, (defined in the body-fixed frame of particle %) is given by [61] where Q* = [Qio, + 9i2J + 9i3*0] » (E.2b) and we have defined the operator „ / d d - d d A _ s V % o % i 5?,-2 % 3 ; The calculation of requires the evaluation of partial derivatives of the form dq ik /O O P O O / " O O _ . / / sT'ih : r r e - ^ T ' ^ + ^ ^ - P ^ ' - P ^ r ^ r , ^ , (E.3a) -co J—oo J—oo where k = 0,1, 2, 3 and the elements of T i f c are T : B = ^ . (E.3b) Appendix E. Evaluation of the Generalized Gaussian Overlap Integral 153 The partial derivative can be written as = -y2ses{ ,co foo rco 2^9['Ii.„+Tj:(i-p){-'-p)]dr:sdrydrx + *kxx J _ 0 0 J _ c o 7_oo ,00 roo yoo r2^s[Tr.rr+Tj:(r-P)^-P)]drzdrydrx + Tlkyy J_oo J_00 J_0O V i - c o i - c o J-oo 7_oo 7-oo •/-co ,oo yco i-oo -s[Tt:Tr+Tf.(v-p)(r-p)]drz(irydrx + 9J 1. I l l r x r z e " T-kXZ j _ c o J_co J - O O ,oo yoo yoo -siTr.rx+Tj-.ir-PHr-Pndrzdrydr^ . l*yz y_oo J-oo -/-oo (E.4) The required integrations are laborious but straightforward. It is useful to choose the integration order such as to simplify the final expressions obtained. For the integrals involving r 2 and rxry the most convenient order of integration is drzdrydrx (i.e., drz followed by dry by drx). This yields /oo ,oo ,00 , / / r2xe-{ar*+hr*+c)drzdrydr: -co J — oo J — oo T T 3 / 2 1 fl b" VaV7a" \ 2 + Aa" //2 / / / / 6 /4a — c (E.5) 7T 3/2 /CO ,oo ,co , / / rxrye-(ar*+br*+cUrzdrydrx = - o o «/ — O O « ^ — O O h\rx) 11 2d a" + a' a" \ 2 ' 4a" //2 / / // s6 /4a - c (E.6) Vad a" where a, 6, c, a', a", b' and 6" are as defined in the previous section. For the integrals with r 2 and r y r z , the most convenient order is drxdrzdry. This gives expressions of the form as the previous integration order, but with the parameters — T- +T-ixx i 3 xx ' (E.7a) Appendix E. Evaluation of the Generalized Gaussian Overlap Integral 154 a s a s b_ s q q _ C2 '-'zz'^'xx u xz q UXX IT^ + Tjl SzzSxx S. 2 ' xz b\ (ry) _ k i r + _ili s v + s ' q q 2 ( S uxzuxy Sx (E.7b) (E.7c) (E.7d) (E.7e) b\ Hi s 5* 2 I -fL(T3xxPx + T3xyPy + T i * z P z ) - ( r J x , P * + TiyzPv + T3zzP* Again, b"/s can be written as the dot product of two vectors, b - = B p , 5 (E.7f) (E.7g) where q q _ C2 '-''ra'-'.sz uxz B„ = 9 9 T.T.UZ SXz Bz - 2 q q _ q2 ^xx^zz uxz q uxx q ^ xy q ^xz T-3 xx Jxy T-X3xz SZx q uzy Szz q uxx q ^xy q uxz T 3yx T Jyy T 1^yz Szx SZy Szz Sxx Sxy q T-X3zx 3zy T 3 zz q uzx SZy Szz (E.7h) (E.7i) (E.7j) Appendix E. Evaluation of the Generalized Gaussian Overlap Integral 155 For the integrals involving r\ and rxrz, the integration order drydrxdrz yields equa-tions with a s a s - = T +T-s -Liyy T J-lyy , C C _ C2 >~>xx>~'yy uxy = 9 ' uyy IT^ + Tjl Q C _ Q2 ' ^xx'^yy ^xy ' 2 1 5 S 5 S = 2 S r 9 <? uxy uyz s yy b-f = 2(^(T3xypx + T3yypy+TJzypz)-(TJxxpx + TJyxpy + TJzxp2) Also, we have where - = B p , s S1XX Q uxy Q '-'xz - 2 p J D x ~ q q c-2 >->xx>~>yy >->Xy Q uyx Q uyy q uzy T 3 XX T-3xy T 3 xz (E.8a) (E.8b) (E.8c) (E.8d) (E.8e) (E.8f) (E.8g) (E.8h) B„ Q Q _ Q2 <Jxx*~>yy ^ x y q ^ XX q ^ xy q Syx Syy Szy Jyx Jyy T 3yz (E.8i) q ^XX q ^ xy SXz - 2 Bz = z q q q2 ^xx^yy >->xy Syx q °yy Szy T-3 zx Tj 3zy T, 3 ZZ (E.8j) Appendix E. Evaluation of the Generalized Gaussian Overlap Integral 156 T h e torque on particle j due to i, r^-, is given by 4 = - \ ^ ( v q i «(»i) - [ Q ~ V Q J U ( ? J ) ] Q , ) • (E.9) If we substitute r' = r — p in E q . ( E . l ) , V Q j u ( i j ) can be evaluated by interchanging the labels i and j and changing the signs of all components of p in all relevant equations. T h e numerical calculations can be simplified by noting that all integrals in the torque integrals contain the factor 7 r 3 / 2 7 2 e 6 " 2 / 4 a " - c " ^(p,q*,qj) = r^r, > yaa'a" which is independent of integration order and choice of origin. It is only necessary to calculate this term once, for example, when the pair energy is evaluated. This means that only parts of the prefactors occuring in these integrals must be recalculated in order to determine all the relevant integrals. This is also true of the j t h particle torques. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 89 a temperature scan is performed on a system.with b/a = 1.45, and p* = 0.792. There is a sharp rise in the polarization around T* = 1.9. Below and above this temperature, < P i > is near-zero. Above T* = 1.9, the system is simply isotropic, as indicated by < P 2 > « 0. However, for reduced temperatures below T* = 1.9, < P 2 >?s 0.8. This suggests that strong columnar correlations are present, and that the loss of polarization might be due to antiferroelectric columns similar to those observed in the dipolar cut-sphere systems [9,56]. The values of b/a, and p* that were chosen for this plot are not unique. Systems with aspect ratios between 1.4 and 1.5 and reduced densities ~ 0.8 all have similar behavior. In systems where the dipole lies along the symmetry axis of the particle, we can describe the local structure of the fluid by the pair correlation functions [9,10,56] g000(r) = < r * ~ r ) > , (5.2.4a) g ™ < r ) = 3 < 6(rj - TJ - r)[3(A,- • h3){fa • h3) - fa • fa] > ^ ( 5 _ 2 > 4 b ) y i i o ( r ) = < Eh3,^ S ( r , r ) [ f a • fa] > _ { 5 2 A c ) g000(r) is the usual radial distribution function. The average interaction energy between a pair of dipoles separated by a distance r is related to g112(r). The function g110(r) gives an indication of parallel or antiparallel dipolar correlations as a function of the intermolecular distance. Also, we can define parallel and perpendicular correlation functions, denoted by <7||(r||) and g±(r±), respectively. If the instantaneous director is given by d , then r|| and rj_ are vectors such that r*|j = (r • d ) d and r - j . = r — r||. The correlation functions are defined as = <B(r||,r| l + ft1|)> " i i v < n 0 ( r | | , r | | + 6ru) > v '
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Computer simulation study of ferroelectric liquid crystal...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Computer simulation study of ferroelectric liquid crystal phases Ayton, Gary Steven Douglas 1996
pdf
Page Metadata
Item Metadata
Title | Computer simulation study of ferroelectric liquid crystal phases |
Creator |
Ayton, Gary Steven Douglas |
Date Issued | 1996 |
Description | Molecular dynamics and Monte Carlo simulations are used to study ferroelectric liquid crystals. The effect of molecular shape on the stability and nature of ferroelectric phases is examined. The effect of molecular flatness on the stability of uniaxial nematic liquid crystal phases is studied. Most liquid crystal molecules are elongated and have a degree of flatness. However, most models used in computer simulation of liquid crystals are axially symmetric. A series of molecular dynamics simulations, specially designed to isolate the effects of molecular flatness on liquid crystal phases, are performed. Molecular shape is approximated with a generalization of the Gaussian overlap model [ B. J. Berne and P. Pechukas, J. Chem. Phys. 56, 4213 (1972)]. The form presented here has been extended to include ellipsoidal particles with non-degenerate semi-axes. It is found that small amounts of molecular biaxiality can drive an isotropic to nematic phase transition. Simulations of randomly frozen and dynamically disordered dipolar soft spheres are used to study ferroelectric ordering in spatially amorphous materials. Systems where the dipole moment has 1, 2, and 3 components are considered. It is found that the 1 component (Ising) model has ferroelectric phases. The systems with 2 and 3 dipolar components form disordered phases at low temperatures. A ferroelectric phase diagram is constructed for oblate molecules with point dipoles embedded along the particle symmetry axes. The role of particle shape on the stability of ferroelectric liquid crystals is examined with molecular dynamics simulation. Molecular shape is modeled with the generalized Gaussian overlap. Ferroelectric phases are found in systems with weak (short-ranged) columnar correlations. Systems with long-ranged columnar order are found to be antiferroelectric. It is found that the stability of ferroelectric phase is very sensitive to details of molecular shape, and exists only in small regions of the phase diagram. A more robust model for a ferroelectric liquid crystal is developed. Monte Carlo calculations are used to examine ferroelectric order in fluids of disc-shape particles with embedded dipoles. The dipoles are uniformly distributed over a circular "patch" of finite size, placed in the central plane of the particle. It is shown that such systems may undergo spontaneous polarization to form a stable ferroelectric discotic nematic phase. |
Extent | 12482777 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-19 |
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.0059610 |
URI | http://hdl.handle.net/2429/6239 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1996-14724X.pdf [ 11.9MB ]
- Metadata
- JSON: 831-1.0059610.json
- JSON-LD: 831-1.0059610-ld.json
- RDF/XML (Pretty): 831-1.0059610-rdf.xml
- RDF/JSON: 831-1.0059610-rdf.json
- Turtle: 831-1.0059610-turtle.txt
- N-Triples: 831-1.0059610-rdf-ntriples.txt
- Original Record: 831-1.0059610-source.json
- Full Text
- 831-1.0059610-fulltext.txt
- Citation
- 831-1.0059610.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-0059610/manifest