Monte Carlo Studies of Phase Transitions in Open Systems by Christopher David Daub B.Sc , The University of Guelph, 1997 M . S c , The University of Guelph, 1999 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S F O R 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 The Faculty of Graduate Studies (The Department of Chemistry) T H E U N I V E R S I T Y OF BRITISH C O L U M B I A November 24, 2004 © Christopher David Daub, 2004 11 Abst rac t This thesis consists of the investigation of two different problems involving phase transitions using grand canonical Monte Carlo simulations. The first problem sought to understand the nature of the condensation phase transi-tion in simple ionic fluid models, and in particular the restricted primitive model. A specific focus of this investigation was the resolution of a dis-crepancy between calculations of the constant volume heat capacity near the condensation critical point in the canonical ensemble versus the grand canonical ensemble. By doing calculations on a fluid of hard spheres with at-tractive power law interactions, we have demonstrated that the near-critical heat capacity calculated in the grand canonical ensemble leads to results which are inconsistent with analytical results for the criticality of these flu-ids. Although the reasons for this discrepancy remain unclear, we believe that it stems from an incomplete consideration of the finite size effects which affect the grand canonical calculation. We then return to the original prob-lem of the criticality of ionic fluids. We have simulated the critical region of a model of charged hard dumbbells, and we have done a mixed field fi-nite size scaling analysis and calculated the canonical heat capacity. We find no qualitative differences between the criticality of the dumbbells and the ionic model, suggesting that the pairing of ions does not affect the critical behaviour. We surmise that the lack of a peak in the canonical heat capacity of ionic fluid models is a consequence of the very slow convergence of this calculation toward the correct results in the thermodynamic limit, and is not a signature of non-Ising criticality. The second problem we investigated is the adsorption of CO2 on MgO crystals. This study was driven by experimental results indicating a phase transition between a low temperature, ordered monolayer and a high tem-perature, disordered monolayer. There is also experimental evidence that this disordered phase has a higher density than the ordered phase. We were able to simulate the formation of an unusual disordered high density phase, however, some discrepancies between the experimental results and our sim-Abstract i i i ulations remain to be resolved. iv Contents Abstract ii Contents iv List of Figures vi List of Tables viii List of Symbols i x Acknowledgements xii 1 Introduction 1 2 Theoretical Background and Simulation Methodology ... 4 2.1 Interparticle Potentials 4 2.1.1 Long-Range Forces 7 2.2 Monte Carlo Simulation 9 2.2.1 Grand Canonical Monte Carlo 13 2.2.2 Histogram Reweighting and Umbrella Sampling . . . . 16 2.3 Theories of Phase Transitions and Critical Phenomena . . . . 17 2.3.1 Critical Exponents and Mean Field Theories 18 2.3.2 Renormalization Group Theory 22 2.3.3 Mixed Field Finite Size Scaling 25 3 Phase Transitions in Coulombic Fluids 28 3.1 Survey of Past Work 29 3.1.1 Experimental Results 29 3.1.2 Analytical Theories 32 3.1.3 Simulation Results 36 3.2 Ensemble Dependence in Monte Carlo Estimates of Cy . . . . 39 Contents v 3.2.1 The Heat Capacity of Models with Pair Potentials u(r) oc • - l / r n 42 3.2.1.1 n=3.1 45 3.2.1.2 n=4 47 3.2.1.3 n=6 48 3.3 Simulations of Ionic Fluid Models 57 3.3.1 MFFSS Analyses 58 3.3.1.1 Restricted Primitive Model 59 f 3.3.1.2 Charged Hard Dumbbells 60 3.3.1.3 2:1 Primitive Model 62 3.3.2 Heat Capacity Measurements 64 3.4 Discussion 68 3.4.1 Ensemble Dependence of Cy 68 3.4.2 The Critical Peak in Cv for Ionic Fluids 74 4 The Adsorp t ion of C 0 2 Gas on M g O Crystals 76 4.1 Survey of Past Work 77 4.1.1 Analytical Treatments of Adsorption and Surface Phases 77 4.1.2 Experimental and Theoretical Work on the Adsorption of C 0 2 on MgO 78 4.2 Model and Simulation Methods 82 4.3 Results 89 4.4 Discussion 106 4.4.1 Comparion of Simulation vs. Experiment 106 4.4.2 The Structure of the High Density Monolayer 109 5 Conclusions and Future W o r k I l l 5.1 Simulations of the Condensation of Coulombic and Related Fluids I l l 5.2 Simulations of C 0 2 Adsorption on a MgO Substrate 114 Bibl iography 116 vi Lis t of Figures 3.1 Cv calculated for the R P M 43 3.2 Cv calculated for the n = 3.1 potential 45 3.3 Configurations of the n = 3.1 system for a box size L — 10/cr . 46 3.4 Cy calculated for the n = 4 potential 47 3.5 Configurations of the n = 4 system for a box size L = 10/CT . . 48 3.6 Cv calculated for the n = 6 potential 49 3.7 Configurations of the n = 6 system for a box size L = 12/a . . 50 3.8 Configurational energy calculations in the canonical ensemble for the n = 6 system 51 3.9 Cv calculations in the canonical ensemble for the n = 6 system. 52 3.10 Plot of the position of the crit ical peak in Cv for the n — 6 system 55 3.11 Plot of the maximum value of Cv for the n = 6 system . . . . 56 3.12 Matching to the Ising universal order parameter distribution for the R P M 59 3.13 Matching to the Ising universal order parameter distr ibution for C H D ' s 61 3.14 Extrapolat ion to L —> oo for T*(L) measured via M F F S S analysis of C H D ' s 63 3.15 Extrapolat ion to L —>• oo for p*c(L) measured via M F F S S anal-ysis of C H D ' s 64 3.16 Matching to the Ising universal order parameter distribution for the 2:1 P M with L = 15 65 3.17 Cv calculated for the R P M , without umbrella sampling . . . . 66 3.18 Cv calculated for the R P M , with umbrella sampling 67 3.19 Cv calculated for the C H D system 68 4.1 The (2y/2 x y/2)R45° and [y/2 x v^)J?45° structures of C 0 2 on M g O 81 4.2 Adsorption isotherms of C 0 2 on M g O at 80 K and 100 K . . . 90 List of Figures vii 4.3 A monolayer formed by C 0 2 on MgO at 80 K 91 4.4 Adsorption isotherms of C 0 2 on MgO at 150 K 92 4.5 Distribution of tilt angle 9 in the low-density monolayer . . . . 94 4.6 Distribution of azimuthal angle <p in the low-density monolayer 95 4.7 A dense monolayer formed by C 0 2 on MgO at 150 K on a square MgO lattice of size L — 12a 96 4.8 A dense monolayer formed by C 0 2 on MgO at 150 K on a square MgO lattice of size L = 24a 97 4.9 The average number of adsorbed C 0 2 molecules as a function of height, (N{z)), at T = 150 K 98 4.10 Distribution of tilt angle 9 in three monolayer coverages at T = 150 K 99 4.11 Distribution of azimuthal angle 0 in three monolayer coverages at T = 150 K 100 4.12 A dense monolayer formed by CO2 on MgO at 150 K on a square MgO lattice of size L = 11a 101 4.13 Monolayer structures of C 0 2 on MgO at 5 K 103 4.14 Monolayer structures of C 0 2 on MgO at 5 K, with the MgO substrate not shown 104 4.15 Distribution of tilt angle 9 in three monolayer coverages at T = 5 K 105 4.16 Distribution of azimuthal angle (f> in various monolayer struc-tures at T = 5 K 106 4.17 The average number of adsorbed C 0 2 molecules as a function of height, (N{z)), at T = 5 K 107 viii List of Tables 2.1 Critical exponents for some common universality classes . . . 23 3.1 The critical parameters of models with u(r) oc — l / r n pair interactions 44 3.2 Finite size scaling fits of the positions of the peak in CV • • • 54 3.3 Finite size scaling fits of the height of the maximum in Cy • • 56 3.4 The critical parameters of ionic fluid models 58 3.5 The size-dependent critical parameters of CHD's 62 4.1 The parameters in the electrostatic model for C 0 2 83 4.2 The potential parameters for the repusion and dispersion in-teractions between CO2 molecules 85 4.3 The potential parameters for the repusion and dispersion in-teractions between CO2 molecules and the MgO surface . . . . 86 4.4 The binding energy of monolayers at T — 150 K 101 4.5 The binding energy of monolayers at T = 5 K 102 ix Lis t of Symbols The same symbol has not been used to identify two different quantities except where absolutely necessary. In most of those cases, if a symbol has two possible meanings, the section with the alternate usage has been noted. We have also neglected to list quantities only mentioned in passing. The average of a quantity Y is always indicated by placing angle brackets, (Y), around the quantity. Reduced dimensionless quantities are always indi-cated with a superscript asterisk (*). A critical quantity is always indicated with a subscript "c". Symbol . Definition Section introduced a Interaction range 2.3 only van der Waals parameter 2.3.1 only Lattice separation 3.2 only Like-ion separation 4.2 A Born-Mayer repulsion parameter 2.1 a Born-Mayer steepness parameter 2.1 Heat capacity critical exponent 2.3.1 Polarizability 4.2 cx Underlying matrix of the Markov chain 2.2 b van der Waals parameter 2.3.1 B Adams' B parameter 2.2.1 P {kBT)-1 2.2 Coexistence curve critical exponent 2.3.1 Cn Dispersion series parameters 2.1 Cv Constant volume heat capacity 2.2 cP Constant pressure heat capacity 3.1.1 List of Symbols x Symbol Definition Section introduced d Number of spatial dimensions 2.3.1 7 Compressibility critical exponent 2.3.1 s Critical isotherm exponent 2.3.1 Creation function 2.2.1 Destruction function 2.2.1 Maximum trial spatial displacement 2.2 Maximum trial angular displacement 2.2 E Electric field 4.2 e Dielectric constant 2.1 Well depth 2.3 and 3.2.1 only eo Permittivity of free space 2.1 V Structure factor critical exponent 2.3.1 Viscosity 3.1.1 only fn Potential damping functions 2.1 h Planck's constant 2.2.1 k Wavevector 2.1.1 kB Boltzmann's constant 2.2 K Inverse Debye length 3.1.2 KT Isothermal compressibility 2.3.1 L Linear system size 2.1.1 A Thermal DeBroglie wavelength 2.2.1 Conductance 3.1.1 only m Particle mass 2.2.1 M Magnetization 2.3.1 M Magnetization operator 2.3.3 Chemical potential 2.2.1 Dipole moment 2.1 N Number of Particles 2.1 V Correlation length critical exponent 2.3.1 List of Symbols xi Symbol Definition Section introduced P Probability or probability distribution 2.2 Pressure 2.3.1 Transition probability from state m to state n 2.2 Azimuthal angle 4.3 ®zzzz Hexadecapole moment 4.2 $ Electrostatic potential 4.2 Q Charge 2.1 Qzz Quadrupole moment 4.2 r Interparticle separation 2.1 r Set of particle coordinates 2.2 r* Position of particle i with respect to the origin 2.1 P Particle density 2.2 P Limiting probability distribution 2.2 s Mixing parameter 2.3.3 a Hard sphere diameter 2.1 t Scaled temperature 2.3.1 T Temperature 2.2 T Number of Monte Carlo trials 2.2 9 Correction to scaling exponent 2.3.2 Til t angle 4.1.2 U,U(r) Potential energy 2.1 u Energy density 2.3.3 Interaction potential between particles i and j 2.1 V Volume 2.2 W Weighting function 2.2.2 W Probability distribution 3.2 £ Correlation length 2.3.1 X(A,B) Fluctuations in quantities A and B 2.2.1 z Height 4.2 Activity 4.2 only Z Configurational integral 2.2 C Discretization parameter 3.2 c Orientation vector 2.2 CO Weighting function 3.2 only Empirical correction to scaling exponent 3.2.1 only X l l Acknowledgements Oedipa wondered whether, at the end of this (if it was supposed to end), she too might not be left with only compiled memories of clues, announcements, intimations, but never the central truth itself, which must somehow each time be too bright for her memory to hold; which must always blaze out, destroying its own message irreversibly, leaving an overexposed blank when the ordinary world came back. In the space of a sip of dandelion wine it came to her that she would never know how many times such a seizure may already have visited, or how to grasp it should it visit again. Perhaps even in this last second -but there was no way to tell. - Thomas Pynchon, The Crying of Lot 49 My six years here in Vancouver at U B C sure haven't been boring! Looking back, I'd have to say if I could do it over again, I'd do some things differently. On the other hand, I've never believed that saying about not having any regrets. The way I see it, if you don't regret some things you've done with your life, then you haven't really been living. Come to think of it, if I could boil down my Vancouver experience to its most concentrated essence, that last sentence comes pretty close. Numerous people have contributed to the completion of this thesis. On the academic side, I'd of course like to thank my supervisor, Dr. Gren Patey. The first time I met him, he launched straight into an enthusiastic discussion of possible projects, and this attitude was enough to make me want to work with him. Since then, his vast knowledge and experience have helped me learn a great deal about science and academia in general. I also owe a great debt to all of the members of the Patey research group whom I have had the pleasure of sharing a lab with over the years. In particular, during my first two years at U B C Dr. Philip Camp was a great mentor to me, and never shirked from answering my naive questions about simulations and fluid physics. More recently, I had the pleasure of interacting with Dr. David Jack and the members of his research group at Concordia University. This A ckn owledgemen ts xiii taught me a great deal about the process of scientific collaboration, both its difficulties and its rewards. Outside of the Ivory Tower, I have enjoyed numerous great experiences and friendships in my time here, so many that I cannot hope to start naming names without forgetting someone, so I won't try! I will, however, thank the lovely and talented Karen Hardie for being such a great partner over the last year or so. She's changed my life in so many ways, it's hard for me to imagine where I'd be without her. I'd also like to thank Bob Hope for many great conversations. Finally, I'd like to thank my parents and my family. They're far away, but I still think about them a lot. After all, if it wasn't for them, I sure wouldn't be here now! Christopher Daub Vancouver, October 2004 Chapter 1 i Introduct ion The idea of using computers to simulate real processes has been with us as long as the computers themselves. Some of the first computer algorithms to be developed were the Metropolis Monte Carlo procedure [1] and John von Neumann's cellular automata [2], both of which took advantage of the computer's abil i ty to do simple mathematics very rapidly. B y the 1960's techniques for the simulation of physical and chemical systems were beginning to be developed, and simulation studies of models of these systems played a significant role in improving our knowledge of the intermolecular forces which hold matter together [3]. As the computer technology has advanced, ever larger and more complex systems have been simulated. Today, systems as complicated as the ion channels in cell membranes are routinely simulated. However, for the most part simulation studies have been done in closed systems, where the number of particles is fixed. This restriction limits the scope of simulation studies to situations where there is no exchange of matter between the system and its environment. In studies of systems near a cri t ical point, for example, fluctuations in density are an important feature, and it may be difficult to understand the crit ical region without incorporating these fluctuations in the simulation. Chemical systems of importance to biology are usually considered as open systems, where particles can be exchanged between the system and its environment. If computer simulations of biological systems are ever to be considered meaningful, they must take the open nature of these systems into account. The framework for doing simulations in open systems has been developed for many years, mainly due to the work of Gibbs in formulating the concept of statistical mechanical ensembles [4], although of course the first computer simulations in the grand canonical ensemble came much later [5]. B y sim-ulating a system in the grand canonical ensemble, which allows for density fluctuations, rather than in the canonical or microcanonical ensemble, which must keep the particle number fixed, the open nature of a chemical system Chapter 1. Introduction 2 may be taken into account in a rather straightforward way. The development of theories of phase transitions and critical phenom-ena such as renormalization group theory [6] have greatly improved our un-derstanding of such processes. Application of these theories to simulations of phase transitions has for the most part been focused on simulations of phase transitions in the canonical ensemble, and techniques such as finite size scaling analysis [7] have proven valuable in relating the results of these simulations to the theories of real phase transitions. What remains unclear is whether similar analyses can still be brought to bear when we do simula-tions in the grand canonical ensemble, without changing the framework of the theory in any major way. The main focus of this thesis is on Monte Carlo simulations in the grand canonical ensemble on systems undergoing phase transitions. We have ap-plied this technique in two separate projects. The first investigation involves the study of simple models of ionic fluids near their vapour-liquid, or con-densation, critical points. By studying this phase transition, we hoped to gain a better understanding of the critical region of ionic fluids, including a final proof of the universality class to which ionic fluids belong. However, in the end we were led to make another investigation into subtle differences be-tween the results of simulations in different statistical mechanical ensembles. These surprising results must now be taken into consideration in reformu-lating theories of finite size scaling which can be applied to grand canonical simulations. The second project is an investigation of adsorption of CO2 gas on a MgO substrate at low temperature and low pressure, using realistic interatomic po-tentials and simulation methods. Although systems similar to this have been studied by Monte Carlo simulation in the canonical ensemble, to our knowl-edge this is the first attempt to do so in the grand canonical ensemble. The use of the grand canonical ensemble is predicated in our case by experimen-tal data suggesting a transition from an ordered monolayer of low density to a higher density disordered structure above a certain temperature. This unusual behaviour, with a transition between a low density ordered phase and a higher density disordered structure, is reminiscent of the unique prop-erties possessed by water, which allows ice to float on liquid water. If we wish to investigate this phase transition, and elucidate the behaviour of the disordered phase, the density must be allowed to fluctuate and hence the use of the grand canonical ensemble is a necessity. This thesis is written in five chapters, beginning with the current intro-Chapter 1. Introduction 3 duction. In Chapter 2, all ofthe theoretical background needed to understand the results of the two research projects will be briefly outlined. Beginning with a discussion of intermolecular potentials, next Monte Carlo simulations in the canonical and grand canonical ensemble will be described. Then, a short review of theories of phase transitions and criticality will be given. In the next chapter, Chapter 3, the project to study condensation phase tran-sitions in fluids will be described. This begins with a review of the history of experimental, analytical and simulation studies of ionic criticality, leading into the observation of a conflict between the work of two different research groups. The resolution of this conflict by our work, the other investigations it inspired, and our return to the initial question of the universality class which ionic solutions belong to are the subject of the remainder of the chapter. Chapter 4 is devoted to the study of adsorption in the CG^/MgO system. The first part of the chapter introduces the general problem of physisorption of gases on cold crystal surfaces, leading into a review of the experimental work which drove our investigation. This is followed by a section on the specification of the model parameters for the system, as well as the means to compare our simulation results with those of the experiments. Finally, the results of our simulations are reported, and compared with the experimental results. The final chapter, Chapter 5, offers conclusions about both projects, speculates on problems which remain unanswered, and makes suggestions for future work. 4 C h a p t e r 2 T h e o r e t i c a l B a c k g r o u n d a n d S i m u l a t i o n M e t h o d o l o g y The general idea of atomistic computer simulation is to set up a particle configuration as a list of coordinates. Then, knowing the nature of the in-teractions between the particles, the potential energy U of the configuration can be calculated. In a classical molecular dynamics simulation, the con-figuration would then be evolved according to the Newtonian equations of motion. By contrast, in a Monte Carlo simulation, the configuration space is explored by attempting trial moves from the initial configuration. Moves are accepted or rejected according to certain criteria which depend on the statistical mechanical ensemble being simulated. In the most widely used en-semble, the canonical ensemble, the criterion depends on whether the energy of the configuration rises or falls after the move. 2.1 Interparticle Potentials A description of the forces acting between atoms and molecules is necessary before any calculation of the energy of a given particle configuration can be attempted. The potential field arising due to these forces will have terms for pairwise additive interactions u(ri,Tj) which depend on r$ and r,-, the positions of particles i and j. In principle, the total potential also contains higher order terms pertaining to many-body interactions. The many-body interactions are sometimes negligible, or can be included in an ad hoc manner by using an effective pair potential which accounts for the many-body inter-actions [3]. The total potential energy for all N particles can then simply be obtained from the sum of pairwise interactions, JV (2.1) i<j Chapter 2. Theoretical Background and Simulation Methodology 5 where the sum is taken over all unique pairs, as noted by the subscript i < j. This avoids double counting of pairs. In general, for electrically neutral particles the pair potential for two par-ticles separated by a distance r = |r| = |r! — r 2 | will consist of a repulsive part for small particle separations, and an attractive part for larger separa-tions. Unless the particles themselves are asymmetric, both of these terms depend only on the interparticle separations, and not on their orientations. If there are electrostatic forces present, either due to ionic charges or due to higher multipole moments, they must be included as well, and these terms may depend on the orientations of the particles. The total pair potential can therefore be expressed as the sum of three terms, «(ri, r 2) = urep{r) + uatt(r) + ueUc(ri, r 2). (2.2) The repulsive part of the potential accounts for the Pauli exclusion prin-ciple which forbids the overlap of the wavefunctions of fermionic matter, either due to the overlap of core electrons for small internuclear separations, or the overlap of the electrons with the nuclear wavefunctions themselves. The mathematical form chosen to model this effect may be as simple as a hard-core model, which excludes any interparticle separations less than a given diameter a. This is equivalent to saying that urep(r) = co for r < a. In other cases, the repulsive part may be modelled by an inverse power law, such as in the well-known Lennard-Jones potential which uses urep(r) oc 1 / r 1 2 . A more realistic model for the repulsive part is the Born-Mayer potential, given by urep(r) = Ae-ar, (2.3) where A and a are parameters related to the maximum of the potential at r = 0 and the slope of the exponential decay, respectively. While there is little theoretical justification for preferring one form of the repulsive potential over another, the Born-Mayer potential does provide a better fit to experimental data and ab initio SCF calculations [8]. The attractive part of the potential arises due to the induced dipole and higher multipole forces caused by the interactions of the electrons surround-ing each particle. For a spherically symmetric particle, the induced dipole forces give rise to a potential uatt(r) = C$/r6, and the higher multipoles give terms of all higher even powers in r, so that in general the attractive potential Chapter 2. Theoretical Background and Simulation Methodology 6 is described by a dispersion series 2n uatt(r) = J2C^r2n- (2-4) i=3 This series may be truncated at different points, depending on the level of accuracy required. The various models for the repulsive part and the attractive part of the intermolecular potential can be combined in different ways, and many of these combinations have proved to be useful for different applications [9]. The combination of a Born-Mayer repulsion term with the first two terms of the dispersion series (C 6 and C 8) is called the Buckingham potential. This potential suffers from the drawback that the Born-Mayer term reaches a maximum Ae~ar —> A as r —» 0, and hence is dominated by the dispersion terms for very small interparticle separations, since C2n/r 2n —> oo for r —> 0. This problem can be avoided by the introduction of damping functions for the dispersion series, which takes on the new form uatt(r) = -^2 n=3 L f2n{r) y2n (2.5) with the damping functions f2n(r) given by 2n , h n ( r ) - l - ^ ^ ' a r - (2-6) The combination of the Born-Mayer repulsion term with the damped disper-sion series is called the Tang-Toennies potential [8]. It is especially useful because the damping functions depend only on the Born-Mayer potential parameters, and so no additional data is required in the formulation of the potential. In the atomistic simulations we consider in this thesis, the particles are al-ways spherical, and so the core interactions are spherically symmetric and do not depend on the orientations. If the particles in question were asymmetric, as might be the case in a coarse-grained simulation of polymer chains, the modelling of directional interactions would require modifying these spheri-cally symmetric potentials to include the angular variations in the interpar-ticle forces. Some model systems will also be considered in this thesis which Chapter 2. Theoretical Background and Simulation Methodology 7 interact via longer ranged attractive pair potentials. Because these simple models are well-understood theoretically, they provide a valuable means of testing simulation techniques for use with more physically meaningful sys-tems. The electrostatic interactions between charges centered on two particles are described by the Coulomb potential, which in the SI system of units is expressed as ucouiir) = (2.7) where e0 = 8.854 x 10~ 1 2 C 2 ( J m ) _ 1 is the permittivity of the vacuum and qi are the charges on the particles. If the charges are the same sign, this is a repulsive interaction; if they are of opposite sign, it is an attractive interaction. There may also be interactions between higher multipoles, which have more complicated functional forms due to the vectorial nature of these interactions. If a permanent dipole moment Pi is placed at the centers of the two particles, there will be terms added to w e / e c ( r 1 , r 2 ) corresponding to the charge-dipole interactions and the dipole-dipole interactions, given by [10] and „ ( v r \ - ( M l ' ^ 3 ( / x 1 - r ) ( ^ 2 - r ) U d d { r i > r 2 ) - " W T " 47re0r« ' ( 2 ' 9 ) respectively. More complicated expressions arise for interactions between higher multipoles [10], but they will not be needed in this thesis. In this thesis models of solutions, particularily ionic solutions, will be considered. Often, it is unnecessary to consider the interactions between the solvent and the solute particles explicitly. Instead, the solvent can simply be described as a continuum with a dielectric constant of e. In this approx-imation, if the particles are in solution, all terms in the pair potential are modified by a factor of 1/e. 2.1.1 Long-Range Forces Any computer simulation must deal with the difficulty of modelling a large, macroscopic system with a model system small enough to be tractable. Tech-niques such as the use of periodic boundary conditions can mitigate against Chapter 2. Theoretical Background and Simulation Methodology 8 the presence of spurious surface effects, however, other artifacts of the sim-ulation procedure require more complex solutions. One of these problems is the accurate calculation of long-range interactions such as the Coulombic interaction between charges. Usually it is sufficient to ignore all interactions of the particle with its images outside of the central simulation box, but for long-range interactions the influence of the periodic images will be significant. A technique called Ewald summation [5,11] has been developed to ac-curately calculate the contribution of the periodic images. Although this technique has been applied to many different simulation geometries, we will only concern ourselves here with the use of cubic simulation boxes. It has been applied to cases including point charges or ions, point dipoles and point quadrupoles, but we will only consider point charges in this thesis. The Coulombic interaction between all ions in the simulation box, as well as all the periodic images of all the ions can be written as 1 1 , ( N N U c o u l = 2 4 ^ ^ ' ^ ^ m ^ u n \ i = l j=l Here n = [nxL, nyL, nzL] is a vector which identifies each image of the central simulation box, of linear dimension L. The prime indicates that for n = 0 the term i = j is not included in the summation, since the interaction of an ion with itself is not included (though its interaction with its own images is). The factor of 1/2 appears in order to cancel the double counting of interactions of ion i with ion j, and ion j with ion i. The Ewald summation technique consists of taking each point charge, and surrounding it with a charge distribution of equal magnitude and opposite sign. This distribution is usually taken to be a Gaussian, with the functional form p?(r) - qiK3exp{-K2r2)/^2. (2.11) The parameter K determines the width of the distribution, while r is the position relative to the center of the distribution. This distribution serves to screen the interaction of the point charges, such that the interactions will now be short-ranged and will converge better. Of course, the influence of this charge distribution must be removed in order to recover the proper potential energy for the original distribution of point charges. This is done by adding a cancelling distribution equal and opposite to the distribution of Equation 2.11, taking the Fourier transform, adding all of these contributions Chapter 2. Theoretical Background and Simulation Methodology 9 in reciprocal space, then transforming it back to real space. The final result of this calculation is given by [5] TT 1 1 V ^ V ^ I/ > T > ' e r / c(K | r i j + n|) 47ren 2 ' / — ' V ^ \TH + n u j=l j=l \ |n|=0 1 u 1 + (1/TTL 3) ^g io i(4 7r 2/A; 2)ea;p(-rc 2/4rv 2)cos(k • r^) k^o ' 1=1 The first term in this equation is the real-space sum of the point charges and their cancelling distributions. Here erfc(a;) is the complementary error function. If K is large enough, this summation converges rapidly, and only the central simulation box need be included (i.e. n = 0). The second term is the summation in reciprocal space, with k = 27rn /L 2 being the wavevectors included in the sum. Enough wavevectors must be included in the summation to ensure that it converges for the system under investigation. The third term is a correction term, due to the fact that the interaction of the cancelling distribution with its own point charge must be explicitly removed from the final result. This result allows for the calculation of the Coulombic potential energy of a large, spherical (if enough wavevectors are included in the calculation) array of simulation cells, each of which contains many charged ions, surrounded by a good conductor with a dielectric constant e' = oo. For a simulation where the charges are surrounded by vacuum (e = 1), an additional term is required to account for the formation of a dipolar layer on the surface, so that ^ C W , e ' = l _ Ucoul,t =oo 6 e 0 £ 3 E (2.13) 2.2 M o n t e C a r l o S i m u l a t i o n Monte Carlo simulation [1,5] is a method used for calculating the equilibrium macroscopic properties of a system given knowledge of the energy of a particle configuration. One of the chief tasks of traditional statistical mechanics is Chapter 2. Theoretical Background and Simulation Methodology 10 to calculate the configurational integral Z, given for a canonical ensemble of states with the same particle number N, volume V, and temperature T by where U(r) is the potential energy of the state specified by the set of particle coordinates r, and j3 = l/kBT, where kB — 1.381 x 1 0 ~ 2 3 J / K is Boltzmann's constant. From the configurational integral, quantities such as the average configurational internal energy (U) can be calculated quite easily. In principle, by simply generating rmax random particle configurations, calculating the energy of each configuration, and weighting the contribution of each configuration according to the Bol tzmann probability distribution, P ( r ) oc e~Pu(r\ one could get an estimate of the configurational integral using However, this would be a very inefficient process, since many ofthe randomly generated configurations would contribute negligibly to Z^VT- For example, configurations wi th significant overlaps between particles (so that \uaU\ « \urep\) would be very unlikely. A much better way of sampling states would be to concentrate on configurations which have some reasonable probability of occuring in the system. Monte Carlo simulations provide a way to explore the phase space in a very efficient manner. They also allow for the calculation of thermodynamic quantities like (U) without the knowledge of the entire configurational in-tegral. In the Metropolis Monte Carlo algorithm, the phase space of the system is traversed by means of a Markov chain, which is a sequence of tr ial configurations where the outcome of each tr ial depends only on the outcome of the t r ia l immediately before. The probability of going from state m to a new state n is given by a transition probability 7 r m n . In order to ensure that the l imi t ing probability distribution p is the correct one, corresponding to, for example, the canonical Bol tzmann distribution PNVT, if is necessary to ensure that the sequence of states generated obeys the principle of detailed balance (also known as microscopic reversibility), that is, the probability of going from state m to state n must be the same as the probability of the reverse, (2.14) max (2.15) T = l mn — Pn^i nm (2.16) Chapter 2. Theoretical Background and Simulation Methodology 11 where p m is the element of the limiting distribution corresponding to state m. The total transition probability for a state m to become any other state must be equal to 1. This places an additional constraint on the transition matrix, since 5>mn = l . (2.17) n A matrix with this property is also called a stochastic matrix. If this relation does not hold, it means that certain states of the system are inaccessible from other states. A system where all the states can be reached from every other state is called an ergodic system. Combining Equations 2.16 and 2.17 and summing over all states, we obtain the relation ^ ^ Prn^mn — ^ ^ Pn^nm — Pn ^ ^nm = Pn (2-18) m m m or, in vector notation, P7T = p (2.19) where p is the limiting distribution which satisfies the eigenvalue equation. As long as the initial distribution and the transition matrix satisfy the con-dition of ergodicity, after sufficient iterations of the eigenvalue equation the proper limiting distribution must be obtained. There remains some freedom in choosing the transition matrix with which to explore the phase space. The scheme suggested by Metropolis considers three cases, 7I"mn — Oimn, pn ^ Pm ^mn — ®mn{pn/'Pm)i Pn *^ Pm n^m Here cx is also called the underlying matrix of the Markov chain. It is a symmetric stochastic matrix, however, non-symmetric a matrices can be constructed which satisfy Equations 2.16 and 2.17. Other prescriptions for constructing a suitable transition matrix 7r also exist, but will not be dis-cussed here. Applying the Metropolis algorithm in an attempt to simulate a fluid sys-tem now requires the specification of the underlying matrix a. This matrix Chapter 2. Theoretical Background and Simulation Methodology 12 must take a state m and choose any neighbouring state n with uniform prob-ability. In Monte Carlo simulations of a canonical ensemble of a fluid system, a reasonable working definition of a neighbouring state is to choose one par-ticle i at random from the configuration, initially centered at r ^ , and move it with equal probability to a new center ri>n such that its maximum dis-placement in each direction x,y,z is A r m a x . Defining the number of points in a cube C of side Armax as Nr, the elements of a. are then O^mn — ^i,n £ C Cirnn = 0, Titn $ C. (2.21) In the simulation, a particle center will be chosen at random, and displaced a random distance between 0 and Armax in each direction. If the particles are molecules and not just spherical atoms, a scheme for rotating the particle randomly must be included as well. Care must be taken to ensure that these rotations are isotropic. One prescription for ensuring that the rotations are uniform in all directions is to take the randomly chosen molecule's orientation vector £, and rotate it within a spherical cap by some angle. Two random angles are generated, one in the interval [0, A6MAX], describing the size of the rotation, and one in the interval [0, 27r], describing the direction within the cap to rotate in. This procedure is only valid for axially symmetric particles. The translations and rotations may be attempted as part of the same move, or they may be done separately as two different trial moves. Once the translation and/or rotation of particle i is done, the move is accepted or rejected according to the probabilities defined by Equation 2.20. The difference in interaction energy between states m and n, AUnm = Un — Um, is calculated by summing over the pair interactions involving particle i before and after the move, N N AUnm = 53«(ry,n) - ^2u(rij>m),i ^ j. (2.22) j=l j=l If the move lowers the energy (A£7 n m < 0), then pn > pm and the move is accepted automatically. This is consistent with Equation 2.20, since now num. = Cimn. On the other hand, if the move increases the energy, then Pn < Pm and the move may be accepted, but only with a probability pn/pm-Again, this is done to stay consistent with the Metropolis scheme, since for this case iTmn — amn(pn/pm). The ratio of probabilities is easily expressed in Chapter 2. Theoretical Background and Simulation Methodology 13 terms of the ratio of Boltzmann factors of the two states, so that n 7,~l P-Pu* p-PUnp-PAUnm Pn _ ^ N V T e _ e e _ „-(9A<y n m (O O Q \ o ~ Z~l e-W™ ~ e-Wn ~ ' 1 ' Hm N V T Thus, the move is accepted with a probability dependent only on the Boltz-mann factor of the difference in energies between the initial and final states. If the move is rejected, then the next state in the Markov chain is a repeat of state m, in accordance with the final line of Equation 2.20. If the simulation is run for a sufficiently large number rmax of trial con-figurations, then after a period of equilibration, the configurational energy of the system can be estimated simply as the statistical average of the con-figurational energies of each configuration, (U)NVT = Sn=T ^he maximum sizes of the translations and rotations, A r m a x and A9max, should be adjusted to ensure that approximately 50% of the moves are accepted. This will ensure that the simulation results converge as quickly as possible. These maxima will depend on the system under study, as well as the ther-modynamic variables like the temperature T and the density p. As a general rule, larger moves can be supported at higher T and lower p, and vice versa. The calculation of the configurational energy as a statistical average also gives a direct route to calculating CV, the heat capacity at constant volume, since r r (du\ - ^ N V T ~ WNVT \dTjv hBT* This calculation can be compared with a direct differentiation of the energy calculated as a function of temperature in order to ascertain if a simulation has been allowed to run long enough [12]. Note that this calculation does not include the kinetic or ideal gas contribution to the heat capacity. If desired, this contribution can be added to the result afterwards, since the kinetic energy K is easy to calculate in simple models. For example, for a collection of N particles without internal degrees of freedom, (K) = (3/2)NkBT and CvMeal = (3/2)7VA;B. 2.2.1 Grand Canonical Monte Carlo In the grand canonical ensemble, the chemical potential p is held fixed while the number of particles fluctuates (T and V remain fixed). Grand canonical Monte Carlo (GCMC) simulations must therefore allow for the possibility Chapter 2. Theoretical Background and Simulation Methodology 14 of particle creations and destructions, in addition to the particle displace-ments already described [13]. For a particle destruction, a particle is selected randomly from the configuration, and the change in the interaction energy accompanying this removal is calculated. The ratio of the probability of the state n after the removal to the probability of the initial state m is Pm V In a creation attempt, a new particle is inserted into the configuration at a random location. The energy change due to this insertion is calculated, and the ratio of probabilities of the states is given by Pm A 3 ( / Y + l ) In both Equations 2.26 and 2.25, N is the number of particles in the initial state m and A is the thermal DeBroglie wavelength, defined by A=J ^ _ (2.27) with m here representing the mass of one particle, and h = 6.6261 x 10 - 3 4 Js being Planck's constant. The extra factors in Equations 2.26 and 2.25 in-volving p, N, A and V which are not a part of Equation 2.23 all arise in order to enforce the criterion of microscopic reversibility of the Markov chain when particle insertions or deletions are considered. Since p, A and V are all constants in the simulation, they may be replaced by a new function, called B, defined as [14,15] B = pp-ln(A3/V). (2.28) Equations 2.25 and 2.26 can be redefined in terms of B as PjL = e-PAUnm-B+\nN _ e - / 3 A D n m ^ 29) Pm and P n _ e-PAUnm+B-ln(N+l) _ e - / 3 A C „ m 2^ 30) Pm respectively, where ADnm and A C n r n are functions with the dimensions of energy which conveniently express the chances of accepting the removal or Chapter 2. Theoretical Background and Simulation Methodology 15 insertion of a particle. For the destruction, if ADnm < 0, the new state is accepted automatically, and the particle is removed from the configuration. If ADnm > 0, the new state is accepted with probability exp(—(3AD n m). For the creation, the same criteria for acceptance or rejection apply to the function ACnm. The GCMC simulations in this thesis were all executed at constant B, V and T. The chief advantage of the use of the quantity B is that it explicitly removes the thermal wavelength A from consideration. This is a quantity related to the kinetic energy of the particles in the simulation, which is not usually of interest. The chemical potential p can always be obtained after the simulation has been run by calculating the value of A for the system under study. As in the simulations in the canonical ensemble, the heat capacity and other response functions can be related to statistical fluctuations in the vari-ables of the GCMC simulation, but will be described by somewhat more complicated formulae due to the fact that the particle number fluctuates as well as the energy. In the case of the constant volume heat capacity, CV, the formula is [4] c^=^{x{u'u)-wn£) (231) where the notation X(A,B) refers to the quantity defined by X(A,B) = ((A-(A))(B-(B))), (2.32) so that, for example, X(U,U) is the mean-square fluctuation in energy and Equation 2.24 can be rewritten as C - X ^ U ) (2 33) <^V,NVT — , rpn • [Z.66) A few other points about GCMC simulations should be made. In or-der to maintain microscopic reversibility, on average there must be as many attempts to remove a particle as there are attempts to create a new parti-cle, unless appropriate weightings are applied to the acceptance criteria. It should also be noted that the technique as described here only applies to at-tempts to create or destroy single particles. Extensions to the method which can be applied to the creation or deletion of pairs of oppositely charged ions will be described in Section 3.2. This basic method also suffers from the problem that the probability of accepting a random creation or destruction, Chapter 2. Theoretical Background and Simulation Methodology 16 especially at low temperature and/or in dense phases, may be minuscule, leading to very low acceptance ratios and slow convergence of the simulation results. Depending on the system under study, some biasing techniques may be useful in improving the acceptance ratios. One such biasing technique will be described in Section 3.2. 2.2.2 Histogram Reweighting and Umbrella Sampling Canonical Monte Carlo simulations are run at a fixed particle number N, volume V, and temperature T. However, it is possible to scale the results of a simulation done at a temperature 7\ so that they apply to a neighbouring temperature T 2 [16]. This is done by noting that the probability Pi{T) of generating microstate i at temperature T is given by Pi(T) = v^-rr-\e~Ui/kBT- <2-34) Therefore, the relative probabilities of generating state i at temperatures Ti and T 2 will be given by PATi) ZNVT(T2) HL(-L.-A.\ ^ 7 7 f T = 7NVT)T2 e*aU n). (2.35) This means that the microstates generated at temperature Ti will also appear in a simulation at temperature T 2 , although with a different weighting u>i, given by wi(T1,T2) = - . (2.36) \ ekB\T2 Tj The denominator in the weighting equation 2.36 ensures that the reweighted probability distribution will remain properly normalized. The reweighted average energy at the temperature T 2 can then be calculated as (U(T2)) = J2W^>TM- (2-37) i=l This technique is known as histogram reweighting. Care must be taken to ensure that the two temperatures Ti and T 2 are sufficiently close, and Chapter 2. Theoretical Background and Simulation Methodology 17 that the simulation is run long enough, so that the reweighted probability distribution still has enough data to be considered an accurate distribution. In this way, a single long simulation at a single temperature can be used to provide a measurement of the average energy over a range of temperatures. Alternately, if very accurate results at a single temperature are desired, the results of simulations at several different temperatures may be recombined via the reweighting of equation 2.36 into a single calculation at one temper-ature. This is especially useful in the calculation of response functions such as the heat capacity, since accurate calculation of the fluctuations in energy ((AU)2) = (U2) — (U)2 depends on being able to sample the rare events at the edge of the probability distribution. This use of histogram reweighting is known as umbrella sampling [16]. Histogram reweighting techniques have also been extended to the grand canonical ensemble [17]. In this case, the probability P{(T, p) of obtaining a microstate i at temperature T and chemical potential p will depend on both the energy Ui and the number of particles Ni, and is given by P(T, p) = _ J _ e - W - ^ ) / * B T _ { 2 3 8 ) ^HVT{J- ) The microstates generated at T\ and pi will appear in a GCMC simulation at T 2 and pi with a weighting given by ekB\T2 TJ kByr2 T J v>i(TuT2,»UK) = 7 — — — (2-39) .7=1 2.3 T h e o r i e s o f P h a s e T r a n s i t i o n s a n d C r i t i c a l P h e n o m e n a Phase transitions have been a subject of experimental study for many years, dating back at least as far as the studies of Andrews on carbon dioxide [18]. Soon after, van der Waals developed his theory of liquid-vapour criticality [19], which was the first modern theory of phase transitions. In the context of this thesis, it is interesting to consider what minimal features the potential model must possess in order for a simulation of particles interacting via the potential to undergo the phase transitions seen in real Chapter 2. Theoretical Background and Simulation Methodology 18 world materials. The dimensionality of the system under investigation is important, but for the time being we shall only consider simulations of three dimensional systems. It is convenient to define a reduced, dimensionless density p* in terms of the particle diameter a as This definition of the reduced density corresponds to a value of p* = 1.0 for particles in a regular, simple cubic lattice. A system of hard spheres of diameter CT, interacting with a pair potential u(r;j) = co for pair separations rij < CT, and with no interactions for larger separations, will undergo a phase transition between a disordered phase with a reduced density p* — 0.9454 and an ordered, solid phase with p* = 1.0409 [20]. In fact, any system interacting only via repulsive interactions can only undergo transitions between solid and fluid phases, including the formation of liquid crystal phases for systems with highly anisotropic interactions. In order to generate condensation phase transitions between a dilute, gaseous phase and the denser liquid phase, it is necessary to include some sort of attractive interaction. Even the simplest such model, the so-called square well, with a pair potential given by (with e being the well depth and a being the range of the attraction) will undergo phase transitions between the vapour and the liquid state. Defining a reduced temperature T* as the ratio of kBT to the potential well depth, for the square well we get T* = kBT/e. The best current estimate of the liquid-vapour critical point of the square well fluid with interaction range a = 1.5a is T* = 1.2179 ± 0.0003 and p* = 0.3067 ± 0.0004 [21]. 2.3.1 Critical Exponents and Mean Field Theories The van der Waals model of matter [19] was designed to extend the concept of a non-interacting, ideal gas to a more realistic model interacting via a po-tential similar to the square well. It adds terms to the macroscopic equation of state of an ideal gas, P = NkBT/V = pkBT, which account for the volume p* = pa3 = (N/V)a3. (2.40) u{rij) ~ 0 0 ) r < ° = —e , CT < r < a = 0 ,r > a (2.41) Chapter 2. Theoretical Background and Simulation Methodology 19 taken up by the particles, as well as the attractive forces between particles. The van der Waals equation of state is where a and b are system dependent parameters which depend on the strength of the attraction between particles and the size of the particles, respectively. This equation can also be expressed as a cubic equation in p, abp3 - ap2 + p{kBT + Pb) - P = 0. (2.43) If the pressure P is expressed as a function of density, P(p), the behaviour of the equation differs depending on the value of T. At higher temperatures, there is a single solution to the cubic equation, and P is a monotonically increasing function of p. The critical point is defined as the point where the three roots of Equation 2.43 merge into a single solution. This point is defined by the conditions that (dP/dV)T=TC — 0 and (d2P'/'dV2)T=Tc — 0, and can be expressed in terms of the van der Waals parameters as kBTc = 8a/27b. If the temperature is lower than the critical temperature, there are three solutions to the equation. This behaviour is unphysical, as it implies that there must be a region where the isothermal compressibility, KT = —^(dV/dp)T, is negative, and this is forbidden thermodynamically [22]. Instead, the presence of multiple solutions to Equation 2.43 is interpreted as being an indication of the onset of phase coexistence. Starting from the gaseous state, if the system is compressed along an isotherm at a temperature below Tc, there will be some density pg at which the liquid will begin to condense out of the gas phase. The density of this liquid is pi. If the system is further compressed, the liquid will continue to condense, while the density of the remaining gas stays at pg, until all of the gas has condensed. Then, the density of the liquid will continue to increase from pi. The densities at which phase coexistence begins, pg and pi, can be deter-mined by Maxwell's equal-area construction. The curve described by these densities at different values of T < Tc is called the coexistence curve. At tem-peratures below Tc but not too far below, it it found that the width of the coexistence curve can be related to the scaled temperature t = (T — Tc)/Tc as Pi-Pg~ \tf (2-44) where /3 is termed a critical exponent. In the van der Waals model, the value of ti is 1/2. Chapter 2. Theoretical Background and Simulation Methodology 20 The behaviour of many other measurable quantities in the near critical region can also be described by critical exponents. The heat capacity at constant volume, Cy, in the critical region is described by the relation cv ~ | t | - a , T ^ r + ~ \t\~°' , T —>• T~ (2.45) with a and a being other critical exponents. In general, the exponent could be different depending on whether Tc is approached from above or below. In fact, the exponents are the same regardless of the direction of approach, so that a = a . For the van der Waals model, it is commonly stated [6] that Cv = Cv,\deai = (3/2)Nkg, so that a — 0 and Cv displays no un-usual behaviour through the phase transformation. This calculation, how-ever, mistakenly assumes that the van der Waals equation of state is correct in the coexistence region of the phase diagram below Tc. A better calculation considers the subcritical region more correctly by doing a Taylor expansion about the critical point and applying the conditions for phase coexistence (i.e. Pg — Pi, etc.) along the coexistence curve. This procedure leads to very different results, given to lowest order by [22] CV ~ Cv,ideal = \NkB (1 + gt + . . . ) ,*< 0 0, O O , (2.46) so that there is non-trivial structure in the heat capacity for T < Tc. More complete derivations of and expressions for this quantity and other critical properties of the van der Waals model are given in Sengers and Sengers [23] and most recently in Kim et al. [24], who extended the Taylor expansion to higher order terms. In any case, because all of these treatments assume that the equation of state must remain analytic at the critical point, the exponent a = 0 in all treatments of the van der Waals model, and though Cy may display discontinuities or cusps in the critical region, it does not diverge. The critical behaviour of the isothermal compressibility, KT, can also be described by a critical exponent. This relation is K r =4(f)T H (i"7 ( 2 - 4 7 ) where 7 is another critical exponent. The van der Waals model gives 7 = 1, and hence, KT diverges at the critical point. The physical interpretation of Chapter 2. Theoretical Background and Simulation Methodology 21 this phenomenon is that near the critical point, the system is extraordinarily sensitive to changes in pressure. A very small increase in P would shift the system towards a higher proportion of the liquid phase, and vice versa. Yet another critical exponent, 5, describes the shape of the critical isotherm: \P - Pc\ ~ \V - Vc\ 5 at T = Tc. (2.48) In the van der Waals case, 5 = 3. Taken together, the four exponents a, j3, 7 and S are termed the thermo-dynamic critical exponents, since they describe the critical behaviour of the macroscopic, thermodynamic variables. There are other exponents which describe the critical behaviour of the spatial correlations between the mi-croscopic particles of the system. One such exponent describes the critical behaviour of the correlation length £, itself defined in terms of the form of the two-point correlation function G(r), which describes the degree of correlation between one randomly selected particle and all of the pairs it can form with the other particles in the system. Far from the critical region, the form of G ( r ) is e - | r | / C C « - | r | ( r f - 1 ) / 2 ^ - 3 ) / 2 (2-49) where d is the number of dimensions. G ( r ) can be related to the isothemal compressibility KT, and it can be shown that if KT diverges, then £ also diverges, leading to the relation e~i*rv (2-50) with v = 1/2 for the van der Waals model. There is one additional ex-ponent, 77, which describes the divergence of the structure factor S(k) = p Jddre~lkvG(r) as the wavenumber k —> 0 at T = Tc. At the critical point, S(k) ~ k~2+v, and 77 = 0 for the van der Waals model. The divergence of the structure factor is observed experimentally as the phenomenon of critical opalescence. G(r) measures the fluctuations in density of the fluid, which are able to scatter light. The intensity of scattered light is proportional to S(k), and so, as it diverges near the critical point, the light is strongly scattered and the critical fluid appears milky or opaque. The preceding treatment briefly summarizes how to consider a van der Waals model of a fluid near criticality by means of critical exponents. Many other types of systems, from magnets to cosmological models of the early Chapter 2. Theoretical Background and Simulation Methodology 22 universe, also undergo phase transformations. Remarkably, many of these other systems can be described by the same critical exponents, provided suitable replacements of variables are made. In the case of the Ising model of a ferromagnet, for example, a very similar treatment can be applied which leads to the same values of the critical exponents, provided the width of the coexistence curve pi — pg is replaced by the net magnetization M and the isothermal compressibility KT is replaced by the isothermal susceptibility XT = (dM/dH)T, where H is the applied magnetic field, which can be related to the pressure of the fluid. In fact, any analysis of critical phenomena which makes the same as-sumptions as those made in the early theories of criticality in ferromagnets and fluids will lead to the same values for the critical exponents. The rea-sons for this were put on solid ground by the development of Landau theory. This showed that all treatments of criticality which assume that the relevant thermodynamic variables remain analytic functions at the critical point must lead to the same values of the critical exponents, for reasons related to the dimensions of the variables in the respective models. A l l of the these theo-ries can be collectively termed mean-field or classical theories of criticality. Landau theory also showed why systems as disparate as fluids and magnets could be described by the same theory, by showing that both systems could be characterized by an order parameter, which would exhibit the same near critical behaviour irrespective of the particular physical system in question. 2.3.2 Renormalization Group Theory Unfortunately, despite the compelling success of Landau theory in showing the connections between the critical behaviour of such wildly different sys-tems, there were two glaring problems, one of an experimental nature, the other theoretical. The experimental problem was that, in many observed sys-tems, the critical exponents were measurably different from those predicted by mean-field theory. The theoretical problem was that when Onsager de-rived the exact solution to the two-dimensional Ising model [25], the critical exponents he arrived at were different from the mean-field exponents. These problems led to a new treatment of critical phenomena called renormaliza-tion group theory [6,26], a much more technically sophisticated theory which diverges from Landau type theories by not assuming that the thermodynamic variables must behave analytically at the critical point. A thorough explanation of the details of renormalization group theory Chapter 2. Theoretical Background and Simulation Methodology 23 Universality Class Exponent: a P 7 6 Mean-field 0 1/2 1 3 2-D Ising [25] 0 1/8 7/4 15 3-D Ising [27] 0.109(4) 0.326(1) 1.240(1) 4.804 3-D Heisenberg [27] -0.12(1) 0.366(3) 1.390(5) 4.798 Universality Class Exponent: v V e Mean-field 1/2 0 2-D Ising 1 1/4 3-D Ising 0.630(1) 0.034(3) 0.504(8) 3-D Heisenberg 0.707(4) 0.035(3) 0.55(1) Table 2.1: Critical exponents for some common universality classes is beyond the scope of this thesis. However, some of the results obtained with the theory must be mentioned. One of these important results is the universality exhibited by near-critical physical systems, whereby the values of the critical exponents are determined by only three parameters. One of these is the number of spatial dimensions. If there were four dimensions, then all of the critical exponents would be equal to those obtained by the mean-field theories. The second factor is the number of degrees of freedom of the order parameter; for example, different critical behaviour is expected in a magnetic system modelled as a lattice of spins of absolute magnitude |S| — 1/2, depending on whether the spins can only take on values ± 1 / 2 (the Ising model) or if they are vectors which can point in all directions (the Heisenberg model). Finally, in some cases the range of the interactions is important. The effect of the interaction range will be explored more fully in the context of models of charged fluids in Section 3.2.1. Physical systems with the same values of these three parameters have the same critical exponents, and are said to belong to the same universality class. The critical exponents for some of the more commonly encountered universality classes are listed in Table 2.1. Renormalization group calculations also showed that the critical expo-nents were not all independent, but were in fact connected by scaling rela-tions. One such relation is the Rushbrooke scaling law, a + 2/3 + 7 = 2. (2.51) Chapter 2. Theoretical Background and Simulation Methodology 24 There is another relation which involves the exponent 5, PS = /3 + >y, (2-52) which we have used in the derivation of 8 in Table 2.1. There are also relations which apply to the correlation length exponents 77 and v, 2-a = du (2.53) and 7 = 1/(2-77). (2.54) Equation 2.53 is termed a hyperscaling relation, since it depends on the number of dimensions d. The hyperscaling relation is not obeyed, however, by the systems interacting via long-range attractive power law potentials which we study in Section 3.2.1 [6]. Another result which must be mentioned is the fact that the critical ex-ponents of mean-field theory describe only the leading order terms of the critical behaviour, and so can only be considered valid very near the critical point. Further away from the critical point, corrections to scaling come into play. The full scaling relations are expressed as a Wegner series [28], for example, for the heat capacity, Cvit) = A\t\~a (1 + B\t\e + . . . ) . (2.55) Here 9 is the correction-to-scaling exponent. Like the other critical expo-nents, it is the same for all systems in a given universality class. Values for 9 are listed with the other critical exponents in Table 2.1. Finally, finite-size effects on the critical behaviour must be mentioned. This is especially important in simulation studies of criticality, since all com-puter simulations are necessarily limited to a finite system. Because critical points, and phase transitions, are characterized by discontinuous changes in thermodynamic quantities, such as the density, true phase transitions can only be detected in the thermodynamic limit, where the number of particles can be considered to be infinite. In a simulation of a finite d-dimensional system of size L = Vlld, any discontinuities seen at a critical point must be replaced by a somewhat rounded peak. As an example, the divergent heat capacity at Tc of a system belonging to the Ising universality class will be replaced with a peak with a maximum at a temperature Tp(L) or a reduced temperature tp(L) somewhat below the true critical temperature Tc. The Chapter 2. Theoretical Background and Simulation Methodology 25 height of the peak Cv,max(L) will become larger with increasing system size. The scaling of peak position and peak height is given by the equations [6] tp(L) oc Z T 1 / V (2.56) and CV,max(tp,L) « L ^ , (2.57) respectively. If corrections to scaling are considered, these equations can have additional terms added to accomodate the corrections far from the critical point. Up until very recently, the correct form of the corrections to scaling for the peak positions was thought to lead to an adjusted form of Equation 2.56 [29], tp(L) oc ZT 1 / V (1 + blT 9'" + cL~llu). (2.58) However, new work has shown that when the contribution of the pressure is mixed into the formulation of the scaling fields (as briefly discussed in Section 2.3.3), the corrections to scaling lead instead,to a relation of the form [30] tp(L) oc ZT 1 / V (1 + b,L-0'u + cL-W"). (2.59) As for the peak heights, the corrections to scaling can be accomodated by including another parameter in Equation 2.57, resulting in an expression Cv,maX(tp, L) oc (L + 6L)a', (2.60) where 5L has been described as an "aid to extrapolation" [31] which accounts for the corrections to scaling in an approximate way. The question of the proper form of the corrections to scaling, and whether they may vary in different ensembles, is open and controversial. Some aspects of this will be discussed further along in the thesis. 2.3.3 Mixed Field Finite Size Scaling The study of systems near their critical points via Monte Carlo simulations was greatly aided by the development of the mixed field finite size scaling (MFFSS) technique, introduced by Wilding and his collaborators [32-34]. It has been relied on quite heavily in the study of Coulombic fluids described in Section 3.3 of this thesis. In all systems of a given universality class, it is possible to find an ordering operator M whose probability distribution PU(M) at the critical point is Chapter 2. Theoretical Background and Simulation Methodology 26 universal. For a model such as the Ising model, this operator is simply the magnetization. In a continuous fluid model, owing to the lack of particle-hole symmetry, the ordering operator is a linear combination of the density p and the energy density u = U/V, where s is a system-dependent mixing parameter which determines the degree of mixing of p and u. In recent work, it has been pointed out that M must also include a factor due to the pressure [21,24,30,35]. The inclusion of this factor leads to many novel effects, including the presence of a Yang-Yang anomaly in the heat capacity [35-38]. Although the MFFSS procedure employed here cannot account for the pressure term, in most cases these contributions are believed to be small and in this thesis we have made use the MFFSS technique in spite of its possible shortcomings. The universal ordering operator distribution PU(M) for a given univer-sality class can be obtained from G C M C simulations at the critical point of a system known to belong to that universality class. It is then possible to test if another system belongs to the same universality class by generating the distribution P(M) via G C M C simulations performed at the system's critical point and comparing it with the universal distribution. In practice, it is very difficult to perform the simulations precisely at the critical point, since the accurate determination of the critical point is usually one of the goals of the simulation. However, by taking advantage of the histogram reweighting techniques described in Section 2.2.2, the simulation can be run near the critical point, at a given value of T and p, and a joint distribution of energy density and particle density P(p, u) is generated. Then, the mixing parameter s is adjusted and the distribution is reweighted to new values of T and p such that the distribution P(M) matches the universal distribution PU(M) as closely as possible. This is most easily accomplished if the probability distribution is normalized, symmetric, and has unit variance. This is done with the aid of the finite size scaling relation for the standard deviation of the magnetization, given by A change of variables x = a^L^u, with a" 1 being a system-dependent pref-actor which guarantees that the distribution will be normalized, will give a distribution with the desired characteristics. M = p — su, (2.61) (2.62) Chapter 2. Theoretical Background and Simulation Methodology 27 The values of T and p which give the best fit to the universal ordering operator distribution are the best possible estimates of the critical temper-ature and chemical potential, Tc and pc. However, these estimates are still dependent on the system size L. Estimates of the true critical parameters in the thermodynamic limit L —» oo can be made by way of the appropriate finite size scaling relations, given by with /3C = l/kBTc. Plotting the size-dependent estimates TC(L) and pc(L) arrived at for several different system sizes versus L~(e+1^u and extrapolating to zero will give an estimate of the true critical parameters. The critical number density and energy density also obey finite size scaling relations, Tc(oo) - TC(L) ~ L - ^ 1 ) / " , / W < » ) - /3cpc(L) ~ L ~ ^ l \ (2.63) (2.64) pc(oo) - pe(L) ~ L-l'-M, uc(oo)-uc(L)~L-«-1'''l (2.65) (2.66) Estimates of the true critical parameters pc(oo) and uc(oo) can then be ob-tained by similar means. 28 C h a p t e r 3 P h a s e T r a n s i t i o n s i n C o u l o m b i c F l u i d s An ionic fluid consists of a fluid where the most important interactions take place between charged species. This may be either a molten salt, or an electrolyte solution where the dominant interactions are between the solvated ions. Our interest in this thesis is in condensation phase transitions occuring in ionic fluids. In a molten salt, this is nothing but a liquid-vapour phase transition similar to that taking place in non-ionic fluids, but at much higher temperatures. In solutions, the transition is between salt-rich and salt-poor regions of the phase diagram. By analogy with the discussion of the Van der Waals model in Section 2.3.1, we can understand what is seen experimentally in these phase transitions. With the temperature held constant below a critical temperature T c , as the mole fraction x of salt is increased past a certain mole fraction xi, a phase with a much higher mole fraction Xh will begin to form. This concentrated phase will remain in eqilibrium with the dilute phase, and the two phases will coexist until the overall mole fraction x > Xh- Then, the dilute phase will disappear and the mole fraction in the remaining concentrated phase will continue to increase. This sort of phase transition is known as a liquid-liquid demixing transition, because if the temperature is lowered below the coexistence curve, the solution will separate ("demix") into two phases, one of lower mole fraction xi and one of higher mole fraction x^. Experiments on ionic criticality usually focus on demixing transitions in electrolyte solutions. In the theoretical and simulation studies, the solvent is not usually treated explicitly and is most often modeled as a background continuum with dielectric constant es. Within this approximation, it is more natural to speak of a condensation phase transition taking place between a dilute vapour phase and a concentrated liquid or solid phase. We will use both of these descriptions interchangeably throughout this thesis. Chapter 3. Phase Transitions in Coulombic Fluids 29 3.1 S u r v e y o f P a s t W o r k The first experimental observations of a demixing phase transition in a liquid ionic mixture were reported in 1903 for a mixture of KI + SO2 [39]. A few other reports of demixing electrolyte solutions appeared in the next fifty years, many of which are mentioned in an article by Friedman [40]. In general, there was very little work done until, beginning in the 1980's, an abundance of work began to appear, both theoretical and experimental, and this continues to this day. The impetus which drove these investigations was the unusual nature of the critical region of some of these demixings. As this thesis is concerned with Monte Carlo simulations of models of ionic fluids, a thorough survey of the experimental results will not be at-tempted, although a brief summary is in order, since the investigations of ionic criticality were originally driven by experimental results. A good re-view of the experimental situation can be found, for example, in an article by Weingartner and Schroer [41]. The experimental summary will be followed by a more detailed summary of the past theoretical work, first considering analytical theories and then the results of simulations. 3.1.1 Experimental Results As already mentioned, molten salts are the most obvious example of an ionic fluid. Unfortunately, studying the liquid-vapour criticality of molten salts is very difficult owing to the high boiling points of these substances. As an example, the critical temperature for NaCl has been estimated to be Tc = 3300 K [42]. In order to facilitate comparisons between different systems, a reduced temperature T* can be defined as before, as the ratio of / c s T to the potential depth. For a symmetrical 1:1 electrolyte with cations and anions having equal diameters a and ionic charges ±q, this gives a reduced temperature T* = ea/Pq2, (3.1) where e is the background dielectric constant and a is the mean diameter of the ions. For molten salts, e is the relative permittivity of the vacuum e = eo-For NaCl, this gives an estimate T* ~ 0.058. The critical mass density of NaCl has been estimated as 0.18 g/cm 3 [42]. Taking the mean ionic diameter a ~ 2.76 A for NaCl , this leads to a reduced critical density, p* ~ 0.08. As can be seen by comparing these values with the critical parameters of non-ionic fluids, such as those for the square well fluid mentioned in Section Chapter 3. Phase Transitions in Coulombic Fluids 30 2.3, the reduced critical parameters of ionic fluids are much lower than in non-ionic fluids. In two component ionic fluids, such as aqueous solutions, e in Equation 3.1 is the dielectric constant of the solvent, es. Assuming the estimated T* for ionic fluids to be approximately correct and universal, this suggests that a monovalent salt dissolved in a solvent of dielectic constant es ~ 5 would have a critical temperature around room tempearture, and so the critical region would be more readily accessible to experiment. Familiar 1:1 ionic aqueous solutions, such as NaCl dissolved in F£20, are not good candidates for study, since if water is the solvent es ~ 80. At room temperature, these solutions exist in a state far above the criti-cal temperature, and are well approximated as dilute solutions. The wa-ter in these solutions would freeze before the critical temperature for the demixing phase transition was reached. Liquid-liquid demixing has been ob-served in a few aqueous systems involving large organic ions, such as tetra-n-pentylammonium bromide, with a critical point at Tc = 409K [43]. How-ever, in these cases the criticality is driven by the hydrophobicity of the alkyl groups, as opposed to the ionic nature of the salts. These solvophobic demixings must be considered separately from demixings driven primarily by electrostatic effects [Coulombic demixings) [44-46]. Like liquid-vapour phase transitions in non-ionic fluids, such as the condensation of noble gases, the criticality of the solvophobic demixings is unambiguously characterized by the critical exponents of the Ising universality class. More promising systems for observing a truly Coulombic demixing in-volve replacing the water with a different polar solvent such as an alcohol or an ether. The first such investigations by Pitzer's group on a solution of tetra-n-butylammonium picrate (Bu4NPic) in 1-chloroheptane [47] quite surprisingly gave evidence that the coexistence curve was best described by a critical exponent (3 = 1/2, the value one would expect from classical, mean-field theory. However, owing to experimental difficulties these measurements were subject to an uncertainty of A T ~ 0.1K, leaving open the possibil-ity that this system exhibits crossover from mean-field to Ising behaviour very near the critical point. They subsequently investigated a solution of triethyl n-hexyl ammonium triethyl n-hexyl borate in diphenyl ether [48], and in this case they were able to measure the temperature much more pre-cisely. The experimental findings were fit better with a coexistence curve exponent j3 equal to the classical value of {3 = 1/2, rather than the Ising value /3 ~ 0.326. In retrospect, the findings on the latter system have been Chapter 3. Phase Transitions in Coulombic Fluids 31 brought into question by conflicting results [49], which may perhaps be due to chemical instabilities in the system [50]. Nevertheless, these findings led to many more experiments being done, mostly on the Bu 4 NPic system in a wide variety of different solvents, mainly long-chain alcohols [51]. The gen-eral result from these investigations was that the criticality was Ising-like, with deviations towards mean-field results as the solvophobic effects were lowered by increasing the size of the solvent molecules. The early experiments by Pitzer's group consisted simply of visual obser-vation of the formation of meniscii signalling the onset of phase separation. Their later experiments measured the difference in the refractive index be-tween the phases, since this difference is proportional to the mole fraction of salt and hence is a good choice for the order parameter [44,48]. Since their work, experiments measuring many different quantities near the critical point have been undertaken. These include measurements of light scattering [52] and turbidity [53], which are especially useful since they provide a route to the measurement of the scattering exponent u, and hence the compressibility exponent 7 via the scaling relation of Equation 2.54. Like the coexistence curve measurements on the same system, the turbidity results on Pitzer's bo-rate/diphenyl ether system are ambiguous, with the first experiment giving a value of 7 = 1.01 ± 0.01, consistent with mean-field criticality, at a range oft between 1 0 _ 1 and 1 0 - 4 [53], but a second experiment by the same group gave conflicting results which supported Ising criticality [49]. Turbidity mea-surements on the Bu 4 NPic/alcohol systems are more consistent. These data indicate a crossover between mean-field criticality far enough away from Tc and Ising criticality as the critical point is approached more closely [54]. Measurements of the viscosity 77 near the critical point can also give in-formation about the universality class. Viscosity data for Bu4NPic in 1-tridecanol are consistent with the turbidity data, supporting Ising criticality with a crossover to mean-field further from Tc [55]. Another experiment di-rectly measured by calorimetry the constant pressure heat capacity Cp of a mixture of ethylammonium nitrate in n-octanol [56]. These measurements showed a critical divergence of the heat capacity which was well described by an exponent a — 0.11, consistent with Ising criticality. The broad picture, then, for Coulombic criticality, would seem to suggest that although the asymptotic behaviour is indeed Ising-like, there is an un-usual crossover to mean-field, classical behaviour at temperatures quite close to the critical point. This is quite different than the critical behaviour of ordinary, non-ionic fluids, which is firmly and unambiguously Ising-like, with Chapter 3. Phase Transitions in Coulombic Fluids 32 crossover to mean-field exponents only seen quite far away from the critical region, if at all. One final interesting set of experiments, namely measurements of the electrical conductance A, should be mentioned. Conductance measurements in the vicinity of Tc on the Bu 4NPic system in 1-tridecanol [57] and 1-chloroheptane [58] both show a minimum in the value of the Walden product of conductance x viscosity, An, at a salt concentration corresponding to a re-duced density of « O.Olp*. This is an indication that in the dilute phase, the salt exists predominantly as neutral ion pairs, which are non-conducting. At very low densities, these pairs dissociate and the solution becomes conduct-ing. At densities above the critical density, the ionic structure is liquid-like, and conducting. Below Tc, the phase separation is between a dilute, non-conducting phase, corresponding to a "gas" of bound ion pairs, and a dense, conducting phase, corresponding to a complex liquid-like structure. As we shall see in the next sections, these results are also seen in analytical theories and in simulation results. The nature of the ion pairing near the critical point has been hypothesized to be a key factor influencing the critical be-haviour in the vicinity of the demixing phase transition in Coulombic ionic liquids [59,60]. 3.1.2 Analytical Theories The starting point for analytical theories of strong electrolyte solutions is Debye-Hiickel (DH) theory [61]. We will focus on the simplest model of an electrolyte solution, called the restricted primitive model (RPM), which consists of N equisized hard spheres of diameter cr, half of which have charge q, and half of which have the opposite charge —q, immersed in a solvent of dielectric constant e. For this restricted case, DH theory shows that the effective pair potential between widely separated ions becomes UDH(T) = Y—> (3-2) 47re 0er where the additional factor of e~KT arises as a consequence of the ionic at-mosphere of oppositely charged ions which builds up around each charge in solution. This will screen the effect of the central charge and lower the ef-fective potential. A measure of the magnitude of this screening is provided by K, which is also called the inverse Debye length. If, instead of the ionic atmosphere, a single charge of opposite sign were placed at a distance 1/re, Chapter 3. Phase Transitions in Coulombic Fluids 33 the effect on the potential would be exactly the same. For the R P M , the Debye length 1/K is given by 1 / W V " . ( 3 3 ) Although the Debye-Huckel theory has been studied for many years, only fairly recently was it established that the theory implies that a phase coex-istence similar to what is observed experimentally in ionic solutions should occur in the R P M [62-64]. In this analytic approach, the pressure P(T, pj, pj) is expressed as a function of the temperature T and the chemical potentials Pj and densities pj of each species j. If there are multiple values of the den-sities or chemical potentials which give the same pressure, then this implies that multiple phases with these different densities are coexisting. The critical point is defined according to the usual conditions, dP\ dVj and d2P dV = 0, (3.4) T,N = 0: (3.5) T,N The equation of state for the pressure of a system of N ions can be expressed as P = NkBT/V + Pxs, (3.6) with the first term being the ideal gas pressure and the second including all other effects. According to D H theory, the excess pressure due to the electrostatic interactions between the ions is described by the equation Pxs = -~^-Jl + Ka- —^ 21n(l + Ka)) . (3.7) Solving the equation of state according to the conditions of Equations 3.4 and 3.5 yields values for the critical parameters [62] T* = 1/16 = 0.0625 and p* = 1/64-7T w 0.005. The value for T* is comparable to the experimental results, but the critical density is too low by an order of magnitude. The reason for this discrepancy is due to the neglect of ionic association in the simplest D H theory. As suggested by Bjerrum not long after Debye and Hiickel's work Chapter 3. Phase Transitions in Coulombic Fluids 34 was published [65], not all ions in an electrolyte solution are dissociated, especially at low temperatures, and certainly in the critical region. Instead, most of the ions form cation-anion pairs. Using a theory of association due to Fuoss [66], McGahay and Tomozawa found that the degree of dissociation in the R P M at the critical point was only 3.2% [63]. Since the associated ions are considered to be neutral and inert in Bjerrum's theory, simply dividing the critical density by the degree of dissociation gives a new value of p* = 0.1594. Soon after, Fisher and Levin used the slightly different treatment of association originally suggested by Bjerrum to find a value of p* = 0.045 [64], quite close to the experimentally observed value (the value of T* is unchanged by the inclusion of ionic association effects). It should also be noted that Friedman and Larsen [67] used a different thermodynamic stability condition to find that a modified Bjerrum theory gave a coexistence curve with a critical point at Tc = 1/16, pc = 0.0237 a full decade before McGahay's work and Fisher's subsequent study were published. Friedman's approach, however, leads to a strange, "spiked" coexistence curve, which Fisher argues misrepresents the Bjerrum theory [46]. The simple Bjerrum-like treatment of ion association is incomplete in its turn due to the fact that it treats the associated ion pairs as completely neu-tral, inert species. This leads to a peculiar, highly skewed coexistence curve, where the gas phase persists up to quite high densities at low temperatures, where there are few free ions. Because the paired ions are ignored in the treatment, there is nothing in the model which can spur the formation of the liquid phase at the correct density. To begin to correct this, Fisher and Levin modelled the ion pairs as dipoles, which interact with the free ions. The inclusion of this effect makes the coexistence curve look much more rea-sonable. When the excluded volume effect of the hard spheres is also taken into account, values of T* = 0.057 and p* = 0.028 are obtained [64]. Further improvements to these DH-based methods are possible. For ex-ample, the influence of the formation of charged trimers, neutral tetramers, and even larger clusters in addition to the ion pairs can be included using a treatment similar to that of Gillan's [68]. However, since these larger clus-ters only form at temperatures and densities far below the critical point, they have little effect on the criticality. Another improvement comes from considering the effects of dipole-dipole interactions between pairs, in addi-tion to the dipole-ion contribution. Although this has been done in several different approximate ways [69,70], the results for the critical temperature do not change appreciably, and the results for the critical density actually Chapter 3. Phase Transitions in Coulombic Fluids 35 become worse (p* — 0.00863 in one treatment [71]). There are other analytic approaches which can predict the shape of the coexistence curve and the position of the liquid-vapour critical point of the R P M . In fact, the first published study to do so antedates the D H based treatments by over a decade. This was a study by Stell et ai, which used a variety of integral equation theories to arrive at different estimates of the coexistence curve and the critical point [72]. Their best estimate of the R P M critical point after considering all of their results was T* = 0.085 and p* = 0.011. Stell and his coworkers have improved on the theory since then, using as the starting point the mean spherical approximation (MSA) to the solution of the Ornstein-Zernike equation for the R P M . As in the DH-based approaches, these advances consist of improving on the M S A by accounting for ion association using approximate schemes. In one approximation, they consider the pairing of the ions, but as in the early attempts by Fisher, they do not explicitly consider the contribution to the thermodynamics of the as-sociated pairs. This approach leads to a critical point at T* = 0.0748 and p* = 0.025 [73], a marked improvement on the first M S A results. Further im-provements to the theory come from considering the influence of the dipolar pairs themselves, with the most accurate approach giving a critical point at T* = 0.0745 and p* = 0.0245. As in the DH-based approaches, there is little improvement once the ionic pairing is accounted for in some reasonable way. In addition, attempts to improve on the MSA using better approximations to the solution of the Ornstein-Zernike equation, such as the hypernetted chain approximation, fail due to numerical instabilities which appear in the solution near the critical region [74]. The best analytical result comes from a recent calculation by Jiang et al. [75] which is based on the MSA, but does not involve solving the Ornstein-Zernike equation. It incorporates an improvement on the M S A called the binding-MSA (BIMSA) in order to ac-count for the influence of the dipolar pairs. This calculation locates the R P M critical point at T* = 0.0525, p* = 0.0640. To summarize the analytical work on the R P M , then, both the Debye-Hiickel based approaches and the integral equation approaches predict a liquid-vapour critical point with reasonable critical parameters, provided the effects of ion association are accounted for. However, further improvements to these theories do not lead to any great improvements in the estimated crit-ical parameters, especially when compared with the results of recent Monte Carlo calculations on the R P M , which we shall examine in Section 3.1.3. The best available Monte Carlo estimates of the critical parameters are [76] Chapter 3. Phase Transitions in Coulombic Fluids 36 T* = 0.04917 ± 0.00002 and p* = 0.080 ± 0.005. Both analytical approaches severely underestimate the critical density, and while the DH-based theories predict a reasonable critical temperature, the integral equation based theories overestimate T* by about 20%. The reasons for the failure of these theories vary from case to case. How-ever, since they are all mean-field theories, one thing they all fail to account for correctly is the non-analyticity ofthe order parameter at the critical point. As such, all these theories predict mean-field critical exponents, including co-existence curves characterized by a critical exponent /3 = 1/2. Therefore, it is not surprising that there are inconsistencies in the critical parameters derived by these theories. To go beyond the level of these theories would require a renormalization group treatment of the R P M , and although there have been some advances in this direction [77,78], a complete R G treatment is still lacking. There are, however, model systems of relevance to the discussion' of R P M criticality which have been successfully described by R G analyses. These are systems interacting via attractive pairwise potentials of the form U(r) oc — r~n. It has been rigorously shown by R G calculations [79-81] that such systems exhibit three regimes of critical behaviour if the number of di-mensions d < 4. If n > 2 + d — r], with 77 the correlation length critical exponent, the critical behaviour is that of the Ising universality class. Con-versely, if n < 3d/2, the critical exponents are mean-field. In between these limits, the critical exponents interpolate between the mean-field and Ising values. 3.1.3 Simulation Results The earliest Monte Carlo simulations to attempt to locate the critical point of the condensation phase transition in the R P M were undertaken by Vorontsov-Veliaminov et al. in 1970 [82]. They made use of the isothermal-isobaric or NPT ensemble in their calculations, in which instead of the volume, the pressure is held constant. This necessitates the inclusion of volume changes as part of the Markov chain to explore the phase space. Their estimate of the critical point was T* = 0.094 and p* w 0.35. However, like many early computational studies, this result is not trustworthy due to a lack of consideration of the necessity of allowing for sufficiently long equilibration and the like. They also failed to properly account for the long-range nature of the Coulombic interactions by using Ewald summation. Chapter 3. Phase Transitions in Coulombic Fluids 37 The next Monte Carlo attack on the problem came from Valleau, who simulated a system of 32 ions in the canonical ensemble at several different values of T* and p* [83]. His estimate of the critical point was T* — 0.07, p*c = 0.07, but this value is invalidated by the small size of the system and the use of a minimum image cutoff in the potential instead of the more accurate Ewald summation. Soon after, a study by Panagiotopoulos [84] arrived at an estimate of T* — 0.056, p*c = 0.04, including Ewald summation and making use of a new simulation technique called Gibbs ensemble simulation, in which two distinct simulations are run at the same temperature, but connected via exchanges of particles and volume between the two simulations. The acceptance criteria for these exchanges are set up in such a way to ensure that the systems remain at equilibrium, with the same chemical potential and the same pressure in both simulation cells. After equilibration, and in the region where phase coexistence is expected, each simulation cell will correspond to one of the coexisting phases. Another study by Caillol et al. [85] arrived at an estimate of T* = 0.057 and p*c = 0.04. This was a simulation which again used the Gibbs ensemble technique, but instead of Ewald summation, the simulation is run on the surface of a four-dimensional hypersphere. This is a more efficient way to accurately account for long-range interactions, with Caillol et al. citing a decrease by a factor of three in the computer time required for the simulation of a typical configuration of several hundered ions [86]. Thus, for several years the conventional wisdom [46] held that the critical point for condensation in the R P M was in the range of T* between 0.052 and 0.056 and p* between 0.023 and 0.035. However, a return to the problem by simulation researchers revealed that the earlier work had not accounted properly for finite size effects. In addition, it had become clear that Gibbs en-semble simulations, while very useful for simulating coexisting phases, could not adequately describe the near-critical region [87]. By running simulations in the grand canonical ensemble, and applying the MFFSS analysis tech-niques of Bruce and Wilding [32] to the problem, reliable estimates of the critical point were finally able to be obtained. The first such estimate was that of Caillol [86], who reported T* = 0.0487 ± 0.0001, p* = 0.065 ± 0.002. Soon after, the same group revised their estimate to T* = 0.0488 ± 0.0002, p* = 0.080 ±0.005 [88]. Panagiotopoulos also returned to the fray, arriving at an estimate of T* = 0.0490 ± 0.0003 and p* = 0.070 ± 0.005 [89]. The group of de Pablo also made an attempt [90], using a technique they developed called hyper-parallel tempering Monte Carlo in conjunction with MFFSS to Chapter 3. Phase Transitions in Coulombic Fluids 38 get an estimate of T* = 0.0492 ± 0.0003 and p* = 0.062 ± 0.005. The mag-nitude of the finite size effects in this study is larger than that seen by the other groups, which explains why their estimate of the critical density is so much lower than the others, but the reason for the larger finite size effects themselves remains unclear. There have been two very recent simulation studies, one by the Caillol group [76] and one by Panagiotopoulos [91], which extended their analyses to larger system sizes. The critical point estimates were T* = 0.04917±0.00002, p* = 0.080±0.005 and Tc* = 0.0489±0.0003, p* = 0.076±0.003, respectively. Despite the fact that some discrepancies remain, perhaps due to the different boundary conditions used (hyperspherical versus Ewald summation), or due to the use of a discretized R P M in the Panagiotopoulos study, both of these results seem very reliable, and any further change to these estimates outside the statistical uncertainty is very unlikely. Besides providing the means to accurately estimate the critical point, the MFFSS analysis also provided the first evidence from simulations about the universality class of the R P M condensation. Because the MFFSS analy-sis requires the fitting of the ordering operator distribution to the universal distribution at the critical point of the appropriate universality class, an as-sumption about the universality class must be made before the analysis can proceed. It is hoped that, if this assumption should prove to be incorrect, then a poor fit to the universal distribution or other internal inconsistencies will appear in the analysis which indicate that the universality class is differ-ent. Camp and Patey have shown that if Ising universality class is assumed for systems which are mean-field, inconsistencies do appear in the MFFSS analysis [60]. However, it is not at all clear that this should always be the case. Bearing that in mind, all of the MFFSS analyses of the R P M critical point were made assuming Ising criticality, and no inconsistencies have been reported. The same cannot be said for another piece of evidence in the debate on R P M criticality. This concerns the behaviour of the constant volume heat capacity, Cy. Valleau and Torrie made the first attempt at this measure-ment [92], making use of their temperature and density scaling Monte Carlo technique to examine Cy at several values of temperature and density close to the critical point, in the canonical ensemble, using several different system sizes. They measured Cy as a function of both temperature and density, and reported that there were no peaks or divergences along both near-critical isochores and isotherms that would indicate Ising criticality. In addition, Chapter 3. Phase Transitions in Coulombic Fluids 39 they demonstrated that for a simulation of a non-ionic fluid interacting via a Lennard-Jones pair potential (so-called "Lennard-Jonesium"), known to belong to the Ising class, there is a peak in Cy expressed as a function of density along a near-critical isotherm, and it gets sharper with increasing system size. Soon after, Luijten, Fisher and Panagiotopoulos published their own heat capacity calculations [93]. They objected to the results of Valleau and Tor-rie on three grounds. First, Valleau and Torrie examined the behaviour of Cy at reduced temperatures T* > T*, when due to the finite-size round-ing, peaks should be expected somewhat below T*. Secondly, the use of the canonical ensemble, rather than the grand canonical, might suppress the density fluctuations near the critical point, and hence the finite size rounding of the divergence would be greater. Third, the much lower critical density of the R P M compared with that of Lennard-Jonesium might be sufficient to explain the lack of a Cy peak in the R P M , and not a difference in the uni-versality class. Luijten et al. [93] supported their assertions with their own calculations of the heat capacity, using a discretized version of the R P M and running their simulations in the grand canonical ensemble. These calcula-tions showed that there was a peak in Cy measured as a function of T* along a near-critical isochore, located somewhat below T* and getting sharper with increasing system size. 3.2 E n s e m b l e D e p e n d e n c e i n M o n t e C a r l o E s t i m a t e s o f Cy The conflict between the results of the studies by Valleau and Torrie, and those of Luijten, Fisher and Panagiotopoulos (hereafter in this section re-ferred to as V T [92] and L F P [93], respectively, for convenience) led us to undertake our own investigation into the behaviour of the near-critical heat capacity of the R P M . These results have been published [94,95]. In their response to L F P [96], V T pointed to the possibility that the use of the discretized version of the R P M by L F P might be the cause of the discrepancies between the two investigations. It is known from analytical [97-99] and simulation [100] studies that the simple lattice gas version of the R P M , with the ions limited to positions on the sites of a cubic lattice, displays quite different behaviour than the continuum R P M . In particular, Chapter 3. Phase Transitions in Coulombic Fluids 40 the lattice gas R P M has a line of phase transitions separating charge-ordered and charge-disordered phases, ending in a tricritical point at a density p*cp = 1/3, and nothing similar is observed in the continuum R P M . The finely discretized lattice R P M studied by L F P is somewhat different, in that the lattice size a is smaller than the ionic diameter a. The fineness of the lattice is characterized by a discretization parameter £ = o/a. Although simulation studies have shown that the liquid-vapour criticality of the continuum R P M is essentially the same as that of the lattice R P M with ( > 3, with only slight differences in the location of the critical point [100], it is conceivable that there might be subtle differences between the models which would manifest as small differences in the critical behaviour of response functions like Cy. To answer this question, we performed G C M C simulations near the crit-ical point for both the finely discretized lattice R P M with £ = 5 and the continuum R P M , for the same simulation box size (expressed as a multi-ple of the particle diameter o) L/o = 10. The critical point of the lat-tice R P M with C = 5 has been estimated as T* = 0.05069 ± 0.00002, p* = 0.0790 ±0.0025 [101]. However, at the time we did our calculations, the critical density in the lattice R P M was thought to be smaller, and we used the earlier estimate of p* = 0.068 (the same estimate used by L F P [93]) as the fixed density for our canonical simulations. The continuum R P M canon-ical simulations were done at a reduced density p* = 0.079, within the range of the estimates of p* for the R P M . The simulations implemented the Ewald summation method, with conducting boundary conditions (e' = oo). Some calculations were repeated with vacuum boundaries (e' = 1), since these were the same boundaries implemented by LFP . There were no important differences in the results obtained with different boundary conditions. In the G C M C simulations, the chemical potential p was tuned to ensure that the critical density was being sampled. Because we used no reweighting techniques in these calculations, it was very difficult for us to extend our G C M C calculations below T*. We did, however, make use of the biasing scheme of Orkoulas and Panagiotopoulos [102] to improve the acceptance ratios for the creation and deletion of ion pairs. In this method, in a creation attempt, a cation is placed randomly, and an associated anion is placed at a distance r from the cation according to a decaying exponential distribution W(r) of the form e —ar W(r) C ,o < r < L/2 (3.8) Chapter 3. Phase Transitions in Coulombic Fluids 41 with a being a positive parameter with units of inverse length which controls the steepness of the decay, and C being the normalization constant, C=- (e-aa - e~La/2) . (3.9) a In this way, the interionic distance r is generated with a frequency W(r)/W0(r) compared with the random placement of the counterion. This must be ac-counted for in the acceptance criteria to ensure that microscopic reversibility is maintained. W0(r) denotes the probability density for the random inser-tion, which is proportional to the area of the spherical shell at radius r, and therefore, 47rr 2 Wo(r) = - j T 3 - , a < r < L / 2 . (3.10) For insertions at a distance greater than L/2, the expression for Wo(r) is quite complicated, as it must allow correctly for insertions into the corners of the simulation box [103]. For this reason, we have chosen to restrict all of our creation and deletion attempts to ionic separations less than L/2. At near-critical densities, the probability of accepting the creation or removal of widely separated ion pairs is infinitesimally small, and our results are unaffected by the restriction to r < L/2. In an attempted pair deletion, out of a number of ions N = 2NP, the first ion is selected at random, regardless of sign. Then, one of the Np counterions is selected with a probability Pi which depends on its distance rj from the first ion according to a weighting ufa), so that P, = -g™-. (3.11) The weighting function we use is an exponential, of the form The acceptance probabilities for insertion and deletion are modified from the expressions of Equations 2.25 and 2.26 to account for the insertion and dele-tion of ion pairs rather than single particles, as well as the biased weightings Chapter 3. Phase Transitions in Coulombic Fluids 42 for selecting or placing counterions. The final expression for the relative prob-abilities of states n (the new state) and m (the initial state) in the Markov chain is, for the removal of an ion pair, Vm~ 6 { V ) [NpPl(r)}m • { 6 A S ) The additional factor of arises from the deletion of a pair as opposed to a single particle. The rest of the additional factors are present to maintain the microscopic reversibility of the Markov chain. The subscripts m and n are reminders that the quantities indicated are calculated with respect to the correct state. For a pair insertion, the relative probabilities are given by PL = e-wme-to ( y ym+m{r)]n Pm \(Np + l)A*J [W(r)/W0(r))m- ^ 1 Applying these biased criteria gave much better acceptance ratios for the particle insertions and deletions. We did not use any pair or cluster displace-ment moves, as the benefits of using such complicated Monte Carlo strategies did not seem to outweigh the time required to implement them in our study. Once sufficient data was collected, Cy was calculated using the appro-priate fluctuation formula for the grand canonical ensemble, given in Equa-tion 2.31. These results are shown in Figure 3.1, along with the results of LFP . A l l results are shown with the kinetic or ideal gas contribution Cvjdeai = (3/2)Af£;£ included. It is clear that we have been able to reproduce (within statistical uncertainty) the results of L F P ; moreover, simulations of the continuum R P M in the grand canonical ensemble lead to very similar values of the heat capacity, including the presence of a peak just below T*. However, when these same systems are simulated in the canonical ensemble, and their heat capacities calculated with Equation 2.33, no peak is evident. In the canonical calculations, Cy was also determined from differentiation of the internal energy U modelled as a third order polynomial in y/T*. Both methods of determining Cy give comparable results. 3.2.1 The Heat Capacity of Models with Pair Potentials u(r) oc —l/rn Apparently, there is some discrepancy arising with the use of different sim-ulation ensembles. The next question is, which ensemble, if any, should be Chapter 3. Phase Transitions in Coulombic Fluids 43 -0.2 0 0.2 0.4 t Figure 3.1: Cv calculated for the R P M . The system sizes are given as a ratio to the hard sphere diameter o. Triangles represent our canonical re-sults, while diamonds indicate our grand canonical results. The dashed line represents the results of L F P [93]. used in the calculation of Cy if the results are to be trusted as a signature of criticality? To answer this question, we did similar simulations of a system of hard spheres interacting via attractive power-law potentials of the form u(r) = -e(o/r)n, r>o (3.15) = oo, r < a, (3.16) where e denotes the well depth and o is the hard sphere diameter. Reduced parameters are defined in the usual way, so that T* = ksT/e and p* = po3. A reduced energy is also defined as U* = U/e. As mentioned above in Section 3.1.2, the universality class which these systems fall into depends on the value of the power-law exponent n. In three dimensions, if n < 4.5, the criticality is classical or mean-field, while if n > 5 — r\ ~ 4.966, the system belongs to the 3-D Ising universality class. The critical behaviour of the heat capacity measured in simulations of these systems should follow from the Chapter 3. Phase Transitions in Coulombic Fluids 44 universality class. For the Ising systems, with a ~ 0.109, there should be a peak corresponding to a divergence rounded by finite size effects. For the mean-field systems, a variety of behaviours are consistent with the fact that the critical exponent a — 0. Either there is no critical feature in the heat capacity, or there could perhaps be a discontinuity in the heat capacity, or a logarithmic divergence. In any case, if the heat capacity is to be used as a signature of criticality, there should be an unambiguous distinction between the behaviours of the near-critical heat capacity in systems belonging to different universality classes. We investigated systems interacting according to power-law potentials with exponents n = 3.1, n = 4 and n = 6. In earlier work in our group [60], Camp and Patey investigated the critical region of these same three models, and used MFFSS analyses to determine the critical points of these systems. These critical parameters are listed in Table 3.1. They also demonstrated that discrepancies appear in the MFFSS analyses of the systems whose uni-versality class differed from the assumed Ising form used in the analysis. In fact, the MFFSS study was in many ways a precursor to the heat capacity study considered in this thesis, since it sought to determine how simula-tion results could be used to discriminate between two different universality classes. Table 3.1: The critical parameters of models with u(r) oc —l/rn pair inter-actions. The numbers in parentheses represent the uncertainty in the last digit. A l l data is taken from Ref. [60]. The pair interactions were cut off for all interparticle separations r > L/2, and a long-range correction ULR was included by assuming that the pair correlation function g(r) —>• 1 for r —>• oo. This gives [60] n 3.1 4 6 11.452(8) 0.247(5) 1.3724(1) 0.2993(1) 0.5972(1) 0.3757(4) (3-17) The constant volume heat capacities were measured as functions of T along isochores as close as possible to the critical isochores, given the finite Chapter 3. Phase Transitions in Coulombic Fluids 45 size of the simulations. In order to facilitate comparisons between different systems, all results are plotted as functions of t — (T* —T*)/T*. Calculations were done in both the canonical and the grand canonical ensembles, and the heat capacities were determined by collecting the appropriate histograms [p(U) for the canonical, p(U, N) for the grand canonical) and using Equations 2.33 and 2.31, respectively. The kinetic or ideal gas contribution to the heat capacity Cv,ideai = (3/2)NkB was added as well. We have published the results of this section. The discussion of ensemble dependence in all three of the power law models we study is reported along with the discussion of RPM criticality in References [94] and [95]. The de-tailed results of our more focused study on the n = 6 system are published in Reference [104]. 3.2.1.1 n=3.1 Figure 3.2: Cy calculated for the n — 3.1 potential. The results for n = 3.1 are shown in Figure 3.2. They provide striking con-firmation of the discrepancy between Cy calculations in the two ensembles. In the canonical calculations, the heat capacity function is rather featureless, Chapter 3. Phase Transitions in Coulombic Fluids 46 Figure 3.3: Configurations of the n — 3.1 system from canonical simulations with a box size L = 10/a. Snapshots are shown for temperatures (from left to right) t = —0.64, t = —0.16, and t = 0.15. There are 250 particles in the simulation box. showing only a slight decrease with increasing temperature which cannot be resolved on the graph. The kinetic contribution to CV is much larger than the configurational part in this system, and in fact the total heat capacity is not far off the results that would be obtained from a system of non-interacting hard spheres, for which Cy = Cy^deai = (3/2)JV&£. This is perhaps not sur-prising, when we consider that for a system with n = 3, the thermodynamic l imit does not exist [6], the energy density U/V diverges and there would be no measurable heat capacity. With n = 3.1, the configurational energy per particle calculated in the canonical ensemble is rather constant over the en-tire range of the simulations, with a value (U*X)NVT/N ~ —15.9. Snapshots of the n — 3.1 system at temperatures above and below the phase transition are displayed in Figure 3.3. In agreement with the behaviour of (U*X)NVT/N, the snapshots do not show much variation with temperature. In the grand canonical simulations, on the other hand, a large peak is seen in Cy with both system sizes studied (L/a = 6 and 10). The peak is sharper and closer to T* with the larger system size, just as a finite-size rounded divergence ought to do. In fact, in the absence of all other information, based on the GCMC results one would be tempted to conclude that the n = 3.1 system belongs to the Ising universality class, contrary to the RG calculations. The presence of this peak is supported by the GCMC results for the configurational energy (U*X), which drops dramatically as T* is approached from above. At T* — 11.0, our lowest temperature GCMC Chapter 3. Phase Transitions in Coulombic Fluids 47 result with Lo = 10, (U*x)^VT/'(N)uVT « -20.6. 3.2.1.2 n = 4 12 10 8 NkB g 4 2 \ ^ L/cr = 6, N V T T ( + L / a = 10, N V T L/cr = 6, G C • T .. \ \ L \ \ i L/cr = 10, G C \ \ H ^ f e ^ A •—o -0.6 -0.4 -0.2 0 0.2 0.4 t Figure 3.4: Cy calculated for the n = 4 potential. The heat capacity measurements for n = 4 are shown in Figure 3.4. The grand canonical calculations again show a large peak, with a maximum somewhat higher than that seen in the n — 3.1 system, but lacking any qualitative differences. The canonical calculations do not show any sign of a critical divergence. There is an increase in the heat capacity at low temperatures, but this does not appear to be related to the liquid-vapour criticality. The weak peak or cusp in Cy at rj « —0.4 in the Ljo = 10 system is more likely related to the liquid-solid phase transition, since the triple point Tt* for n = 4 is located at T* ~ 0.83, corresponding to t ~ -0.395 [105]. In partial support of this, configurational snapshots of the system at temperatures below t = —0.4 show that the system is associating into large, solid-like clusters at these low temperatures, yet they still are amorphous and there is no sign of an order-disorder transition taking place. It is possible that at even lower temperatures, there will be signs of an ordering transition Chapter 3. Phase Transitions in Coulombic Fluids 48 Figure 3.5: Configurations of the n — 4 system from canonical simulations with a box size L = 10/cr. Snapshots are shown for temperatures (from left to right) t = -0.63, t = -0 .41 , and * = -0.10. There are 300 particles in the simulation box. to a solid phase, but at the same time there is no reason to assume that ordering must occur even at T = 0 in a system at the critical density. Sample configurations of the n = 4 system at three different temperatures are shown in Figure 3.5. The results for n = 4 are interesting as well because of their connection to the question of the criticality ofthe RPM. Some authors have speculated that the longest range effective interactions in the RPM are power-law attractions with n = 4, caused by the interactions of dipolar pairs with free ions. We will have more to say on this subject in Section 3.3. 3.2.1.3 n=6 Cv results for the n = 6 system are displayed in Figure 3.6. Once again, other than an increased maximum in the peak heights as compared to the n = 3.1 and n = 4 systems, there is little change in the GCMC results. There are, however, significant differences in the canonical ensemble results. There now appears to be a peak below T*, in agreement with the RG analysis and the expectations of finite size scaling. The Cv peak in the canonical calculation is smaller in magnitude and has a maximum at temperatures lower than the corresponding peak in the GCMC calculation. I t seems clear that the canonical calculations of Cy give a better indication of the universality class than do the GCMC simulations, since the GCMC simulations erroneously suggest Ising criticality for systems which are clearly mean-field. Chapter 3. Phase Transitions in Coulombic Fluids 49 -0.6 - 0 . 4 - 0 . 2 0 0.2 0.4 0.6 t Figure 3.6: Cy calculated for the n — 6 potential. There remains, however, the possibility that the peak seen in the canoni-cal ensemble calculations is not a critical peak. It is well known that in other related systems, effects such as ion pairing can cause peaks in Cy which are not related to criticality, but appear near enough to the critical point that they could be mistaken for critical peaks [59]. There is also the possibility that these peaks are signatures of the freezing phase transition, like the cusp seen in the n = 4 results. T£ in the n = 6 system is located at T* ~ 0.5 or t ~ —0.16 [105], which is near where the Cy peak appears in the smaller systems studied. There is also a sharp peak near t = —0.4, which is certainly the signature of an ordering transition between the liquid and an ordered, crystal-like structure. This is borne out by the energy measurements, which show a marked discontinuity near t — —0.4, as well as by snapshots of the system at t < —0.4, which show the system in a crystalline, hexagonally close packed arrangement. Of course, this is not an accurate representation of the solid phase, since the density is much lower in our simulations. Snapshots of the n = 6 system at three different temperatures are displayed in Figure 3.7. To verify that the broad peak in Cy just below T* is indeed a critical divergence rounded by finite size effects, a more detailed investigation of Chapter 3. Phase Transitions in Coulombic Fluids 50 Figure 3.7: Configurations of the n = 6 system from canonical simulation with a box size L = 12/o. Snapshots are shown for temperatures (from left to right) t = -0.50, t = -0.14, and t = 0.09. There are 649 particles in the simulation box. the n = 6 system was undertaken. If the peak is a true critical peak, then the peak should sharpen, with a higher maximum closer to the true T*, as the system size is increased, in agreement with finite size scaling theory and Equations 2.56 and 2.57. Cy and (U*,x) were calculated in the canonical ensemble in a broad temperature range along the critical isochore for system sizes L/o — 6, 10, 12, 17 and 22. The heat capacity was calculated both directly from the fluctuation formula (Equation 2.33) and from differentiation of a [4,4] Pade approximant fit to the calculated configurational energy per particle as a function of T, (U*X)(T)/N. The [I,J] Pade approximant used has the form (U;X)(T)/N = ^ , (3.18) 1 + Ybs* where x = y/T* and a* and bj are the fit parameters. Although a function of this form cannot fit a true divergence, it has been found to be adequate for fitting the finite size rounded peaks measured in simulations [60,92]. The energy fits are shown in Figure 3.8. The Fade approximant fit was made to the energies calculated at temperatures within the liquid regime of each system, and so the range over which the fits were done varies with system size. The constant volume heat capacity found from the fluctuation Chapter 3. Phase Transitions in Coulombic Fluids 51 (Uex) N -3 - 4 . 5 - 0 . 6 0 . 4 - 0 . 2 t 0 . 2 Figure 3.8: Configurational energy calculations and Pade approximant fits for the n = 6 system in the canonical ensemble. System sizes shown are as follows: diamonds, Ljo = 6; squares, Ljo = 10; stars, L/o = 12; triangles, L/a = 17; inverted triangles, Ljo = 22. formula, as well as the differentiation of the Pade approximant fit to the configurational energy, are displayed in Figure 3.9. In order to pinpoint the location of the Cy critical peak as precisely as possible, the fits considered only the measurements below T* = 0.6. The Pade approximant form is unable to correctly incorporate both the high temperature results and the critical region, as can be seen by the poor fitting of the Cy results in the high temperature region. Nevertheless, the quality of the Pade approximant fits to (U*x) is very good, with the fits passing right through the individual points in all cases. It is very difficult to estimate the error in a Pade approximant fit in any quantitative way, unless one simply repeats the fit several times on different sets of data and estimates the statistical uncertainty. Because this would have been prohibitively time consuming, we have been forced to rely on a more ad hoc means of estimating the uncertainties in the fits. We attempted to fit different forms of the Pade approximant, with a different number of Chapter 3. Phase Transitions in Coulombic Fluids 52 NkB 40 30 20 10 0 • oo-Q-o 0 0 0 0.6 - 0 . 4 0.2 0 0.2 t Figure 3.9: Constant volume heat capacity calculations and Pade approxi-mant fits for the n = 6 system in the canonical ensemble. The symbols are the same as in Figure 3.8. For clarity, the plots have been vertically offset by 5 units for each increase in system size. The error bars represent the uncertainty in the location and the height of the maximum in Cy. terms in the power series expansions. Although there is little difference between the [4,4] and [5,5] fits, the Cy peak position does change somewhat, with differences of « 1% in the estimate of T*(L). We incorporate this in our estimate of the error. The location of the Cy peak has been found as accurately as possible, and does not seem to be unduly sensitive to the fitting procedure, provided the fitting is restricted solely to the critical region. The error bars in Figure 3.9 give an approximate idea of the uncertainty in the peak positions and peak heights, but due to the complex nature of the fitting problem, they should not be interpreted as being true statistical uncertainties. The location of the critical peak in Cy is subject to a finite size scaling re-lation given by Equation 2.56, or with corrections to scaling, either Equation 2.58 or Equation 2.59. By fitting the peak positions to these equations, and extrapolating as L —> oo, estimates of the critical parameters, and perhaps Chapter 3. Phase Transitions in Coulombic Fluids 53 the critical exponents themselves, can be obtained. Equation 2.56 can be re-expressed as T ( r ( L ) - r ; ( c o ) = C l L - 1 ^ (3.19) and a nonlinear least-squares fitting can be attempted to find the optimal values of Tc*(oo), cx and v. A direct fit to Equation 3.19 gives an estimate of v K, 0.86 which is somewhat different from the accepted Ising value v « 0.630. However, we believe that the poor estimate of v is simply a manifestation of the fact that the more parameters one wishes to fit, the more data is required, and we have data points at only 5 different values of L. This is also reflected in large error estimates for C\ in the free fits to v. Since we are unable to obtain a reliable value for the exponent v, we instead simply fix it at the Ising value in the rest of our fits. By fixing the value of v, a somewhat more robust fit can be obtained to Equation 3.19, which again does not include the corrections to scaling. If the corrections to scaling are considered, the fit to the results at smaller system sizes should improve. However, there are difficulties in attempting to fit data for a small number of different system sizes to the equations which account for corrections to scaling, again due to the lack of data points. Therefore, following Ferrenberg and Landau [29], we have chosen to attempt a fit to a function of the form T*{L) - rc*(oo) = ClL-1/u{l + c2L-u). (3.20) This equation folds all of the corrections to scaling into a single exponent, LO. As discussed in [29], all of the hypothetical correction terms are numerically close (1/u « 1.587, 0/v « 0.80, and 2/3/V « 1.035), so it may be hoped that the corrections to scaling could be accounted for in some average way by a single correction term. We are in fact unable to obtain a meaningful value for the correction to scaling exponent in our fits, due to a lack of data over enough different system sizes. However, by arbitrarily fixing the correction to scaling exponent to = 1, and fitting for c\, ci and T*(oo), we can fit our peak position data for all system sizes. The results of the best fits to Equation 3.19 and to Equation 3.20 are shown in Table 3.2. The raw data and the fits are also shown in Figure 3.10, with the results for T*(L) found from the MFFSS analysis [60] and the corresponding finite size scaling fit also shown for comparative purposes. The magnitude of the finite size effects is much more severe in the Cy extrapolation compared with the MFFSS analysis. The extrapolated Tc*(oo) arrived at by the Cy route also differs from the Chapter 3. Phase Transitions in Coulombic Fluids 54 MFFSS result by 2 to 4%, but given the uncertainty in the peak positions and the lack of simulation data for larger systems, this is not surprising. A fit to Equation 3.20 which fixes the extrapolated T*(oo) to the MFFSS result also passes through the peak positions, once the uncertainty is taken into account. Range Tc(oo) Cl V c 2 CO L > 6 0.566(5) -7(3) 0.52(6) L > 6 0.594(12) -1.7(3) 0.86(10) L > 6 0.575(2) -3.8(1) 0.630+ L > 6 0.568(5) -3.3(2) 0.630+ L > 6 0.581(4) -5.0(4) 0.630+ -1.9(3) 1+ L > 6 0.5972+ -6.6(4) 0.630+ -2.7(3) 1+ Table 3.2: Finite size scaling fits of the positions of the peak in Cy to Equa-tion 3.19 or 3.20. The numbers in parentheses represent the statistical un-certainty in the last significant figure (s). A + indicates that the parameter was fixed at the given value. We have also attempted to fit the Cy maximum peak heights to a finite size scaling relation. From Equation 2.60, we have obtained an expression for the finite size scaling of the inverse of the heat capacity peaks, given as NkB C V,max c3(L + 8L)-a^. (3.21) The inverse of the peak height should approach zero as the size of the system increases towards the thermodynamic limit. The advantage of fitting the peak heights in this way is that it leads to a direct estimate of the heat capacity critical exponent a. Recalling that in classical criticality, a = 0, while in the Ising systems, a « 0.11, it should be easier to distinguish whether an exponent is zero or positive, compared with, say, distinguishing whether v = 0.5 or v & 0.63. As in our fitting to the peak positions, we have considered several different fittings of the inverse peak heights. We have attempted fits with and without the correction to scaling term 5L, including or excluding the smallest system size Ljo = 6, and with and without some exponents fixed at their Ising values. The fitted parameters are listed in Table 3.3, and the fits and the raw data are shown in Figure 3.11. The estimates of ajv derived from these fits are close to the Ising value (a/u ~ 0.17), Chapter 3. Phase Transitions in Coulombic Fluids 55 Figure 3.10: Plot of the position of the critical peak in Cy measured for the n = 6 system, as well as the finite size scaling fits. The parameters of the fitted lines are those listed in the given rows of Table 3.2. The MFFSS results are from Ref. [60]. especially when the correction to scaling term is included in the fit, and a fit with a/v fixed at the Ising value describes the data quite well within the uncertainties. A l l of the peak height fits are convincingly not consistent with the mean-field value of a = 0, in which case no predictive relation would exist between the system size and the height of the Cy peak. Chapter 3. Phase Transitions in Coulombic Fluids 56 Range c 3 8L ajv L > 6 0.149(9) ot 0.24(2) L > 6 0.172(6) ot 0.298(17) L > 6 0.1222(98) -3.27(76) 0.189(28) L > 6 0.1174(8) -3.65(12) 0.175+ Table 3.3: Finite size scaling fits of the height of the measured maximum in Cy for a system of particles interacting via a pair potential U(r) oc r - 6 . The numbers in parentheses represent the statistical uncertainty in the last significant figure (s). A t indicates that the parameter was fixed at the given value. 0 0.025 0.05 0.075 0.1 0.125 0.15 (a/L) Figure 3.11: Plot of the maximum value of Cy measured in simulations of particles interacting via a pair potential u(r) oc r~ 6 , as well as the finite size scaling fits. The parameters of the fitted lines are listed in the given rows of Table 3.3. Chapter 3. Phase Transitions in Coulombic Fluids 57 3.3 S i m u l a t i o n s o f I o n i c F l u i d M o d e l s Having shown that the behaviour of Cy measured in the grand canonical ensemble cannot be considered a valid signature of the universality class, it is clear that the peak in the RPM heat capacity detected by Luijten et al. [93] is not a reliable measurement. On the other hand, the lack of a divergence in the canonical calculations of Valleau and Torrie [92], as well as our own calculations [94], is still puzzling, especially when compared with the MFFSS analyses, which show no deviations from Ising criticality. Using these studies as our starting point, we undertook our own investigations into the criticality of ionic fluid models. The interactions of charges and electrical multipoles may be easily charac-terized in terms of the power of the exponent n in the power law u(r) oc r~ n describing the way the interactions decay with interparticle separations r. The interactions of bare charges have n = 1, while charge-dipole interactions have n = 2, or n = 4 once the averaging of the possible orientations of the dipole is taken into account. Two dipoles interact according to a power law with n = 3, which becomes n = 6 when orientational averaging is accounted for [60]. We have seen in Section 3.1.2 that in the dilute phase, the RPM exists al-most entirely as ion pairs, which may be approximated as dipolar dumbbells. In the condensed phase, there are more free ions. Were it not for the possible breakdown of ion pairing near the critical point, there would be no reason to expect anything other than Ising criticality, since the longest ranged interac-tions would have n = 6. Even if ions are present, one would expect that the Debye screening should screen the ionic and ion-dipole interactions enough to make the effective interactions decay exponentially, which would make the criticality Ising-like in accordance with most of the experimental results. However, this may not be the case. It may be possible that if only a few free ions are present, the ion-dipole interactions would not be screened as effec-tively, and these interactions would be sufficiently long-ranged to influence the observed critical behaviour. To resolve the question of whether a breakdown of ion pairing in the RPM near the critical point might be sufficient to influence the universality class of the system (at least as seen in simulations), we undertook a study of a system of charged hard dumbbells (CHD's). These consist of ionic hard spheres as in the RPM, but each ion pair is constrained to remain bound and in contact. In this way, the ion-ion and ion-dipole interactions are explicitly removed, Chapter 3. Phase Transitions in Coulombic Fluids 58 so that all that remain are the dipole-dipole interactions, which with n = 6, are certainly in the Ising universality class. We have done a MFFSS analysis of the CHD system, and calculated the canonical heat capacity. Shelley and Patey have done simulations of CHD's showing that the liquid-vapour coexistence curve is quite similar to that of the RPM [106], but they were unable to determine the critical parameters with high accuracy and also did not measure the heat capacity. Our new CHD results have recently been published [95,107]. For comparison with our CHD simulations, and to check our results against those of other researchers, we also simulated the ordinary ionic RPM, and did our own MFFSS analysis and calculated the canonical heat capacity for larger systems than had been studied previously. We also ran simulations of an asymmetric primitive model, where the cations had a charge twice as large as the anions (2:1 PM). The results of all of these calculations will be described in this section, first dealing with the MFFSS analyses for all three systems, then proceeding to the heat capacity results. 3.3.1 M F F S S Analyses MFFSS analyses of the ionic RPM, the CHD's, and the 2:1 PM were under-taken according to the prescription described in Section 2.3.3. The critical parameters obtained via this approach for all these systems are shown in Table 3.4. The same definitions for the reduced parameters T* and p* (see Equations 3.1 and 2.40) can be used in all three systems, noting that in the 2:1 model we have used the charge of the anion in the definition of T*, not the doubly charged cation. System T; p*c RPM 0.0493(2) 0.079(5) CHD 0.04911(3) 0.101(3) 2:1 PM 0.093(2) 0.096(8) Table 3.4: The critical parameters of ionic fluid models obtained from a MFFSS analysis. The numbers in parentheses represent the uncertainty in the last digit. Chapter 3. Phase Transitions in Coulombic Fluids 59 3.3.1.1 Restricted Primitive Model Because the MFFSS analysis of the R P M had been done previously by several others [76,89-91], our main intent in doing our own MFFSS study was to verify our methods and results against those done by others. Therefore, we did not seek to complete simulations in a broad range of system sizes. The ionic R P M was simulated only at system sizes L/a = 15 and 17. The best fittings of the order parameter distribution to the universal Ising distribution are shown in Figure 3.12. With only two system sizes, we were unable to do a proper analysis of the finite size scaling of the critical parameters. However, since these effects are very small (see, for example, Figure 3.10), a reasonable estimate of the L —> co critical parameters can nevertheless be obtained. Our estimates for the critical parameters (Table 3.4) compare well with previously published results (T* = 0.04917 ± 0.00002, p*c = 0.080 ± 0.005 [76]). - 2 - 1 0 1 2 x = a^L^' v{M - (M)) Figure 3.12: Matching to the Ising universal order parameter distribution for the R P M . Chapter 3. Phase Transitions in Coulombic Fluids 60 3.3.1.2 Charged H a r d Dumbbel ls No MFFSS study of the CHD model has been done previously, and so we endeavoured to do a full MFFSS analysis, including a proper extrapolation of the finite size results to the infinite system limit. The system of CHD's was simulated by GCMC at four different system sizes, L/o = 10, 13, 15 and 17. We used the same procedure described in Section 3.2 for the RPM simulations, except that each cation and anion pair is now constrained to be a dumbbell bound together at contact. Because of this, we did not need to apply any distance biasing in our Monte Carlo procedure in order to achieve efficient acceptance rates. In fact, in a way the CHD simulations can be thought of as a limiting case of the distance biased R P M simulations, where only minimally separated ions are considered for insertion or removal. The spatial translation moves were changed to consist of translations of or rotations about the center of mass of each dumbbell. After allowing for equilibration, the values of the density and the total internal energy were output from the simulation every 500 simulation steps, and stored in histograms for the MFFSS analysis. The best fittings to the Ising distribution for sizes Ljo — 10 and 17 are shown in Figure 3.13. It is worth noting the lack of data on the low density side of the distribution for the Ljo = 10 simulation, which is caused by pronounced finite size effects in the vapour phase. Because the critical density in ionic systems is much lower than in non-ionic systems, the importance of these effects is magnified in the ionic systems. These low density effects are the main source of statistical uncertainty in the determination of the size dependent critical parameters, especially at smaller system sizes. The fitting of the data to the Ising order parameter distribution is necessarily a process subject to uncertainties, since there is no unambiguous way to fit the probability histrograms obtained in the simulations to the Ising form. Various criteria have been discussed, including setting the peak heights to be equal, or similarily equating the area under each half of the distribution [21]. In our analysis, we have attempted to be as systematic as possible by tuning the values of the MFFSS parameters T*(L), (3cpc(L), and s until the minimum of the distribution was centered at x = 0 and the peak heights were as near as possible to the peaks of the Ising distribution. Once the best fit was obtained, the value of p*c(L) for each distribution was obtained directly from the reweighting procedure. The uncertainty in p*(L) was estimated by repeating the reweighting at the upper and lower limits of T*(L), and looking at the effect of this change on Chapter 3. Phase Transitions in Coulombic Fluids 61 the measured critical density. - 2 - 1 0 1 2 x = a^L^M - (M)) Figure 3.13: Matching to the Ising universal order parameter distribution for CHD's. The values of the size dependent critical parameters are listed in Table 3.5. These values were then fit to Equations 2.63 and 2.65 to extrapolate to T*(oo) and p*(oo), respectively. These extrapolations are shown in Figures 3.14 and 3.15. It is apparent from these plots that the finite size effects on the MFFSS procedure are extremely small, which makes the extrapolation to L —> oo a somewhat questionable procedure at best. Regardless of these difficulties, the results for the bulk critical parameters obtained via various fitting procedures are not very different from what would be found from a simple average of the results at all four system sizes. The results we have obtained for the bulk critical parameters of CHD's (Table 3.4) agree well with the results found by other workers who performed simulations on finely discretized versions of various tethered dimer models, where the ions were paired but the separation of the ions was not necessarily fixed to the ion diameter a [108,109]. It also agrees fairly well with the critical point found by Jiang et ai, who used their BIMSA theory to calculate the critical point of CHD's as well as the R P M [75]. They locate the critical Chapter 3. Phase Transitions in Coulombic Fluids 62 Ljo number of configurations in P(p,u) / T I * C Pi 10 4.80 x 10 6 0.04902(5) 0.107(1) 13 8.32 x 10 6 0.04910(5) 0.1050(6) 15 11.52 x 10 6 0.04910(3) 0.1059(3) 17 15.36 x 10 6 0.04907(3) 0.1038(3) Table 3.5: The cri t ical parameters of C H D ' s obtained from fitting G C M C simulation results to the Ising distribution. The numbers in parentheses represent the uncertainty in the last digit. point of C H D ' s at T* = 0.0627, p* = 0.0471. The effect on the crit ical parameters of constraining the ions to remain paired is quite small . The cri t ical temperature changes only very slightly ( « 0.1%), while the crit ical density rises by about 25%. The rise in the crit ical density can be easily rationalized if we consider that the constraint on the interionic separations simply makes the ion pairs in the C H D simulation effectively smaller than the ion pairs in the R P M simulations. This trend is not observed by the work of Jiang et a/., who instead see a rise in the crit ical temperature of about 20% and a lowering of p*c by about 25% compared with their R P M . 3.3.1.3 2:1 Primitive Model Simulations of the 2:1 primitive model were undertaken in response to in-teresting analytical work, which predicted the presence of a lower crit ical solution temperature as well as three phase coexistence in three component ionic systems where, in addition to the singly charged cations and anions, a third component of higher charge is introduced as well [110]. Al though this project was not pursued very far, it st i l l provided the incentive to do simu-lations of the asymmetric 2:1 electrolyte. Camp and Patey did a simulation study of clustering effects in asymmetric electrolytes up to 4:1, where they found that in the dilute phase, the formation of neutral pairs in the R P M is replaced by the formation of neutral clusters [111]. The formation of these larger clusters poses additional challenges to successful G C M C simulation. The biasing technique for creations and deletions described in Section 3.2 had to be altered to accomodate the creation and deletion of neutral triplets, instead of pairs. Because of the additional particle, the acceptance ratios for the insertions and removals were approximately an order of magnitude lower Chapter 3. Phase Transitions in Coulombic Fluids 63 0.0492 0.04915 0.049 0.04895 0.001 0.002 0.003 0.004 Figure 3.14: Extrapolation to L —> oo for T*(L) measured via MFFSS anal-ysis of CHD's than in the RPM, even with the biasing. Because of this, we were unable to extend our calculations to more asymmetric cases. The 2:1 PM system was simulated for a single system size, Ljo — 15, and the best fitting of the order parameter distribution is shown in Figure 3.16. As in the standard RPM case, without completing simulations for more system sizes, no proper analysis of the finite size effects is possible. However, assuming the finite size effects on the measured critical parameters to be as small as they are in the other MFFSS analyses of Coulombic systems [76], a good estimate of the critical parameters (see Table 3.4) is still possible. The critical temperature rises by an approximate factor of 2, which follows simply from the definition of the reduced temperature; the critical temperature of a symmetric 2:2 electrolyte would be exactly four times as high as that of the standard 1:1 RPM. The results of our simulations agree well with a more detailed simulation study of criticality in charge- and size-asymmetric electrolytes [112]. Chapter 3. Phase Transitions in Coulombic Fluids 64 Pc(L) 0.04 Figure 3.15: Extrapolation to L —> oo for p*(L) measured via MFFSS anal-ysis of CHD's 3.3.2 Heat Capacity Measurements We now turn to the measurement of the constant volume heat capacity along the critical isochore for both the ionic RPM and the CHD system. Although we have already discussed a measurement of Cy in simulations of the RPM with box size Ljo = 10 in our discussion of the ensemble dependence of Cy, here we extend our measurements to much larger systems. In agreement with Luijten et al. [93], we measure the heat capacity below the critical temperature, in order to search for the shifted peaks which are the signature of true divergences. However, having demonstrated that Cy measured in the grand canonical ensemble cannot be considered a reliable indication of the universality class in Section 3.2.1, we follow Valleau and Torrie [92] in making our calculations within the canonical ensemble. We have measured the heat capacity of a system of size Ljo = 17 for both the R P M and the CHD model. As in the study of the u(r) oc r~n systems, we have calculated the heat capacity using both the direct fluctuation formula, and as a derivative of the [5,5] Pade approximant fit to the average internal Chapter 3. Phase Transitions in Coulombic Fluids 65 Figure 3.16: Matching to the Ising universal order parameter distribution for the 2:1 PM with L = 15. configurational energy, (U*x). In addition, we have done umbrella sampling calculations in order to achieve optimal accuracy in the energy histograms, as well as to provide additional data points at temperatures between those at which the simulations were actually run. For the RPM, we have also been able to estimate the uncertainties by splitting the data into two parts and calculating the heat capacity using the fluctuation formula separately for each data point in each of the split data sets. The results for the RPM from the direct calculations, without umbrella sampling, are displayed in Figure 3.17, where we also include the Ljo = 10 results from Figure 3.1 for comparison. The measurements utilizing umbrella sampling are shown in Figure 3.18. The corresponding results obtained for the CHD model are shown in Figure 3.19. Though there is a significant amount of noise in the heat capacity results, even with the use of umbrella sampling techniques, the consistency of the results for the separate data sets and with and without umbrella sampling, as well as the Pade fits, indicate that these heat capacity results are quite reliable. The results for the RPM show that there is a very broad peak in Chapter 3. Phase Transitions in Coulombic Fluids 66 8 Flue, set 1 Flue., set 2 Flue., total Pade, total L/a = 10 0 • * 5 • • • -0.3 -0.2 -0.1 0 0.1 0.2 t Figure 3.17: Cy calculated for the RPM system in the canonical ensemble without umbrella sampling, L/a = 17, p* = 0.079. "Flue." denotes Cy measured by the fluctuation formula. "Pade" denotes Cy measured as the derivative of the Pade approximant fit to (U*x). For comparison, we have included our results for Ljo = 10. All results include the translational contribution to the heat capacity (Cyjrans/NkB = 3/2) Cy(T) measured below Tc. To verify whether or not this peak is a very weak signature of a critical divergence or not, we would require very accurate data over a broad range of system sizes. Although our results for L/a = 10 are lower than the L/a = 17 system, and may be consistent with what one would expect of a critical divergence, we stress that the computational challenges presented by the ionic models make the kind of systematic study we made on the n = 6 system much more difficult on the RPM, and we have not made such an attempt. We note, however, that heat capacity measurements by Caillol et. al. [76] in the RPM up to L = 34 also do not show the behaviour one would expect of a divergent heat capacity, showing a narrowing of the peaks as the system size is increased, but no systematic rise in the peak height or approach of the peak position towards T*(oo). Our calculations of Cy in the CHD system do not extend as far down Chapter 3. Phase Transitions in Coulombic Fluids 67 8 7.5 7 Cj 6.5 6 5.5 SP A 2? o • • F lue , set 1 0 Flue , set 2 A Flue., total — Pade, total - 0 . 3 - 0 . 2 5 - 0 . 2 - 0 .15 - 0 . 1 - 0 . 0 5 0 t Figure 3.18: Cy calculated for the R P M system in the canonical ensemble with umbrella sampling, Ljo = 17, p = 0.079. "Flue." denotes Cy measured by the fluctuation formula. "Pade" denotes Cy measured as the derivative of the Pade approximant fit to (U*x). A l l results include the translational contribution to the heat capacity (Cy^rans/Nks = 3/2) in temperature as our R P M calculations, and so we cannot confirm whether the same broad peak occurs in the CHD fluid. When one accounts for the differing kinetic contributions to Cy in the models, the heat capacities appear nearly identical in the temperature range we have simulated. The contrast between the C H D result and the results shown earlier for the n = 6 system (see Figure 3.9) is striking, since the longest range effective interactions in the CHD's must be n = 6 as well. There is no sign of a peak in the CHD calculations anywhere near the location of the peak in the n = 6 calculations. Since the C H D model must be an Ising system, it seems that the lack of a peak in Cy in this model, at least, must be due to reasons other than the universality class of the system. Chapter 3. Phase Transitions in Coulombic Fluids 68 8.5 - 0 . 2 -0 .15 - 0 . 1 - 0 . 0 5 t Figure 3.19: Cy calculated for the CHD system in the canonical ensemble, L/a = 17, p* = 0.1034. The dagger (t) indicates those results obtained with umbrella sampling. "Flue." denotes Cy measured by the fluctuation formula. "Pade" denotes Cy measured as the derivative of the Pade approximant fit to (U*x). All results include the translational and rotational contributions to the heat capacity (Cyitrans+rot/NkB = 5/2) 3.4 Discussion The discussion of the results of our investigations can be brought into better focus by posing two questions. First, why are the results of calculations of the heat capacity so radically different in calculations with different statistical ensembles? Second, if the models for ionic fluids are in fact Ising systems, why is the critical peak in Cy so difficult to detect? 3.4.1 Ensemble Dependence of Cy We would first like to stress that this ensemble discrepancy is emphatically not a numerical problem due only to a lack of precision or a lack of data. This should be clear from an inspection of Figure 3.1, which shows that the Chapter 3. Phase Transitions in Coulombic Fluids 69 results obtained earlier by Luijten et al. [93] are completely reproducible. Our results for Cy calculated in the grand canonical ensemble for the power law potentials are similarily reproducible. The observed discrepancy is sys-tematic, having to do with the finite size effects which are inherent in the results of Monte Carlo simulations. It is usually taken for granted that all statistical mechanical ensembles are equivalent in the thermodynamic limit. For example, if one takes a grand canonical ensemble of systems with fluctuating particle numbers N, the size of the density fluctuations will vary as 1/y/N, and if N ~ 102 3 or so, then there can be no difference between a calculation in this grand canonical ensemble and a canonical ensemble calculation with N fixed, since the density fluctuations are insignificant in the macroscopic system. Any quantity determined as an average over the states in the ensemble will have the same value in the thermodynamic limit, regardless of the constraints on the ensemble. In considering computer simulations of criticality, this simple conclusion is invalidated by the fact that, as we have mentioned elsewhere in this the-sis, computer simulations must be done on finite systems with rather small numbers of particles compared to a real, bulk system one might study in the laboratory. Although techniques like periodic boundaries can be applied to mitigate against these finite size effects, they can never be completely removed. One method of approaching simulations of the critical region is to rec-ognize that no single simulation can ever be considered an accurate repre-sentation of critical behaviour, but if the simulation is repeated for several different system sizes, then the scaling behaviour of the response functions and other quantities as the system size changes can be used to extrapolate towards the thermodynamic limit. This is the basis ofthe MFFSS techniques of Bruce and Wilding, and also lies behind much of the work of Fisher and his coworkers. However, we believe that our results demonstrate that the scaling hypotheses which underlie much of the work on simulations of the critical region must be carefully re-examined. In our case, we have studied the behaviour of the constant volume heat capacity, Cy, near the condensation critical point of several model fluids. We recall that Cy can be computed from the fluctuation formulae appropriate Chapter 3. Phase Transitions in Coulombic Fluids 70 for the statistical ensemble, i.e. Cv,NVT — X(U,U) kBT2 (3.22) in the canonical ensemble, and kBT2 1 X(U,U) X(N, U)2 X(N, N) (3.23) in the grand canonical ensemble. In the canonical ensemble, there is a single term in the equation, X(U, U), which is the mean-square fluctuation in the internal energy U. This quantity diverges as the system size increases, but weakly, with an exponent a/is, in agreement with Equation 2.57. In the grand canonical ensemble, however, there are three separate terms in the expression, and all of them diverge strongly. For example, in the thermodynamic limit [113], where as we have seen, KT is the isothermal compressibility, and in three dimensions it diverges as T —> Tc with an exponent 7 = 1 in classical systems or 7 = 1.24 in Ising systems. In a finite system the divergence of KT should be replaced by a peak which increases in height as the system size increases, such that Since the exponent 7/V is much larger than the a/v exponent for the heat capacity, we should expect this peak to be much more pronounced than the peak in Cv- X(U,U) must also diverge more strongly compared to the canonical case, since in the grand canonical ensemble the energy fluctuations will be very strongly tied to the density fluctuations. This explains how a sharp peak might appear in the n = 3.1 grand canonical calculation, despite the fact that X(U,U) ~ 0 in the canonical case. In the thermodynamic limit, the three terms in Equation 3.23 must combine in some particular way such that they become equal to the value of X(U, U) in the canonical ensemble. This should occur for sufficiently large systems, where the density fluctuations do not extend over the entire simulation box. Another way of looking at the situation is to think about the physical picture of what is happening in a Monte Carlo simulation. In a canonical X(N, N) (N2)kBTKT V (3.24) ^T,max oc L?lv. (3.25) Chapter 3. Phase Transitions in Coulombic Fluids 71 ensemble, one takes a fixed number of particles N at temperature T confined to a volume V, and allows them to attain an equilibrium in some thermody-namic state. If this simulation is attempted inside the coexistence region of the phase diagram, as all of our simulations along the critical isochore with T < Tc are, then the results cannot be considered an accurate representation of reality unless the system is large enough to contain domains consisting of each phase, as would be seen in an experiment. This can only be achieved with very large systems, well out of the reach of contemporary simulations of three dimensional continuum models. In our small simulations, only the barest vestiges of the critical energy fluctuations remain. This explains why the finite size effects on the location and size of the heat capacity peak in the n = 6 system (see Figures 3.10 and 3.11) are so large. By contrast, in the grand canonical simulations, the particle number fluc-tuates while the chemical potential p is kept constant. If the simulation is run on the critical isochore, where p = pc and (N)/V = pc, then in the two-phase region, even a relatively small system can flucutate strongly in density and both phases will be sampled, provided the simulation runs long enough. It is for this reason that the grand canonical ensemble is often said to be a more natural ensemble to use to study phase transitions where the density is different in each phase [93]. Unfortunately, our results suggest that the standard analysis of the grand canonical simulations, which inherently assumes that the fluctuation formulae such as Equation 3.23 apply far from the thermodynamic limit, erroneously magnifies the true critical fluctuations. Of course, in a very large grand canonical simulation, the fluctuations in N would be insignificant compared to the system size, both phases would coexist in a single simulation box as they would in a large canonical simulation, and there would be no difference between the canonical and the grand canonical results. Although we have focused on the calculation of Cv in grand canonical simulations near the critical point, there is no reason to believe that our concerns are not applicable to the calculation of other response functions via simulation, since like CV, these quantities can all be related to second and higher moments of the probability distribution function p(N, U). For exam-ple, one quantity often examined to provide information on the universality class is the dimensionless fourth-order magnetization cumulant ratio Q(T, L) Chapter 3. Phase Transitions in Coulombic Fluids 72 (closely related to the Binder cumulant ratio [114]), defined as Q{T,L) <A4,L> ' (3.26) This function takes on different forms depending on the universality class, and by measuring the ratio over a range of different temperatures and system sizes, it can be used to distinguish different critical behaviours. This analysis has been used successfully to confirm the expected universality class of several models, including spin models with algebraically decaying interactions [115— 117], as well as in work in our group on fluid models [60]. However, in the light of the current study, it is unclear whether the finite size effects, which are considered in the derivation of the different, universality class dependent functional forms for Q(T, L), are being dealt with properly. Since (M) = p — su = (N)/V — s(U)/V, the presence of strong divergences in the higher moments of the probability distribution might render the accepted expressions for Q(T, L) invalid for finite systems. The same critique can be applied to much of the approach used by Fisher et al. to the determination of universality classes via simulation. In this approach, many higher moments of the probability distribution are collected and their behaviour analyzed to arrive at unbiased estimates of the critical exponents, in contrast to the Bruce-Wilding MFFSS analysis, which relies on the assumption of a form for the order parameter distribution which depends on the universality class [31,35]. Unfortunately, these higher moments might be subject to finite size effects in the same way that the grand canonical heat capacity calculations are. Given these concerns, it is perhaps surprising that the results obtained by the Fisher methods are quite good. In one study, they were able to confirm that the square well fluid has critical exponents ft and 7 equal to the Ising values, and definitely excluded the Heisenberg and self-avoiding walk universality classes, which have very similar values of the critical exponents [21]. Our speculation is that in the case of quantities which depend on the larger critical exponents, such as ft and 7, the rate of convergence of the finite size results is much more rapid than in the case of the heat capacity, with its weak exponent a. This is also consistent with the fact that the MFFSS analysis is less susceptible to finite size effects, since the order parameter distribution scales with an exponent ft/v. Another factor that improves the Fisher analysis is the use of higher order moments which have maxima above Chapter 3. Phase Transitions in Coulombic Fluids 73 the two-phase region, where T > T c . This would mitigate against there being large finite size effects. There is one aspect of the Fisher group's work which we strongly feel is suspect, due to its connections to heat capacity calculations. This is the work on the detection of the Yang-Yang anomaly. The Yang-Yang anomaly was first hypothesized by Yang and Yang [36], who noted that starting from the Gibbs-Duhem relation, This relation implies that if the heat capacity diverges at the critical point, then either or both of the terms on the right hand side must diverge. While it is well known that the pressure derivative diverges, the behaviour of the chemical potential term is less certain. On the one hand, theories of lattice gases all feature a chemical potential which remains analytic through the critical point, which would prevent a diverging second derivative [118,119]. On the other hand, Yang and Yang suggested that real gases should have diverging contributions from both the pressure and chemical potential second derivatives. The strength of the divergence of (d2p/dT2)v has been dubbed the "Yang-Yang anomaly". One of the consequences of a Yang-Yang anomaly is that the formulation ofthe scaling fields in the Bruce-Wilding method must be extended to include a pressure term [35]. Fisher and his collaborators have undertaken several studies of both ex-perimental data [38] and simulations of fluid models, including the hard-core square well and the R P M [21,24], which indicate that there is a significant Yang-Yang anomaly in real fluids and simulations of continuum fluid models. While we do not wish to question the existence in principle of a Yang-Yang anomaly, or the necessity of including a pressure term in the formulation of the scaling fields, all of these simulation studies rely on calculations of quan-tities related to the constant-volume heat capacity, and they are all done in the grand canonical ensemble. This fact makes the conclusion [21] that there exists a large Yang-Yang anomaly premature, at best. Until a simulation in the canonical ensemble can be done which demonstrates the anomaly, no cer-tain determination of the size or existence of the Yang-Yang anomaly can be SdT - VdP + Ndp = 0, (3.27) it is easy to arrive at the following relation for the heat capacity: (3.28) Chapter 3. Phase Transitions in Coulombic Fluids 74 made. There are also many ambiguities in the experimental situation [120], which must be dealt with before firm conclusions can be reached. We believe that there is a need for a more careful examination of the Fisher method for detemining universality class. A l l of the analyses made by Fisher and coworkers have been made on models which belong to the Ising universality class. We would like to see a Fisher-type analysis of a model, like the n = 3.1 model, which displays classical criticality. If this case could be analyzed in the same way as the R P M and the Ising model, and the results were consistent with classical criticality, we would be much more confident in their techniques. 3.4.2 The Critical Peak in Cv for Ionic Fluids There is no doubt that Cv diverges near the liquid-vapour critical point of both real gases, and lattice gas models and Ising models. This divergence has been well established experimentally for over 40 years [121], and in canonical simulations of "Lennard-Jonesium" [92]. In one and two dimensions, while no true Cv divergence exists since a = 0, the critical behavior of Cv has been well studied, in analytical work on exactly soluble square lattice Ising models [122], as well as in canonical ensemble simulations of adsorbed ni-trogen gas on sodium chloride [123]. These results are fully consistent with the expectations, i.e., a logarithmic divergence can be detected which gets sharper with system size. However, although there is no other evidence to suggest that ionic fluids are not members of the Ising universality class, no canonical ensemble simulation suggests that there is a critical peak in the heat capacity. Although it is still possible that ionic fluids are not members of the Ising universality class, we believe that the evidence from Bruce-Wilding MFFSS analyses and the Fisher group's analyses strongly suggests that they are indeed Ising systems. Far more likely, in our opinion, is that accurate calcu-lation of the heat capacity via simulation simply requires much larger system sizes than have been studied to date. This is consistent with our heat capac-ity calculations in the n — 6 model, which only approach the infinite system results very slowly, and with far larger finite-size effects than in the MFFSS analysis. In the R P M and other ionic fluids, detection and finite size scaling analysis of the critical divergence in Cv, if it exists, would seem to require system sizes L/o » 34, well out of the range of current simulation work, but perhaps becoming within reach in a decade or so. Chapter 3. Phase Transitions in Coulombic Fluids 75 Finally, we note that the original question we sought to answer with the CHD simulations, namely, whether the critical behaviour of ionic fluids might depend on the breakdown of ion pairing near T c , remains only partially an-swered. Although we have not detected any difference between the criticality of the R P M and the CHD's, it is entirely possible that any differences in Cy would only begin to appear in simulations on much larger systems. 76 Chapter 4 T h e A d s o r p t i o n of CO2 Gas on M g O Crysta ls Adsorption refers to a process whereby particles of one species become trapped on the surface of another species. The adsorption of gases is a problem of relevance to many industrial processes, especially those involving catalysis, which often include reactions between gases (or other molecular substances) and the surfaces of condensed phases. Adsorption processes may be divided into two main classes, distinguished by the forces involved in the interaction between the adsorbate and the substrate. In chemisorption, the adsorbate forms a chemical bond with the substrate. The electronic configurations of the adsorbate and the substrate will undergo large changes accompanying the formation of the bond, and there may be reconstruction on the first one or few atomic layers of the substrate. In physisorption, on the other hand, the adsorbate and the substrate interact via weaker forces, such as van der Waals interactions, which do not change the electronic structure. Chemisorption has been studied extensively in the laboratory, and is of more interest from a practical viewpoint, since it is implicated in many cat-alytic processes of industrial importance. Much less is known about ph-ysisorption, due to its lesser practical impact as well as the experimental difficulties in studying weakly interacting systems, which include a need for ultra-high vacuum (UHV), low temperature, and very clean substrate sur-faces. On the other hand, physisorption is much easier to study theoretically than is chemisorption. This is due to the lack of chemical bond formation, and surface reconstruction is less of a concern due to the weakness of the van der Waals interactions compared to the covalent or ionic interactions within the surface. These two factors allow for the use of approximations which greatly simplify the simulation and modelling of the physisorbed system. In this chapter, an attempt to model a particular system of a gas ad-sorbed on to a crystal surface, that of C 0 2 adsorbed on the MgO (100) surface, via G C M C simulations shall be described. Starting with a brief Chapter 4. The Adsorption of C02 Gas on MgO Crystals 77 survey of analytical attempts to understand adsorption processes in general, we will continue with a brief survey of the experimental and theoretical re-sults on the CO^/MgO system which motivated our simulation study. This will be followed with a description of the potential model and the simulation techniques employed, as well as the means whereby the simulations may be compared with experimental data. Then the results of the simulations, both in the canonical and the grand canonical ensemble, will be shown, and finally, some discussion and conclusions and ideas for future work will be presented. 4.1 Survey of Past Work 4.1.1 Analytical Treatments of Adsorption and Surface Phases The qualitative features describing the physisorption of gas molecules by a crystal surface may best be understood by first considering the gas on its own. Because of the weakness of the van der Waals interactions between the gas molecules, solidification takes place at quite low temperatures. As long as the temperature is less than the triple point temperature Tt, then there is no distinct liquid phase, and the substance undergoes sublimation from gas to solid at some pressure, likely well below atmospheric pressure. When the crystal surface is introduced to this low temperature, low pres-sure gas, as long as the van der Waals attraction between the gas and the surface is stronger than the self-interaction of the gas itself, the gas molecules will interact with and adsorb on the surface. In fact, at sufficiently low pres-sure and temperature adsorption will take place on the surface regardless of the relative strength of the gas-gas and gas-surface interactions. There will likely be a narrow range of pressure (below the pressure where the gas simply solidifies or liquefies) where a monolayer of gas atoms forms on the crystal surface. As the pressure rises, further layers may be adsorbed, and finally, a bulk layer of solid (or liquid, if T > Tt) will form on the surface. In a theoretical sense, the monolayer formed on the crystal surface may best be described as a quasi-two dimensional, surface phase of some sort. If the lateral interactions between the gas molecules in the monolayer are quite weak compared to the interaction with the surface, this is best thought of as a two dimensional lattice gas. If there are large lateral interactions, though, then it is more like a two dimensional liquid or crystal. The details Chapter 4. The Adsorption of C02 Gas on MgO Crystals 78 of the geometry of the adsorbed monolayer depend on the relative size and lattice spacings of the gas and substrate molecules, and on the details of their interactions. If the gas and the substrate molecules have similar sizes, and similar equilibrium geometries, the system is described as epitaxial, and will generally be easier to understand, since there is a one-to-one correspon-dence between the substrate lattice and the available binding sites. In a non-epitaxial system, the mismatch between the preferred geometries of the gas and the substrate can lead to much more complicated and interesting behaviour on the surface. More generally, one may have to distinguish be-tween commensurate phases, where the geometry of the adsorbed monolayer is based in some simple way on the geometry of the substrate, and incom-mensurate phases, where there is no simple relation between the geometry of the monolayer and that of the substrate. We have already seen how the use of Ising models of spin systems is rel-evant to understanding fluid criticality, once it is recognized that an order parameter can be defined in both systems which will behave the same way near the critical point. In the study of physisorption in systems with com-peting interactions, similar insights have been furnished by the study of Ising systems with anisotropic interactions. These systems show the same sorts of complex behaviour seen in physisorbed systems, such as the phenomenon of the "devil's staircase", where the competing interactions can lead to a sequence of extremely closely separated phases of increasingly longer period-icities. An extensive review of these theories and models has been prepared by Bak [124]. Another useful survey of the possible structures which can form on various substrate lattices is the pair of articles by Domany et al. [125,126]. 4.1.2 Experimental and Theoretical Work on the Adsorption of C O 2 on MgO The motivation for our work on modelling the adsorption of CO2 gas on MgO lies in several interesting experiments done in the last two decades on the system. Chief among these is the work of Suzanne et al, who made low en-ergy electron diffraction (LEED) studies of C 0 2 adsorbed on MgO [127,128]. In their study, they cleave MgO crystals in UHV (below 10 - 1 0 Torr) before introducing a low pressure (4x 10 - 9 Torr) of C 0 2 . By measuring the decrease in the intensity of the MgO LEED pattern, they observe the formation of a monolayer after several minutes of exposure to the gas at 60 and 80 K. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 79 They also observe the replacement of the MgO L E E D pattern with a new pattern, and they attribute this to the formation of an ordered monolayer of C02- This monolayer has a (2\/2 x y/2)R45° unit cell, commensurate with the MgO (100) surface, with two perpendicular glide planes, indicating that the CO2 molecules lie flat on the surface. They propose that the monolayer consists of flat C 0 2 molecules, with each oxygen atom bound equally to a M g 2 + binding site, in a "herringbone" pattern matching the observed unit cell. Further, they observe that the L E E D pattern disappears when the tem-perature is raised higher than 93 K, indicating that the ordered monolayer becomes disordered at higher temperatures. The pattern does not reap-pear when the temperature is lowered back down below 93 K . Based on the L E E D pattern, they calculate the areal density (i.e. the surface area occu-pied by one adsorbed molecule) of the ordered CO2 monolayer to be 17.7 A2/molecule, larger than the areal density of the densest (111) plane of bulk CO2, that is 13.45 A2/molecule [129], indicating that the C 0 2 molecules in the monolayer are quite widely spaced. The adsorbed monolayer is denser in both the (2\/2 x y/2)R45° C 0 2 monolayer which forms on graphite, with an areal density calculated to be about 15 A2/molecule [130], and in the (2 x 1) monolayer of C 0 2 on NaCl, which has a measured areal density of 15.85 A2/molecule [131]. Suzanne et al. also report on volumetric isotherm measurements at 156 K of C 0 2 adsorbed on MgO powder done as part of a thesis [132], citing an areal density of 11.9 A2/molecule. These results suggest the intriguing possibility that at temperatures above 93 K , the ordered commensurate structure is replaced by a disordered struc-ture with a higher density. This is a very unusual situation, since denser structures are usually more ordered due to the restrictions on motion caused by the close packing. One well-known exception to this tendency is the be-haviour of liquid water, with a maximum density of 999.97 kg m - 3 , which is significantly denser than ice (917 kg m~3) [133]. In water, the breaking of the hydrogen bonds which form in solid ice allows a slightly denser dis-ordered liquid structure to form. Something similar could occur in the CO2 monolayer, whereby the flat CO2 molecules might, upon heating, tip up off the surface and rotate by 45° so that the carbon atom can interact favorably with the O 2 - anion, leaving additional binding sites which could allow the formation of a higher density monolayer. There are other experimental and theoretical studies on the C 0 2 / M g O system which furnish additional information. Investigations by Heidberg Chapter 4. The Adsorption of C02 Gas on MgO Crystals 80 et al. [134-136] using Fourier transform infrared spectroscopy (FTIR) show three separate infrared absorbtion bands for the adsorbed monolayer. Fur-ther analysis of the F T I R results shows that the CO2 monolayer consists of two sets of translationally inequivalent molecules, which are bound to en-ergetically equivalent binding sites, supporting the herringbone structure of the L E E D study. However, analysis of the FTIR spectra taken with s and p polarized light shows that, rather than lying flat, the C 0 2 molecules tilt off the surface with an angle 9 — 27° ± 10°. There is a single glide plane, on the short side of the unit cell, rather than the two glide planes predicted by the L E E D results. This would suggest that, instead of balancing between two binding sites, the C 0 2 molecules choose a preferred binding site, and one oxygen atom will be more tightly bound to a M g 2 + site than the other, even if it is still blocking further gas from adsorbing. This makes the sce-nario suggested above for the formation of the higher density layer even more reasonable. Another interesting piece of experimental data on the system is from laser induced thermal desorption (LITD). This experiment provides information on the binding energies of the adsorbed molecules. The energy required to desorb the molecules varied from 9.4 kcal/mol at low coverage, down to 7.1 kcal/mol at full coverage, in this case corresponding to a monolayer with only half the density ofthe (2\/2 x y/2)R45° monolayer [137]. This conflict in the reported density of the monolayer makes the LITD results questionable, as does the finding that the binding energy of the complete monolayer is less than that of the single molecule. This would indicate that the lateral interactions between adsorbed gas atoms were unfavorable, hence lowering the total binding energy. One possible explanation for the large low coverage result is that this is the binding energy of C 0 2 molecules which adsorb on the defects between steps of the MgO lattice. These molecules would be bound closer to the crystal surface, thus having a larger binding energy [138]. The volumetric isotherm measurements mentioned above [132] also give an estimate of the monolayer binding energy, in this case 11.0 kcal/mol. One final study of adsorbed monolayers of C 0 2 on MgO deserves atten-tion. In a theoretical study, Girardet et al. [138] calculated binding energies using a model similar to the one we develop in the next section based on pair-wise atom-atom potentials. The Girardet study has a Lennard-Jones 6-12 potential acting between all of the gas atoms, and between the gas atoms and the substrate atoms (the interactions within the substrate itself are ignored). The electrostatics are dealt with by placing a point quadrupole at the center Chapter 4. The Adsorption of C02 Gas on MgO Crystals 81 of each CO2 molecule which interacts with the charges in the substrate. The results for a single adsorbed gas molecule have it lying flat, balanced between two M g 2 + binding sites, as in the LEED investigations. The binding energy of a single molecule at T = 0 was found to be 7.9 kcal/mol. The binding energy of the (2\/2 x y/2)RAb° monolayer was calculated to be 10.9 kcal/mol. This structure is found to be much more stable than either the (\/2 x \/2)RA§° phase with all C 0 2 molecules parallel, or the (2\/2 x 2y/2)R4b°, "cross-hatched" phase. Snapshots of the (2\/2 x \/2)R4h° herringbone structure and the (\/2 x y/2)R4h° parallel structure are shown together in Figure 4.1. Figure 4.1: The (2^2 x v/2)#45° and {y/2 x ^)R45° structures of C 0 2 on MgO. These snapshots are from canonical simulations at T = 5 K on a square MgO lattice of size L = 8a These results show some consistency in the observation of the herringbone pattern in the adsorbed monolayer. However, the observed ti l t in the FTIR study and the curious LITD results suggest that more work must be done to explain what is occurring in this system. The possibility of observing an unusual disordered phase with a higher density than that of the ordered monolayer is also tantalizing. These are the chief motivations behind our own simulation study. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 82 4.2 Model and Simulation Methods The general concepts underlying the construction of the potential models used in this study have already been described in Section 2.1, so here we give only the details specific to the current calculation of the energy of CO2 interacting with a MgO (100) surface. In contrast with the work described in Chapter 3, in this investigation we cannot make use of a simple model, because such a model would be unable to adequately capture the details of the interactions. For example, modelling the problem as a collection of point quadrupoles (representing the C02's) bound to a two-dimensional lattice has been done [139], but in our study, modelling the electrostatics as a single molecular quadrupole is inadequate for describing the interactions with the very strong surface potential, since the surface field felt by one end of the C 0 2 molecule is significantly different from that experienced by the other end. The formulation of the accurate interaction potential between an adsor-bate and a crystal surface is a rather complicated problem, since it requires the input of a large amount of experimental and calculated data. A complete description of the problem appears in References [140] and [141]. In addi-tion, it should be noted that the model parameters used in this study were verified and provided by Dr. Abdul-Wahab Sallabi of Dr. David Jack's re-search group. For these reasons, only the most relevant aspects of the model specification will be described. The geometry of C 0 2 is linear, with the carbon between the two oxygen atoms. The CO bond-length is 1.163 A [129]. The MgO crystal has a face-centered-cubic (fee) structure, with a lattice constant of 4.20 A, and the (100) face has a (2 x 2) unit cell, corresponding to a like-ion separation on the (100) face of a = 2.97 A [127]. It is helpful to separate the problem into the interactions between the gas molecules themselves, and between the gas and the surface. The interactions within the surface itself are ignored. As already discussed in Section 2.1, the core repulsion and attraction terms can be separated from the purely electrostatic terms. The electrostatic contributions to the interaction energy are computed after deciding on a particular model for the charges and electric multipoles within CO2 and MgO, and matching this model with available experimental data and results of ab initio calculations. In CO2, the molecular symmetry is such that the leading terms in the multipole expansion are the quadrupole Chapter 4. The Adsorption of C02 Gas on MgO Crystals 83 moment Qzz and the hexadecapole moment $ z z z z (all other components of these tensors can either be expressed in terms of the qiven quantities or are zero due to the molecular symmetry). As explained above, simply placing a quadrupole in the center of the carbon atom is inadequate for our investi-gation, owing to the large field gradients above the MgO surface. We have chosen to model the CO2 molecules as a collection of three charges centered on each of the three atoms, as well as two point dipoles of equal magnitude but opposite sign, centered on the oxygen atoms, and both pointing towards the central carbon atom. The values of the charges and dipoles are chosen to reproduce the experimentally determined quadrupole moment, Qzz = —4.3 Debye A [142], and the calculated hexadecapole moment, $ z z z z = —1.68 De-bye A 3 [143]. These values are listed in Table 4.1. The polarizabilities are also listed in this table, but since they are determined in conjunction with the dispersion forces, their derivation will be dealt with below. Atom Q M a \ \ ot± C 5.4402 0 1.6123 0.7763 0 -2.7201 0.65741 1.2186 0.5868 Table 4.1: The parameters in the electrostatic model for CO2. Units are as follows: q, electron charges e = 1.602 x 1 0 - 1 9 Coulombs; /J, Debye; ay and QX, A 3 The calculation of the electrostatic contribution to the binding energy within the CO2 gas, then, requires only the summing up of all the interac-tions of the charges and dipoles in the adsorbed molecules according to the formulae of Section 2.1. For the gas-surface interactions, we can take advan-tage of the symmetry of the crystal lattice in order to make our calculation as efficient as possible. The electrostatic potential existing above the (100) face of an fee ionic crystal can be written as [144] [cos(27nr/a) + cos(27ry/a)], (4.1) where a = 2.97 A is the like-ion separation on the MgO lattice and qs is the ionic charge (qs = 2e for MgO). The coordinate system of the vector (x, y, z) denoting the position above the surface is centered on the site of a surface M g 2 + cation, with the z axis perpendicular to the surface and the x and y axes pointing towards neighbouring cations. Then, the surface electric field 1 / x 4<7S $(x,y,z) = — a 1 + e - ^ Chapter 4. The Adsorption of C02 Gas on MgO Crystals 84 created by this electric potential can be computed by E(x,y,z) = -V*{x,y,z), (4.2) and the energy of interaction of a CO2 molecule with this surface field can be calculated as 3 ^coul = ^ (EXi(x + EyXy + EZi(2)}, (4.3) i=i where the index i loops over the three atoms in each molecule and £ = [Cx, Cy, C] is the orientation vector of the CO2 molecule. Another electrostatic effect we consider is the polarization of the gas molecules by the surface field (the self-polarization of the gas caused by the field generated by the gas itself is not included in our calculation). The polarizabilities of each of the atoms both parallel, ct\\, and perpendicular, a±, to the C 0 2 molecular axis are calculated from knowledge of the molecular polarizability of the CO2 molecule. We use experimental estimates of the molecular polarizabilities of a>\\ = 4.05 A3 and a± = 1.95 A3 [9]. The atomic polarizabilities are apportioned among the atoms by optimizing the structure and binding energy of bulk CO2 using the computer program W M I N [145]. This was done previously for another project by members of Jack's group [146], and we will simply use their results. These values are listed with the rest of the electrostatic parameters in Table 4.1. The contribution to the energy due to the polarization of a single CO2 molecule can then be expressed as 3 tfpoi = - E K ^ i + El + El) + (am - «L,i)(s*c. + + EZiQ2], 1=1 (4.4) where the symbols have the same meaning as in Equation 4.3. For the core repusion and dispersion terms, we use either the Buckingham or the Tang-Toennies potential. For the gas-gas interactions, the Buckingham potential is used, and we use only CQ and Cs terms in the dispersion series. Including C10 or higher terms leads to difficulties unless the Tang-Toennies form of the potential including damping is used. The repulsion parameters for the gas-gas interactions were determined along with the polarizabilities by Hu et al. [146] using W M I N [145]. The values of the gas-gas repulsion and dispersion parameters are listed in Table 4.2. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 85 Interaction A a C76 C8 C-C C-O 0-0 19959.052 3.094059 438.7693 3334.3695 19784.555 3.287635 331.7595 2519.9482 19670.626 3.507049 250.8479 1904.4496 Table 4.2: The potential parameters for the repusion and dispersion interac-tions between CO2 molecules. Units are as follows: A, kcal/mol; a, A - 1 ; C$, A 6kcal/mol; C 8 , A 8kcal/mol. The calculation of the gas-surface interactions can be completed more efficiently by taking advantage of the symmetry of the crystal lattice. Instead of summing up the interactions between the gas molecules and the surface atoms individually each time the energy is calculated, by expressing the surface potential as a two-dimensional Fourier series, the interactions can be computed for different heights above the surface just once and then stored in a lookup table. The general procedure is outlined in the book by Steele [147]. In our case, the potential parameters appropriate for the calculation of the interaction between CO and MgO were already known from a previous project done by members of the Jack group [123]. This potential was com-puted by summing up the interactions of successively deeper layers of the MgO crystal with CO until the energies converged, and were tabulated for heights above the MgO surface 0 < z < 8.0 A and separated by Az — 0.001 A. Between successive values of z, the energy was calculated as a linear in-terpolation between the results at the nearest two values of z. Of course, the atom-atom potentials for calculating the interaction of CO2 with MgO will differ slightly from those which describe the C O / M g O interactions, and ideally, we would have repeated the gas-surface potential calculation for this system. Unfortunately, these calculations are lengthy and complicated, re-quiring mainframe computing power to sum over many layers of the MgO crystal. We bypass the problem by assuming that all C02/MgO interactions below the first surface layer cannot be distinguished from the interactions of the same atoms when they are found in the C O / M g O system. Since the most important differences will be in the interactions with the topmost layer, we compute these interactions directly by summing over all the gas and surface atoms. In the case of the gas-surface potential, we extend the dispersion series out to Cw and use the Tang-Toennies potential which includes the damping functions on the dispersion series. The potential parameters for the repulsion and dispersion interactions for all atom-atom pairs under con-Chapter 4. The Adsorption of C02 Gas on MgO Crystals 86 sideration are listed in Table 4.3. They are taken from both ab initio and experimental data [148]. The repulsion parameters are subject to a certain amount of adjustment in order that reasonable values of the binding energy can be obtained. Interaction A a c 6 C8 Cio C - M g 2 + c - o 2 -0 - M g 2 + o-o 2 -133838.56 3229.529 161762.92 2865.57 4.68165 2.42954 5.13953 2.54731 45.9063 585.3665 34.6978 442.4436 • 241.6571 5965.3749 182.6542 4508.8704 1558.352 74470.396 1177.865 56287.722 Table 4.3: The potential parameters for the repusion and dispersion inter-actions between CO2 molecules and the MgO surface. Units are identical to those used in Table 4.2. The units of Cw are A 1 0kcal/mol. Because we only have used the proper, Tang-Toennies form for the gas-surface core interactions, our other pair potential models break down for small interparticle separations. In fact, even using the proper Tang-Toennies potential including the damping of the attractive interactions, the electro-static interactions also must be damped, or they too will dominate the Born-Mayer term as the particle separation r —>• 0. These breakdowns are especially troublesome for G C M C simulations, since a newly created CO2 molecule might randomly be placed at the same location as a previously ex-isting molecule. To avoid all of these problems, we have simply introduced a small hard-sphere radius r#s = 1.5 A and no molecules may be created or moved such that any interparticle separation r < rHs results. This dis-tance is significantly smaller than any interparticle separation which might be expected to occur naturally in this system. Now that the binding energy can be calculated to our satisfaction, the rest of the Monte Carlo procedure can be built around it. We utilize periodic boundary conditions in the x and y directions, with each side of the.periodic box having an equal length L. Different simulations can be run with different lattice sizes, where L must always be a multiple of the lattice constant a in order to ensure that the periodicity of the lattice is properly maintained. In the z direction, the adsorbed molecules are confined to heights THS < Z < 20.0 A, sufficient to fit five or six full layers of adsorbed C 0 2 . No part of the molecule is allowed outside this slab. We have not employed Ewald summation in this investigation, since the van der Waals and quadrupolar interactions (the longest ranged electrostatic Chapter 4. The Adsorption of C02 Gas on MgO Crystals 87 interaction) between CO2 molecules decay sufficiently rapidly that the in-teractions with the periodic images need not be considered. Instead of the Ewald summation we have instituted a minimum image convention which considers only contributions from atoms separated by less than a cutoff ra-dius rcut. This cutoff radius must be less than half of the size of the periodic square. If this restriction is not obeyed, then spurious anisotropies will de-velop due to the fact that contributions from the larger distances can only come from particles located in the corners of the periodic box, where the minimum image distance r > L/2. Unless otherwise specified, we use a cut-off radius rcut = 13.5 A. This is sufficiently large to capture the significant interactions in the adsorbed layer. Although we did implement a long-range correction to account for the dispersion forces from particles separated by more than rcut, the effects of including this are negligible. We have done Monte Carlo simulations of the C0 2 /MgO system in both the canonical and the grand canonical emsembles. The canonical calcula-tions are useful for calculating the binding energy of a single molecule, or of a monolayer with a fixed number of particles. In our canonical calculations, a random displacement of the center of mass of a randomly selected CO2 molecule is attempted with probability 5/12, and a random rotation about the center of mass is attempted with the same probability. The maximum amplitude of these displacements is tuned to ensure acceptance probabilities of about 50%. In addition, it is important in our study to allow for large lat-tice moves, where the adsorbed CO2 molecules may be rotated azimuthally by 45 or 90° and moved over to different lattice sites. Although the accep-tance probabilities of these moves will always be low, without these moves the CO2 molecules will always remain fixed on a given lattice site. Two kinds of large moves are implemented; a rotation of 90° with a shift by a/2 in both x and y directions, and a rotation of 45° with a shift by a/4 in both x and y. The first lattice move corresponds to a shift from a parallel orientation to the herringbone, while the other is close to the sort of rearrangement we expect in the formation of the denser monolayer. These lattice moves are each attempted with probability 1/12. In our investigation of the possible transition to a higher density mono-layer structure, allowing fluctuations in the number of particles in the sim-ulation is important. The GCMC simulations provide a powerful method for doing this, but comparing the simulation results with experimental data such as adsorption isotherms requires some additional work. In the grand canonical ensemble, simulations are run at fixed volume V, temperature T, Chapter 4. The Adsorption of C02 Gas on MgO Crystals 88 and chemical potential p. The density of the system will be determined by the mean particle number, (N). In typical adsorption isotherm measure-ments, T and V are fixed, and the pressure P is varied while the amount of adsorbed material is measured. What is needed, then, to compare the two is a way to convert the chemical potential at which the simulation is run into a pressure corresponding to the state of the gas above the adsorbed layer. The gas above the surface is in equilibrium with the adsorbed molecules, so we can write that A W = Pgas, (4.5) where psm{ and p g a s are the chemical potentials of the surface layer and of the gas, respectively. /_i s u rf simply corresponds with the chemical potential of the G C M C simulation. To relate the pressure of the gas to its chemical potential, we follow the procedure outlined in Hil l [4]. Beginning with the definition of the activity, z, we can express the pressure as an infinite series in the activity, p °° i £ T = S > 0 V , (4-7) i > l where the bj(T) are cluster integrals, corresponding to the interactions of increasingly large clusters containing j molecules. In particular, &i (T) = 1. It is known that this series converges, and as long as the pressure is low enough, then we can assume the clustering is unimportant and Equation 4.7 reduces to = z. In our simulations, the pressure is quite low, never exceeding 10 Torr, and this approximation should be adequate. Combining this expression with Equation 4.6, we obtain the expression p g a s = kBT\n l — j . (4.8) Of course, as already mentioned, we choose to run our G C M C simulations at constant B, V and T, where B has been defined in Section 2.2.1. Using the definition of B, and Equation 4.8, Equation 4.5 may be rewritten as B=in(£?)- <4'9' Chapter 4. The Adsorption of C02 Gas on MgO Crystals 89 With the aid of Equation 4.9, the B value at which the simulation is run may easily be translated into the equivalent experimental pressure P. The advantage of introducing B is that the thermal wavelength A is explicitly removed from consideration, and need not be considered at any point in the simulation (see also Section 2.2.1). At each simulation step, the creation, deletion, or moving of a CO2 molecule is attempted with equal probability. If a particle move is attempted, the kind of move is selected with the same probability as it is in the canonical simulations (i.e. translation and rotation with probability 5/12 each, two lattice moves with probability 1/12 each). 4.3 Results Simulations in the grand canonical ensemble have been performed at temper-atures of 80 K, 100 K and 150 K. Results of simulations at 80 K and 100 K are shown in Figure 4.2. We have done simulations both starting from the empty surface, and starting from a complete monolayer, in order to get an idea of the hysteresis in the system. Our results at these temperatures do show the formation of a stable monolayer at the experimental pressure. On the other hand, although the monolayer formed in GCMC simulations starting from the empty state does tend to have a slightly higher surface density, there is no sign of any sharp transition to a denser monolayer at 100 K, contrary to our expectations. This may not be a very accurate picture, however, as the very low acceptance ratios for the creation and deletion of CO2 molecules indicate that the simulations approach a proper equilibrium state only very slowly at these low temperatures. This is a common problem in GCMC simulations at low temperatures. Another indication of these convergence issues includes the fact that the monolayers formed at low temperature consist of multiple domains, and contain many defects in the surface structure. Hysteresis is also important in the system, as evidenced by the contrast between simulations run from different initial states; if the perfect monolayer is formed, it is clearly a very stable configuration. We also attempted tak-ing a configuration resulting from a GCMC simulation at a given pressure or chemical potential, altering the pressure, and observing whether the cov-erage changed. At 100 K, even in very long simulations, the defect-ridden monolayer maintains its initial state. A typical snapshot of the monolayer formed in a GCMC simulation at 80 K is shown in Figure 4.3. Despite the Chapter 4. The Adsorption of C02 Gas on MgO Crystals 90 log(P e x p / l Torr) 2 . 5 2 (AO 1.5 1 0 . 5 -©-T = 80 T = 80 T = 100 T = 100 T = 100 K Kt 24a i i • * 4*A--10 -8 -6 log(P/l Torr) -2 Figure 4.2: Adsorption isotherms of CO2 on MgO at 80 K and 100 K from GCMC simulations. The MgO lattice size is L = 12a, unless otherwise specified. Isotherms which had the low density monolayer as the initial state are indicated by a dagger(+). The coverage is expressed as a fraction of the low density monolayer coverage on the given lattice of linear dimension n (Nid = n 2/2). The experimental pressure P e x p = 4 x 10~9 for the LEED experiments [127] is indicated by the vertical line. Note that the lines joining simulation results are shown only to guide the eye. Formation of bulk solid CO2 occurs at the pressures shown where the coverage increases rapidly past the maximum vertical range. convergence and hysteresis issues, the proposed mechanism for forming the denser monolayer does seem to be active, with several of the adsorbed CO2 molecules tilted and rotated to only cover one binding site. We also have observed bilayer and multilayer formation at these temper-atures. However, the range of pressure where the multilayers are stable is so narrow that it is difficult to see anything other than a fully occupied simu-lation box at higher pressures, corresponding to the formation of bulk solid C02- The pressure where bulk formation begins corresponds with the point on the isotherms where the coverage increases rapidly past the maximum Chapter 4. The Adsorption of C02 Gas on MgO Crystals 91 vertical range of Figure 4.2. Figure 4.3: A monolayer formed by C 0 2 on MgO. This snapshot is from a GCMC simulation at T = 80 K and PCo2 = 3.15 x 10~7 Torr on a square MgO lattice of si'/e L = 12a By running simulations at 150 K, the acceptance ratios for particle cre-ation and destruction are much better, hysteresis causes less problems, and we can be much more confident that the results are properly converged and represent a realistic representation of the monolayer structures. The results of these simulations on square lattices of varied sizes are shown in Figure 4.4. Although the pressures for monolayer formation are much higher in these simulations than in the low temperature experiments, they match quite well with the pressures for monolayer formation measured in the experiments of Trabelsi, who found that monolayers of C 0 2 form on MgO powders at T = 156 K and 0.2 < P < 15 Torr [149]. The simulations at 150 K tend to have a broad region of pressure where the C 0 2 monolayer forms in a struc-ture with every C 0 2 molecule covering two M g2 + binding sites, but with a significant t i l t angle, so that one end of each adsorbed molecule is much Chapter 4. The Adsorption of C02 Gas on MgO Crystals 92 log(P/l Torr) Figure 4.4: Adsorption isotherms of C 0 2 on MgO at 150 K from GCMC simulations. Isotherms which had the low density monolayer as the initial state are indicated by a dagger(+). The coverage is expressed as a fraction of the low density monolayer coverage on the given lattice of linear dimension n {N\d = n2/2). Note that the lines joining simulation results are shown only to guide the eye. more tightly bound than the other. The density of the simulated monolayer matches the LEED results, although there are still some discrepancies in the structural details. First, the large tilt angle does not agree with the observa-tion of two glide planes in the LEED experiment. Second, in our simulations the (y/2 x y/2)JR45° parallel structure appears to be more stable than the (2y/2 x y/2)R45° herringbone structure, even at low temperature (we de-scribe these results below, see Table 4.5). Within our model, we find that these structures are very similar in energy, and suitable refinements in the model might make the herringbone structure come out with a higher bind-ing energy. In any case, small differences in the energies of the two possible ordered low-density monolayers should not have a large effect on the forma-tion of a denser monolayer. The low density monolayers have a coverage of Nm — n<2:/2, where n = L/a is the number of binding sites along one linear Chapter 4. The Adsorption of C02 Gas on MgO Crystals 93 dimension of the simulation box. We have generated angular distributions of both the vertical tilt angle 6, being the angle between the CO2 molecular axis and the surface, and the azimuthal angle (j), defined as the angle between the projection of the molecu-lar axis into the xy-pl&ne and the a;-axis. These distributions were generated from canonical simulations, starting from a (y/2 x \/2)R45° configuration. The distribution of 9 for the low density monolayers is displayed in Figure 4.5, while the distribution of 0 is shown in Figure 4.6. The average tilt angle at all temperatures is just under 45°. At T = 80 K and 100 K , the azimuthal angle varies in a narrow range about 7r/2, indicating that all of the orienta-tions of the adsorbed molecules remain the same throughout the simulation. At T = 150 K , however, the molecules are distributed more or less evenly about four different orientations, indicating that the molecular axis of each CO2 does oscillate between the two possible directions on each axis, as well as often switching between the parallel and the herringbone configurations. This provides further reassurance regarding the possibly erroneous higher binding energy of the parallel configuration. It is apparent that at the higher temperature the influence of entropy ensures that neither of the low den-sity symmetric configurations is strongly preferred, and so the fact that the parallel configuration is the ground state in our model should not affect the results at 150 K . The more interesting results at T = 150 K occur at somewhat higher pres-sures/chemical potentials. Here, we have observed the formation of mono-layers with additional adsorbed molecules. As we hypothesised earlier, these additional molecules are adsorbed when a CO2 molecule lifts up off one of its binding sites and becomes tilted nearly vertically. Once one molecule has tipped up, the adsorbed molecules on either side tend to follow suit, until the full diagonal on the MgO (100) surface has tipped up and exposed a binding site. This allows the neighbouring row of molecules to move over, and frees up another diagonal somewhere on the MgO (100) lattice for more molecules to adsorb. Thus, another diagonal across the MgO lattice becomes filled with vertical CO2 molecules adsorbed to single binding sites, raising the total number of molecules in the monolayer by n, so that the total coverage on a lattice of even size is now Nm = n2/2 + n. A typical snapshot of this structure produced via G C M C simulation is shown in Figure 4.7. In simula-tions on a very large lattice, it was also possible to see a structure with four diagonal rows of singly bound C 0 2 ' s and a coverage of Nm = n2/2 + 2n. A snapshot of this configuration is shown in Figure 4.8. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 94 0.016 P{9) 6 (rad) Figure 4.5: Distribution of tilt angle 9 in the low density monolayer at three temperatures. The distributions are from canonical Monte Carlo simulations on a lattice of linear size L = 24a. There are often defects in the diagonals formed, especially in the larger lattices. One common defect which can be seen in Figure 4.8 is the shifting over by one lattice site of the singly bound diagonal. It may be more ac-curate in these cases to speak of two adsorbed C 0 2 molecules sharing three binding sites and alternating between them over the course of the simulation. It should be noted that even in the presence of the defects, the coverage nev-ertheless follows the formulae given above with surprising consistency. This is best shown in a graph of the average density of adsorbed molecules as a function of height, (N(z)}. A plot of this function shows that the first Nm C 0 2 molecules are bound to the surface with the carbon atom at a height just under 3 A above the MgO lattice, with additional molecules occasionally being created in a second layer approximately 6 A above the substrate. A plot of (N(z)) for each of the observed coverages on the L = 24a lattice is shown in Figure 4.9. Orientational distributions of the same structures are Chapter 4. The Adsorption of C02 Gas on MgO Crystals 95 0.035 0.03 0.025 0.02 T = 80 K T = 100 K T = 150 K P(4>) 0.015 h 0.01 h 0.005 4> (rad) 3JT 2 2TT Figure 4.6: Distribution of azimuthal angle </> in the low density monolayer at three temperatures. The distributions are from canonical Monte Carlo simulations on a lattice of linear size L — 24a. plotted in Figure 4.10 for the distribution of 9 and in Figure 4.11 for the distribution of (p. The angular distributions clearly show that a significant proportion of the adsorbed molecules have tilted upwards with much higher tilt angle 9 « 75°, and have rotated their azimuthal orientation so that they are oriented along the line x = y, midway between the x and y axes. Bilayer and multilayer formation are also seen at even higher pressures and chemical potentials, but they are usually more disordered and stable only in a very small pressure range, if at all, and it is difficult to make any conclusions about them. Still, we have observed the formation of bilayers with an equal coverage of 60 molecules in each layer in G C M C simulations on the MgO lattice of size L — 10a. The particles in the second layer occupy the spaces between the diagonal rows of the first layer, but the higher degree of disorder present in the second layer obscures any further details about its structure. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 96 Figure 4.7: A dense monolayer formed by CO2 on MgO. This snapshot is from a GCMC simulation at T = 150 K and Pco2 = 4.7 Torr on a square MgO lattice of size L = 12a. There are 84 molecules in the monolayer. The highest monolayer coverage we have been able to observe via GCMC simulation is a coverage of 60 C 0 2 molecules on the L = 10 lattice. This corresponds to a fractional coverage of 60/50 = 1.2. The cited experimental areal density of 11.9 A2/molecule equates to a coverage of approximately 74 CO2 molecules on the L — 10a lattice, or a fractional coverage of 1.49, well larger than any of our results. This indicates that despite the successful for-mation of high density monolayers, there are still some serious discrepancies between our simulation results and the experimental data. The dependence of the coverage on the linear size of the lattice indicates that some severe finite size effects are at play in this system. It appears that the possible structures of the high density monolayer depend on how well they agree with the periodicity of the lattice. This is seen most strongly when a simulation is made with a MgO lattice where n is an odd number. In this situation, the regular low density monolayer structures do not match with the periodicity Chapter 4. The Adsorption of C02 Gas on MgO Crystals 9 7 Figure 4.8: A dense monolayer formed by CO2 on MgO. This snapshot is from a GCMC simulation at T = 150 K and Pco2 = 6.3 Torr on a square MgO lattice of size L = 24a. There are 336 molecules in the monolayer. of the lattice, and are unstable. Instead, a structure with a single diagonal of tipped, singly bound CO2 molecules and a total coverage of N = (n2 + n)/2 is favoured. This structure is shown in Figure 4.12. The relation between the stability of the higher density structure and the lattice size can also be detected in the adsorption isotherms on lattices of different dimensions. It is apparent that the simulations on a square lattice of size L — 10a more easily go to the denser structure than do the simulations on a lattice of size L = 12a, suggesting that the most optimal dense structure matches well with the periodicity of the L = 10a lattice. To investigate the high density monolayers more thoroughly, we have cal-culated the binding energies of the various different monolayers in order to determine which structures are most stable. First, we did canonical simu-lations of all of the structures observed during our GCMC simulations at a temperature of 150 K. These results are shown in Table 4.4. Second, we Chapter 4. The Adsorption of C02 Gas on MgO Crystals 98 400 350 300 250 200 </V(z)> 150 100 50 h P = 1.39 Torr P = 3.79 Torr P = 6.71 Torr TV = 288 N — 312 N = 336 Figure 4.9: The average number of adsorbed CO2 molecules as a function of height, (N(z)), collected from GCMC simulations at T = 150 K and the indicated pressures. The number of molecules in each of the observed monolayer structures on the L = 24a lattice is also indicated for comparison. did a more systematic study of the possible ideal, ordered structures which might form at low temperatures, where the entropy should not have a sig-nificant effect. These simulations were done in the canonical ensemble, on the smallest possible lattices which could accomodate a given structure, as well as have r c u t large enough that the binding energy would be accurately calculated. The smallest lattice we have used for simulations is a square lattice with L = 8a, which means that r c u t > 4a = 11.88 A. We used this value of r c u t in all of the low temperature simulations in order to facilitate comparisions between different lattice sizes. The results of the T = 5 K simulations are shown in Table 4.5, and snap-shots of the most interesting structures are shown in Figure 4.13. In order to make the differences between the structures clearer, the same configurations are shown with the underlying MgO lattice omitted in Figure 4.14. In ad-Chapter 4. The Adsorption of C02 Gas on MgO Crystals 99 P{6) 0.009 1 Nm = 288 0.008 - IT Nm Nm = 312 . = 336 0.007 0.006 -0.005 l/ J -0.004 0.003 0.002 \ • 0.001 0 e (rad) Figure 4.10: Distribution of t i l t angle 9 in three monolayer coverages at T = 150 K. The distributions are from canonical Monte Carlo simulations on a lattice of linear size L = 24a. dition to the symmetric low-density configurations, we also simulated a low density structure with one diagonal in the orthogonal orientation (structure 1 in Table 4.5) on the L = 8a lattice. Moving on to the high density cases, based on the GCMC simulations we concentrated on the L = 10a lattice with 60 adsorbed CO2 molecules, where the high density monolayer is char-acterized by the presence of a singly bound high occupancy diagonal row, alternating with two rows in the low density structure. These structures cor-respond with structures 2, 3, and 4. We also briefly considered a structure (structure 5) with a high density row alternating with a single low density row, because the areal density in this configuration of 13.23 A 2/molecule is a reasonable match with the experimental report of 11.9 A 2/molecule [132]. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 100 0.007 0.006 Nm Nm Nm = 288 = 312 = 336 /-0.005 •\ ) 0.004 - l 1 / /M" / f '\' / ^ I 0.003 // \ 0.002 // \ If r 0.001 0 '•..\\ if i ji V\ ^ ' V N " " ' , ~ ' ' ' • " ^ ' x j'N 4> (rad) Figure 4.11: Distribution of azimuthal angle (j) in three monolayer coverages at T = 150 K . The distributions are from canonical Monte Carlo simulations on a lattice of linear size L = 24a. Note that the distribution in the interval from <j) = Tr to 4> = 2K is identical to the distribution shown, and so is not shown. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 101 Figure 4.12: A dense monolayer formed by C 0 2 on MgO. This snapshot is from a GCMC simulation at T = 150 K and PCo2 = 1.5 Torr on a square MgO lattice of size L = 11a. There are 66 molecules in the monolayer. n L/a Nm Binding Energy (kcal/mol) 10 50 6.214(2) 60 6.232(1) 11 66 6.236(3) 12 72 6.195(2) 84 6.226(4) 20 200 6.164(2) 220 6.229(8) 24 288 6.132(3) 312 6.222(2) 336 6.238(3) Table 4.4: The binding energy of monolayers at T = 150 K, found via canoni-cal ensemble simulation. The numbers in parentheses represent the estimated uncertainty in the last digit. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 102 n = L/a Structure Nm Binding Energy (kcal/mol) 8 (y/2 x V2)R45° 32 7.132 (2y/2 X y/2)R45° 32 7.081 1 32 7.103 10 (V2 x ^ # 4 5 ° 50 7.133 2 60 7.130 3 60 6.944 4 60 7.084 12 5 96 6.946 Table 4.5: The binding energy of monolayers at T = 5 K, found via canonical ensemble simulation. The numbered structures are displayed in Figure 4.13. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 103 Figure 4.13: Monolayer structures of C 0 2 on MgO from canonical Monte Carlo simulations at 5 K. The structures correspond sequentially with the numbered binding energy results in Table 4.5. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 104 Figure 4.14: Monolayer structures of CO2 on MgO from canonical Monte Carlo simulations at 5 K. The structures correspond sequentially with the numbered binding energy results in Table 4.5. The MgO substrate has been omitted for clarity. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 105 Even within the restriction to the L = 10a lattice, there are many possible configurations of molecules which can be imagined, and a systematic investi-gation of all of the possible structures would be extremely time consuming. In order to ensure that the configurations we have selected are indeed the most likely of several "nearby" states of similar energies, the configurations were first annealed at a temperature of 30 or 60 K . Angular distributions of some of the configurations are also plotted, with Figure 4.15 giving the distri-bution of 9 and Figure 4.16 giving the distribution of <p. Finally, (N(z)), the average number of adsorbed molecules as a function of height, of the same configurations is displayed in Figure 4.17. 0.05 0.04 0.02 0.01 o J, (2\/2 x v/2)/J<15° I \ (N/2 x y/2)IU5° - Structmc 2 Structure 3 - Jj • • • • Structure 4 - II A. i V /V 1 // \.. -// ••. i •' % •• a •' \ :. ii •' y f • i\ i : 'A •. / .• >.\ / .• >;••. / .• >.\ 6 (rad) Figure 4.15: Distribution of tilt angle 9 in various monolayer structures, generated from canonical simulations at T = 5 K . The distributions shown correspond to the given energies listed in Table 4.5. Chapter 4. The Adsorption of C02 Gas on MgO Crystals 106 0.05 0.04 0.03 P(4>) 0.02 0.01 0 1 i i ( 2 \ / 2 x V2)H45' 1 x v/2)/4i° Structure 2 Structure 3 - s ill • • • • Structure 4 I i i I! ii i i i\i \ ii; iii I • •si •'i' ijii |i J! li Ij 5 :i it) • i i i i i i i ) •' li ii j; 4 i' i ; i ; i '. ;i :i i ji i i S ii 1 « ii 1 . Ill II i is \ i i 4 \ \ : 1 i s c ' . l i 4> (rad) Figure 4.16: Distribution of azimuthal angle 4> in various monolayer struc-tures, generated from canonical simulations at T = 5 K . The distributions shown correspond to the given energies listed in Table 4.5. 4.4 Discussion Several aspects of the results of our simulations of the C 0 2 / M g O system should be discussed. We need to determine how our results compare with the experimental picture, and then, we can discuss the various monolayer structures. 4.4.1 Comparion of Simulation vs. Experiment Our determinations of the monolayer binding energy generally come out to « 6.2 kcal/mol at T = 150 K , rising to « 7.1 kcal/mol at T = 5 K. These results are somewhat lower than the prior experimental and theoreti-cal determinations, which we recall from the discussion of Section 4.1.2 vary from about 7.1 kcal/mol for the somewhat questionable LITD results [137], Chapter 4. The Adsorption of C02 Gas on MgO Crystals 107 10 h z(A) (v/5 x </2)R4&° Structure 2 Structure 3 Structure 4 3.5 Figure 4.17: The average number of adsorbed CO2 molecules as a function of height, (N(z)), in various monolayer structures, collected from canonical simulations at T = 5 K . to 11.0 kcal/mol for the determination based on the volumetric adsorption isotherms [132]. This discrepancy is a cause for concern, however, it should be pointed out that it is easy to justify changing many of the model pa-rameters we have used, and slightly different values can lead to significant changes in the binding energy and the structure of the adsorbed monolayer. In the early stages of this project, we experimented with different models and values of the model parameters before settling on the current model and parameters. In one prototype model, the binding energy was much too small, and though monolayers were stable in canonical simulations, when we extended to the grand canonical case the monolayer would either evaporate, or, at high pressures, it would condense into the bulk solid. In another at-tempt, with different parameters and lacking the Cg term in the CO2-CO2 interaction, the binding energy was somewhat higher, more closely matching the experimental data. In addition, the adsorbed CO2 molecules lay flat on Chapter 4. The Adsorption of C02 Gas on MgO Crystals 108 the surface equally sharing two M g 2 + binding sites, in agreement with the L E E D results. Unfortunately, with this parameter set it was impossible to observe any higher density monolayers. It seems that in order to undergo the transition to the higher density structures, there must be a significant tilt angle already present, or else the tilt-and-rotate mechanism for the liberation of new binding sites cannot operate. The model and parameters we have settled on seem to be a good com-promise. Though our calculations may underestimate the binding energy, models with higher binding energies do not form the interesting high den-sity structures which drove this investigation in the first place. It is also possible that the experimental data is wrong, and severely overestimates the binding energy. The fact that we do see the formation of monolayers in pres-sure ranges consistent with both the adsorption isotherms and the L E E D experiments suggests that our model is reasonably accurate. There are other issues regarding the L E E D experiments which require some explanation. In Ref. [127], Suzanne et al. reported that they observed monolayer formation at T = 60 K with a CO2 gas pressure P e x p = 4 x 10~9 Torr. They also give the vapour pressure of bulk CO2 at that temperature as P v a p = 3.1 x 1 0 - 1 3 Torr. Since P e x p > PVap) it is difficult to imagine how a monolayer might form under these conditions, instead of forming bulk CO2. We find that the very close agreement between the binding energies of the various low- and high-coverage monolayers makes it probable that a sig-nificant number of defects and grain boundaries occur in the monolayer at all temperatures. Below T — 93 K , there are likely large domains of molecules in the ordered (2\/2 x v / 2)P45° structure, which produce the measured L E E D pattern, but there should also be large domains containing defects, and if a distinct high density phase exists, there may also be rows of singly-bound tipped molecules at all temperatures. The sharp disappearance of the L E E D pattern above T = 93 K is consistent with the formation of the high density monolayer, as it seems that there are enough defects in the singly bound diagonals to obliterate any long range order. The finite size dependence of the high density monolayers makes it impossible for us to say anything quan-titative about the preferred density of the high density monolayer without extending our simulations to a MgO lattice large enough to accomodate a very large number of coexisting singly bound diagonals, so that any effects dependent on the periodicity of the lattice would disappear. Chapter 4. The Adsorption of C 0 2 Gas on MgO Crystals 109 4.4.2 The Structure of the High Density Monolayer Leaving aside for a moment the important issues regarding the model specifi-cation and the experimental work, it is apparent that within our model, high density monolayers are quite stable. In fact, the high density state on the L = 10a lattice comes very close to being the ground state, with a difference of only 0.03 kcal/mol between the most stable low density and high density structures at T = 5 K (Table 4.5). In addition, at higher temperatures it seems that the low density monolayer disorders slightly more than the higher density layer, and this leads to a large difference in the binding energies at T = 150 K favouring the high density monolayer, with the results on the largest lattice showing a difference of « 0.1 kcal/mol. It is apparent even from the G C M C simulations that when a diagonal row of C 0 2 molecules tips up and goes to the singly bound configuration, the diagonal rows of molecules on either side of it find it very favourable to take on orthogonal orientations. In addition, the doubly-bound rows tilt towards the singly-bound row, in this configuration: IJnjP. This is consistent with the preferred orientation of point quadrupoles. In two dimensions, a lattice of point quadrupoles will ideally form T-shaped structures, with the axis of one quadrupole pointing at the center of another [10,150]. The (2\/2 x y/2)RAb° flat monolayer is a consequence of this interaction, except for the shift required to maintain commensurability with the lattice. If the molecules are tilted, then the quadrupolar interaction is not as favourable. However, assuming for the moment that the adsorbed molecules must be tilted, the high density configuration shown, with the molecular axes of neighbouring C 0 2 ' s nearly orthogonal, seems to be the most favourable. The configuration of three diagonal rows of high density C 0 2 molecules in the structure shown above is quite robust, and occurs in all of the G C M C simulations. Of course, many structures can be imagined which are consis-tent with the preferred configuration. On the L = 10a square lattice, three different orientations of the two low density rows between each high density row must be considered. These structures correspond with structures 2, 3, and 4 in Table 4.5 and Figure 4.13. As already mentioned, the presence of the high density diagonal row makes the orthogonal orientation of the neigh-bouring low density rows much more strongly favoured. This is reflected in the significantly lower binding energy calculated for structure 3 compared Chapter 4. The Adsorption of C02 Gas on MgO Crystals 110 with structure 2. In fact, structure 3 has nearly an identical binding energy per molecule to the higher density structure 5, which we have not been able to produce in any of our G C M C simulations. Given the fairly substantial energy difference- between the herringbone and parallel configurations, it is initially surprising that structure 4, with alternating double rows of parallel C02's, does not have a lower binding energy than structure 1. A closer look shows that the parallel rows must tilt in opposite directions in order to preserve the strong quadrupolar interactions with the high density rows. Apparently, the energy loss required to tilt parallel rows in opposite directions cannot compensate for the energy gained by taking on the parallel configuration instead of the herringbone. Another way to view the energy differences between the structures is simply to compare the tilt angles and heights above the surface, shown in Figures 4.15 and 4.17. The distribution of 9 shows that among the low density configurations, the parallel configuration has the lower tilt angle and the largest binding energy. The same rule seems to apply to the high density configurations; structure 2 has the lowest tilt angle (in the doubly bound diagonal rows) of any of the structures studied, and it has the highest binding energy. Finally, our result for the higher density structure 5, where one singly bound row alternates with one doubly bound row, shows that this structure is significantly less stable than the less dense structures. In fact, when this structure is simulated in the canonical ensemble on the L = 12a lattice at the higher temperature of 150 K , it is unstable. At this temperature one of the singly bound rows lifts off the surface, these CO2 molecules bind in a second layer, and the monolayer relaxes quickly to the (Nm) = 84 structure. This makes the very high coverage reported by Trabelsi [132] seem too high. It is possible that this value refers to a measurement taken very close to the onset of bilayer formation, so that in fact, it is a measurement of the density of the monolayer plus the partially covered bilayer. Further work, both experimental and theoretical, is obviously required to make a more accurate determination of the density in the high-coverage monolayer. I l l Chapter 5 Conclusions and Future Work 5 . 1 Simulations of the Condensation of Coulombic and Related Fluids Over the course of the work done on simulating the condensation of ionic and other fluids, the focus changed from attempting to determine the universality class of the R P M to attempting to understand why calculations of the heat capacity lead to vastly different results when calculated in different ways. As such, the most important conclusions of our work are quite different from anything that could have been predicted at the outset. Our research into heat capacity calculations was driven by the need to understand why two different research groups, one working in the canonical ensemble [92] and one working in the grand canonical ensemble [93], came up with such qualitatively different results in their calculations of Cv- Af-ter considering several other possibilities, we determined that the difference between the two calculations arises due to the use of the different ensembles. Having explained this discrepancy, we proceeded to examine which of these calculations, if any, can be considered trustworthy. We did this by calculating the heat capacity of models of hard-sphere fluids with attractive interactions u(r) oc r~n. Unlike the R P M , for which no good analytical theory exists which can predict the critical behaviour, renormalization group theory provides a good understanding of the critical behaviour of hard-sphere fluids with attractive interactions. By calculating Cv in the u(r) oc r~n fluids near the condensation critical point, in both the canonical and grand canonical ensembles, and comparing the results with the predictions of the renormalization group theory, we have shown unambiguously that Cy calculations in the grand canonical ensemble are not consistent with the analytical predictions. By contrast, the canonical calculations of Cy , while being very sensitive to finite size effects, are consis-tent with all available analytical predictions. A l l of our calculations of CV in Chapter 5. Conclusions and Future Work 112 both ensembles are reproducible within statistical uncertainties, and so this is not simply a numerical issue. The reasons for the ensemble discrepancy are still unclear, however, we can offer some ideas. It appears that the finite size effects which affect all computer simulations of bulk, macroscopic systems manifest in different ways depending on the ensemble which is used in the simulations. These differences are especially important in simulations near a critical point, where density fluctuations become very large. These fluctuations lead to a divergence of the isothermal compressibility KT at the critical point in the thermodynamic limit, regardless of the universality class of the system. The challenge for theorists now must be to carefully re-examine the be-haviours of the various thermodynamic response functions near a critical point, determine how they are affected by finite size effects, and hence, de-cide whether G C M C simulations can provide trustworthy representations of the behaviour of these response functions near a critical point. This improved understanding of finite size effects is necessary in order to take advantage of the great utility of the grand canonical ensemble in simulating phase transi-tions. The canonical ensemble is easier to work with, but as we have seen in the calculations of Cy in the n = 6 system, finite size effects are large and it is difficult to extrapolate canonical results to the thermodynamic limit. This is especially true of heat capacity calculations, which diverge only very weakly at the critical point. A recent theoretical attempt to understand finite size effects in the microcanonical ensemble is encouraging [151], suggesting that a similar attempt on the grand canonical ensemble would be worthwhile. A fruitful avenue for further research would be a large-scale canonical ensemble simulation of a lattice gas at fixed density, or equivalently, a fixed magnetization simulation of an Ising model where only spin exchanges are attempted, rather than spin flips. Simulations of this simple model could be done to an exhaustive degree of accuracy, and the analysis of the data would be beneficial to theoretical attempts to improve the understanding of finite size effects in simulations. In a similar vein, the heat capacity of a spin model with classical critical exponents should be calculated using G C M C . It might be possible to extend these calculations to systems sufficiently large that the equivalence of ensembles in the thermodynamic limit could be seen. A simulation study of the critical behavior of a continuous Ising model in one and two dimensions with varied interaction ranges has been done [152], but unfortunately, their results for Cy appear to be wrong, since they use the canonical fluctuation formula for Cy, while their simulation allows fluctua-Chapter 5. Conclusions and Future Work 113 tions in the magnetization. The only verifiable analytical results for finite size effects on the heat capacity are calculations on square Ising lattices, which can be solved exactly [122,153]. This calculation is, however, also a grand canonical result. It would be interesting to repeat this analytical calculation in the canonical ensemble with the total magnetization kept constant. It may even be worth simulating the square Ising lattice, in order to determine how simulation results in both canonical and grand canonical ensembles differ from corresponding exact expressions for Cy on finite Ising spin systems. It might also be worth studying a simple model such as a van der Waals model analytically, and determine exactly how finite size effects enter into the expressions for the heat capacity. Starting from analytical expressions for the canonical and grand canonical partition functions, it is relatively easy to obtain expressions for the heat capacity in the thermodynamic limit. However, when finite size effects are incorporated explicitly, the mathematics become much more challenging. An alternate approach is to start from a Taylor expansion around the critical point, as we briefly discussed in Section 2.3.1. Again, this has been done in the thermodynamic limit [23,24], but it is unclear how finite size effects might be included explicitly. Finally, we were able to return to the question we asked at the outset of this investigation: what is the universality class ofthe R P M and other related ionic fluid models? We did G C M C simulations of the R P M and the 2:1 asymmetric P M , as well as a model consisting of dumbbells of associated ion pairs, and analysed them using MFFSS. This technique gave good estimates ofthe critical parameters for all ofthe models, which are fully consistent with those arrived at by other workers. There is no evidence from MFFSS that any of the ionic fluids models considered belong to anything other than the Ising universality class. We also calculated Cy in the canonical ensemble for the R P M and the dumbbells, and found no convincing evidence of a critical peak in either model. We believe the most likely reason for the lack of this peak is simply the slow approach towards the thermodynamic limit of heat capacity calculations with small systems. The fact that the critical density in the R P M is a full order of magnitude lower than it is in the n = 6 model may be the only reason no peak has been seen yet in the ionic fluids, even though simulations of quite large systems have been done. As for the influence of ion pairing on the criticality of the ionic fluids, the only important difference we have been able to discern between the ionic R P M and the dumbbell system is that the dumbells have a slightly lower critical density, which we have attributed to the fact that the dumbbells, since Chapter 5. Conclusions and Future Work 114 they cannot leave contact with each other, are effectively smaller particles. We have not seen any evidence for a difference in universality class caused by the absence or presence of free ions, as had been speculated. However, we stress that our conclusions about the difficulty of observing Ising-like critical behaviour in Cy still apply in this case. It remains possible that simulations on larger ionic systems may show differences in Cy based on whether the ion pairs break down near the phase transition or not. 5.2 Simulations of C O 2 Adsorption on a MgO Substrate Crudely speaking, our simulation results do agree with the experimental picture, in that an adsorbed monolayer forms and is stable at simulation temperatures and pressures comparable with the experimental temperatures and pressures. In addition, the coverage of the low density, ordered monolayer is consistent with that reported in the L E E D experiments. The adsorbed CO2 molecules each cover two M g 2 + binding sites, with a large tilt angle. However, there are still some discrepancies between our results and the L E E D data with respect to the preferred structure of the ordered monolayer. We have also been able to simulate the formation of higher density mono-layers at elevated temperatures. These structures are formed when adsorbed CO2 molecules tilt up and rotate, so that additional binding sites are made available for more molecules to adsorb. The singly bound molecules form in diagonals on the MgO (100) surface. Typically, there are kinks in these diagonals which effectively destroy any long-range order which might exist. Because of severe finite size effects, it is difficult to draw definite conclusions about the true structure of the high density phases in the thermodynamic limit. Nevertheless, we can conclude with confidence that these high density monolayers have binding energies comparable with the low density structures, and they may occur in this and in similar systems. The simulation of surface adsorption with realistic models is a relatively young field, and much work of a foundational nature remains to be done. The models used for the interparticle potentials contain a large number of parameters, many of which are interdependent. Only with further work can the potential models be refined to the point where we can be absolutely confident in the results we obtain via simulation. Chapter 5. Conclusions and Future Work 115 This is the first attempt we are aware of to simulate physisorption of polyatomic molecules on a crystal surface in the grand canonical ensemble with a realistic, continuum model, rather than a lattice or some other simple model. There are many improvements which should be made before making future attempts on similar systems. As only one suggestion, biasing the cre-ation and deletion attempts to account for the preferred binding sites of the CO2 molecules on the MgO surface would greatly improve the convergence rate of the G C M C calculations. It would be interesting to apply G C M C techniques to other, related sys-tems. One example would be the N20/MgO system, which has a similar surface structure and vibrational spectrum to the C02 /MgO system [136]. It would be another candidate for competing low- and high-density monolayer structures. Another tantalizing case is the C O / M g O system. Although this system has been simulated extensively using canonical Monte Carlo simula-tions by the Jack group [123], it would be interesting to do G C M C simulations on the system, due to the possibility of observing many different high-order commensurate phases of different coverages (the "devil's staircase" [124]). G C M C would allow us to directly test the relative stability of surface phases of different densities. Bibliography 116 [1] N . Metropolis and S. Ulam. J. Am. Stat. Ass., 44:335, 1949. [2] J. von Neumann. Theory and organization of complicated automata. In A . W. Burks, editor, Theory of Self-Reproducing Automata [by] John von Neumann, chapter 1. University of Illinois Press, Urbana, 1966. Based on lectures delivered at the University of Illinois in December 1949. [3] J. S. Rowlinson. Cohesion: A Scientific History of Intermolecular Forces. Cambridge University Press, Cambridge, 2002. [4] T. L. Hi l l . Statistical Mechanics: Principles and Selected Applications. Dover Publications, Inc., New York, 1987. [5] M . P. Allen and D. J. Tildesley. Computer simulation of liquids. Claren-don Press, Oxford, 1987. [6] N . Goldenfeld. Lectures on Phase Transitions and the Renormalization Group. Addison-Wesley, Reading, 1992. [7] M . N . Barber. Finite size scaling. In C. Domb and J. L. Lebowitz, edi-tors, Phase Transitions and Critical Phenomena, volume 8, chapter 2. Academic Press, London, 1983. [8] K . T. Tang and J. P. Toennies. J. Chem. Phys., 80:3726, 1984. [9] J. 0 . Hirschfelder, C. F. Curtiss, and R. B. Bird. Molecular Theory of Gases and Liquids. John Wiley and Sons, New York, 1964. [10] C. G. Gray and K . E. Gubbins. Theory of Molecular Fluids, Volume 1. Clarendon Press, Oxford, 1984. [11] P. Ewald. Ann. Phys., 64:253, 1921. Bibliography 117 [12] M . Plischke and B. Bergersen. Equiibrium Statistical Mechanics, 2nd Edition. World Scientific, Singapore, 1994. [13] G. E. Norman and V . S. Filinov. High Temp. (USSR), 7:216, 1969. [14] D. J. Adams. Mol Phys., 28:1241, 1974. [15] D. J. Adams. Mol. Phys., 29:307, 1975. [16] J. P. Valleau and G. M . Torrie. A guide to monte carlo for statisti-cal mechanics: Byways. In B. J. Berne, editor, Statistical Mechanics A: Equilibrium Techniques, Modern Theoretical Chemistry, volume 5, chapter 5. Plenum, New York, 1977. [17] A . M . Ferrenberg and R. H. Swendsen. Phys. Rev. Lett, 61:2635, 1988. [18] T. Andrews. Phil. Trans. Roy. Soc, 159:575, 1869. [19] J. D. van der Waals. On the continuity of the gas and liquid states. PhD thesis, Leiden, 1873. English translation by J. S. Rowlinson, in Studies in statistical mechanics, Amsterdam, 1988, v. 14. [20] P. M . Chaikin and T. C. Lubensky. Principles of Condensed Matter Physics. Cambridge University Press, Cambridge, 2000. [21] G. Orkoulas, M . E. Fisher, and A. Z. Panagiotopoulos. Phys. Rev. E, 63:051507, 2001. [22] H. E. Stanley. Introduction to Phase Transitions and Critical Phenom-ena. Oxford University Press, Oxford, 1971. [23] J. V . Sengers and J. M . H. Levelt Sengers. Critical phenomena in classical fluids. In C. A . Croxton, editor, Progress in Liquid Physics, chapter 4. John Wiley and Sons, Chichester, 1978. [24] Y . C. K i m , M . E. Fisher, and G. Orkoulas. Phys. Rev. E, 67:061506, 2003. [25] L. Onsager. Phys. Rev., 62:117, 1944. [26] M . E. Fisher. Rev. Mod. Phys., 70:653, 1998. [27] R. Guida and J. Zinn-Justin. J. Phys. A: Math. Gen., 31:8103, 1998. Bibliography 118 [28] F. G. Wegner. Phys. Rev. B, 5:4529, 1972. [29] A . M . Ferrenberg and D. P. Landau. Phys. Rev. B, 44:5081, 1991. [30] Y . C. K im, M . E. Fisher, and G. Orkoulas. Phys. Rev. E, 68:041506, 2003. [31] G. Orkoulas, A . Z. Panagiotopoulos, and M . E. Fisher. Phys. Rev. E, 61:5930, 2000. [32] N . B. Wilding and A. D. Bruce. J. Phys.: Condens. Matter, 4:3087, 1992. [33] N . B. Wilding and M . Muller. J. Chem. Phys., 102:2562, 1995. [34] N . B. Wilding. Phys. Rev. E, 52:602, 1995. [35] Y . C. K i m and M . E. Fisher. J. Phys. Chem. B, 108:6750, 2004. [36] C. N . Yang and C. P. Yang. Phys. Rev. Lett, 13:303, 1964. [37] M . E. Fisher and G. Orkoulas. Phys. Rev. Lett, 85:696, 2000. [38] G. Orkoulas, M . E. Fisher, and C. Ustiin. J. Chem. Phys., 113:7530, 2000. [39] P. Walden and M . Centnerschwer. Z. Phys. Chem., 42:432, 1903. [40] H. L. Friedman. J. Phys. Chem., 66:1595, 1962. [41] H. Weingartner and W. Schroer. Adv. Chem. Phys., 116:1, 2001. [42] Y . Guissani and B. Guillot. J. Chem. Phys., 101:490, 1994. [43] M . L. Japas and J. M . H. Levelt Sengers. J. Phys. Chem., 94:5361, 1990. [44] K . S. Pitzer. Acc. Chem. Res., 23:333, 1990. [45] H. Weingartner and W. Schroer. J. Mol. Liq., 65/66:107, 1995. [46] M . E. Fisher. J. Stat. Phys., 75:1, 1994. Bibliography 119 [47] K. S. Pitzer, M. C. P. De Lima, and D. R. Schreiber. J. Phys. Chem., 89:1854, 1985. [48] R. R. Singh and K. S. Pitzer. J. Chem. Phys., 92:6775, 1990. [49] S. Wiegand, M. E. Briggs, J. M. H. Levelt Sengers, M. Kleemeier, and W. Schroer. J. Chem. Phys., 109:9038, 1998. [50] H. L. Bianchi and M. L. Japas. J. Chem. Phys., 115:10472, 2001. [51] M. Kleemeier, S. Wiegand, W. Schroer, and H. Weingartner. J. Chem. Phys., 110:3085, 1999. [52] H. Weingartner, S. Wiegand, and W. Schroer. J. Chem. Phys., 96:848, 1991. [53] K. C. Zhang, M. E. Briggs, R. W. Gammon, and J. M. H.Levelt Sen-gers. J. Chem. Phys., 97:8692, 1992. [54] T. Narayanan and K. S. Pitzer. J. Chem. Phys., 102:8118, 1995. [55] M. Kleemeier, S. Wiegand, T. Derr, V. Weiss, W. Schroer, and H. We-ingartner. Ber. Bunsenges. Phys. Chem., 100:27, 1996. [56] T. Heimburg, S. Z. Mirzaev, and U. Kaatze. Phys. Rev. E, 62:4963, 2000. [57] H. Weingartner, T. Merkel, U. Maurer, J. P. Conzen, Ff. Glassbrenner, and S. Kashammer. Ber. Bunsenges. Phys. Chem., 95:1579, 1991. [58] D. R. Schreiber, M. C. P. De Lima, and K. S. Pitzer. J. Phys. Chem., 91:4087, 1987. [59] P. J. Camp and G. N. Patey. Phys. Rev. E, 60:1063, 1999. [60] P. J. Camp and G. N. Patey. J. Chem. Phys., 114:399, 2001. [61] P. W. Debye and E. Hiickel. Phys. Z., 24:185, 1923. [62] V. McGahay and M. Tomozawa. J. Non-Cryst. Solids, 109:27, 1989. [63] V. McGahay and M. Tomozawa. J. Chem. Phys., 97:2609, 1992. Bibliography 120 [64] M . E. Fisher and Y . Levin. Phys. Rev. Lett, 71:3826, 1993. [65] N . Bjerrum. Kgl. Dan. Vidensk. Selsk. Mat-fys. Medd., 7:1, 1926. [66] R. M . Fuoss. J. Am. Chem. Soc., 80:5059, 1958. [67] H. L. Friedman and B. Larsen. J. Chem. Phys., 70:92, 1979. [68] M . J. Gillan. Mol. Phys., 49:421, 1983. [69] B. Guillot and Y . Guissani. Mol. Phys., 87:37, 1996. [70] V . C. Weiss and W. Schroer. J. Chem. Phys., 108:7747, 1998. [71] H. Weingartner, V . C. Weiss, and W. Schroer. J. Chem. Phys., 113:762, 2000. [72] G. Stell, K . C. Wu, and B. Larsen. Phys. Rev. Lett, 37:1369, 1976. [73] Y . Zhou, S. Yeh, and G. Stell. J. Chem. Phys., 102:5785, 1995. [74] J. S. Hoye, E. Lomba, and G. Stell. Mol. Phys., 75:1217, 1992. [75] J. W. Jiang, L. Blum, and O. Bernard. Mol. Phys., 99:1765, 2001. [76] J. M . Caillol, D. Levesque, and J. J. Weis. J. Chem. Phys., 116:10794, 2002. [77] A . Ciach and G. Stell. J. Mol. Liq., 87:253, 2000. [78] O. V . Patsahan and I . M . Mryglod. J. Phys.: Condens. Matter, 16:L235, 2004. [79] M . E. Fisher, S. Ma, and B. G. Nickel. Phys. Rev. Lett, 29:917, 1972. [80] J. Sak. Phys. Rev. B, 8:281, 1973. [81] G. Stell. Phys. Rev. B, 8:1271, 1973. [82] B. P. Chasovskikh and P. N . Vorontsov-Vel'yaminov. High Temp. (USSR), 14:174, 1976. [83] J. P. Valleau. J. Chem. Phys., 95:584, 1991. Bibliography 121 [84] A . Z. Panagiotopoulos. Fluid Phase EquiL, 76:97, 1992. [85] J. M . Caillol. J. Chem. Phys., 101:2161, 1994. [86] J. M . Caillol, D. Levesque, and J. J. Weis. Phys. Rev. Lett, 77:4039, 1996. [87] J. P. Valleau. J. Chem. Phys., 108:2962, 1998. [88] J. M . Caillol, D. Levesque, and J. J. Weis. J. Chem. Phys., 107:1565, 1997. [89] G. Orkoulas and A. Z. Panagiotopoulos. J. Chem. Phys., 110:1581, 1999. [90] Q. Yan and J. J. de Pablo. J. Chem. Phys., 111:9509, 1999. [91] A. Z. Panagiotopoulos. J. Chem. Phys., 116:3007, 2002. [92] J. P. Valleau and G. Torrie. J. Chem. Phys., 108:5169, 1998. [93] E. Luijten, M . E. Fisher, and A. Z. Panagiotopoulos. J. Chem. Phys., 114:5468, 2001. [94] C. D. Daub, P. J. Camp, and G. N . Patey. J. Chem. Phys., 118:4164, 2003. [95] P. J. Camp, C. D. Daub, and G. N . Patey. Liquid-vapour criticality in coulombic and related fluids: What can be learned from computer simulations? In D. Henderson and M . F. Holovko, editors, Ionic Soft Matter: Modern trends in theory and applications. Kluwer Academic, Dordrecht, 2004. in press. [96] J . P. Valleau and G . Torrie. J. Chem. Phys., 117:3305, 2002. [97] A . Ciach and G. Stell. J. Chem. Phys., 114:382, 2001. [98] A. Ciach and G. Stell. J. Chem. Phys., 114:3617, 2001. [99] V . Kobelev, A . B. Kolomeisky, and M . E. Fisher. J. Chem. Phys., 116:7589, 2001. Bibliography 122 [100] A. Z. Panagiotopoulos and S. K. Kumar. Phys. Rev. Lett, 83:2981, 1999. [101] E. Luijten, M. E. Fisher, and A. Z. Panagiotopoulos. Phys. Rev. Lett, 88:185701, 2002. [102] G. Orkoulas and A. Z. Panagiotopoulos. J. Chem. Phys., 101:1452, 1994. [103] D. N. Theodorou and U. W. Suter. J. Chem. Phys., 82:955, 1985. [104] C. D. Daub, P. J. Camp, and G. N. Patey. J. Chem. Phys. in press. [105] P. J. Camp. Phys. Rev. E, 67:011503, 2003. [106] J. C. Shelley and G. N. Patey. J. Chem. Phys., 103:8299, 1995. [107] C. D. Daub, G. N. Patey, and P. J. Camp. J. Chem. Phys., 119:7952, 2003. [108] J. M. Romero-Enrique, G. Orkoulas, A. Z. Panagiotopoulos, and M. E. Fisher. Phys. Rev. Lett, 85:4558, 2000. [109] S. Moghaddam and A. Z. Panagiotopoulos. J. Chem. Phys., 118:7556, 2003. [110] A. G. Moreira and R. R. Netz. Eur. Phys. J. D, 13:61, 2001. [Ill] P. J. Camp and G. N. Patey. J. Chem. Phys., 111:9000, 1999. [112] A. Z. Panagiotopoulos and M .E. Fisher. Phys. Rev. Lett, 88:045701, 2002. [113] J.-P. Hansen and I. R. McDonald. Theory of Simple Liquids. Academic Press, London, 1986. [114] K. Binder. Introduction. In K. Binder, editor, Topics in Applied Physics: The Monte Carlo Method in Condensed Matter Physics, vol-ume 71, chapter 1. Springer-Verlag, Berlin, 1992. [115] H. W. J. Blote, E. Luijten, and J. R. Heringa. J. Phys. A, 28:6289, 1995. Bibliography 123 [116] E. Luijten and H. W. J. Blote. Phys. Rev. B, 56:8945, 1997. [117] E. Luijten and H. W. J. Blote. Phys. Rev. Lett., 89:025703, 2002. [118] B. Widom and J. S. Rowlinson. J. Chem. Phys, 52:1670, 1970. [119] N. D. Mermin. Phys. Rev. Lett, 26:957, 1971. [120] A. Kostrowicka Wyczalkowska, M. A. Anisimov, J. V. Sengers, and Y. C. Kim. J. Chem. Phys., 116:4202, 2002. [121] M. I. Bagatskii, A. V. Voronel', and V. G. Gusak. Sov. Phys. JETP, 16:517, 1963. [122] A. E. Ferdinand and M. E. Fisher. Phys. Rev., 185:832, 1969. [123] A. K. Sallabi and D. B. Jack. Phys. Rev. B, 62:R4841, 2000. [124] P. Bak. Rep. Prog. Phys., 45:587, 1982. [125] E. Domany, M. Schick, J. S. Walker, and R. B. Griffiths. Phys. Rev. B, 18:2209, 1978. [126] E. Domany and M. Schick. Phys. Rev. B, 20:3828, 1979. [127] J. Suzanne, V. Panella, D. Ferry, and M. Sidoumou. Surf. Sci. Lett, 293:L912, 1993. [128] V. Panella, J. Suzanne, P. N. M. Hoang, and C. Girardet. J. Phys. I France, 4:905, 1994. [129] W. H. Keesom and J. W. L. Kohler. Physica, 1:167, 1934. [130] L. W. Bruch. J. Chem. Phys., 79:3148, 1982. [131] J. Heidberg, E. Kampshoff, R. Kuhnemuth, and O. Schonekas. Surf. Sci, 251/252:314, 1991. [132] M. Trabelsi. PhD thesis, University of Aix-Marseille 2, France, 1991. [133] R. C. Weast and M. J. Astle. CRC Handbook of Chemistry and Physics, 60th Edition. CRC Press, Boca Raton, 1980. [134] J. Heidberg and D. Meine. Ber. Bunsenges. Phys. Chem., 97:211, 1993. Bibliography 124 [135] J. Heidberg, D. Meine, and B. Redlich. J. Electron Spectrosc. Relat. Phenom., 64:599, 1993. [136] J. Heidberg and B. Redlich. Surf. Sci., 368:140, 1996. [137] D. L. Meixner, D. A. Arthur, and S. M. George. Surf. Sci., 261:141, 1992. [138] C. Girardet, S. Picaud, and P. N . M. Hoang. Europhys. Lett, 25:131, 1994. [139] V. M. Rozenbaum and S. H. Lin. J. Chem. Phys., 112:9083, 2000. [140] J. C. Polanyi, R. J. Williams, and S. F. O'Shea. J. Chem. Phys., 94:978, 1991. [141] V. J. Barclay, D. B. Jack, J. C. Polanyi, and Y. Zeiri. J. Chem. Phys., 97:9458, 1992. [142] A. D. Buckingham and R. L. Disch. Proc. R. Soc. A, 273:27, 1963. [143] C. S. Murthy, S. F. O'Shea, and I. R. McDonald. Mol. Phys., 50:531, 1983. [144] J. E. Lennard-Jones and B. M. Dent. Trans. Faraday Soc, 24:92, 1928. [145] J. Heidberg, E. Kampshoff, R. Kiihnemuth, and O. Schonekas. Surf. Sci., 272:306, 1992. [146] W. Hu, M. A. Saberi, A. Jakalian, and D. B. Jack. J. Chem. Phys., 106:2547, 1997. [147] W. A. Steele. The Interaction of Gases with Solid Surfaces. Pergamon Press, Oxford, 1974. [148] A. K. Sallabi. PhD thesis, Concordia University, Montreal, 2001. [149] M. Trabelsi. private communication, 2004. [150] V. E. Klymenko and V. M. Rozenbaum. J. Chem. Phys., 110:5978, 1999. [151] A. D. Bruce and N . B. Wilding. Phys. Rev. E, 60:3748, 1999. Bibliography 125 [152] E. Bayong and H . T. Diep. Phys. Rev. B, 59:11919, 1999. [153] M. E. Fisher and A. E. Ferdinand. Phys. Rev. Lett, 19:169, 1967.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Monte Carlo studies of phase transitions in open systems
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Monte Carlo studies of phase transitions in open systems Daub, Christopher David 2004
pdf
Page Metadata
Item Metadata
Title | Monte Carlo studies of phase transitions in open systems |
Creator |
Daub, Christopher David |
Date Issued | 2004 |
Description | This thesis consists of the investigation of two different problems involving phase transitions using grand canonical Monte Carlo simulations. The first problem sought to understand the nature of the condensation phase transition in simple ionic fluid models, and in particular the restricted primitive model. A specific focus of this investigation was the resolution of a discrepancy between calculations of the constant volume heat capacity near the condensation critical point in the canonical ensemble versus the grand canonical ensemble. By doing calculations on a fluid of hard spheres with attractive power law interactions, we have demonstrated that the near-critical heat capacity calculated in the grand canonical ensemble leads to results which are inconsistent with analytical results for the criticality of these fluids. Although the reasons for this discrepancy remain unclear, we believe that it stems from an incomplete consideration of the finite size effects which affect the grand canonical calculation. We then return to the original problem of the criticality of ionic fluids. We have simulated the critical region of a model of charged hard dumbbells, and we have done a mixed field finite size scaling analysis and calculated the canonical heat capacity. We find no qualitative differences between the criticality of the dumbbells and the ionic model, suggesting that the pairing of ions does not affect the critical behaviour. We surmise that the lack of a peak in the canonical heat capacity of ionic fluid models is a consequence of the very slow convergence of this calculation toward the correct results in the thermodynamic limit, and is not a signature of non-Ising criticality. The second problem we investigated is the adsorption of CO₂ on MgO crystals. This study was driven by experimental results indicating a phase transition between a low temperature, ordered monolayer and a high temperature, disordered monolayer. There is also experimental evidence that this disordered phase has a higher density than the ordered phase. We were able to simulate the formation of an unusual disordered high density phase, however, some discrepancies between the experimental results and our simulations remain to be resolved. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-24 |
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. |
IsShownAt | 10.14288/1.0059473 |
URI | http://hdl.handle.net/2429/17348 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2005-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2005-994643.pdf [ 11.55MB ]
- Metadata
- JSON: 831-1.0059473.json
- JSON-LD: 831-1.0059473-ld.json
- RDF/XML (Pretty): 831-1.0059473-rdf.xml
- RDF/JSON: 831-1.0059473-rdf.json
- Turtle: 831-1.0059473-turtle.txt
- N-Triples: 831-1.0059473-rdf-ntriples.txt
- Original Record: 831-1.0059473-source.json
- Full Text
- 831-1.0059473-fulltext.txt
- Citation
- 831-1.0059473.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-0059473/manifest