Transactions of the Materials Research Society of Japan 30[3] 909-912 (2005)

Atomistic Studies on High-Temperature Elasticity and
Auxetic Behavior in Si0 2 Polymorphs
Hajime Kimizuka*•t, Shigenobu Ogata+•* and Yoji Shibutani*•+• Department of Mechanical Engineering, Graduate School of Engineering, Osaka University,
2-1 Yamadaoka, Suita, Osaka 565-0871, Japan
Fax:+81-6-6879-4121, e-mail: kimizuka@comec.mech.eng.osaka-u.ac.jp
t Engineering Technology Division, The Japan Research Institute, Limited, 16 Ichiban-cho, Chiyoda-ku,
Tokyo 102-0082, Japan:1 Science and Technology Center for Atoms, Molecules and Ions Control, Graduate School of Engineering,Osaka University, 2-1 Yamadaoka, Suita, Osaka 565-0871, Japan

We have investigated the topological structure of three-dimensional network glass, such as
vitreous Si02, and its effect on the elastic properties with a molecular-dynamics (MD) approach.
The topological analysis based on the graph theory is employed to characterize disordered networks in the model systems of vitreous Si0 2 , and thus the results concerning the structure ofn-membered rings are presented. The adiabatic elastic constants of vitreous Si0 2 are evaluatedat finite temperatures, with the fluctuation formula for the internal stress tensor using the MD
method. The quantitative relation between the statistics of rings in vitreous Si0 2 and the elastic properties are studied. In addition, the corresponding MD results for the high-temperatureelasticity and auxeticity of crystalline Si0 2, a- and ,8-cristobalite, are presented for comparison.
Keywords: negative Poisson's ratio, elastic constant, molecular dynamics, cristobalite, vitreous Si02

1 INTRODUCTIONof the principal network glass-forming systems, has
a topologically-disordered network structure and exhibits an amorphous state.
It is believed that many important macroscopic
properties of materials, including elastic and mechanical properties, depend on the underlying microscopic
structures. The main objective of the present paper is to provide results of the molecular-dynamics(MD) simulations that can be used to investigate
the microscopic structures of cristobalite and vitreous Si02, and their effect on the macroscopic elastic
properties, in particular, the high-temperature elasticity and auxeticity. Using the MD configurations,
we also investigate the medium-range order structureof the amorphous state in vitreous Si02, such as then-membered rings and their distribution.

One of the most important properties of vitreous
Si02 is its extremely low coefficient of thermal expansion. This makes the material particularly usefulfor optical flats, mirrors, furnace windows and critical
optical applications which require minimum sensitivity to thermal changes. It is particularly known thatdevitrification and particle generation due to the nucleation of cristobalite are limiting factors in the hightemperature performance of vitreous Si02. In orderto properly assess the mechanical degradation of vitreous Si02 under high-temperature conditions, it isimportant to understand the elastic behavior of bothcristobalite and vitreous Si02 over a wide temperature range.
Cristobalite, a low density polymorph of Si02,
shows an anomalous elastic property of a negative
Poisson's ratio [1, 2]. Negative Poisson's ratio materials undergo lateral contraction under longitudinal
compression, and also lateral expansion under longitudinal tension, so-called "auxetic" behavior [3]. In
cristobalite, the negative Poisson's ratio appears overthe entire temperature range of 300-1800 K, where
the bulk modulus is extremely low compared withthe shear modulus [4]. Markedly different mechanicalproperties are found between the two silica, includingthat cristobalite is auxetic whereas the vitreous Si02
is not.
Silica exists in a variety of topologically-different
crystalline forms, such as quartz, cristobalite,
tridymite, coesite, etc. The elementary unit in thesepolymorphs at ambient pressure is an almost idealSi(0 1; 2)4 tetrahedron. Their structures are made up
of different arrangements of corner-sharing tetrahedral elementary units. Vitreous Si02, which is one

2 COMPUTATIONAL METHOD
To obtain the model system of vitreous Si02, equilibrium molecular dynamics (MD) simulations havebeen performed within the Nose-Hoover [5, 6] andParrinello-Rahman [7] isothermal-isobaric ensemble.
We start from the model system of a-cristobalite [9],and heat the system gradually from 300 K to 5000 K
at the rate of 50 K/ps, after which the system is equilibrated over 80 ps. The integration time step in all
MD simulations is 2 fs. The amorphous state is obtained from a well-thermalized molten state at 5000
K by quenching to room temperature 300 K, at therate of -50 K/ps. In the quenched state, atoms only
vibrate around their equilibrium positions, and no diffusion occurs on the time scale of the MD runs. The
fluctuation formula for the internal stress tensor [8] is
adopted to obtain a complete set of adiabatic elastic

