PRESSURE INDUCED PHASE TRANSFORMATIONS IN SILICA


TAHIR CAGIN, ERSAN DEMIRALP, and WILLIAM A. GODDARD, III

Materials and Process Simulation Center,
Beckmann Institute, (139-74),
California Institute of Technology, Pasadena, CA 91125

ABSTRACT

Silica, SiO 2, is one of the most widely studied substance, and it has some complex and unusual properties. We have used a recently developed 2-body interaction force field [1] to study the structural phase transformations in silica under various pressure loading conditions. The specific transformations we studied are alpha-quartz to stishovite, coesite to stishovite and fused glass to stishovite-like dense, a dominantly six-coordinated glassy phase. Molecular dynamics simulations are performed under constant loading rates ranging from 0.1 GPa/ps to 1.0 GPa/ps, pressures upto 100 GPa and at temperatures 300, 500, and 700 K. We observe the crystal to crystal transformations to occur reconstructively, whereas it occurs in a smooth and displacive manner from glass to a stishovite-like phase confirming earlier conjectures.[2] To elucidate the shock loading experiments, we studied the dependence of transition pressure on the loading rate and the temperature. To assess the hysterisis effect we also studied the unloading behavior of each transformation.

INTRODUCTION

Despite its simple chemical formula, SiO2, silica has some of the most complex, interesting and unusual properties of any material. It is ubiquitous in the earths crust, and silicon and oxygen are abundant elements in the solar system. The manner in which silica responds to pressure and stresses also touches many fields of study, from geology to various areas of technology, from making simple tools to glass fibers to micro- and nano-electronic device design, processing and manufacturing. Also due to rich variety of its polymorphs and amorphous states, it also serves as a model system for studies of high pressure structures and structural phase transformations.

Over the last decade increasing number of researchers from different fields of science attempted to investigate various aspects of silica, structure, melting, phase behavior, energetics using molecular mechanics, molecular dynamics and monte carlo techniques. If silica is exposed to shock wave, various low density silica polymorphs and its glass transform to high density stishovite along the Hugoniot isentrope. Several years ago Stolper and Ahrens[2] pointed out a difference in the alpha-quartz to stishovite vs. glass to stishovite transformation, "this transformation is more displacive than reconsructive". In this study, we aim at shedding some light into this sound conjecture by analizing the microstructures of the model systems along the transformation path.

Central to all classical computational approaches is the representation of atomic interactions. We recently developed a simple 2-body interaction potential with long range electrostatic interactions are handled via variable charge model whereas a short range morse term supplies the combination of valence and van der Waals interactions.[1] The partial charges are determined by the charge equilibration method of Rappe and Goddard.[3] We have succesfully applied this potential to characterize silica polymorphs and especially the microstructure of glassy form. [4]

Our objective is to present a systematic study on crystal-crystal, glass- crystal transformations observed in silica and assess the role of each factor, such as temperature, pressure loading rate on the transition pressure. Also using the microscopic level information to investigate the details of transformation to relate to the experiments, which are mostly the shock wave experiments at these pressures, and loading rates. Furthermore, to develop an understanding of behavior of various silica polymorphs and glasses subjected to high pressures.

In the Theory and Method section, we briefly describe the interaction potentials used in this study, and procedures used in molecular mechanics and molecular dynamics simulations. In the Results section, we present and discuss the simulation results we obtained.

THEORY AND METHOD

Interaction Potentials for Silica

Silica can exist in several crystalline polymorphs as well as in the amorphous state. A number of interatomic potential functions have been developed which can, with reasonable accuracy, reproduce the structural characteristics of one or more of the silica polymorphs. For studying silica most widely used forms of interaction potentials have a 2-body Born-Mayer-Huggins[4-8] term along with an attractive dispersion term, i.e. an exp-6 potential. However they may slightly differ in handling the electrostatic interaction, screened coulomb (erfc -complementary error function) to formal charge models. Some model potentials also include a Stillinger-Weber type 3-body term. [9]

We developed the new pure 2-body interaction potential for simulation of all phases of silica, glasses and melts and extended it to alimuna as well. This new MS-Q interatomic potential differs from the previous potentials by allowing the charges to adjust self-consistently to changes in the environment. These charges are determined via the Charge Equilibration (Qeq) method of Rappe and Goddard.

A Morse type function is assumed for short range bonding forces. Thus the MS-Q two body interatomic potential function has the form:

U(Rij) = qi qj/ Rij + Do [ exp { gamma ( 1-Rij / Ro ) } - 2 exp { gamma ( 1-Rij / Ro )/2 } ]

