Transformations in the medium-range order of fused silica under high pressure

L lian P. D vila,1,2,* Maria-Jos Caturla,2,3 Alison Kubota,2 Babak Sadigh,2

Tom s D az de la Rubia,2 James F. Shackelford,1 Subhash H. Risbud,1

and Stephen H. Garofalini 4

1

Dept. of Chemical Engineering and Materials Science, University of California Davis, Davis, CA 95616

2

Lawrence Livermore National Laboratory, Livermore, CA 94550

3

Dept. Fisica Aplicada, Universidad de Alicante, Alicante, Spain, E-03690

4

Department of Ceramic and Materials Engineering, Rutgers University, Piscataway, NJ 08854

(Dated: September 23, 2003)

Molecular dynamics simulations of fused silica at shock pressures reproduce

the experimental equation of state of this material and explain its characteristic

shape. We demonstrate that shock waves modify the medium-range order of this

amorphous system, producing changes that are only clearly revealed by its ring size

distribution. The ring size distribution remains practically unchanged during elastic

compression but varies continuously after the transition to the plastic regime.

PACS numbers: 61.43.-j, 64.30.+t, 02.70.Ns, 78.55.Qr

High-pressure experiments reveal an unusual equation of state (EOS) of

amorphous silica, displaying an initial increase in compressibility with pressure and

irreversible densification [1-2]. Shock loading and neutron irradiation [3-5] can also

result in higher density silica glass. Under hydrostatic loads at ambient

temperatures, the compaction begins at 8-10 GPa, below which the compression

remains elastic [6-7]. Samples recovered from shock compressions at 10-16 GPa

show evidence of permanent densification and a transition to a high-density form at

about 16 GPa [8]. Experimental studies have revealed large (20-30%) permanent

densifications under high pressures [9-10].

Electronic address: abqkzp@r.postjobfree.com

1

The densification of amorphous silica has been studied by many groups over

several decades [2-3, 11-12]. However, only in the 1990s was it suggested that the

high-density forms are unique in their structure [7,13]. Structural changes were

revealed by in situ vibrational spectroscopy measurements [14-16] and X-ray

diffraction [17]. The increase in density was attributed to the rotation of the silica

tetrahedra. X-ray [18] and neutron [9] diffraction, infrared and Raman spectroscopy

revealed a decrease in the average Si-O-Si angle and a small increase in the Si-O

bond distance [13]. Moreover, X-ray diffraction measurements [17] have shown that

the characteristic first sharp diffraction peak (FSDP) associated with the medium-

range order (MRO) of these materials is significantly reduced upon densification.

The structural origin of the MRO and the dynamic behavior of silica glass at high

pressures have been studied extensively using molecular dynamics (MD)

simulations [9, 19-22] but a completely satisfactory understanding has not yet been

achieved. Recent papers have shown how the medium-range order can be

manipulated experimentally [23-24], resulting in materials of scientific and

technological interest, motivating further understanding of the features that

determine this ordering.

We have used MD simulations with empirical potentials to study the

behavior of fused silica under high pressure. The initial amorphous structure was

generated using the interatomic potential reported in Webb and Garofalini [25].

We present two different simulation schemes for studying the structural

transformations in fused silica under high pressures, and show that they give similar

P-V relationships. In one case, a shock wave is generated in a slab containing

240,000 atoms (two free surfaces and periodic boundary conditions in the lateral

directions) at 300 K, through a piston, consisting of a few designated atoms at one

2

end of the cell set to a constant velocity. This method (noted as MI) was previously

applied to the study of shock propagation in crystalline materials [26]. The shock

simulation is performed typically until the shock wave reaches the free surface at the

other end of the simulation box. With this method the propagation of the shock

wave can be followed and its velocity measured. The results can therefore be

directly compared to flyer plate experiments.

The second simulation scheme (or MII) involves small periodic supercells

containing 1,536 atoms at 300 K, that are compressed by reducing the volume of the

supercell and allow for relaxations. The temperature in these simulations was kept

constant at 300 K. A comparable method has been applied to study compression in

vitreous silica [21] using minimization techniques.

The initial fused silica structure was obtained through a MD melt-quench

technique similarly used in previous studies [25,27]. Periodic cells in the cubic

-cristobalite structure were melted at 7000 K for 25 ps (with a MD time step of 0.5

fs). Thereafter, they were quenched down to room temperature by a series of steps.

Our starting configuration has a density of 2.205 g/cc, that is, equivalent to the

experimental value for fused silica.

Using these starting configurations we applied methods I and II described

above to compute the pressure-volume relationship. In the case of MI, different

initial velocities of the piston were used to obtain different shock velocities. The

final volume and pressure was then calculated after the shock reached the back

surface. For the MII simulations, the volume was reduced sequentially in decrements