909

Atomistic Studies on High- Temperature Elasticity and Auxetic Behavior in Si02 Polymorphs

910

constants at finite temperatures. In order to checkthe validity of statistical sampling, the averages are
separately taken every 40 ps during a 0.64 ns run.
We have employed the nonempirical pairwise potential by van Beest, Kramer, and van Santen [10](BKS), derived from the cluster calculations of potential energy surfaces using the ab initio method.
This potential has been found to adequately describe
a wide scope of properties of Si02 through numerous previous studies. Here, to eliminate spurious effect, driven by the dispersion terms which overwhelm
all the other terms at very short distances, the original BKS interatomic potential is modified by adding
a Lennard-Jones (24-6) contribution [11]. The functional form of the modified version of BKS potential(BKSm) is expressed as....."~.
....

i\
I

1•

I

,""'•\_,',.\,)--------- _:t~~_.lf\-------- __:''·\

a-cristobalite\,..•·~, .. ••,.,.,"''"',,'", ... ~·· • • • • ., ... •• ..

P-cristobalite........... . ... -·-,···-......- ...····--····-·····

liquid······----------/..~-...•....._.......... ~--------·-····-----·-·-··- ..-------···

vitreous silica (3)0.0 0.1 0.2 0.3 0.4 0.5 0.6 0. 7 0.8 0.9 1.0
lnteratomic distance [nm]where u(rab) is the interaction energy between the
two atoms a and b (a,b E {Si,O}), rab is the interatomic separation, and (za, Zb, aab, bab, Cab) areparameters fitted to Hartree-Fock calculations for
H4Si04 molecules as basic building blocks of various phases of Si02. The force-field parameters andthe adjustment parameters (Eab and CTab) used in thiswork are given in Refs. [10, 11]. The Coulomb interactions are calculated using the Ewald method.
The number of atoms in the system is 576 (containing 48 unit cells), and we use periodic boundary conditions. To estimate the finite-size effects,
the MD calculations are also performed on the 4608and 36864-atom systems (576x2 3 and 576x43 , respectively).

3 RESULTS AND DISCUSSION
3.1

Structure of vitreous Si0 2

To characterize the structures of amorphous statesfrom the MD results, the radial distribution functions (RDFs) for 0-0 distances in crystalline, liquid, and vitreous Si02 were calculated (Fig. 1). In
Fig. 1, we observe that the first-nearest 0-0 distanceremains essentially unchanged during the heating andquenching of the system (between 300-5000 K), andthe rigid Si(0 1; 2)4 tetrahedral units in Si02 structure
are maintained over the phase transformations. However, the peaks, which indicate medium-range order
beyond the nearest neighbors, become much broader
and even decay in the amorphous state. This indicates that oxygen sites are less ordered in vitreous
Si02, and Si(0 1 ; 2)4 tetrahedral orientations (specified by corner oxygen positions) widely distribute dueto the destruction of the long-range periodicity. The
profiles for the three vitreous Si02 models of different system sizes (from 576 to 36864 atoms) look sosimilar that they are hardly distinguished each other.
In the network-forming structure such as vitreous
Si02, it is important to investigate the character anddynamics of the underlying microscopic structures. It
is possible to describe the nature of the connectivity