where qi are the charges and Rij is the distance between atoms i and j. Do, Ro and gamma are adjustable parameters.

The charges for the Si atoms in the optimized silica polymorphs varied between 1.22 and 1.4. This compares with potential derived charges and Mulliken charges obtained from a Si(OH)4 cluster ab initio Hartree Fock calculation using a 6-31G** basis set of 1.528 and 1.493 respectively. Thus, the Qeq charges appear reasonable. The Qeq charges for O vary from -0.61 to -0.70.

In our model glass simulations with this potentials, we found the pair distribution functions, g(r), and angular distribution functions for O-Si-O and Si-O-Si are in good agreement with experiment.[4]

Model Systems and Simulation Procedure

In our studies we employed both molecular mechanics (for optimization of the structure and lattice parameters/volume, i.e. objective function is the enthalpy) and variable cell Parrinello-Rahman-Nos\' e molecular dynamics. Model systems are supercells of unit alpha-quartz, coesite and stishovite, and glasses which were obtained from melt through lengthy quench procedure with sizes ranging 576 atom to 640 atoms. Energy, force and stresses in the simulations described below all employed Ewald summation method.

In molecular mechanics simulations we did not impose any crystal symmetry constraints. We used an external hydrostatic pressure, and the pressure on the model systems were varied by 1 GPa intervals upto 10 GPa, 2 GPa intervals upto 30 GPa, and 5 GPa intervals beyond 30 GPa. In order to explore the release ``isentrope'' we also re-traced the loading curve backwards for stishovite and alpha-quartz cases only. Termination criteria for optimization was 0.1 kcal/mol/A. The results of these are presented in the next section.

Molecular dynamics simulations are performed under constant temperature condition, in contrast to shock wave experiments. Furthermore, the imposed pressure load or profile was isotropic rather than uniaxial. In these simulations 1 femtosecond is used as the integration step. In the course of simulations the pressure is incremented with a specific period ( step function compressions) to give approximate pressure loading rates 0.1, 0.25, 0.5, and 1.0 GPa/ps. The end points of each pressure segment were further continued for additional 25 ps subjected to the same pressure in order to study the micro structure and obtain averages.

RESULTS

Molecular Mechanics simulations are used to obtain 0 K P-rho ``isotherms'' for alpha-quartz, coesite, stishovite and glass. The results are plotted in Fig. 1.a Since we used the second derivative matrix, Hessian, to obtain the higher order mechanical properties; using a smaller size (maximum 2x2x2 supercell) models, we also calculated the Birch-Murnaghan equation of state parameters for the crystalline phases only. Bulk modulus, B, and its pressure derivative, B', at zero pressure. For alpha-quartz, B=36.4 (37) GPa, for coesite, B=85.3 (89) GPa and for stishovite B=290.7 (298) GPa, values in parantheses are the experimental values.

Figure 1. Molecular Mechanics simulation of Equation of State for silica polymorphs and a model silica glass: a) Pressure vs. Density for alpha-quartz , coesite, stishovite and glass models. b) Loading and release isotherms (at 0 K) for quartz and stishovite

The pressure-density relationships reveal sharp transitions from alpha-quartz to stishovite (with a number of defects) at 24 GPa, and at 20 GPa for coesite, whereas transition in glass is gradual, as the pressure is increased the 5- 6- coordinated local structures slowly emerges. Even at a pressure as high as 100 GPa, the ratio of 6-coordinated atoms are less than that of structures obtained from quartz and coesite. When the loading curves retraced for stishovite and alpha-quartz and the transformation to stishovite from alpha-quartz is found to be not reversible. Resulting structure has a density of 3.7 g/cm3. From shock wave release isentropes similar observations are made.[10]

We studied the rate and temperature dependence of the phase transformation using isothermal-isotension molecular dynamics, to mimic a non-equilibrium process where in the course of simulation the system is subjected to increasing isotropic pressure as step functions. The time interval between steps are varied to investigate the loading rate dependence of pressure loading. The rates are varied from 0.1 GPa to 1.0 GPa per picoseconds. The simulations for 1 GPa/ps rate were repeated at two higher temperatures 500 and 700 K, to investigate the temperature dependence of the transition preassure in the case of alpha-quartz. These temperatures are considerably lower than the temperatures observed in shock wave experiments. The results of pressure loading rate and temperature effects are displayed in Fig. 2.a and 2.b respectively.

Figure 2. Temperature and pressure loading rate dependence of transformation pressure for alpha-quartz to stishovite; a) Temperature dependence, pressure loading rate is 1.0 GPa/ps, b) Pressure loading rate dependence, T = 300 K