of 2% and pressure was calculated upon reaching equilibrium, typically after 50 ps.

Subsequent compression simulations were performed down to 0.58 V/Vo, where Vo is

the starting volume.

3

Shock Experiments (Sugiura et al)

1

I

Shock Simulations (M)

II

Compression Simulations (M )

0.9

0.8

V/V0

0.7

0.6

0.5

0 5 10 15 20 25

Pressure (GPa)

FIG. 1. Pressure-volume behavior calculated from MD

simulations compared to experiments [4].

Figure 1 depicts the resulting EOS curve in the form of pressure versus

relative volume, V/Vo. This figure shows the experimental results (filled circles)

obtained by Sugiura et al. [4] using flyer plate experiments, the results from shock

simulations (crossed squares) and the results from continuous compression

calculations (open circles). The agreement between these MD simulations and the

shock experiments is remarkable, particularly considering that the interatomic

potential used was fitted only to ambient pressure and temperature conditions. It is

also interesting to note the similar results between the two different types of

calculations, revealing the slow dynamics and large relaxation times associated to

this system. The high strain rates of our compression simulations imply that the

conditions are closer to those found in shock experiments than in quasi-static

experiments [6-7].

We can clearly identify two regions in the EOS curve in Fig. 1 associated with

1) linear compressions up to 0.82 V/Vo or elastic behavior and 2) asymmetric

compressions starting from 0.82 V/Vo or plastic response. It is important to point out

4

that irreversible compaction of the glass is found in the range 9-10 GPa [6-7]. The

plastic region occurs in the regime between 10-23 GPa, as shown in Fig. 1, which

falls in the range of pressures known as transformation region [6]. We find that

our simulations can accurately predict the densification of fused silica at various

shock pressures. For the pressure ranges studied here values of compression up to

42% are reached during shock loading. This densification however is not

permanent. We have performed relaxation of the dynamically loaded structure for

the highest densification (42%). Upon relaxation the lattice expands to a final

permanent densification of 20% in agreement with experiments [4-5,10].

In order to understand the structural changes occurring as the pressure

increases, and to correlate them with the calculated EOS, we have performed a

detailed analysis of the structure factor, the pair distribution function, the ring size

distribution and the bond angle distribution of the system. In particular, the

structure factor can be measured experimentally through X-ray or neutron

diffraction and can provide important information about the medium-range order of

amorphous materials.

5

Neutron Expts (Inamura et al)

Unshocked (or V/Vo = 1.00)

1.25 Km/s (or V/Vo = 0.82)

2 1.5 Km/s (or V/Vo = 0.70)

2.5 Km/s (or V/Vo = 0.583)

FSDP

1.5

S(q)

1

0.5

0

0 1 2 3 4 5 6 7

-1

q (A )

FIG. 2. Structure factor from MI calculations. Experimental

results for zero pressure are also included for comparison.

Unshocked (or V/Vo = 1.00)

1.25 Km/s (or V/Vo = 0.82)

2.5

1.5 Km/s (or V/Vo = 0.70)

2.5 Km/s (or V/Vo = 0.583)

2

S SiSi (q)

1.5

1

0.5

0

0 1 2 3 4

-1

q (A )

FIG.3. Partial structure factor for Si-Si from shock simulations

(MI).

6

One of the interesting features of silica glass is the presence of a characteristic

FSDP as reported in experiments [9,17,28]. Figure 2 shows the structure factor for

different pressures calculated from our MI calculations. Similar results are obtained

from MII. In the same figure we include the experimentally measured S(q) at 0

pressure obtained by Inamura et al. [28] for comparison. Although there is

experimentally available data for the S(q) of silica at higher pressures, these have

been obtained after recovery of the sample or from diamond anvil cell experiments

[9,28]. Therefore direct comparison with the simulations presented here would not

1

be appropriate. The FSDP position occurs in the range 1-2 and, as observed

experimentally [9,17,28], it decreases in intensity and shifts towards larger q values

as the pressure increases. Figure 2 also denotes that short-range order prevails with

increasing pressure in this glass as noted by the third peak in the structure factor

calculations up to about 23 GPa as seen in Fig. 1. The origin of this FSDP is still

controversial and many different models have been proposed, such as the one by

Elliott [29] based on an interstitial-void ordering. In our simulations we observe a

continuous decrease in the intensity of the FSDP with pressure, however no

significant changes are observed at the point of transition between elastic and plastic

behavior. The same behavior is observed in the pair correlation function. Although

changes in the height and width of the peaks are observed with pressure,

comparable to previous studies [9,22], none of these changes reveal clearly the

transition between elastic and plastic regimes, as seen in the EOS.

The total structure factor shown in Figure 2 is a linear combination of the Si-

Si, Si-O, and O-O partial structure factors (PSF). It is thus instructive to study the