Fig. 1: Radial distribution functions (RDFs) of
crystalline, liquid, and vitreous Si02 for 0-0 distance. Solid curves denote the RDFs for the(1) 576-atom, (2) 4608-atom, and (3) 36864-atomsystems of vitreous Si0 2, respectively. The RDFs
of crystalline (a- and ;'J-cristobalite) and liquid
Si0 2 are represented by dashed curves.
of the elementary units, beyond the nearest neighbor, in terms of n-membered rings. Fig. 2 shows then-membered ring structures of different sizes, wheren = 4 and 6. They are defined as the closed paths containing n Si-0 segments, and black and gray spheres
represent silicon and oxygen atoms, respectively. The
medium-range order structure in crystalline and disordered networks can be characterized in terms of
their ring structures and distributions derived witha shortest-path analysis (see, for example, Ref. [12]).
Rings defined by the method are a convenient way
to analyze Si02-type covalent structures. If a silicon
atom has four oxygen nearest neighbors, then there
are six 0-Si-0 paths emanating from each Si atom.
Shortest-path analysis involves taking every pair ofbonds on each silicon atom in turn and finding thesmallest size ring. For a four fold-coordinated network such as Si02 at ambient pressure containing Ns;
silicon atoms, there are 6Ns; rings. An n-memberedring therefore contains 2n Si-0 bonds. Using thisring analysis, one finds that cristobalite consists entirely of six-membered rings, whereas quartz has sixmembered and eight-membered ring.
Table I shows the statistics of rings of sizes from 2
to 9. The listed values exhibit the population of each
ring, normalized by the number of silicon atoms ( Ns;)in the system. Note that the three systems of vitreous
Si02 have almost the same ring distribution. In allcases, the distribution is peaked at a six-memberedring, and nearly symmetric around the value. Also,
the population of five-membered rings is slightly but

Hajime Kimizuka et al.

Transactions of the Materials Research Society of Japan 30[3] 909-912 (2005)

Table I: Statistics of n-membered rings per Si(0 1; 2 )4 tetrahedron in vitreous Si02.
Ring size (n)a- and ,6-cristobalite
vitreous Si02 (N = 576)
vitreous Si02 (N == 4608)
vitreous Si02 (N == 36864)

2

3

4

5

0.05
0.04
0.05

0.05
0.09
0.14

0.68
0.69
0.85

2.06
1.78
1.90

{a)

Fig. 2: Examples of n-membered rings for (a) n =4 and (b) n = 6. Black and gray spheres represent silicon and oxygen atoms, respectively. Dotted lines within the rings indicate the Si(0 1 ; 2) 4
tetrahedral sub-units.
clearly larger in the 576-atom system, compared withthe others. It suggests that the relative populationsof rings are not well described in a smaller system,
due to the limitation of size.
3.2

Elastic properties of cristobalite and vitreous Si02

The adiabatic elastic stiffness constants ( O;j) ofvitreous Si02 were evaluated using the MD simulations with the BKSm potential. The obtained O;jvalues are given in Table II, for the systems of 576,
4608, and 36864 atoms. We have included the experimental values (1] for comparison, as well as theresults of previous MD calculations (4] for a- and ,6cristobalite using the original BKS potential.
In the calculations of elastic constants of vitreous
Si02, we observe the three finite elements, 0 11 , 044,
and 012· The three components can be reduced totwo based on the equation 044 == ~(On- 0 12), in
the isotropic media such as vitreous Si02. In Table II, we confirm that this isotropic solid criterion
is well satisfied for the three vitreous Si02 models.
Also, the bulk component On is much larger than
the shear component 044 in vitreous Si0 2 , as in ordinary solids, whereas On is smaller or almost equal
to 044 in a- and ,6-cristobalite. This fact suggeststhat the peculiar effect of microscopic structures on
the macroscopic elasticity exits in the crystalline Si02
form, cristobalite.
In ,6-cristobalite, it is noteworthy that the 0 11 exhibits extremely close value to the 044, in contrast
to the Cauchy relation. This result is unusual in
the ordinal cubic crystals, whereas the well-known
Born stability criteria for cubic crystals are well sat-

6
6
2.17
2.16
2.18

7

8

9

1.27
1.34
1.34

0.35
0.40
0.39

0.03
0.04
0.03

isfied [4]. Thus we predict that the elasticity of {3cristobalite can be described by only two quantities;
however, it clearly differs from the isotropic media
such as vitreous Si0 2 , as in Table II.
The bulk modulus (K) and shear modulus (G)for an isotropic aggregate of cristobalite and vitreous Si02 were evaluated from the O;j 's by means of
the Voigt-Reuss-Hill average. Table Ill shows K, G,
and the Poisson's ratios (v) of the single-phased aggregate of cristobalite and vitreous Si02. Here, thesign of a Poisson's ratio is determined by the relationship between K and G, as in the following formula:
V== (1/2)(3K- 2G)/(3K +G).
In vitreous Si0 2, the bulk modulus is extremely
high compared with the shear modulus, and thus the
Poisson's ratio exhibits a positive value as in ordinary
solids. Also, while the shear modulus exhibits almost
the same value in the three amorphous systems, thebulk modulus is clearly different between the 576atom system and the others. Since these three systems are found to have almost the same density in
the MD calculations (p == 2.69, 2.68, and 2.71 gjcm 2
for the 576-, 4608-, and 36864-atom systems, respectively), it is predicted that this result is due to thedifference in the medium-range order structure, especially in the population of five-membered rings. It issuggested that the large population of rings smaller
than size 6 brings relatively high incompressibility in
Si02-type covalent structures, due to the limitationof the cooperative motions of Si(0 1; 2)4 tetrahedra.
In a-cristobalite, the negative Poisson's ratio results from the typical framework structure of acristobalite which consists of uniform six-memberedrings [1, 4]. These rings have an inverted characteristic similar to re-entrant honeycombs, the typicalmodel which induces a negative Poisson ratio [13, 14].
Here, the bulk modulus is extremely low comparedwith the shear modulus, as in Table III. We predictthat the softening of the bulk modulus in cristobalite
occurs due to the cooperative motions of Si(01;2)4
tetrahedra, whereas the anomalous properties are localized due to the destruction of long-range periodicity in vitreous Si02.

4 SUMMARY
This paper presents the molecular-dynamics (MD)results concerning the structural and elastic properties of cristobalite and vitreous Si02, the typical polymorphs of Si02. Using the equilibrium MD methodwith a stress-fluctuation formula, we have successfully
evaluated the adiabatic elastic constants of (a,(3)cristobalite and vitreous Si02 at finite temperatures.
The nature of connectivity of the elementary units

911

912

Atomistic Studies on High-Temperature Elasticity and Auxetic Behavior in Si02 Polymorphs

Table II: Adiabatic elastic constants (in GPa) of cristobalite and vitreous Si0 2 .

Cn
c33
c44

c6a
c12
c13
l.(Cn- C12)
aReference [1].
bReference [4].

a-cristobalite
Expt.a BKSb
59.4
64.5
42.4
37.9
67.2
69.5
25.7
27.6
3.8
6.5
-4.4
-0.7

,8-cristobalite
BKSb
49.6

vitreous Si02 (300 K)
N=576 N = 4608 N = 36864
113.0
135.3
120.4

49.7

40.1

41.8

38.8

5.7

44.4

36.7

36.0

45.5

38.2

42.2

Table III: Isotropic bulk modulus (Kin GPa), shear modulus (Gin GPa),and Poisson's ratio (v) of cristobalite and vitreous Si02.
a-cristobalite
Expt.a BKSb
K
16.37
18.8
G
39.05
39.4
V
-0.163 -0.12
aReference [1].
bReference [4].

,8-cristobalite
BKSb
20.3
35.8
-0.05

beyond the nearest neighbor is characterized in termsof the ring distribution by shortest-path analysis, andthe results can be associated with the medium-rangeorder structure of amorphous state in vitreous Si02.
It is suggested that the microscopic structures invitreous Si02 affect the elastic properties directly.
Also, we have investigated the characteristic hightemperature elasticity of cristobalite, using the MD
simulations based on the BKS pairwise potentials.
We confirm that cristobalite exhibits a negative Poisson's ratio over the wide temperature range of 3001800 K, whereas vitreous Si02 exhibits the positive
value at 300 K. To this end, our simulations provide
the microscopic insight into the intimate relation be·tween the elastic behavior and the nature of connectivity of the elementary Si(0 1; 2)4 units, in the Si02
polymorphs.

vitreous Si02 (300 K)
N = 576 N = 4608 N = 36864
74.3
61.9
63.9
41.9
40.2
39.8
0.26
0.23
0.24

[5] S. Nose, Mol. Phys. 52, 255 (1984).
[6] W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
[7]. M. Parrinello and A. Rahman, J. Appl. Phys.
52, 7182 (1981).
[8] J. R. Ray, M. 0. Moody, and A. Rahman, Phys.
Rev. B 32, 733 (1985).
[9] J. J. Pluth, J. V. Smith, and J. Jr. Faber, J.
Appl. Phys. 57, 1045 (1985).

[10] B. W. H. van Beest, G. J. Kramer, and R. A. vanSanten, Phys. Rev. Lett. 64, 1955 (1990).
[11] Y. Guissani and B. Guillot, J. Chem. Phys. 104,
7633 (1996).

(12] J. P. Rino, I. Ebbsjo, R. K. Kalia, A. Nakano,and P. Vashishta, Phys. Rev. B 47, 3053 (1993).
[13] R. F. Almgren, ]. Elasticity, 15 427 (1985).

ACKNOWLEDGEMENTS
H.K. would like to thank Hideo Kaburaki and
Yoshiaki Kogure for their support and encouragements during the research.

REFERENCES[1] A.

Yeganeh-Haeri, D. J. Weidner,
J. B. Parise, Science 257, 650 (1992).

and[2] N. R. Keskar and J. R. Ohelikowsky, Nature(London) 358, 222 (1992).

[3] K. E. Evans, M. A. Nkansah, I. J. Hutchinson
and S. 0. Rogers, Nature (London) 353, 124(1991).
[4] H. Kimizuka, H. Kaburaki, and Y. Kogure,Phys. Rev. Lett. 84, 5548 (2000).

[14] L. J. Gibson and M. F. Ashby, Cellular Solids:
Structure & Properties (Pergamon, Oxford,
1988).
(Received December 23, 2004; Accepted April 28, 2005)