The temperature dependence of transformation pressure is stronger than the rate dependence, at least for the rates we have utilized in this study. The phase transition pressure changes by about 3 GPa for a temperature increase of 400 K. Experimental Hugoniot data show evidence of transition starting at about 14 GPa for 500 K [10], which is in agreement with our result, 15-16 GPa at 500 K.

Finally, to check the conjecture of Stolper and Ahrens, we analyzed the trajectories for microstructural evolution. alpha-quartz or glass to stishovite transition can easily be characterized by the variation of the distribution of O-Si-O angles in the structure. This angle is 100 % tetrahedral for quartz, dominantly tetrahedral, but not all, for glass structure, however in stishovite this angle is octahedral, i.e. 90 degrees. In Figure 3.a and 3.b we present the variation of the O-Si-O angles for quartz and silica glass. We see an abrupt change in the distribution for quartz whereas the change in glass is continuous even at pressures upto 120 GPa. The angle distribution for stishovite like structure obtained from the model glass is almost twice broader at the octahedral angle value. As a note, the peaks at 45 and 135 degrees are due to the silicon bonds to oxygens on the bisectors of 90 degree O-Si-O angles. We found that the transformations from alpha-quartz to stishovite are reconstructive whereas the transformation in the glass is displacive, with the peaks of the angle distributions progressively approaching to the 6-fold coordinated values which correspond to stishovite.

Figure 3. Microscopic analysis of pressure induced transformations using O-Si-O angle distribution as a precursor: a) alpha-quartz, b) model silica glass

CONCLUSIONS

We have used a recently developed interaction potential to study the pressure induced structural transformations in silica, using molecular mechanics and molecular dynamics techniques. The effects of hydrostatic pressure, hydrostatic pressure loading rates and the temperature are systematically investigated. We identified that the transition pressure is more strongly dependent on the temperature at the loading rates we used in our simulation. This is plausible since the crystal to crystal transformation is require an activation energy which is more readily available at higher temperatures. We found good agreement with the shock wave transition pressure data at 500 K. The transition pressure will be further depressed to 7-9 GPa range as in some shock experiments where measured shock temperatures are closer to 2000 K.[11] As a caution these computational experiments based on described molecular dynamics procedure by no means mimicking the conditions, and thermodynamics of a shock experiment. However, it clarifies the progression of phase transitions under hydrostatic pressure and isothermal conditions. With our simulations we confirmed the validity of the conjecture by Stolper and Ahrens for the glass to stishovite transformation by investigating the microstructure of the model systems as a function of increasing pressure load. We found that the transformations from alpha-quartz to stishovite are reconstructive whereas the transformation in the glass is taking place in a displacive manner, with the peaks of the angle distributions progressively approaching to the 6-fold coordinated values as in stishovite.

ACKNOWLEDGEMENT

This research was supported by the NSF (CHE 95-22179 and ASC 92-17368). The facilities of the MSC are also supported by grants from DOE-AICD, Chevron Petroleum Technology, Asahi Chemical, Chevron Chemical Co., Owens Corning, Aramco, Seiko-Epson, Avery Dennison, Chevron Refinery Technology, and the Beckman Institute.

REFERENCES

  1. E. Demiralp, T. Cagin, W.A. Goddard, III, MRS Spring 1997 Meeting, unpublished.
  2. E.M. Stolper and T.J. Ahrens, Geophys. Res. Let. 14, 1231 (1987).
  3. A. K. Rappe and W.A. Goddard, III, J. Phys. Chem. 89, 8900 (1990).
  4. E. Demiralp, T. Cagin, W.A. Goddard, III, MRS Spring 1997 Meeting, unpublished.
  5. S. Tsuneyuki, M. Tsukada, H, Aoki, Y. Matsui, Phys. Rev. Lett. 61, 869 (1988).
  6. A.B. Belonoshko, L.S. Dubrovinsky, Geoch. Cos. A. 59, 1883 (1995).
  7. Nakano-A Kalia-RK Vashishta-P, Phys. Rev. Lett. 73, 2336 (1994).
  8. B.W.H. van Beest, G.J. Kramer and R.A. van Santen Phys. Rev. Lett. 64, 1955 (1990)
  9. S.H. Garofalini, J. Am. Ceram. Soc. 67, 133 (1984).
  10. J. W. Swegle, J. Appl. Phys. 68, 1563 (1990).
  11. R. F. Grieve, F. Langenhorst, D. Stoffler, Meteor. Planet. Sci. 31, 6 (1996).