evolution of these PSF with pressure. We find that the Si-O and O-O short-range

order is only slightly changed with pressure, however the pre-peak in both cases

follow the same trend as in the total structure factor, i.e. moving out to larger q-

7

values with diminishing intensity. A more interesting and unique behavior is

observed in the Si-Si PSF, where the pre-peak and the first peak are merged together

and eventually only one peak is observed, as shown in Figure 3 for the case of MI.

The origin of this behavior can be traced to changing Si-Si short-range order with

pressure. The real-space Si-Si pair-correlation function reveals that the average

nearest-neighbor distance between Si atoms decreases with pressure, in agreement

with previous calculations [9,18,29], while the corresponding first peak in the Si-Si

pair-correlation function reduces in magnitude and broadens significantly into a

structure with almost double-peak character. From this analysis, however, we still

can not extract a clear correlation between the changes in the structure factor and the

EOS.

Most interesting are those changes revealed by the ring size distribution,

which characterizes the medium-range order of these amorphous materials. Rings

are defined as the shortest closed loop of Si-O bonds. Fused silica consists of a

distribution centered on six-member rings. Figure 4 shows the number of rings

(normalized to their value at zero pressure) as a function of densification for the two

types of simulations presented here. Note the ring evolution (a) for MI and (b) for

MII and for different sizes of rings. The ring distribution was calculated using the

shortest-path criterion [26], assuming that stable bonds are formed at distances

smaller than 2 .

8

6

Normalized Number of Rings

1 0.9 0.8 0.7 0.6 0.5

V/Vo

300

Normalized Number of Rings

1 0.9 0.8 0.7 0.6 0.5

V/Vo

FIG.4. Normalized number of rings as a function of volume in

I

(a) shock simulations (M ) and (b) continuous compression

simulations (MII).

9

These simulations show that the ring size distribution does not change from 1

to 0.82 V/Vo, consistent with the linearity of the EOS curve in this same region.

Therefore, in this regime, no structural changes are occurring, and the compression

is elastic. Above 18% compression, large changes in the ring distributions in both

simulations are observed. These changes correlate well with the EOS curve between

0.82 and 0.64 V/Vo. There is a continuous decrease in the number of rings of size 5

and 6 during the plastic regime and an increase in the number of both smaller and

larger rings. Qualitatively both simulations show the same behavior. However,

dynamic simulations show a decrease in rings of size 6 at lower pressures than for

the case of continuous compression. This discrepancy could be due to size and

temperature differences between the two simulations. The temperatures reached in

method M I are higher than in MII. This effect would be particularly noticeable at the

highest pressures, but would not change qualitatively the trends in the EOS.

Although significant changes in the structure are observed during increasing

pressure, the average ring size for silica glass is still 6, since it preserves the

tetrahedral structure.

Figures 4 (a and b) also show that the smaller (3- and 4-member) rings

increase in total number upon compression in the region between 0.82 and 0.64 V/Vo.

These small rings are very close to planar rings and present vibration frequencies

that can be detected by Raman spectroscopy. Raman measurements of compressed

silica glass have shown enhancement of the D1 and D2 lines associated to rings of

size 4 and 3 respectively [15], in agreement with the results of our model. Moreover,

the number of rings of size 3 and 4 increases with pressure up to about 0.625 V/Vo at

which point it saturates for the case of rings of size 3, while the number of 4-member

rings seems to decrease. Recent experiments by Krol et al. have shown similar

behavior in the Raman spectra of laser irradiated silica glasses [16].

10

The bond angle distribution of both O-Si-O and Si-O-Si also change

continuously with pressure, even in the elastic regime. This implies that although

the distribution and number of rings is not changing in the elastic regime, their

geometry must be changing, as expected. The O-Si-O bond angle distribution

broadens with pressure while keeping a value centered around 109 to 108 for the

highest pressures. The Si-O-Si bond angle distribution shifts towards lower values as

pressure increases, in agreement with experimental observations [17]. At the highest

pressures this distribution presents two broad peaks, as shown earlier in simulations

by Vashishta et al. [20]. It is interesting to note that this change from a single peak to

a double peak seems to occur in the elastic to plastic regime, although it is hard to

identify clearly the point of transition.

The above analysis shows that the behavior of the EOS can be explained by

the changes observed in the ring size distribution and that these changes are not

clearly reflected in other characteristic parameters, such as the intensity of the FSDP

or its width. Both of them decrease continuously with pressure even during the

elastic regime, where we find no variations in the ring size distribution. Therefore,

the origin of this peak cannot be related only to the total number of rings present or

its distribution.

In summary, we have investigated the behavior of an amorphous system,

fused silica, under pressure using atomistic simulations. Our MD simulations of

both dynamic shock loading and continuous compression at constant temperature

reproduce the equation of state of flyer plate experiments. The transition between

