C O M P U T E R SIMULATION STUDY OF FERROELECTRIC LIQUID CRYSTAL PHASES By Gary Steven Douglas A y t o n B. Sc. (Chemistry) U . B . C. 1990 A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T OF T H E R E Q U I R E M E N T S 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 F A C U L T Y 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 degree at the thesis in partial fulfilment of University of British Columbia, I agree freely available for reference copying of department this or publication of and study. thesis for scholarly by this his or her Department of The University of British Columbia Vancouver, Canada requirements that the I further agree purposes representatives. may be It thesis for financial gain shall not permission. DE-6 (2/88) the that advanced Library shall make it by the understood be an permission for extensive granted is for allowed that without head of my copying or my written Abstract Molecular dynamics and Monte Carlo simulations are used to study ferroelectric liquid crystals. T h e effect of molecular shape on the stability and nature of ferroelectric phases is examined. T h e 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 . C h e m . Phys. 56, 4213 (1972)]. T h e 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. T h e 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. T h e 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. n Systems with long-ranged 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. T h e 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 E q u i l i b r i u m Statistical Mechanics 7 2.2 Intermolecular Potentials 10 2.3 Computer Simulation Techniques 12 2.3.1 Periodic Boundary Conditions 13 2.3.2 T h e Treatment of Long-Ranged Potentials Under Periodic Bound- 2.3.3 2.3.4 3 ary Conditions: T h e Ewald Sum for Dipolar Systems 13 T h e Monte Carlo Method and Importance Sampling 16 Molecular Dynamics Method 20 A Simple Model Potential for Liquid Crystal Molecules 23 3.1 25 Generalized Gaussian Overlap Integral iv 3.2 3.3 M D Simulation Results: T h e Effect of Molecular Biaxiality on Liquid Crystal Phases 29 3.2.1 34 Molecular Biaxiality and the Uniaxial Nematic Phase Summary 45 4 Ferroelectric Phases of Amorphous Dipolar Systems 4.1 46 M o d e l and Simulation Details 48 4.2 Results for Randomly Frozen Systems 49 4.3 Decoupled M D Simulation Technique 53 4.4 A M e a n Field Theory: Isolating the Effect of Long-Ranged Interactions 4.5 Summary . 63 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 T h e Cut-Sphere M o d e l 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 E Quaternion Parameters and the Orientation of Rigid Bodies 144 Evaluation of the Generalized Gaussian Overlap Integral 149 E. l Evaluation of the Forces and Torques for the Gaussian Overlap Integral vi . 152 List of Figures 3.1 A sketch showing the molecular axes (x',y',z') a 2 3.2 and the principal radii o i , and <T which characterize the ellipsoid model 27 3 (|r(i*) - r(0)| ) for 2a 2 x : 2a 2 : 2a 3 = 2.75 : 3 : 9 (dash line) and 2a x : 2a 2 : 2<7 = 1.5 : 3 : 9 (dotted line). These systems are isotropic and uniaxial 3 nematic phases, respectively 3.3 32 Order parameters as a function of particle thickness at fixed reduced temperature and density. T h e squares are the uniaxial order parameter QQ Q and the circles are the biaxial order parameter Ql . 2 T h e error bars repre- sent one standard deviation 3.4 35 Snapshot of an instantaneous configuration for an isotropic phase with 2ax : 2 a : 2<r = 2.875 : 3 : 9 2 3.5 36 3 Snapshot of an instantaneous configuration for a uniaxial nematic phase with 2<7i : 2<r : 2cr = 2.125 : 3 : 9. Here, the director is along the z' axis 2 3 (up and down on the page) 3.6 37 Snapshots of instantaneous configurations for a biaxial nematic phase with 2(7i : 2<r : 2cr = 0.5 : 3 : 9. Here, the directors are along x' and z'. T h e 2 3 snapshots are taken along the axes perpendicular to the z' director. 3.7 . . . <jf (r*) for biaxial 2<7i : 2cr : 2cr = 1 : 3 : 9 (solid), uniaxial nematic 2a 000 2 3 x 38 : 2<r : 2cr = 2 : 3 : 9 (short dash), and isotropic 2o\ : 2 a : 2<r = 3 : 3 : 9 2 3 2 (long dash) systems 3 40 vii 3.8 < P (x' (0) 2 1 • x' (r*)) 2 > for biaxial 2a : 2a x 2 = 1 : 3 : 9 (solid), : 2a 3 uniaxial nematic 2<7i : 2<r : 2cr = 2 : 3 : 9 (short dash), and isotropic 2 2(7! : 2a 2 3.9 3 : 2<r = 3 : 3 : 9 (long dash) systems 41 3 < x'(t*) • x'(0) > for a system with 2a x : 2<r : 2<r = 2 : 3 : 9 (solid), 2 3 and a system with 2o\ : 2cr : 2cr = 1.5 : 3 : 9 (short dash). T h e inset is 2 3 In < x(t*) •x'(0)> 4.1 42 T h e radial distribution function g(r*) for soft spheres at p* = 0.8 and T* = 10.5 4.2 T h e polarization < P > at T^ 51 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 4.3 T h e root-mean-square spin length, < S 52 r m s >, for the iV = 256 frozen X Y Z model 4.4 54 < 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 4.5 < P 58 > 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 4.6 < P > vs T * (rotational) for the Ising model. 59 Results are shown for dynamically random systems with m * = 1 (squares) and m* = 5 (triangles) and for the randomly frozen case (circles). T h e heat capacities per particle, C /N, v are plotted vs T* (rotational) in the inset vm 60 4.7 T h e mass dependence of the disordered to ferroelectric transition temperature Tp (rotational). T h e squares and triangles are results for the and XYZ models, respectively. T h e values of T F of the heat capacity, C /N, were obtained from plots vs T * (rotational) and a typical example is v shown in the inset. XY T h e error bars represent estimated uncertainties in the peak position 4.8 62 < P > vs T* from the reaction field, as found from a simple mean field theory. T h e short dash, solid and long dash lines are for the XYZ, and Ising models, respectively. XY, T h e circles are the frozen TV = 256 Ising results 4.9 67 T h e reaction field contribution to the energy per particle < U^ > /TV for F the Ising model, as calculated from simple mean field theory. T h e inset is the reduced heat capacity per particle C /N 69 v 4.10 T h e reaction field contribution to the energy < U^ > /TV for the model, as calculated from simple mean field theory. T h e inset is the re- F XYZ duced heat capacity per particle C /N 70 v 4.11 T h e short and long-range (reaction field) contributions to the energy per particle < U* > /TV for the dynamically decoupled XYZ model. solid line corresponds to < U* > /TV, the short dash is < U^p the long dash is < Ug T > /TV. The > /TV and T h e squares, triangles and circles are for m* = 1, 5 and 10, respectively 4.12 72 T h e short and long-range (reaction field) contributions to the energy per particle < U* > /TV for the dynamically decoupled XY line corresponds to < the long dash is < U$ T U* > > /TV. model. T h e solid /TV, the short dash is < U^ F > /TV and T h e squares, triangles and circles are for m* = 1, 5 and 10, respectively 73 ix 4.13 T h e 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, long dash is < Ug > /N. T the short dash is < UR F > /N and the T h e squares, and triangles are for m * = 1 and m* = 5, respectively 5.1 74 A sketch of a tetragonal I lattice. T h e circles are "top" views of columns of particles. T h e 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. T h e circles represent columns of particles. There is no interdigitation between columns of particles 5.3 Dipolar energy per particle U^, /N for different lattices of oblate dipo- D lar particles. 80 Triangles are for an antiferroelectric columnar lattice, and squares are a ferroelectric tetragonal I lattice 5.4 T* — 1.9 results for the < P 2 82 > 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 5.5 86 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 5.6 b/a = 1.45 p* = 0.792 < P 2 87 > (circles) and polarization < P i > (triangles) as a function of reduced temperature, T* 5.7 g (r*) 000 (solid), g (r*) 112 88 (short dash) and g (r*) no (long dash) for a p* = 0.792, b/a = 1.45, T* = 1.35 system 5.8 0||(r[j) (solid) and g (rl) ± (dash) for a p* = 0.792, b/a = 1.45, T* = system 91 1.35 93 x 5.9 A computer "snapshot" of an b/a = 3, p* = 0.792, T* = 2.5 system. The phase is antiferroelectric columnar 5.10 g (r*) 000 (solid), g {r*) 94 (short dash) and g (r*) 112 lw (long dash) for a / = 0.792, b/a = 1.45, T* = 1.9 system 5.11 flf||(rfj) (solid) and g±(rl) 95 (dash) for a p* = 0.792, b/a = 1.45, T* = 1.9 system 5.12 g (r*) 000 96 (solid), # 112 ( r * ) , (short dash) and <? (r*) (long dash) for & p* = no 0.792, 6/a = 1.45, T* = 2.5 system. 5.13 <7 (rf|) (solid) and g (r* ) N ± 97 (dash) for a/)* = 0.792, b/a = 1.45, T* = 2.5 ± system 5.14 98 A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 1.35 system. T h e phase is antiferroelectric columnar 5.15 100 A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 1.9 system. T h e phase is ferroelectric 5.16 101 A computer "snapshot" of an b/a = 1.45, p* = 0.792, T* = 2.5 system. T h e phase is isotropic 5.17 102 T h e reduced energy per particle, < U* > /N as a function of reduced temperature, T*. T h e inset is the excess contribution to the reduced heat capacity per particle, < C v 5.18 T * = 1.9 , p* = 0.792, < P 2 > /N >, (circles) and polarization < P 103 x > (triangles) as a function of aspect ratio, b/a = width/height 5.19 A rough phase diagram for dipolar, p* — 3, oblate ellipsoids. 105 T h e re- duced density is p* = 0.792. T h e width of the ferroelectric phase is only estimated. T h e phases are I (isotropic), F E L(ferroelectric liquid crystal), A F E C L (antiferroelectric columnar fluid), and S (solid) 106 6.1 T h e coordinate system used to evaluate the interaction potential. The relevant vectors defined in the text are shown. T h e inset is a cross sectional view of a cut-sphere with point dipoles dispersed over a circular patch of diameter D' 6.2 Ill T h e mean square displacement, ( ( R ( r ) - R ( 0 ) ) ) , for T* = 0.02 (short 2 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 6.3 114 T h e temperature dependence of the polarization < P > and reduced constant volume heat capacity per particle C /N v D/L = 10, D'/L = 9 and p* = 0.525. (inset) for a system with T h e error bars represent one estimated standard deviation 6.4 "Snapshot" of a configuration for a system with D/L T h e system with D'/L 6.5 phase. Note that this configuration consists of loose polarized columns polarization in this system g ( *) s r = 10 and p* = 0.525. = 5 is shown in the unpolarized ( T * = 0.0275) with nearly equal numbers up and down. 6.7 = 10 and p* = 0.525. = 9 is shown in the unpolarized (T* = 0.050) phase. 117 "Snapshot" of a configuration for a system with D/L T h e system with D'/L = 10 and p* = 0.525. = 9 is shown in the ferroelectric (T* = 0.025) phase. 116 "Snapshot" of a configuration for a system with D/L T h e system with D'/L 6.6 115 O n average, there is no net 118 f ° T* = 0.025. No long-range smectic order is detected in this system. 120 r xn 6.8 T h e variation of the reduced electrostatic potential energy, L u (12)/p?, 3 es as one particle is slid over another with their dipoles directed in the same direction and their surfaces in contact. T h e inset shows the same curves on a magnified scale. Curves are shown for single point dipoles (D'/L (long dash) and for D'/L 6.9 squares are the total energy < U* > /N, F > /N, 0) = 4 (short dash), 5 (dotted) and 9 (solid). . . . T h e contributions to the reduced energy per particle, < U* > /N. contribution < U^ ~ 122 The the triangles are the reaction field and the circles are the short-range structural contribution (everything else) < Us T > /N. T h e reduced density is p* = 0.525 123 A.l T h e 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 xin Table of Symbols Symbol Description b/a T h e aspect ratio of an oblate ellipsoid. b Shape tensor. B Rotation matrix. Cv T h e constant volume heat capacity. D' T h e diameter of the charge distribution embedded in the cutsphere. d T h e director unit vector. g(r) T h e radial distribution function. <7||(r||), <7|_(r|_) Parallel and perpendicular pair correlation functions. &B T h e Boltzmann constant. K Kinetic energy. L Box length. m,m' T h e mass and reference mass of the particle. M T h e instantaneous moment of the sample, n A n integer lattice vector . < P > Average polarization. < Pi > T h e average first rank order parameter. < P T h e average second rank order parameter. 2 > QNVT T h e partition function for the canonical ensemble. Ql Uniaxial nematic order parameter. Q22 Biaxial nematic order parameter. . 0 xiv qo, qi, ?2, 93 Quaternion parameters. R T h e reaction field. 5 Gaussian overlap "hardness" parameter. T h e root-mean-square spin length. r T h e absolute temperature. T Orientation-shape tensor. u(12) A pair potential. [/ Instantaneous configurational energy. <U> IN T h e average energy per particle. < U > IN T h e average reaction field contribution to the energy per par- RF ticle. < U ST > IN T h e average structural contribution to the energy per particle. V Volume. ZNVT T h e configurational integral. a T h e convergence parameter for the Ewald sum . a\ Moment of inertia parameter. 1/fcBT V Gradient operator. £ T h e reference unit of energy. e T h e dielectric constant . A* T h e 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 T h e reference unit of distance. 0 i 5 0"2> o"3 Ellipsoid principle radii. xvi 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, m y 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. H e 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 D r . Michel Gingras for the opportunity to work with h i m , and for his infinitely long and valuable comments. Although I a m quite sure that Michel and I agreed on possibly nothing, I found our collaboration to be challenging, and very rewarding. O u r 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. D r . D a n Berard saved m y life. In m y formative years as a theoretical chemist I met m y nemesis, the computer. If it were not for D a n , I would have gone mad. I don't know how many times I said "hey D a n , how do I..." D a n 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 D a n . However, I still hate computers. xvn Finally, 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 w i t h his M O V I E program. A picture is worth a 1000 words. A n d as I have found out i n the course of w r i t i n g m y 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, C l a u d i o , R i c h a r d , M a r k , John, L i a m , and Frank, it's been fun under the p a l m tree. xvin Chapter 1 Introduction W h a t 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. W h e n 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. B u t 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. B u t what is grease? O r 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 throughout, 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. O f 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. T h e 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. T h e fluid phase has short-ranged positional and orientational order. A t 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. A t 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. T h e 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. T h i s 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 G e r m a n 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. T h e cloudy blue fluid was similar to crystals in that it interacted with polarized light, but it also had fluid-like flow properties. A r o u n d 1916, B o r n [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. T h e 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. T h e 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]. T h e y 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. F r o m 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. T h e results of these simple hard models suggest that the anisotropic "shape" of liquid crystal molecules plays an important role i n 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. T h e 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. T h e unique properties of liquid crystals, namely their solid-like orientational order, but fluid-like flow properties, generates very interesting behavior i n 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 i n 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 i n 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 i n the solid phase since the particles cannot rotate, and not easy in the isotropic phase since the aligned state is not stable. to isotropic phase. Once the field is turned off, the ordered fluid would return 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 i n 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. 5 Introduction A ferroelectric liquid crystal would consist of molecules with permanent dipole moments 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. T h e closest candidate is the chiral smectic C * phase, where the particles 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 crystal phase has not been observed experimentally. W h y ? There are two possible reasons. Firstly, such systems do not exist. Secondly, none has been found. A t least theoretically, both of these statements have been challenged. Recently, computer simulations of dipolar hard and soft spheres have been found to exhibit a stable ferroelectric nematic phase [9,10]. T h e orientational order comes from purely dipole-dipole interactions. F r o m a theoretical 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 i n 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? els other than dipolar spheres be found that have ferroelectric Can mod- liquid crystal phases? Is Chapter 1. Introduction there a possibility of constructing 6 a ferroelectric liquid crystal in a laboratory? Com- 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 mechanical 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 amorphous 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. T h e 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. Chapter 2 Statistical Mechanics a n d C o m p u t e r Simulation of F l u i d s Computer simulation of fluids is one application of equilibrium statistical mechanics. In this Chapter, we will discuss the necessary tools to study equilibrium fluids by computer simulation. E q u i l i b r i u m Statistical Mechanics 2.1 In a classical fluid, the Hamiltonian for a system of N particles can be written as H (^p) = K (p) + U (q), N N N (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 = ( q i , q , ...q./v) and the corresponding set 2 of conjugate momenta is p = ( p i , p , •••PN)2 T h e 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/Q r and Q r is large. Here, T is the absolute temperature is the characteristic temperature of rotation [13]. O f 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]. T h e time evolution of is determined by Hamilton's equations [11,12,16] , as given by <9H/v q, = - g j j - , 7 (2.1.2a) Chapter 2. Statistical Mechanics and Computer Simulation Pi where, i refers to the i th of Fluids 8 (2.1.2b) = 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. O f 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. of 12iV variables. In the general case, iJ/v(q, p) is a function These variables are simply coordinates i n 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, (2.1.3) However, as an alternative to calculating mechanical observables as time averages, Gibbs suggested removing the time dependence and instead, calculating ensemble averages. 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. E a c h member of the ensemble has a corresponding phase point, and the ensemble forms a "cloud" of points i n phase space. T h e 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), than others. ical ensemble. then all phase points are accessible, but some are more probable We will only consider the case of constant (NVT), known as the canon- T h e 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, Chapter 2. Statistical Mechanics and Computer Simulation where dc\ — dqidc[2---d(\N, and dp = dpidp2-..dpN- of Fluids 9 T h a t is, if the number of m e m - bers of the ensemble is S, the number i n the small volume element dqdp is given by SPNVT{<\-> p)dqdp. T h e function PNVT{^ p) is the probability density for the canonical ensemble. Mechanical observables can then be calculated from J W ( q , p) NVT(q, <W>= P (2.1.4) p)dqdp, where PNVT(<1, p) can be written as [12] ( ^\ PNVT(q,p) = n 1 ' 1 e x P(-^Af(q,p))/fcBT) ^ • (2.1.5) Here, QNVT is the partition function and has the form, J dqdpexv(-H (q,p))/k T), QNVT = -^j^ where k B N (2.1.6) B is the Boltzmann constant. T h e 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.) T h e factor of of the indistinguishability of the particles, and the factor of h~ xN takes account is from quantum cor- respondence [13]. T h e 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. T h i s gives exp(-U (r,n))/k T) ~ N PNVT = B , (Z.l.t) ZJNVT where the set of generalized coordinates q = ( q i , q , ...qyv) has been split into a set that 2 describe particle positions, r = ( r i , r , • ••YN) and a set that describe particle orientations 2 $7 = ( f i i , f i , . . . O N ) . T h e term ZNVT is the configurational integral, and is defined as 2 Z N V T = J J drdnexp(-U (r,n))/k T). N B (2.1.8) Chapter 2. 2.2 Statistical Mechanics and Computer Simulation of Fluids 10 Intermolecular Potentials In real molecular systems, the configurational energy U arises from all molecular interactions. In general, for TV particles it can be written as [14] = £«(»)+ £5>(*i) + E E E «&•*) + ..., i i j>i i j>i (2-2.i) k>j>i where the sums include distinct pairs, triplets and higher terms. T h e first term arises from any external fields. T h e second term is from pairwise interactions and usually dominates the total energy. T h e third and higher terms come from many-body interactions. Threebody 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 where u rep rep (12) + 1 ^ ( 1 2 ) , (2.2.2) ( 1 2 ) 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. T h e short-range repulsion originates from the overlap of the charge clouds of the two molecules, and is formally a quantum mechanical problem [16,17]. T h e 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. T h e 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. T h e next class of simple model repulsive potentials are continuous functional forms, written as Uct (12) = e / ( - ) , (2.2.4) where, a is the unit of length and is roughly the distance of closest approach for the two particles. T h e intermolecular distance is r = | r | , and e is the fundamental unit energy. For spherical particles, a simple example is the "soft" sphere potential, u (12) = ( - ) . r (2.2.5) 1 2 ss £ T h e 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. T h e 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. T h e it]r(12) contribution is usually a sum of terms, representing dispersion and multipolar 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 incorporates dispersion interactions is the Lennard-Jones (LJ) potential written as u (12) = As ((^) LJ 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. T h e interaction energy of two point dipoles, p,\ and /* , separated by r is 2 given by the well known expression, u {\2) DD = /x -V^(r) 2 3(A*I • r ) ( > 2 • r) p, • p x where 2 (2.2.7) <f>(r) is the potential at // , due to p,\. A simple derivation of the pair potential 2 between two point dipoles is given i n Appendix A . 2.3 C o m p u t e r Simulation Techniques There are two basic techniques used to simulate classical fluids. These are the Monte Carlo ( M C ) and Molecular Dynamics ( M D ) methods discussed below. Very briefly, i n Monte Carlo simulations, configurations are generated according to PNVT-, so that averages can be calculated from M i < >=uE ^ W W (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, t=i 1V1 where W . is the instantaneous value at time t. f (2.3.2) Chapter 2. 2.3.1 Statistical Mechanics and Computer Simulation of Fluids 13 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. T h e 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. T h e 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 ~ 1 0 23 system. T o 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. T h e small system is contained in a "box", or "cell", and the cell is repeated i n a lattice. A s 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 longranged periodic images is negligible. However, if the potential is not short-ranged as i n 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 i n the next section. 2.3.2 The Treatment of Long-Ranged Potentials Under Periodic Boundary Conditions: The Ewald Sum for Dipolar Systems T h e 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. T h e 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/r n 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 U N OO = oT, 2 n N (2.3.3) E E " ( l + ^l r n t=i i=i where L is the length of the cell, and r = r ; — r;. T h e vector n designates an integer 3 simple cubic lattice point with n = {n L, x n L, y n L) z and n , n , r y n z integers. T h e 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 i n 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). T h e lattice sum is conditionally convergent [25], meaning that the value of U will depend on the order that the periodic images are summed. T h e 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. shell becomes "smoother". A s the radius of the sphere becomes larger, the 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. T h e energy is now written as U(s), and USH can be found from 1 CO ' EI>(|r 2- (2.3.4) + nI|)/( ) 3 t= l j = l U SH = 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 , /*3, ••• MAT}, in a cubic box of length L under periodic 2 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 where v^^(\% e) ^EIX D (12; '), B C (2-3.6) e is an effective pair potential and has the form PBC, uDD = \ 12; e) = - ( / i x • V)(/x • V ) * ( r ) 2 4TT + 4^ a 2 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. T h e 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 T h e 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. T h e 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 2.3.3 • p>2- T h e Monte Carlo M e t h o d and Importance Sampling T h e 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 JJWe-P drd£l u <W> = (2.3.9) (2.3.10) 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), Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 17 and for each one, calculating W ( r , fi) and weighting it by e $ ( & ) . T h e estimate of u v < W > is then 1 T M < W > w-e- < pu (2-3.11) ^ = M However, calculating ensemble averages with this crude method of Monte Carlo integration is impractical, since a good estimate of ZNVT large set of configurations. would require generating a V E R Y 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). figurations are selected from PNVT, If con- then ensemble averages can be calculated from E q . (2.3.1). Now the problem is how to sample configurations from the distribution, PNVTT h e method used to sample PNVT was originally developed by Metropolis et al. [28] and involves constructing a Markov chain of configurations. T h e 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 , C 3 . . . . C M } such that the probability of finding the system in a state n at 2 some point c; in the sequence depends only on the state m that the system occupied at C j _ i - T h 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 PNVTGenerating new configurations in the Markov chain is done by choosing an appropriate transition probability, ir , mn defined as the probability of going from the state m to n in one "step". T h e n the probability that the system is in a state n, given by p = n 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. T h e configuration associated with the n th state is written as ( r ,fi ).This equation can be conveniently written i n the matrix n n 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( k,^k)- TT r 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 ir . First we impose the constraint of mn "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 TTmn ""mm = = OC p mn n p m Pn < Pm = «mn(^) 1 ^ 771 ^ TI m^U (2.3.15) Z / n ^ r a ^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). W e 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. T h e move is either Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 19 accepted, and the new configuration is the next i n 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 25r max centered at i*™. T h i s corresponds to a stochastic matrix with elements Qt-mn — if MR if r " is within the cube, and a = 0 if r™ is outside the cube. MR is simply m n all the (finite) positions available within the cube. T h e change in the energy of the system is calculated by considering the energy of the i th 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 p /pm, n in this case given by Pn _ Z = exp(-/W). Pm N V T exp(-/3U )dr n Z exp(-pU )dr NVT m (2.3.17) T h e volume elements dr have been explicitly included, however they are identical i n the m and n states, and therefore cancel. T o 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 [ 1 1 , 1 4 ] . 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 dSl ? = am 0?d6? dip? d<p?. 7 (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 p Zy n N p + 8U) sin Q\ exp(-/3U T m sin Of Z exp(-f3U ) m NVT e x = 1 m P(-W^f • (-- ) 2 3 19 In this form, random orientational displacements have been done by small changes i n the Euler angles themselves, 8(f>, 89, and 8ip. However, if random displacements are i n cos Of rather than Of, this amounts to taking displacements i n a new variable £, such that £f = cos Of, and d£f = — sin0™dcos Of. T h e 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 Molecular Dynamics Method Molecular Dynamics solves the equations of motion for a system of N molecules interacting through a potential U as i n 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 conserved. T h a t is, the Hamiltonian is a constant of the motion. If the particles interact with a continuous potential, the force on the i th Fj = —V£/j, where Ui is the i th particle is F,- and can be found from 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 T h e total force F j , acting on the i th the i th particle due to j is f(ij) particle is then F j = YljLi,i^j f ( u ) where the force on = —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). that is used is described i n Appendix E . T h e method 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) canonical (NVT) ensemble. If we wish to sample the ensemble, the temperature T rather than the total energy E must be constant. T h e equations of motion for a classical particle i n 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- T h e constraint of constant kinetic temperature can be written as y E | r | l 2 - ^ which assumes the equipartition of energy. B T = 0, (2.3.21) From this, an equation constraining the accelerations is N m £ r , - • £ • = (). j=i This equation is simply the time derivative of E q . (2.3.21). (2.3.22) Following Evans [32], the constrained accelerations r , , in terms of the unconstrained Newtonian accelerations, c r'i" = F , / m can be written as f,- c = r-7-A-. m (2.3.23) Chapter 2. Statistical Mechanics and Computer Simulation of Fluids 22 T h e constraint on the accelerations is then N mVr,--(f,- -A^) = u 0, (2.3.24) implying A = t?= l F : (2.3.25) Here, F,- is the total force acting on particle i. T h e constrained acceleration on the particle can be written as *• = M - §fe!*-) F ( ( 2 - 3 2 6 ) For particle rotations, the equation of constraint is if;i:a;,-u; --{iV^r t where a>- is the angular velocity if the i th 4 = 0, (2.3.27) particle in the body fixed axis frame, and I is the diagonal moment of inertia tensor. T h e number of rotational degrees of freedom / is 2 for a linear particle and 3 for a non-linear particle. T h e constraint on the angular acceleration is « ? = <J?_Au;V, (2.3.28) X ) I : w , - ( w - A u ; ) = 0, t=i (2.3.29) N i i A = tit* .""*'. 1 E,-=i I: (2.3.30) «tWi 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. Chapter 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 Molecules L i q u i d crystal molecules are generally elongated, and have a degree of flatness due to the presence of rigid aromatic rings [4]. T h e 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 . T h e 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. T h i s 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. T h e use of axially symmetric particles as a simple model for liquid crystals neglects the effects of particle flatness. introducing particle flatness It would be interesting to know if, and to what extent, 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. T h e 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. T h e potential decays as a scaled Gaussian where it's m a x i m u m 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". O f 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. T h e repulsive form presented here was specifically chosen to isolate the effects of particle shape on liquid crystal phases. T h e strength of the Gaussian overlap model lies in its ability to represent molecular shape and i n it's computational convenience i n M D simulations. T h e calculation of the forces and torques necessary i n 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. T h e 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 i n the formation of uniaxial phases. As the particles are made flatter still, a uniaxial to biaxial transition is also observed. T h e generalized Gaussian overlap model will be used in later chapters to study the effect of molecular shape on ferroelectric liquid crystal phases. 3.1 Generalized Gaussian Overlap 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) = e / s 7 2 s T : ( -P)( -P) . r r (3.1.1) T h e parameter s determines the m a x i m u m of the Gaussian when r = p . W h e n r lies on the surface of the ellipsoid defined by 2T:(r-p)(r-p) = l , (3.1.2) then j(r)= . 7 T h e value of 7 will be chosen after the overlap integral is calculated. (3.1-3) T h e 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'bB , (3.1.4a) Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 26 where b is the tensor / l/2o-l 0 0 V l/2a o \ 0 (3.1.4b) 0 2 o i/2a 3 2 ; T h e parameters <Ti, a , and <7 are the principle radii when the ellipsoid is cast in standard 2 3 form (see Fig. 3.1). B is a rotation matrix (its transpose is denoted B*) which expressed in quaternions has the form [14] ( 2(<?i<?2 + 9o9s) ql + q\ - q\ - ql 2(^143 + qq) 0 0 2(9293 - 9o9i) 2 2 ^ 2(9293 + 9091) ql - q\ + ql-ql 2(9192 - qoqz) \ 2(5^3 - q q ) (3.1.4c) 9o - 9? - 9 ! + 9l ) where (3.1.4d) 90 = cos ^0cos ^(<f) + ip) , 1 1 q = sin -0 cos x (3.1.4e) -((f) - ip) , 1 1 q = sin -6 sin-(<f> - xp) , (3.1.4f) 2 93 = cos ^ s i n ^(<f> + ip) . (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. T h e overlap integral can be written as co /•co rco / -CO J — OO J — OO r _ , m / w ,. Chapter 3. A Simple Model Potential for Liquid Crystal Figure 3.1: A sketch showing the molecular axes (x',y',z') and <r which characterize the ellipsoid model. 3 Molecules 27 and the principal radii <7i, <T 2 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,) = ^ V I ' / T ( + T ^3/ [i (BB/4A-C):pp] V ' I + ? ( 3 L 6 ) 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 must have the 2 units of energy per unit volume. If we consider identical ellipsoids and choose 2V2i(-f/ (2|b|) / 7 = 1 4 4 , (3.1.7) then the interaction energy, u(p,q -,qj), is given by 8 u(p, q,-, q,-) = 4 e / | b , " + b j , | e y/\Ti + *[i+(BB/4A-C) C T ] - ( 3 L 8 ) Tj\ For this particular choice of 7 the energy is 4ee when the ellipsoids are aligned and the s intermolecular distance is zero. This is true regardless of the size or aspect ratios of the ellipsoids. T h e 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 «- (p,q,-,q,-) = _ t V l ' ^ / ^/|T,- + T | 4 £ / b + b e S 2 [i (BB/4A-C):pp] ; + ( 3 L 9 ) i such that the total pair potential is given by u (p, q,., q.) = 4 7' + |b ' Ti l b j l + T,-| (V[i+(BB/4A-C): p] _ P v e S /2[i (BB/4^-C):pp]\ _ + 1 { 3 A l Q ) Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 29 This attractive term is simply a negative Gaussian with a m a x i m u m of e ' . T h e factor 3 2 of four scales the well depth so that the m i n i m u m for two spherical Gaussians is —e, independent of the scale factor, and comparable with the well depth of a Lennard-Jones potential. O f 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. T h e 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. T h e 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 system 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 T * = ksT/s, temperature (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* = Pa /e 3 and the reduced t i m e i * = t^Je/o m, 2 where m is the mass of a particle. Suitable moments of inertia about the molecular axis (x', y', z') (see F i g . 3.1) ), can be calculated by assuming a mass density of the form £ e s / 2 s b : r r ? I 'a', a 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 y'y' ai(a\ + erf) , (3.2.1a) ai(a + al) (3.2.1b) 2 (3.2.1c) where Oi = C 2 \ / 2 V C r C r C r ( 7 r / s ) / 2 In M D simulations it is advantageous 1 2 3 3 / 2 . (3.2.2) to have similar relaxation rates for translations and reorientations. T h a t is, the decay of velocity-velocity and angular-velocity-angularvelocity autocorrelation functions should be similar [14]. T h i s can be accomplished by adjusting (. Also, as usual, it is convenient to define reduced moments of inertia of the form I*, , = a jma . Of course, the equilibrium properties of the classical fluid will 2 Iii a a 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. T h e 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, Ql 0 and Ql , defined as [40,41] 2 (3.2.3a) Q l 2 2 = ( - ( 1 + cos 0) cos 2(f) cos 2ip — cos 9 sin 2(f) sin 2ip) 2 (3.2.3b) 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. QQ is 0 the usual uniaxial order parameter and (except for finite size effects) will be zero in the Chapter 3. A Simple Model Potential isotropic phase. Q for Liquid Crystal Molecules 31 is a biaxial order parameter which will be zero in isotropic and 22 uniaxial phases but non-zero in a biaxial phase. Following Allen [40], we calculate the order parameters by first forming the ordering matrices Q , XX Q , YY and Q Qlp 2 Z with elements defined by ^ E ^ V ^ - M . = (and similar expressions for Q YY 0 = 1,2,3, (3.2.4) and Q p) where X;, y,- and z; are unit vectors along the z a three symmetry axes of particle i. We next determine the largest eigenvalue associated with each of the three matrices. T h e 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. T h e second largest and smallest eigenvalues are used to identify the molecular y and x axes, respectively. T h e associated are then orthogonalized together with Z in order construct the Y and X eigenvectors axis of the laboratory fixed frame. In terms of these axes, the order parameters have the cartesian forms Q 2 Q 2 22 = -(x • Q XX 00 •x + Y •q yy = (Z • C T • Z) , •Y - x • q yy (3.2.5a) •x - Y •q xx • 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 . T h e effects of molecular flatness on liquid crystal phase transitions have been studied by carrying out M D simulations of selected repulsive Gaussian ellipsoids. T h e length 2<7i was varied from 0.5 to 3 while the remaining dimensions were held fixed. T h e a2 and o"3 radii were chosen such that the ellipsoid of revolution obtained when <7i = a 2 Chapter 3. A Simple Model Potential 30 for Liquid Crystal 32 Molecules 1 111111 1 11111111111111. ' ' ' / / ,1 / / / / / / / A C\2 — / - // // // J=L v 10 / / / s / J — 1 1 0 ~jr y 1\ 1 i 1 1 1 1 1 0 2 4 6 8 1 1 1 l l l 1 l 1 1 1 1 . . . 1 1 11 10 12 t* Figure 3.2: ( | r ( f ) - r ( 0 ) | ) for 2o 2 Y 2cri : 2<7 : 2a 2 3 : 2a 2 = 1.5 : 3 : 9 (dotted line). nematic phases, respectively. : 2a 3 = 2.75 : 3 : 9 (dash line) and These systems are isotropic and uniaxial Chapter 3. A Simple Model Potential has 2<7i : 2cr : 2<r = 3 : 3 : 9 . 2 3 for Liquid Crystal Molecules 33 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. T h u s , 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 value of p* used would correspond to p/p cp o 2 and a , the = 0.8/^/2 — 0.5657, where p cp is the close l5 3 packing density. Most calculations were carried out with 216 particles in a non-cubic simulation cell starting from a body-centered orthorhombic lattice. ranging from 1 : 3 : 9 to 2 : 3 : 9 , For systems with 2&i : 2cr : 2<r 2 3 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 2a x < 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 F i g . 3.3, the order parameters QQ and Q\ 0 thickness 2o\. 2 are plotted as a function of the particle Again, it should be noted that for the points shown in this plot, the lengths 2(j and 2cr are held fixed (at 3 and 9, respectively) as are the reduced temperature and 2 3 Chapter 3. density. A Simple Model Potential for Liquid Crystal Molecules 34 Thus F i g . 3.3 displays only the influence of changing width; as 2 c i 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 F i g . 3.1) along the director. This clearly shows that molecular effect even upon the formation of a uniaxial phase. biaxiality aligned can have a very strong It is evident that as the particles become flatter the tendency to form a liquid crystal phase is significantly increased. At ~ 1.2, Q\ 2<Ti increases sharply indicating the formation of a biaxial phase. 2 2<7i ~ 1.1, Also, at 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. respect to F i g . 3.3, it should be noted that Q\\ 2 With 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 1 11 - ^1 0.8 — 1 T 1 II * for Liquid 111 S T Crystal 1 1 1 1 1 M * * at QfJ0> 0.6 - 1i i i | i i i i j_ Q22 " X 0.4 - r - i 0.2 h — 1 * 0 1111111111111111111 0.5 35 Molecules 1 1.5 2 1 2.5 1 1 1 1 1 1 1 T 3 2(7! Figure 3.3: Order parameters as a function of particle thickness at fixed reduced temperature and density. T h e squares are the uniaxial order parameter QQ and the circles 0 are the biaxial order parameter Q\ . 2 T h e error bars represent one standard deviation. Chapter 3. A Simple Model Potential Figure 3.4: 2o x : 2o 2 for Liquid Crystal : 2a Snapshot of an instantaneous 3 = 2.875 : 3 : 9. Molecules 36 configuration for an isotropic phase w i t h Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 37 Figure 3.5: Snapshot of an instantaneous configuration for a uniaxial nematic phase w i t h 2ai : 2a : 2a = 2.125 : 3 : 9. Here, the director is along the z' axis (up and down on the page). 2 3 Chapter 3. A Simple Model Potential for Liquid Crystal Molecules 38 Figure 3.6: Snapshots of instantaneous configurations for a biaxial nematic phase w i t h 2<7i : 2<r : 2<r = 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. 2 3 Chapter 3. A Simple Model Potential for Liquid (2<Ti : 2<72 : 2<7 = 1 : 3 : 9 ) 3 (2ci : 2a : 2<7 = 3 : 3 : 9 ) 2 3 Crystal Molecules 39 uniaxial, (2<7i : 2<7 : 2<r = 2 : 3 : 9 ) 3 2 phases. and isotropic 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 F i g . 3.6. However, in the uniaxial nematic, there is a peak i n 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] <P (x' (0)-x' (r))>, 2 where P (x) 2 1 (3.2.6) 2 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 P 2 shows < P (x[(0) 2 as a function of the intermolecular distance. • x' (r*)) > for biaxial, 2<7i : 2cr 2 2 : 2<T = 1 : 3 : 9 , 3 3 2 3.8 uniaxial nematic 2<7i : 2<r : 2<T = 2 : 3 : 9 , and isotropic 2<7i : 2<r : 2cr = 3 : 3 : 9 systems. 2 Figure 3 W e 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 like) [12], then reorientation occurs through many small uncorrelated steps. of decay also indicates that the particles are not "spinning" about the z (Debye- T h i s type 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. F i g . 3.9 shows < x (t*) • x (0) > and l n ( < x (t*) • x (0) >) for two uniaxial nematic Chapter 3. A Simple Model Potential Figure 3.7: <7° (r*) for biaxial 2o 00 2<7i : 2<r : 2o 2 3 dash) systems. x for Liquid : 2cr 2 • 2cr 3 Crystal 40 Molecules = 1 : 3 : 9 (solid), uniaxial nematic = 2 : 3 : 9 (short dash), and isotropic 2o\ : 2<7 : 2u 2 3 = 3 : 3 : 9 (long Figure 3.8: < P (x' (0) 2 • x' (r*)) > for biaxial 2<7i : 2a 1 nematic 2<7i : 2cr : 2a 2 2 (long dash) systems. 3 2 : 2a 3 = 1 : 3 : 9 (solid), uniaxial = 2 : 3 : 9 (short dash), and isotropic 2o~\ : 2o~ : 2<r = 3 : 3 : 9 2 3 Figure 3.9: < x'(t*) • x'(0) system with 2a\ : 2a 2 > for a system with 2<TI : 2a : 2a 2 3 : 2a 3 = 2:3:9 (solid), and a = 1.5 : 3 : 9 (short dash). T h e 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 : 2cr — 2 : 3 : 9 , 3 and the other is near the uniaxial to biaxial nematic transition with 2o~i : 2o~ : 2er — 1.5 : 3 : 9. If the reorientational correlation is given by 3 2 the simple exponential < x'(t*) • x'(0) >= e x p [ — f / r * ] , where r * is the relaxation time, then for the 2a t : 2cr • 2cr = 2 : 3 : 9 (solid line) and 2o 2 3 (dash line) systems, r * « 3.5 and r * x : 2a 2 : 2cr = 1.5 : 3 : 9 3 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<r but this l5 structure decays quickly. T h e 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 i n 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 conditions, the transition will happen at a specific density. A t 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. A t 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 corresponding 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. T h i s 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. A s 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. the particles about their corresponding z T h e orientation of 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. T h e results of these simulations indicate that particles with aspect ratios around 2<7i : 2a 2 3 : 9 have sufficient molecular biaxiality to form a stable nematic phase. : 2a 3 = 2 : Chapter 3. 3.3 A Simple Model Potential for Liquid Crystal Molecules 45 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 temperature 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 i n a uniaxial nematic phase where the director is associated with a different molecular axis. T h a t 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 i n 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. A s 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. T h e 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. T h e 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]. T h e existence of the ferroelectric nematic phase was believed to be due to these specific short-ranged tetragonal I correlations [9,10]. T h e translational mobility of the fluid allowed specific short-ranged correlations to build up, and the system spontaneously 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. T h e y 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). T h e 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 However, experimental systems of near spherical Fe 0\ 3 magnetic solvent [46-48] do not exhibit ordered phases. interaction is magnetostatic. 47 particles in a frozen non- In these systems the dipolar T h e particles are translationally frozen in random locations, but above a certain temperature they can rotate freely. A s the system is cooled, the particles do not orientationally order, but rotationally freeze in random directions. T h e random spatial structure is believed to be responsible for the disordered low temperature phase [48]. T h e random environment about a particle creates a strong random field, and the direction of the field is different for different particles. A t low temperatures, the moments of the particles will align with the field, resulting in an orientationally disordered phase. T h e experimental results seem to contradict Zhang and W i d o m ' s theory, although it is possible that the experimental dipole density was too low, or impurities in the sample were present. T h e role of short-ranged spatial correlations on ferroelectric phase formation is still not well understood. T h e 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. T h e second section concerns ferroelectric ordering in amorphous frozen dipolar systems. T h e third section introduces a new simulation technique devised to study orientational order in spatially random media. T h e 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 4.1 Phases of Amorphous Dipolar 48 Systems M o d e l a n d Simulation Details. We consider dipolar soft spheres defined by the pair potential u(12) = u (12) +u (l2), gB (4.1.1) DD where u (12) is the soft sphere potential given i n E q . (2.2.5) and ss u (12) DD = - 3 ( > i • r ) ( / i • r ) / r + m • fi /r (4.1.2) 3 5 2 2 is the dipole-dipole interaction as previously given by E q . (2.2.7). Under periodic boundary conditions, the dipole-dipole pair potential is replaced by the effective given by E q . (2.3.7). interaction 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. T h e Ising system consists of dipoles in a 3-dimensional space that can either point in the -\-z or —z direction. T h e X Y model has dipoles i n a 3-dimensional space that are constrained to rotate i n 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). T h e existence of a ferroelectric phase can be detected by calculating the average polarization per particle < P > defined as <P>=1//V(|f> d|), r (4.1.3) Chapter 4. Ferroelectric Phases of Amorphous Dipolar 49 Systems 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 i n the system. In ferroelectric systems, < P > is nonzero , while in disordered phases < P > is near zero. There is a possibility that the system could "orientationally freeze" at low temperatures. Orientationally frozen, or "spin glass" phases can be detected by calculating the root-mean-square spin length < S rms < s rms >= > [49] given by iv- E(< •»-•>•< < >)] > 1 m (- ) 12 / 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 < S rms < S rms > non-zero will be ferroelectric. Systems with < P 0 but > non-zero will not be polarized, but the particles are orientationally frozen. W h e n both < P > and < S rms > are near zero, the system is rotationally "fluid-like". Dipolar soft-sphere systems can be characterized by the reduced density, p* = the reduced temperature, T* = k T/e, B and the reduced dipole moment, p* = Na /V, s (p /ecr ) / . 2 3 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 Results for R a n d o m l y F r o z e n Systems We have performed simulations where the amorphous spatial structure was obtained from a "random" configuration of soft spheres. T h e 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. A t 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]. T h e 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. T h e configuration of soft spheres determines the positions for the point dipoles. Initially, 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. T h e polarization values that are plotted were obtained at the lowest temperature where equilibrium could be achieved, denoted as T ^ . T h e values of T ^ for the Ising, X Y and X Y Z models are 10.0, 4.0 and n n 3.5 respectively. Below these temperatures, simulations that were started from perfectly aligned and random states did not converge to the same result. could possibly be obtained with longer simulations. Lower values of T ^ i n However, the temperatures that were attained were within the range that Zhang and W i d o m predicted a ferroelectric phase. T h e Ising model at T ^ i n = 10.0 is almost completely polarized and shows little or no system size dependence. In F i g . 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 i n 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 F i g . 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 i—i—i—r Dipolar 51 Systems i—i—r i—i—r i—r i J J 1.5 tad 1 0.5 0 Figure 4.1: T* = 10.5. j 0 i i j L L I I L L 1 T h e radial distribution function g(r*) for soft spheres at p* = 0.8 and Chapter 4. Ferroelectric 1.2 fl 1 I Phases of Amorphous I I I I I I Dipolar I I 52 Systems I I I I I r I I i i i i i i i i i i i 0.8 • r-t . s "X 0.6 P-i V 0.4 0.2 0 L- 0 i I i i I i 0.2 i i 0.4 0.6 100/N 0.8 1 Figure 4.2: T h e polarization < P > at T^ vs 100/N (circles), X Y (squares) and X Y Z (triangles). Results are included for 108 ( X Y Z only), n 256, 500 and 864 particles. for the randomly frozen Ising Chapter 4. at T ^ Ferroelectric Phases of Amorphous Dipolar 53 Systems = 3.5 show significant polarization for the N — 108 and TV = 256 systems. i n However, the polarization decreases monotonously with increasing system size. Indeed, it appears to approach zero i n the thermodynamic limit. T h e 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. T h e 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 ^ seen i n F i g . 4.3, where the root-mean-square spin length < S rms reduced temperature T*. T h e non-zero value of < S rms n = 3.5 can be > is plotted versus the > at T^ suggests that the n system is orientationally frozen at this reduced temperature. Furthermore, < S rms ^min s n o w s little or no system size dependence (eg. for N = 256, < for N = 864, < S rms S r m s > at > ~ 0.37 and 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 calculations 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. T h e Ising model has a ferroelectric phase, however spontaneous polarization was observed at T* w 25, which is much lower than the transition temperature predicted by Zhang and W i d o m . 4.3 Decoupled M D 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 Figure 4.3: T h e root-mean-square spin length, < S model. rms Systems 54 >, for the N = 256 frozen X Y Z Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 55 this would result in very long simulations. A s 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. T h e 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. T h e structure of the substrate is not affected by the dipolar interactions. T h e "equilibrium" state of the dipoles, however, will depend on the underlying motion of the substrate. This model is similar in spirit to those used i n recent studies of non-equilibrium phase transitions i n 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 shortranged 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. T h e implementation is straightforward. T h e force between two particles is f(12) V u (12), r ss and the torque between two particles comes from the dipole-dipole interaction. T h e spatial structure of the soft sphere fluid is then only determined by the soft sphere part of the pair potential. T h e 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 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. A t 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. T h e 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'a ) 8t 2 1/2 = 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. T h e results of F i g . 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. T h e dynamically decoupled X Y model behaves much like the X Y Z system. In F i g . 4.5, we see that the three different masses, m * = 1, 5, 10 spontaneously polarize with transition temperatures at T* « 6, 2, and 1.8. Again, the transition temperature decreases with increasing mass. T h e 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. T h e substrate evolved as i n 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. T h e behavior of the dynamically decoupled Ising model is different from both the X Y and X Y Z models. F i g . 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 C v = Cv/NkB, obtained by numerically differentiating the average dipolar energy with respect to the rotational temperature are also shown in F i g . 4.6. T h e randomly frozen and m* = 5 results are very similar and indicate a phase transition at T* « 25. T h e dependence of the ferroelectric transition temperature, denoted as T , F on particle mass m* for the X Y and X Y Z models is shown in F i g . 4.7. T h e 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. A s the mass increases the transition temperature drops for both the X Y and X Y Z models. A s 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. 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. Chapter 4. Figure 4.5: Ferroelectric Phases of Amorphous Dipolar Systems 59 < 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. Chapter 4. Ferroelectric Phases of Amorphous Dipolar 60 Systems "i—i—i—i—I—i—i—i—i—I—i—i—i—i—I—r 1 - 0.8 A V 0.8 0.6 \0.4 I III * 0.6 E- u t- h 0 4 I IIII| ^ ^ | II M 0.2 0.2 - 0 1111 m 11111111 0 10 20 30 40 r 0 i o I II i i i i 10 i i i i i 20 i J I I I L 30 Figure 4.6: < P > vs T* (rotational) for the Ising model. Results are shown for dynamically random systems w i t h m* = 1 (squares) and m* = 5 (triangles) and for the randomly frozen case (circles). The heat capacities per particle, C /N, are plotted vs T* (rotational) i n the inset. v 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 temperature 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. T h e results for the randomly frozen and dynamically decoupled systems can be interpreted 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. T h e reaction field favors ferroelectric order. It arises from the long-range nature of the dipolar interactions and is independent of the spatial correlations. T h e 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. T h i s 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 . T h e 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. T h e reaction field dominates in the Ising system and gives rise to a ferroelectric state. T h e behavior of the dynamically disordered systems can also be understood. If the Chapter 4. Ferroelectric i i Phases of Amorphous i 62 Systems i i i i i i i i i i i i i i i i i i P I I f | I I I Hj I I I I | I I I LJ 10 2.5 8 - 2 ? 1.5 O 6 - E- I" / =- > 1 -i i i 0.5 1 E-i Dipolar Ii i i i Ii i i i Ii i 1.5 2 2.5 i- 3...--" 4 - 0 i 0 i 0.2 0.4 0.6 0.8 (1/rn) 1 Figure 4.7: T h e mass dependence of the disordered to ferroelectric transition temperature T F (rotational). T h e squares and triangles are results for the X Y and X Y Z models, respectively. T h e values of T-p were obtained from plots of the heat capacity, C /N, v T* (rotational) and a typical example is shown in the inset. estimated uncertainties in the peak position. vs T h e error bars represent Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 63 translational diffusion of the substrate is sufficiently fast compared to the dipolar reorientational time, then a ferroelectric phase is possible. T h e rapid movement of the particles limits the extent that spatially dependent frustrating correlations can build up, allowing R to prevail over E . A s 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 w i t h decreasing mass. A s the substrate particles become sufficiently light, the transition temperature is determined only by the reaction field and hence becomes independent of mass. In fact, m* = 1 gives essentially this l i m i t i n g behavior and reducing the mass further has little 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 i n a spherical cavity interacts only w i t h the reaction field. We w i l l begin w i t h a brief summary of the details of the reaction field. A more thorough discussion can be found i n A p p e n d i x 3. From this simple theory, we can calculate the reaction field contribution to the total energy per particle. T h i s theoretical result w i l l be compared to the reaction field contribution for the same model (Ising, X Y and X Y Z ) from the simulations. T h e results from the simulations gives quantitative information on the role of long-ranged and short-ranged interactions on the stability of ferroelectric phases. T h e reaction field R at a point dipole i n a sphere of radius r , is due to the polarization c of the dielectric beyond the sphere. We w i l l restrict ourselves to the case that the dielectric constant of the surroundings, e = oo. If there are N dipoles i n the sphere, and the total moment of the sphere is M = f>, t=i (4.4.1) Chapter 4. Ferroelectric Phases of Amorphous Dipolar 64 Systems then the reaction field is given by [14,52] R = ^fvl, (4.4.2) '"c where the direction of reaction field is along the total moment of the sample. A point dipole //, w i t h i n the sphere w i l l interact w i t h the reaction field, and its interaction energy is given by [14] u = -^p -K. t (4.4.3) t Clearly, the low energy state is found when the dipole aligns w i t h the reaction field. T h e factor of a half takes into account the work of polarizing the surrounding continuum [51,52]. T h e total moment M can be written as (4.4.4) M = Md, where M = |M|, and d is a unit vector i n the direction of M. If the instantaneous polarization P is given by TV E M«: • (4.4.5) D 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 w i t h i n a sphere of radius r , then c N = (4.4.7) where p is the number density. T h e total moment can be written as M = ^J^Pd. (4.4.8) R = ^Pd, (4.4.9) T h e reaction field is then Chapter 4. and the i Ferroelectric Phases of Amorphous Dipolar Systems 65 particle's instantaneous energy is u = % -^-fii-d. (4.4.10) For the N dipoles within the sphere, the total energy is U — J2iLi ii u 2ir then N U = -- p P\Z(ir<l- (4-4.11) 2 P i=i 6 \i p, • d is written as cos (9, then the average value of P for the X Y Z model in the i t canonical ensemble can be written as < P >_ K P(~^ ) JQ exp(—Bu\) sin 9id9i C O S0 1 e x U l s i n 9 l d d (4 4 12) l W i t h the mean field approximation, we replace the instantaneous polarization P with the average polarization < P > resulting i n = < > _ ^cos-^expfl^/V < P > cosfli)sinMfli J ^ e x p ^ / V < P > cos0i)sinM0i ( A 4 ^ ' This can easily be evaluated, to obtain < P >= ^^I'-'V'. (4-4.14) a^e — e ) 37 x where x = 27T/V < P > (4.4.15) 3k T B For the X Y model, the average polarization < P > is n ffcOsMxptf p/3p 2 %exp(^p8p <P> COS 9 )de < P > cos 8 ^ 2 1 i ' 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 i th particle can be written as u i + 27T = -—pp 2 <P> (4.4.17a) Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 66 or 2TT = + ^-pp 2 (4.4.17b) <P>, depending on whether the dipole is directed along ( —) or against (+) R. For a system of TV particles, the configurational integral is N £exp(-/3[/) = n(exp(-/3«,-+) + exp(-/3«,-_)), j=i (4.4.18) t=i where N */ = £ > , - . (4.4.19) t=i T h e number of possible states for an Ising system of iV particles is easily shown to be 2. N T h e average polarization can be written as < f > 4 ^ f f . ^ ; ) ^ » , N ( 4 ,. 2 0 ) E , = i exp(-/9t/) which is written in the form exp(-/3u+)+exp(-/3u_)' 1 J 2ir = ^fPV 2 <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. T h e 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. T h e driving force for this transition is simply the reaction field due to the polarization of the surrounding dielectric. Chapter 4. Figure 4.8: Ferroelectric Phases of Amorphous Dipolar Systems 67 < P > vs T* from the reaction field, as found from a simple mean field theory. T h e short dash, solid and long dash lines are for the X Y Z , X Y , and Ising models, respectively. T h e circles are the frozen TV = 256 Ising results. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 68 T h e 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 <U >/N -2TT = RF T h e heat capacity per particle Cv/N -pp (4.4.22) <P>' 1 can be found from d<U > IN RF Cv/N dT 47T PP <P> 2 d<P> (4.4.23) dT For the X Y Z model we have, d< P > <P> 1-X T dT (4.4.24) 3k T :i-x) B 2TTPH 2 w here x (e 2 X x + e~ ) - ((e x 2 x (e 2 - e -x\2 x — x (4.4.25) e~ ) x 2 For the Ising case, a similar expression is found but with X given by (e - (e + x X x ^-x\2 (4.4.26) e~ ) x 2 T h e 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 < U RF > /N as a function of reduced temperature for the Ising and X Y Z models, respectively. T h e 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 C v = Cv/N kg is continuous across the transition, indicating that the transition is not first order. O f 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 69 Systems Figure 4.9: T h e reaction field contribution to the energy per particle < UR F > /N for the Ising model, as calculated from simple mean field theory. T h e inset is the reduced heat capacity per particle Cy/N. Chapter 4. Ferroelectric Phases of Amorphous Dipolar 70 Systems Figure 4.10: T h e reaction field contribution to the energy < U RF > /N for the X Y Z model, as calculated from simple mean field theory. T h e inset is the reduced heat capacity per particle C v /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. T h e reaction field con- tribution to the energy can easily be calculated in computer simulations by using E q . (4.4.22). A s in the previous section, where we "conveniently" separated the local field F J ' / into E / oca oco ( = R + E , we can separate the total energy per particle into a reaction field contribution < URF > and the rest, < UST >•> < U >=< U R F such that > + < 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, < U RF > /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, < U RF > jN in disordered phases. dominates i n ferroelectric phases, and < U ST In fact, i n ferroelectric phases, < U ST > /N > /N dominates 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. E v e n 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. A s the mass is increased, the structural contribution dominates at lower and lower temperatures. A s the specific details Chapter 4. Ferroelectric Phases of Amorphous Dipolar 72 Systems 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. T h e solid line corre- sponds to < U* > /JV, the short dash is < UR > /N and the long dash is < U* F T h e squares, triangles and circles are for m * = 1, 5 and 10, respectively. ST > jN. Chapter 4. Ferroelectric Phases of Amorphous Dipolar Systems 73 Figure 4.12: T h e short and long-range (reaction field) contributions to the energy per particle < U* > /N for the dynamically decoupled X Y model. T h e solid line corresponds to < U* > /N, the short dash is < U RF > /N and the long dash is < U* > /N. T h e squares, triangles and circles are for m* = 1, 5 and 10, respectively. ST Chapter 4. Ferroelectric Phases of Amorphous Dipolar 74 Systems Figure 4.13: T h e short and long-range (reaction field) contributions to the energy per particle < U* > /N for the dynamically decoupled Ising model. sponds to < U* > /TV, the short dash is < U RF T h e solid line corre- > /N and the long dash is < Ugr > T h e squares, and triangles are for m * = 1 and m* = 5, respectively. /N. 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. Summary 4.5 F r o m 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. W e 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. F r o m 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 reaction 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 D i p o l a r D i s c o t i c P a r t i c l e s Liquid crystal phases of disc-shaped or discotic molecules have been observed experimentally [53,54] and with computer simulation [23,42]. T h e 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]. T h e shape of the particle was a cut-sphere, as described i n 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~(D k T) 3 B u/^J~(D k T) 3 B < 0.125 exhibited unpolarized columns. However, systems with ~ 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. T h e 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]. A t 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 76 m 0.1 Chapter 5. Ferroelectric do not. Phases of Dipolar Discotic Particles 77 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. T h e 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. T h e 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 A s a preliminary investigation, we have constructed different lattices of dipolar particles. T h e lattices were chosen to coincide with the observed local structure in the ferroelectric nematic and the antiferroelectric columnar phases. T h e 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]. T h e antiferroelectric columnar phase has columns arranged with hexagonal symmetry. T h e 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 B o t h 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) where / is the length of the cell. = (n/,m/), T h e 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]). T h e column in the center is shifted by cr/2 i n the z direction. T h e arrangement is shown in F i g . 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 orienting 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 interdigitation, of columns. A sketch of the lattice is shown i n F i g . 5.2. In these calculations, lattices were "prepared" at the same reduced density, p* = pa /\/2 3 = I/A/2, and dipole moment p* = (p /ea ) 2 3 1/2 = 3. T h e 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 headto-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. T h e circles are "top" views of columns of particles. T h e 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. T h e 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 m i n d 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. reference, a dipole with LI* = energy of Upp/N « —21.6 For 3 in an infinitely long column of dipoles will have an [44]. For the spherical case, the tetragonal I lattice has the lowest energy per particle with U£, /N « —24.0. T h e interdigitation between columns D 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. Up /N D For spherical particles, the energy of the antiferroelectric columnar lattice is s=s —22.4. T h e 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. A t 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) /N D 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. T h e antiferroelectric behavior observed in fluids of dipolar cut-spheres [56,57] can be Chapter 5. Ferroelectric _ 1 1 Phases of Dipolar 1 11 Discotic 82 Particles 1 1 .i...L.i..LJ 1 1 1 -22 Q * Q — t - — i 23 - i # — * * " i i i 1i i i I n , 1 1.2 1.4 i i 1.6 i I I I ~ 1.8 b/a Figure 5.3: Dipolar energy per particle UE> /N D for different lattices of oblate dipolar particles. Triangles are for an antiferroelectric columnar lattice, and squares are a ferroelectric tetragonal I lattice. Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 83 explained simply in terms of entropy. T h e 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. F l u i d systems with these intermediate aspect ratios could have characteristics of both lattices. 5.2 F e r r o e l e c t r i c Phases of S i m p l e D i p o l a r O b l a t e P a r t i c l e s In this section we present a partial phase diagram for dipolar oblate particles. In particular, 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 i m u l a t i o n Details Dipolar oblate ellipsoids can be easily characterized by choosing the fundamental unit of length a to be the height of the particle a. strength p* = (p /ea ) / , 2 3 1 2 We then can define the reduced dipolar and reduced density p* = pab . 2 T h e pair potential for this system can be written as u{\2) where u rep = u rep (12) + u D C (12), (5.2.1) ( 1 2 ) 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). simulations UDD(12) is replaced by the effective pair potential, u^p(12; O f course, in e'), given by E q . (2.3.7) and the parameters for the Ewald sum are as given in Chapter 2. T h e reduced Chapter 5. Ferroelectric timestep St* = Sterna je 2 Phases of Dipolar Discotic Particles 84 = 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 < P 2 > order parameters [14,59]. In models where the dipole is along the symmetry axis of the molecule, the instantaneous value of the P 2 order parameter is calculated from the largest eigenvalue of the ordering matrix, Q, with elements Q AP w here p = 1/N £ 1/2(3^4 - M, (5-2-2) is the a component of the unit vector given by fa = (p) l a ^t -. T h e eigenvector 1 4 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 Pi = 1/N (5.2.3) W h e n the dipole is along the symmetry axis of the particle, P i is identical to the polarization P (Eq. (4.1.3)) in ferroelectric phases. and < P 2 For ferroelectric nematics, both < P i > > will be non-zero. For antiferroelectric columnar phases, < P i > will be near- zero, while < P 2 > will be non-zero . T h e 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 ) lattice. T h e 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. T h e 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 > did not exhibit N dependence. We note 2 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 Results We begin discussing the results with Fig. 5.4, where the < P > order parameter versus 2 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. T h e results include systems with b/a = 1.2, 1.3 and 1.5. A l l three systems have < P > ^ 0 below 2 p* « 0.7 and it begins to increase at p* w 0.8. However, the < P > for the b/a = 1.5 2 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. T h e < P i > order parameter versus p* for the same T * , p* and aspect ratios is shown in F i g . 5.5. Each system has an isotropic to ferroelectric phase transition at p* « 0.8. T h e b/a = 1.5 system shows a sharper transition, with < P i > rising to around 0.85. T h i s 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 F i g . 5.6 where Chapter 5. Ferroelectric 1 I I Phases of Dipolar I I I I I Discotic I i 86 Particles i i r i i i i 0.8 - I A I - DTO.6 v 0.4 0.2 0 0.5 I I::::::::::::::* I I I I I I I I 0.7* 0.6 1 I I I 0.8 I I I 0.9 P Figure 5.4: T* = 1.9 results for the < P > order parameter as a function of reduced 2 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 1 Phases of Dipolar I I I I I I Discotic 87 Particles I I I I I I I I I I I I ' l l I I I I I I I I I 0.8 A pTO.6 v 0.4 0.2 0 0.5 I 0.6 0.7* 0.8 0.9 P 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 Figure 5.6: b/a = 1.45 p* = 0.792 < P 2 Particles 88 > (circles) and polarization < Pi > (triangles) as a function of reduced temperature, T*. Chapter 5. Ferroelectric Phases of Dipolar Discotic 89 Particles a temperature scan is performed on a system with b/a = 1.45, and p* = 0.792. There is a sharp rise i n 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, < P 2 > ~ 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 i n the dipolar cutsphere 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] = g (r) 000 112, s _ 3< E,-,j,,?j 8{TJ g^{r) g (r) 000 < = < ^ g ( F j ~ ' ~ r r ) , > - rt - r)[3Q,- • v^jfa ^ ^ - ^ r ' (5.2.4a) • r,-,-) - fa • fa) > ^ > , (5.2.4c) 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 g (r). T h e function <? (r) 112 110 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||. T h e correlation functions are defined as = " l l v <"(;ii-ni + fr| )> l < n ( r , r + Sr ) > 0 N fl u (5.2.5, ' K Chapter 5. Ferroelectric Phases of Dipolar Discotic ( \ < n(r r u Particles 90 + 8r ) > L L < «o(r±,r + 8r ) ± > ± In the parallel case, (n) is the average number of particles in a disc-shape shell perpendicular to d at r|| of width Sr\\ and radius [(L/2) 2 the simulation cell. — r |] 2 1//2 , where L is the length of 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 both cases (n ) 0 — r ] / . In 2 1 2 is the average number of particles one would expect i n the relevant shell calculated using the angle-averaged pair distribution function. Positional order parallel or perpendicular to the director will appear i n 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 i n 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 In F i g . 5.7, g 0 0 0 , g 110 and g 112 L/2. are shown for a system with p* = 0.792, b/a = 1.45, T * = 1.35. T h e sharp peaks i n all three functions at r* « 2.1 and 3.2 indicates columnlike structures where the dipoles are aligned head-to-tail. T h e parallel and perpendicular correlation functions for the same system are shown in F i g . 5.8. T h e peaks i n g\\(r^) at rj| w 1, 2, 3 and the strong peak in g±(r]_) at r]_ = 0 indicates columnar ordering. Peaks i n g± at r* ~ 1.5 and 3 indicate that the columns are arranged with hexagonal x symmetry. T h e peaks in <7||(rj|) at r*| « 1, 2, 3 arise from head-to-tail particles only, indicating 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 i n 5||(^|j) at rj| « 1, 1.5, 2, 2.5 and 3. T h e peaks i n 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 I 0 Figure 5.7: g (r*) 000 Phases of Dipolar Discotic I I I I III I I I I I I I I 1 Particles II I I I I I I I 2 I 3 (solid), <? (r*) (short dash) and g (r*) 112 b/a = 1.45, T* = 1.35 system. 91 lw I I I I I 4 (long dash) for a p* = 0.792, Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 92 a snapshot of a typical antiferroelectric columnar phase is shown in F i g . 5.9. Here, the columns are well defined, and there is strong hexagonal symmetry of the columns. There is no average polarization. T h e pair correlation functions <7 (r*), g (r*) and g (r*) 110 000 for the ferroelectric phase 112 at T* = 1.9, p* = 0.792 and b/a = 1.45 are shown in Fig. 5.10. T h e 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 g (r*) iw decays to a non-zero constant, indicating that there is long-ranged ferroelectric order. T h e functions and <7x( I)) r a r e shown in F i g . 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. T h e 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, g (r*), 000 g (r*), no 5.12 and <7||(r-|j) and fl'i(r^) are presented in Fig. 5.13. correlations in g (r*), 000 g (r*), 110 and g (r*) 112 and g (r*) 112 are shown in F i g . Weak head-to-tail dipole-dipole are found at r* 1, but the structure is very short-ranged. A s expected in isotropic systems, ^(r-jj) and g±(r]_) are uniform. T h e 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. F i g . 5.14 shows a typical antiferroelectric columnar, (b/a — 1.45, p* = 0.792, T* = 1.35) configuration. T h e 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. phase is antiferroelectric columnar. The Chapter 5. Ferroelectric Figure 5.10: # 000 (r*) Phases of Dipolar Discotic (solid), g {r*) 112 p* = 0.792, b/a = 1.45, T* = 1.9 system. (short dash) Particles 95 and g (r*) lw (long dash) for a Chapter 5. Figure 5.11: system. Ferroelectric Phases of Dipolar Discotic g\\(r*\) (solid) and g±(r* ) ± Particles (dash) for a p* = 0.792, b/a 96 = 1.45, T* = 1.9 Chapter 5. Ferroelectric Phases of Dipolar Discotic I I l 0 Figure 5.12: g (r*) 000 l l l I I I I l I I I I Particles I I I 97 I I I I 1 2 (solid), g (r*), 112 p* = 0.792, b/a = 1.45, T* = 2.5 system. 3 (short dash) 4 and g (r*) 110 (long dash) for a Chapter 5. Ferroelectric Phases of Dipolar Discotic Particles 98 i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—i—r 5-H 1.5 J_ OO oo 1 0.5 j 0 i i i I i system. i i I i i i i I i L 1 Figure 5.13: flr (r[j) (solid) and g±(r* ) N i ± (dash) for a p * = 0.792, b/a = 1.45, T* = 2.5 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 < P 2 This contrasts with the antiferroelectric columnar phase, where < P 2 > is large. > ~ 0. Further evidence of a ferroelectric phase transition for the system with b/a — 1.45 p* = 0.792 is given in F i g . 5.17, where the average energy per particle < U* > /N a function of the reduced temperature T*. < C v > /N, is shown as T h e inset is the heat capacity per particle, 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. T h e first isotropic to ferroelectric phase transition involves a significant change in local structure and energy, and is most likely first order. T h e 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. A s the particle is made more oblate, both order parameters rise. A t b/a « 1.45, just as < P 2 the ferroelectric phase is observed. > is rising, A s 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). A t b/a ~ 1.45, the columnar correlations are strong enough for < P 2 > 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. T h e breakdown of the antiferroelectric columnar phase results i n 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. T h e 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. T h e phase is isotropic. Chapter 5. Ferroelectric 6 rr Phases of Dipolar i i i i i -10 - 10 i i i i i r I L j r •mi 0 i I V i 103 9 - 8—ao L<c>/N A i Particles -1 i i i i i i i i i i i i i i 30 * i Discotic I I -12 - I 1 I 1I 1I 1 1 1 1.5 , 2 I I 2.5 / I T -14 16 J 1 I I I I I L 1.5 Figure 5.17: T h e reduced energy per particle, < U* > jN J 2.5 as a function of reduced temperature, T*. T h e inset is the excess contribution to the reduced heat capacity per particle, < C /N. v > 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 diagram". In F i g . 5.19, a T* versus b/a phase diagram is sketched for p* = 3 and p* = 0.792. T h e 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. T h e 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. A t 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. T h e 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 i n 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. systems with b/a For example, = 2 have the antiferroelectric columnar to ferroelectric transition near Chapter 5. Ferroelectric 1 1 Phases of Dipolar 1 1 1 1 11 1i Discotic | i * 0.8 — i i I I | i * 0.6 - i i , i i | i i i | i 1 ' / P*\ i 105 Particles • - : » * % * A P?0A _ v : 0.2 0 <p > 4 2 h i 1 i 1 i i 1 i i 1 i i i i 1 i i i 1 i i i 1 i 1.2 1.4 1.6 1.8 2 b/a Figure 5.18: T * = 1.9 , / = 0.792, < P >, (circles) and polarization < P i > (triangles) 2 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. T h e reduced density is p* = 0.792. T h e width of the ferroelectric phase is only estimated. T h e 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. T h e "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 i n 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 i n small regions of the phase diagram for specific p, T and b/a. Chapter 6 A More Robust Model for a Ferroelectric Liquid Crystal F r o m 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 antiferroelectric 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. T h e 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. T h e "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. T h i s 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. T h e 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]. A t higher densities this model possesses columnar and other spatially ordered phases. T h e 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 F i g . 6.1. T h e dipoles are directed along the symmetry axis and cover a circular "patch" of diameter D'. If the dipole density in the plane is p , es dipole moment of a particle p is given by p = p A, es where A = TTD' /4 2 then the total is the area of the dipolar patch. For this model the pair potential has the form (6.1.1) where u (12) is the hard cut-sphere (cs) interaction and u (12) is the total electrostatic cs es Chapter 6. A More Robust Model for a Ferroelectric (es) potential. Liquid Crystal 110 T h e total electrostatic pair interaction can be obtained by integration over the dipolar distributions of a pair of particles. Explicitly, one has 3 J (12) = pi J[-3(AI • where the integrations are over the areas A\ vectors are shown in F i g . 6.1. patches, s = s 2 — S i , s = s/s and A 2 Sl 2 (6.1.2) of the dipolar patches and the various 2 a , 3 2 In particular, A i the dipoles of particles 1 and 2, Si and s fa}/s d ds s)(/x • s ) + A i • n d A 2 a r e u m 1 : vectors directed along are position vectors of points within the dipolar and s = | s | . At long-range, u (12) clearly must obey the limiting expression es p, /r «es(12) ~ - 3 ( / i i • r)(fi2 • r ) / r + /* • 5 3 x 2 , (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. T h a t is, the symmetry axes of all cut-spheres must be parallel to each other and the dipoles may point only up or down. the angle between p,\ For this case the pair potential is a function of r and of (or p, ) and r , hence only a two dimensional table is required. 2 Therefore, for computational purposes the total electrostatic contribution to the pair potential is given by •u (12) es = u abie(12) + u ^ ( 1 2 ; e ' ) , (6.1.4) c t where «tabie(12) = u (12)- f-3(/ii-r)(/i -r)/r + e s 2 5 A*i-A*2/r l 3 . (6.1.5) Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 111 Figure 6.1: T h e coordinate system used to evaluate the interaction potential. T h e relevant vectors defined in the text are shown. T h e 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 U DE> (12; e' = oo), which is the effective dipole-dipole interaction under periodic boundary C conditions with a conducting surrounding continuum ( E q . (2.3.6)). T h e short-range electrostatic interactions due to the dipolar patch are contained in u bie(12). ta 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* = pLD = 0.525. 2 system the second-rank order parameter < P 2 Furthermore, for the fully rotational > 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. O f course, for the parallel case < P 2 > = 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]. T h e phase behavior of the model was explored by carrying out Monte Carlo ( M C ) calculations i n the canonical ensemble. Calculations were performed at p* = 0.525 varying the reduced temperature T* = ksTL /p . 3 2 In most calculations 288 particles were used, but several tests with 540 particles gave very similar results. T y p i c a l 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 displacement, ((R(r) — R(0)) ), where R ( r ) is the position vector of a particle and r represents 2 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 F i g . 6.2, it is clear that i n 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)) ) becomes constant as r 2 increases. T h e presence of ferroelectric order was detected by calculating the average polarization per particle < P > defined in E q . (4.1.3). O f 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 i n size from D'/L 9. For values of D'/L was formed. D'/L = 4 to ^ 7, spontaneous polarization occurred and a ferroelectric phase For smaller values of D' ferroelectric transitions were not observed. For = 9, the temperature dependence of < P > and of the reduced constant volume heat capacity C v = Cv/Nks are shown i n F i g . 6.3. T h e 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 DebyesA . - 2 Likewise, for L = 10A, the dipole moment is 37.0 Debyes and p es e s = 0.011 = 0.0058 DebyesA . - 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: T h e mean square displacement, ( ( R ( r ) - R ( 0 ) ) ) , for T* = 0.02 (short dash), 2 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 1 1 1 1 * 1 1 1 1 1 1 1 i 0.8 A i i i i i 111111 i j ii 111111 ii I \> 115 Crystal 2.5 F2 F- 0.6 - 1.5 P-i V Liquid I 0.4 - i ii I i i i I i i i I i i i I i i i"l -\ 0 2 4 6 8 10 T* ^ 2 0.2 0 i 0 i i i i 2 i i I 4 6 10 T* 2 i i i i i i i 8 Figure 6.3: T h e temperature dependence of the polarization < P > and reduced constant volume heat capacity per particle C /N v (inset) for a system with D/L — 10, D'/L and p* = 0.525. T h e error bars represent one estimated standard deviation. =9 Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 116 Figure 6.4: "Snapshot" of a configuration for a system w i t h D/L = 10 and p* = 0.525. The system w i t h D'/L = 9 is shown i n the ferroelectric (T* = 0.025) phase. Figure 6.5: "Snapshot" of a configuration for a system w i t h D/L = 10 and p* = 0.525. The system w i t h D'/L = 9 is shown i n the unpolarized (T* = 0.050) phase. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 118 Figure 6.6: "Snapshot" of a configuration for a system w i t h D/L = 10 and p* = 0.525. T h e system w i t h D'/L = 5 is shown i n the unpolarized (T* = 0.0275) phase. Note that this configuration consists of loose polarized columns w i t h nearly equal numbers up and down. O n average, there is no net polarization i n 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 g (r ) s z [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. gives the probability of finding a particle j with an r g (r ) s z z 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 Sr , and its volume z is K i i d e r = 7rL c>r /4, then g (r ) 2 y 2 n s z can be written as ,.(r, + l / M r . ) = <"('y ' . +fr )> (6.2.1) P ^cylinder Here, n(r , r + Sr ) is the instantaneous number of particles found in the range z z z r ,r +6r , z z z and the average is over the number of particles JV, and the number of M C cycles. Figure 6.7 shows g { * ) i " r s z o r the ferroelectric T* = 0.025, p* = .525 system. T h e 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 F i g . 6.6 that at lower temperatures the D'/L = 5 system is developing strong columnar correlations. A t 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 F i g . 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 i n previous simulations of ferroelectric fluids of dipolar spheres [10]. Chapter 6. A More Robust Model for a Ferroelectric Liquid Crystal 120 Figure 6.7: g (r* ) for T* = 0.025. N o long-range smectic order is detected in this system. s z Chapter 6. A More Robust Model for a Ferroelectric Liquid 121 Crystal 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 L u (12) 3 es /p 2 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 (D'/L ~ 0) and for D'/L = 10. Curves are shown for single point dipoles = 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 m i n i m u m 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. T h e ferroelectric transition for the D' = 9 system can be explained by examining the long and short-range contributions to the energy per particle. A s in Chapter 4, we can conveniently separate the total energy < U > into a reaction field contribution < URF > and the rest < UST > given by < u >=< U RF > + < U ST > • (6.2.2) T h e long-ranged reaction field contribution < URF >, depends on the total dipole moment of a particle's charge distribution. T h e 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 / u ( 1 2 ) / u , as 3 2 e s / one particle is slid over another with their dipoles directed in the same direction and their surfaces i n contact. T h e inset shows the same curves on a magnified scale. Curves are shown for single point dipoles (D'/L 5 (dotted) and 9 (solid). ~ 0) (long dash) and for D'/L = 4 (short dash), Chapter 6. A More Robust Model for a Ferroelectric i i i i ii i i i Liquid 123 Crystal i i r i 0 A -2 * V -4 i i ii • i 4 10 T* i i i I I - I I 6 2 Figure 6.9: T h e contributions to the reduced energy per particle, < U* > jN. The squares are the total energy < U* > /N, the triangles are the reaction field contribution > /N, and the circles are the short-range structural contribution (everything else) < U RF < U* ST > IN. T h e reduced density is p* = 0.525. Chapter 6. A More Robust Model for a Ferroelectric contribution < U RF > /N, Liquid Crystal and the structural contribution < Ug T > /N 124 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. T h e total reduced energy per particle < U* > /N negative as the reduced temperature is decreased. becomes more T h e 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. T h e 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 Summary We have described a model discotic nematic liquid crystal which forms a stable ferroelectric phase. T h e 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 columnar ordering. This avoids the antiferroelectric columnar phase previously observed [56] for cut-spheres with single point dipoles. T h e 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. T h e 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 nondegenerate semi-axes. T h e 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. T h e results offer an answer to the long standing question of whether or not spatially random dipolar systems could have ferroelectric phases. T h e 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 regions. 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. T h e 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. T h e 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 stability 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. T h e 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. W h e n the particles align in the liquid crystal phase, the increase in translational 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. T h e possible existence of ferroelectric phases in spatially disordered dipolar materials 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 materials [46-48]. Systems with one (Ising), two ( X Y ) and three ( X Y Z ) dipolar components Chapter 7. Summary and 127 Conclusions 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. T h e results strongly indicate that the systems are unpolarized in the thermodynamic limit. F r o m dynamically decoupled simulations, it was found that specific short-ranged spatial 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. T h e Ising model has a ferroelectric transition 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 embedded 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. T h e stability of the ferroelectric phase depended on details of particle shape and state parameters. Systems 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 negligible, 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. T h e 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. T h e 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 shortranged 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. Bibliography [1] P. J . Collings, Liquid Crystals (Princeton University Press,New Jersey,1990). [2] M . B o r n , Sitz. Phys. M a t h . 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, M o l . 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. 6 9 , 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 . 4 6 , 7783 (1992). [11] Goldstein, H . Classical Mechanics 2 edition (Addison-Wesley,1980). nd [12] J . P. Hansen, I. R. M c D o n a l d , Theory of Simple Liquids (Academic Press, London, 1986). [13] D . A . M c Q u a r r i e , Statistical Mechanics (Harper & Row, New York, 1976). [14] M . P. Allen and D . J . Tildesely, Computer Simulation of Liquids (Clarendon, O x - 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 Oxford, 1984). 129 Fluids V o l . 1 (Clarendon, Bibliography [17] 130 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, A d v . C h e m . 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, C h e m . Phys. Lett. 105 no. 3, 277 (1984); J . W . Perram and M . S. Wertheim, J . C o m p . Phys. 58 409 (1985). [21] G . J . Zarragoicoechea, D . Levesque, and J . J . Weis, M o l . Phys. 75, [22] D . Frenkel, H . N . W . Lekkerkerker, and A . Stroobants, Nature 332, [23] J . A . C . Veerman and D. Frenkel, Phys. Rev. A 45, 998 (1992). 822 (1988). 5632 (1992). [24] H . Azzouz, J . M . Caillol, D. Levesque, and J . J . Weis, J . C h e m . Phys. 93, 3520 (1990). [25] S.W. de Leeuw, J . W . Perram, and E . R . Smith, A n n . Rev. Phys. C h e m . 37, (1986); Proc. R. Soc. 245 A373, 27 (1980); A388 (1983). [26] D . J . Adams and I. R. M c D o n a l d , M o l . Phys. 3 2 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 . C h e m . Phys. 21, 1087 (1953). [29] H . F . Davis, A . D . Snider Introduction Bacon) to Vector Analysis 4 th edition ( A l l y n and 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, [35] J . G . Gay and B . J . Berne, J . Chem. Phys. 74, 4213 (1972). 3316 (1981). [36] D . J . A d a m s , G . R . Luckhurst, and R . W . Phippen, M o l . Phys. 61, 1575 (1987). Bibliography 131 [37] E . de Miguel and L . F . R u l l , M o l . Phys. 74, 405 (1991). [38] E . de Miguel, L . F . R u l l , M . K . Chalam, and K . E . Gubbins, M o l . Phys. 74, 405 (1991). [39] E . de Miguel, L . F . Rull, M . K . C h a l a m , and K . E . Gubbins, and F . van Swol, M o l . Phys. 72, 593 (1991). [40] M . P. A l l e n , Liquid Crystals 8, 499 (1990). [41] J . P. Straley, Phys. Rev. A . 10, 1881 (1974). [42] D . Frenkel, M o l . Phys. 60, 1 (1987). [43] J . Vieillard-Baron, J . C h e m . 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 . K h a n , 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 G u o , D . H . Ryan, M . J . Zuckerman, and M a r t i n 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 2 ND ed. (Elsevier Scientific P u b - 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, M o l . Phys. 75, 989 (1992). [56] G . J . Zarragoicoechea, D . Levesque, and J . J . Weis M o l . 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, M o l . Phys. 73, 439 (1991). [60] D . J . Evans, M o l . Phys. 34, 317 (1977). [61] M . P. Allen, M o l . Phys. 52, 717 (1984). Appendix A The The 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 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. T o begin the derivation we must first define the gradient operator, W)/ = Jiml(^(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 /. T h e negative charge, for convenience, is at the origin. T h e situation is shown in F i g . A . l . T h e potential at some point r , due to q is <^+(r ) = q/r+, + r + + — | r | . Likewise, the potential at r + to the negative charge at the origin is <^-(r_) = — g/r_, where r The + due + + 1 = r_ and r_ = |r_|. potential at r, where r = r _ , due to these point charges is <Kr) = 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 , + #r) = - ( M O " M r - ) ) . As ( + ( A - A - 3 4 ) ) / becomes small, the definition of the gradient operator can be employed to give #r) = - l - V # r ) , 133 (A.5) Appendix A. The Pair Potential for Point 134 Dipoles 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, u (l2) = q(-<f>(r) + <^(r')), DD (A.6) where r' = r + l . (A.7) Again we can apply the gradient operator to obtain the well known result, u (12) DD = <?1- V ^ r ) , (A.8) and if the limit is taken, / —> 0, q —> oo and q\ = p, we can write u (l2) DD = A*-V^(r) Mi'A*2 3(/xi • r)(f* • r) 2 (A.9) endix A. The Pair Potential for Point Dipoles Figure A . l : T h e geometry that is used for the point dipole pair potential. Appendix 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 surrounded by an "infinite" sphere of periodic images. point dipoles. T h e cell is 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. E a c h spherical shell is composed of cubic cells arranged to form a "jagged" sphere. A s 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* = limiEEE'(^ • V ) ( s ^ u z ,-=i j=i n M i • V)r^-r, l T l r where r = rj — r and n are reduced by the boxlength, L. 4 (B.l) n T h e 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|. T h e 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 136 r + n , (B.2) Appendix B. The Ewald Sum for Dipolar Systems 137 where |r + n | is the reduced lattice distance. T h e dipole strength, LI is not reduced by the box length. T h e sum i n 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). T h e energy is now written as U(s), and USH c U SH a n be found from (B.3) = HmU(s). s—>-0 For dipolar systems, the explicit form of U(s) is given by = lim i £ £ J > , • V ) ( ^ . • V ) j — ^ - -U(s)L» '=1 i=i n' where e~ sn 1 (B.4) 1 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 i n [25]. Following Deleeuw, we consider the substitution 1 1 (B.5) °° x- e-\ f dx. f 1/2 T+n x ypK Jo |r + n | T h e energy can then be written as s^O Z^TT i -i + = 1 j = N , n 1 JO N H ^ E E E ( * • )(^' • ) V s^O Zy/TT i = 1 j = 1 n V L -oo e_Sn2 , *- e1/2 *dx, lr+n{2 (B.6) Ja 2 where the integral has been split by introducing the parameter a . T h e sum over n on 2 the first term can be extended to include n = 0 as r —> 0 if this additional contribution is subtracted off. -U(s)L* = l i m -L= £ £ J > , - • V ) f o • V ) e ' 0 t=l j = l n ' 2 / a -W -W*dx X e Appendix B. The Ewald Sum for Dipolar u n V ^ j=i j = i Z V 138 Systems 0 J —1 —1 I J T h e last term in E q . (B.7) can be easily evaluated, giving l i m {p •V) f • V)(p t 3 * x-V'e-W'dx = 4 3 Jo T-*0,]-H (B.8) - 4^, so that the energy expression is 1 -U{s)L* = N N lim — = ^ + E(tt •V ) ( ^ . •V)c-» £ h ^ E E E ( * •v)(^,• s u 1 ^ V ^ i=ij=i ^ , V « v ) n 71 f 7 / 2 x - 2 x-^e-l'+nl^x / ^ - 1 2 - ^ ! 2 - ^ 0 3 + Now consider the following identity, e -|r+n| y 2 = -ir n /y+2*in-r^ ^3/2 n e V 2 (B.10) 2 n and the manipulation —sn — |r + n | x 2 2 = —sn — r x — 2n • rx — n x 2 2 = —(s + x)n = —(5 2 — 2xn • r — x + x 5 2 2 + x ) n — 2xn • r — x r sxr S+ X s+ X . . , 2xn • r x r , - ( s + x)(n H H ^) s+ x ( + a;) sxr X - ( s + x ) 'n H r 5 + x s + x 2 = sxr 2 2 V n 2; 5 2 = 2 s+ x 2 (B.ll) Appendix B. The Ewald Sum for Dipolar Systems 139 We can now make the substitution to give the fourier space term, N 1 -U(s)L = 3 N /-co lim^EEE(MrV)(MrVr 7 n V i=l j=l n ' ^ V 2 | r + n | 2 ^ 01 + - f — 27rm-r ^ - r 3 1 + 4u a w 2 3 (B.12) ^ E ^ ' iv/vr . 3 = 1 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 _ L( 2y + ...) 2 S O T sur + s Jo = I du 2 feu ' = -^ + ^ r 5 3 - r f-u / ^ 1 2 2 1 / 2 3 2 + 0(3)) S + a 2 - ^ + ^r 2 or 3/2 3 + o(s). (B.i3) a We can now write the energy expression as N 1 -t/(,)I 3 = N J V r V V + J V lim-i=EEE(MrV)(M V)e-7 8 =ij=i i—l 7 n — ' /-co Q 2 ^ V ^ d s Q 1 ^ S ^ ' W h e n the gradient operator acts on 7(1 + ^ ) ' (B 14) _ 1 / ' ) the result is zero, since it is independent 2 of r. T h e limit as s —> 0 is also zero. For the term ^-(1 + ^ ) ~ / , the gradient operators 3 2 Appendix B. The Ewald Sum for Dipolar 140 Systems gives (M,-V)(, V)^(l V + -ij)- 3/2 = i)-"V< • |(1 + f>,- (B.15) T h e limit can then be taken, resulting in N V N r. 2 i = l 7 = 1 1 * 4TT N 4EEf^-^ TV 1 -c7(,)Z (- ) B 16 »=ij=i T h i s gives, 3 • J V N J /•CO V lim-^EEE(MrV)(M V)e-7 = r V i=i j = i n ' 2 0 1 +x dx s 1 " Z + " 4* i=l j = l — r f 0 1 4 ^ V , (B.17) T h e limit s —> 0 can now be taken. We note two results relating to the first, and second integrals respectively, i r . - v f x- e~ 2 2 dx erfc(aa) ^ ^ , = -r^e- v2n2/x (B.i8) v 2 n 2 / a (B.19) \ T h e expression for the total energy of the central cell can be written as - w 3 = iEEEh'VK.rV)"'* ' " lr »=lj=l ' n 1 _£ N 1 1 + nl 1 1 ^ A47T N 1^ , 4 ^ 3 ^ Q 3 ' ' (B.20) Appendix B. The Ewald Sum for Dipolar Systems 141 where the pair potential u(12) is u() 12 = w w . . v- 1 ^ e r f c ( a | £ + n|) E-(^- )(M2-V)V 7 n 2irin-r/L —K n /a i + 2 p r 7r E - ( ^ i - v )' (/ V/f ,~ - v•)/ 4vr + 3 2 n^O £ 3 4u a - — 2 u , a 7 V2 ^ 1 3 N L 2 2 2 ^ 2 e n ,^ 3 . . n (B.21) \ 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) = -( * -V)(^ -V)*(r) A 1 2 47r + 4u a 2 ^ erfcfal-+ nh n' / T 3L5"'^-Sb' *w =jE r;- , " + 1 3 IL + n K l 1 ^ —£—• ^in-r/L n#0 , (B - K n 2 n 2 / a 2 ' . 22) <-> B2 3 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). T h e term ^JA^I ' A*2 i important and arises because we are considering a macroscopic s 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. Appendix C The Reaction Field In this Appendix we consider the situation where a large, macroscopic sphere of dipolar material is placed within a continuum with dielectric constant e . T h e radius of the sphere is a, the instantaneous moment of the sphere is M , and the moment is s located at the center of the sphere. T h e 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] R _ 2( '-l)M, - (27TT)^- e ( 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 = sM, where M is the moment of one of the periodic cells. Also, the volume of the sphere can be written as V =fa = 3 sL , 3 (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. T h a t 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 . T h e work required to increase M ( A ) to M ( A + dX) against the reaction field R ( A ) is = dU(d\) and — A M • RdX, (CA) the total work for the charging process A = 0 —• 1 is the reaction field contribution, URF, given by U = RF — f A M • RdA 1 Jo = -1M-R. (C.5) Here R F refers to the reaction field. We see that although the polarization of the dielectric results i n 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 The = - I M . ^ ^ - ^ - M . R F (2e + 1) 3L 2 R t 3 (C.6) ' y work done i n 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 "»< > = -WTT)m*-»- ' 12 The (C 7) total energy of the central cell will be the sum of its configurational energy, and the reaction field contribution, given by ^ (12; ') c where u(12) is given by E q . (2.3.7). e = u(\2) + u(\2) RF (C.8) Appendix Quaternion D Parameters 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 D.l [14,11]. 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*AB , 144 (D.la) Appendix D. Quaternion 145 Parameters 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 — cos <p sin ip — sin (/> cos 9 cos i/> ^ sin <p cos ip + cos <f> cos 9 sin ip — sin <p sin ^ + sin <p sin 0 sin 9 sin ip cos <f> cos 5 cos ip — cos <p sin 0 sin 0 cos ?/> cos 9 j (D.lb) Likewise any vector p i n the body fixed axis frame can be expressed i n 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 i n 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 F i g . 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), and the vector q = q i -f q j x 2 + qk 3 (D.le) such that q = nsin($/2), (D.lf) the rotation can be written as p' = p(?o - 9i - 4 - ql) + 2 c i(q • P) + P * <w. 2 0 (D.ig) Appendix D. Quaternion Parameters Figure D . l : Clockwise rotation of the vector p to p ' about the unit vector n . 146 Appendix D. Quaternion 147 Parameters T h e elements of the rotation matrix can written in terms of q and q by comparing 0 Eqs. D . l c and D . l d and writing the components of p ' in terms of p , go and q , to obtain ^ 9o + ?i -<&-<& 2(qiq 2 y % i ? 2 + 9o9s) - qq) 0 2 ( ^ 3 - qq) 0 ql - q\ + q\ - q\ 3 2(^1^3 + 9092) 2(q q 2 3 ^ + q q) 0 t (D.lh) g - q\ - q\ + q\ j 2(q q -q ) 2 3 2 2 oqi where qo = cos ^6>cos ^(<f) + xp) fD.li] qi = sin ^6 cos ^ (D-lj) q = sin 2 - xp) , i(9 sin ^(<f> - xp) , g = cos^6>sini(<^ + xp) . 3 T h e terms {go, 9 i , 92,93} a r e (D.Ik) (D.ll) known as quaternion parameters [14,29,60] and are con- strained such that 9o + 9 + 9 + 93 = l 2 1 2 2 2 2 (D.2) Briefly, a quaternion, Q , can be written as the combination of a scalar and a vector, [61] Q = (g ,q), (D.3a) 0 q = gii + g j + g k2 3 (D.3b) T h e product of two quaternions can be defined as Q R = (g r 0 0 - q • r, g r + r q + q x r) 0 0 (D.3c) Appendix D. Quaternion 148 Parameters T h e square of the magnitude of a quaternion is | Q | , where 2 |Q| = Q Q , (D.3d) Q = (?o,-q). (D.3e) 2 and Q is the conjugate of Q , given by Appendix E E v a l u a t i o n of the 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 I n t e g r a l T h e overlap integral of two three dimensional Gaussians can be written as f O O /* O O I I -sW,:VV TyAv-v) -p)} J — OO J — OO OO i e / -OO + {r ^ drzdrydrx ( E 1 ) 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 r z yields - arl br,+c) f°° e { + drz ^ ^ / A a - c = ) ? (E2a yja J-co where - = T a - h = 2(T r lxz + T (r x 3xz - = T P x 3ZZ lyz ZZ 3xz x y \fa j ~{a - ( a rV; +S& r 6 Sc W ) ee + Jyz - p) - T p) y y - p ) + T (r x Collecting terms i n descending powers of r t l (E.2b) + T (r y ~ 2[T (r 2 s Uz - p) + T r x , + T lzz y+ 3yz y Jzz . - p )]p y z z , (E.2c) (E.2d) yields d r y = = _ * 6*/4a'-c' . ( E . 3 a ) v'aa aa' J-oo where a Syy ^zz S y z S;zz S h ' ~ ^ S + ^ S 149 S , (E.3b) (E.3c) Appendix E. Evaluation of the Generalized b\ Gaussian Overlap 9 2 S. Integral 150 9 yz (E.3d) xz u u Sz Z b\ Hi 5 2{^(T = c s + r ,ft, + T iw J x z P x = p) Jzz + T^p, + T - (T z 3xyPx )) 3zyPz , (E.3e) T p -2Tj (r -p )p 2 jyy y xy x x y (r,- r + r ^ r - , . - ^ ) - T Iz E - JyzPy T p )' J22 z <~>zz ( Jxz( - 2 T - r P*) - iyzPy)Pz + TjzzPl T (E.3f) > and the tensor S is defined by (E.3g) S = T i + Tj . Finally, terms are collected in descending powers of r to obtain x / / \/aa. i e - ( a r'+b r , + c ) d * = I ^ / i / 4 a -c (E.4a) II \/aa, a, J-co where <s _ 15 c _ C2 ' ^>yy^>zz ^> Q (E.4b) yz = B p, (E.4c) = C : pp . (E.4d) T h e components of the vector B are Tj J xx Q Q ^>yyJzz _ Q2 ^yx Tj Tj J xy Jxz yy yz l^yz Sy Z S Z (E.5a) Appendix E. Evaluation of the Generalized Gaussian T T Ixy -2 B„ Overlap q yz q q zx °zy u TB z q T3 zz 3yz q _ C2 O Szz Ti 3xz x yx u q q yy °yz S S u Sx Z (E.5b) Syz °yy >-> 151 Tn •> zy yy J q OyylJzz Integral z y (E.5c) z z T h e elements of the symmetric tensor C are (S Tj yz = T- zz S z(SyyS 3 xx Z yz yy J — yy S z(SyyS Z ] Z T S 3 zz (Sy Tj Z — J - j xz xy 3xz M Z — Tj yz z z (S y y S yz z z Tj" S ) Jzz zz — zz xy yz xy — S yz z z (S yz y y — z z — Sz Tj )(S Tj Z yy S (SyyS yz ZZ ZZ >zy) z z T S ) (E.6c) ) Ti Ti Jxz Jyz Szz T- z — S S ) yz T (E.6e) T 3zz Szz (E.6f) S z z lzy) T T- (E.6d) 3 zz 3 xz yz z — Tj yz — S — z z S ) zz S 2 yz — S Tj )(S Tj zz zy (E.6b) Jzz S Tj ) — S Tj )(S Tj S {Sy Tj Cyz — zz zz xz ~ 3yz Z Szz(SyyS (Sy*Tj c. 2 yy — Sy ) Z yz T- rp 2 S Tj ) zz (E.6a) Z >->zz z yz {S Tj — 1 T = C ^xy 2 xy — Sy ) zz (S Tj - rp 2 3x S Tj ) xz 3yz z z Appendix E.l E. Evaluation of the Generalized Gaussian Overlap Integral E v a l u a t i o n o f t h e F o r c e s a n d T o r q u e s for t h e G a u s s i a n O v e r l a p 152 Integral 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) T h e torque on particle i (at the origin) due to j (at p ) , T\-, (defined i n 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 %i V%o T h e calculation of d A d 5?,- 2 % 3 _ s ; requires the evaluation of partial derivatives of the form OO POO /"OO _ . / / sT' / -co J—oo J—oo :rre-^T'^+^^-P^'-P^r^r,^ , ih dq ik - where k = 0,1, 2, 3 and the elements of T T : B i f c (E.3a) are = ^ . (E.3b) Appendix E. Evaluation of the Generalized Gaussian Overlap 153 Integral T h e partial derivative can be written as = -y se { 2 s ,co *kxx J _ rco ^ ' .„ T :(i-p){-'-p)]dr dr dr oo 2 f J_ 0 0 ,00 roo yoo Tl yy J_ J_ J_ k oo 9[ Ii + j :s y + x 7_oo c o 00 ^ T . r+Tj:(r-P)^-P)]dr drydr r2 s[ r r z + x V 0O i - c o i - c o J-oo 7_oo 7-oo •/-co 1 " yco i-oo ,oo Ill 9J . T-kXZ j _ c J_co o ,oo r x z r - [T :Tr+Tf.(v-p)(r-p)] ir dr s e t drz( y + x J - O O yoo yoo -siTr.rx+Tj-.ir-PHr-Pndrzdrydr^ . (E.4) *yz y_oo J-oo -/-oo l T h e required integrations are laborious but straightforward. It is useful to choose the integration order such as to simplify the final expressions obtained. involving r 2 and r r x followed by dr y the most convenient order of integration is dr dr dr y z by dr ). x ,oo ,00 , / / / -co J — oo J — oo TT / 3 1 2 / / -oo «/ — O O 3/2 6 Aa" + +hr +c) z //2 b" ,co / y dr z // : // —c /4a (E.5) , r r e-( * * Ur dr dr x y ar +br +c z y = x \r ) x 11 //2 + s 6 // // /4a - c a' a" \ 2 ' 4a" where a, 6, c, a', a", b' and 6" are as defined i n the previous section. with r and r r , the most convenient order is dr dr dr . 2 y (i.e., «^—OO h 2d a" {ar x \ 2 ,oo CO a" x r e- * * dr dr dr 2 fl VaV7a" Vad y This yields oo 7T For the integrals z x z y (E.6) For the integrals This gives expressions of the form as the previous integration order, but with the parameters — Tixx +Ti 3 xx ' (E.7a) Appendix E. Evaluation of the Generalized Gaussian q a q Overlap 154 Integral _ C2 '-'zz'^'xx u xz (E.7b) q s XX U IT^ + Tjl a s _ki b_ s s b\ (ry) (E.7c) 2 ' SzzSxx S. xz + _ili r v (E.7d) s ' q q + 2 (S xz xy u (E.7e) u S x b\ Hi s 5* 2 I -f ( 3xxPx L + 3xyPy T + T T i*zPz) - ( J x , P * + iyzPv r T + 3zzP* T (E.7f) Again, b"/s can be written as the dot product of two vectors, b - = B p, (E.7g) 5 where q q '-''ra'-'.sz _ C2 xz q xx ^ xy ^xz Jxy T3xz T3 xx u Sx Z q xx u B„ = T 9 9 X q Szz zy q q ^xy yy J Szx Sy Sxx Sxy Z xz u T ^yz z q q ^xx^zz _ q2 xz u T3zx X q zx u 3zy Sy Z (E.7i) 1 Szz q -2 B (E.7h) X u T 3yx Sz T.T.UZ q q u (E.7j) T 3 zz Szz Appendix E. Evaluation of the Generalized Gaussian For the integrals involving r\ and r r , x Overlap 155 Integral the integration order dr dr dr z y x yields equa- z tions with - = T +T-Liyy T J- s (E.8a) , lyy C C _ C2 >~>xx>~'yy a s xy u 9 = (E.8b) ' yy u IT^ + Tjl a s S '21 5 = 2 S -f b = 5 (E.8d) yz (E.8e) S 9 <? xy u r u syy + T p +T p )-(T p 2(^(T p 3xy (E.8c) Q C _ Q2 ' ^xx'^yy ^xy x 3yy y Jzy z Jxx + T p x Jyx y + T p) Jzx 2 (E.8f) Also, we have - = B p, (E.8g) s where S XX 1 -2 p J D x ~ Q yx q q c-2 >->xx>~>yy >-> y T 3 XX u X q ^ B„ XX Syx Q Q _ Q2 <Jxx*~>yy ^ x q z = (E.8h) T 3 xz q ^ xy q Syy Szy q2 >-> yy J T (E.8i) 3yz Sz Syx q °yy Szy T3 zx Tj 3zy T, 3 ZZ -2 q q ^xx^yy T3xy q zy u q ^ xy ^XX z Q yy u Q '-'xz y Jyx B Q xy u X xy (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 «(»i) - [Q~V q i Q J U( J)]Q,) • ? (E.9) If we substitute r' = r — p in E q . ( E . l ) , V Q u ( i j ) can be evaluated by interchanging the j 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 3/2 2 6" /4a"-c" 2 7 r ^(p,q*,qj) = 7 e 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 i n 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 a temperature scan is performed on a system.with Particles 89 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. < P 2 Above T* = 1.9, the system is simply isotropic, as indicated by > « 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 cutsphere 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] g (r) 000 g ™< ) r 3< = y g (r) 000 i i o ( = < r *~ r ) 6(rj - TJ - r)[3(A,- • h ){fa 3 r ) = < E ,^ h3 S ( r , r ) [ f a , > (5.2.4a) • h ) - fa • fa] > ^ 3 • fa] > _ ( 5 { _ 2 > 4 b ) 5 2 A c ) 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 g (r). 112 T h e function g (r) 110 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||. T h e correlation functions are defined as = " i i v <B(r||,r| + ft |)> l 1 < n ( r | | , r | | + 6r ) 0 u > 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 |
File Format | 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 |
Graduation Date | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | 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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0059610/manifest