elastic and plastic behavior is correlated to changes in the ring size distribution of

this amorphous system, and is not clearly reflected by other parameters such as pair

correlation function, structure factor or bond angle distributions. Therefore, these

11

changes must be governed by the preserved connectivity of these network structures

revealed by a decrease in 5- and 6-member rings at high pressure that compensates

with an increase of both smaller and larger rings.

Acknowledgements

This work was performed under the auspices of the U. S. Department of

Energy, Lawrence Livermore National Laboratory under Contract No. W-7405-Eng-

48. Funding was provided by the Department of Energy Office of Basic Energy

Sciences. One of the authors (MJC) wants to thank the Spanish MCYT for support

through the Ramon y Cajal program. Thanks to E. Bringa for useful discussions.

12

References

[1]P. W. Bridgman and I. Simon, J. Appl. Phys. 24, 405 (1953).

[2]R. Roy and H. M. Cohen, Nature (London) 190, 798 (1961).

[3]W. Primak, Compacted States of Vitreous Silica-Studies of Radiation Effects in Solids

(Gordon and Breach, New York, 1975).

[4]H. Sugiura, K. Kondo, and A. Sawaoka, J. Appl. Phys. 52, 3375 (1981).

[5]H. Sugiura, R. Ikeda, K. Kondo, and T. Yamadaya, J. Appl. Phys. 81, 1651 (1997).

[6]C. S. Zha, R. J. Hemley, H. K. Mao, T. S. Duffy, and C. Meade, Phys. Rev. B 50,

13105 (1994).

[7]A. Polian and M. Grimsditch, Phys. Rev. B 41, 6086 (1990).

[8]J. Wackerle, J. Appl. Phys. 33, 922 (1962).

[9]S. Susman, K. J. Volin, D. L. Price, M. Grimsditch, J. P. Rino, R. K. Kalia, P.

Vashishta, G. Gwanmesia, Y. Wang, and R. C. Liebermann, Phys. Rev. B 43, 1194

(1991).

[10]J. Wong, D. Haupt, J. H. Kinney, J. Ferreira, I. Hutcheon, S. Demos, and M. R.

Kozlowski, Proc. SPIE 4347, 466 (2001).

[11]R. Bruckner, J. Non-Cryst. Solids 5, 123 (1970).

[12]A. G. Revesz, J. Non-Cryst. Solids 7, 77 (1972).

[13]R. J. Hemley, C. T. Prewitt, and K. J. Kingma, Reviews in Mineralogy, Silica:

Physical Behavior, Geochemistry, and Materials Applications (Mineralogical Society of

America, Washington D.C., 1994), Vol. 29, Chap. 2.

[14]Q. Williams and R. Jeanloz, Science 239, 902 (1988).

[15]S. G. Demos, L. Sheehan, and M. R. Kozlowski, in Laser Applications in

Microelectronic and Optoelectronic Applications V, Proc. SPIE 3933, 316 (2000).

[16]J. W. Chan, T. Huser, S. Risbud, and D. M. Krol, Optics Letters 26, 1726 (2001).

[17]C. Meade, R. J. Hemley, and H. K. Mao, Phys. Rev. Lett. 69, 1387 (1992).

[18]R. Couty and G. Sabatier, J. Chim. Phys. 75, 843 (1978).

[19]L. V. Woodcock, C. A. Angel, and P. A. Cheeseman, J. Chem. Phys. 65, 1565

(1976).

[20]W. Jin, R. K. Kalia, P. Vashishta, and J. P. Rino, Phys. Rev. Lett. 71, 3146 (1993).

[21]D. J. Lacks, Phys. Rev. Lett. 84, 4629 (2000).

[22]J. P. Rino, G. Gutierrez, I. Ebbsjo, R. K. Kalia, and P. Vashishta, Mat. Res. Soc.

Symp. Proc. Vol. 408, 333 (1996).

[23]P. S. Salmon, Nature 1, 87 (2002).

[24]J. D. Martin, S. J. Goettler, N. Fosse, and L. Iton, Nature 419, 381 (2002).

[25]E. B. Webb II and S. H. Garofalini, J. Chem. Phys. 101, 10101 (1994).

[26]A. B. Belonoshko, Science, 275, 955 (1997).

[27]A. Kubota, M.-J. Caturla, J. S. Stolken, and M. D. Feit, Opt. Express 8, 611 (2001).

[28]Y. Inamura, M. Arai, M. Nakamura, T. Otomo, N. Kitamura, S. M. Bennington,

A. C. Hannon, and U. Buchenau, J. Non-Cryst. Solids 293, 389 (2001).

[29]S. R. Elliott, Nature 354, 445 (1991); S. R. Elliott, Phys. Rev. B 51, 8599 (1995).

13

Contact this candidate