Molecular dynamics study of the influence of mobile cations on the reconstruction of an irradiated silicate glass

Journal of Non-Crystalline Solids 330 (2003) 106–121 www.elsevier.com/locate/jnoncrysol

Molecular dynamics study of the influence of mobile cations on the reconstruction of an irradiated silicate glass J.-M. Delaye *, D. Ghaleb CEA, Centre de la Vall
ee du Rh^one, DEN/DIEC/SESC/LCLT, BP17171, 30207 Bagnols-sur-C
eze cedex, France Received 21 October 2002; received in revised form 26 May 2003

Abstract The effect of ballistic collisions in a simplified nuclear glass was investigated by molecular dynamics. Systematic results were obtained in the 0–16 keV energy range. Following a damage peak, reconstruction of the glass structure was observed in terms of the overall degree of polymerization. The reconstruction was facilitated by the presence of mobile cations. The dynamics of restoration of the SiO4 tetrahedrons during displacement cascades and during the formation of a structure from a random configuration can be fitted to curves corresponding to the same analytic formula. This similarity allowed us to examine the influence of mobile cations (alkalis or alkaline earths) on the formation of SiO4 entities in smaller systems formed from random configurations. The formation rate of SiO4 tetrahedrons accelerates with the Na2 O or CaO percentage to reach a maximum rate above about 10 mol% Na2 O or CaO. This threshold corresponds to the disappearance in the glass structure of zones comprising a central bridging oxygen surrounded by two first-neighbor bridging oxygen rings. No mixed alkali effect was observed in systems containing both Na and Cs because the formation of SiO4 does not require long-range diffusion of mobile cations.
2003 Published by Elsevier B.V. PACS: 61.43B; 61.80; 61.43F

1. Introduction As a containment matrix for radioactive elements, nuclear glass is subjected to the effects of b– c irradiation and a [1–3] decay. These irradiation phenomena result in electronic excitation and ballistic collisions. b–c irradiation produces elec-

*

Corresponding author. Tel.: +33-4 66 79 1794; fax: +33-4 66 79 6620. E-mail address: [email protected] (J.-M. Delaye). 0022-3093/$ - see front matter
2003 Published by Elsevier B.V. doi:10.1016/j.jnoncrysol.2003.08.062

tronic excitation, while a decay produces both electronic excitation by the emission of an a particle with an energy of about 4 MeV and ballistic collisions by the emission of a recoil nucleus with an energy between 70 and 100 keV. Numerous experiments have been carried out to observe the effect of different types of irradiation in glass [4–9]. Electronic excitation generates point defects that may assume a wide variety of forms. Several studies have demonstrated the limited impact of ballistic collisions on the glass structure. Many phenomena including swelling or

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

contraction, changes in macroscopic properties, or the appearance of bubbles, remain partially unexplained. The effects of irradiation vary significantly depending on the glass composition. In response to ballistic collisions arising from neutron irradiation, pure silica becomes increasingly dense [10]. Alkali silicate glasses react differently to electronic excitation. Expansion of the structure is followed by condensation [4]. The intensity of the expansion and of the succeeding condensation depends on the molar composition of the glass. The structural origins of these phenomena have not yet been satisfactorily explained. A survey of volume changes in aluminoborosilicate glass samples doped with radioactive elements shows that the intensity of the phenomenon is highly dependent on the glass composition [5,11]. In order to understand the effects of ballistic collisions alone and the influence of different chemical species on glass behavior, we began simulating displacement cascades in simplified nuclear glasses several years ago [12]. The cascade sequence is always the same. A highly disordered liquid zone is created during the first picosecond and progressively reconstructed during the next 5–10 ps. At low energies (less than 3 keV) the structure may be only partially restored leaving the glass depolymerized compared with its initial state, while at energies above 3 keV and up to 7 keV (the maximum energy simulated) the final material is only slightly depolymerized, if at all. The cascades were simulated in glass containing both network formers and modifiers. During the cascades, differences were observed between the displacement dynamics of the formers and modifiers. At the end of the cascade, a cloud of displaced modifiers surrounded the compact zone containing the displaced formers and oxygen atoms. Moreover, modifiers continue to be displaced longer than formers and oxygen atoms, suggesting that the modifiers have a major role in the dynamics of glass reconstruction after the thermal peak. We directly investigated the influence of the mobile cations on the formation rate of basic glass entities by studying the emergence of the liquid structure of an oxide in random atomic configu-

107

rations and measuring the rate at which the entities were formed. After showing that the method is dynamically representative of what actually occurs during the reconstruction of a glass following a displacement cascade, we were able to study the influence of various mobile cations on the restoration of the glass structure. The effect on irradiated glass reconstruction of mixing several alkalis is an interesting subject. A mixture of several alkalis is known to have a limiting effect on the diffusion coefficients of each alkali [13] (the well known mixed alkali effect). This phenomenon has not yet been fully elucidated, but could affect displacement cascades. Our objective was to ascertain the effect on the rates of formation of the basic entities of a glass containing two mobile cations. The morphology of displacement cascades simulated by molecular dynamics is broken down into several zones in which the disturbances range from mild to very severe. The heterogeneity of the structure of oxide glass (both at atomic scale due to the chemical diversity of the constituent species and at nanometric scale due to density or chemical fluctuations) thus affects its irradiation resistance. Increasingly numerous studies have demonstrated the structural heterogeneity of oxide glasses arising from the juxtaposition of zones of different density and local chemical composition [14–16]. These independent experiments converge on a value between 1 and 2 nm to describe the mean size of these heterogeneities. Observations of structures prepared by molecular dynamics also reveal heterogeneous distributions of alkalis and network formers. Ballistic collisions have already been widely studied by molecular dynamics in various types of materials [17]. The energy limit depends on the available computing power. The results of the 16 keV cascade shown here required nearly 1800 h of computer time using a Fujitsu VPP 5000/15 processor. This paper is separated in two parts. The first part is dedicated to study the structural glass behavior under auto irradiation. We show the dynamics of the different ionic species. The second part proposes an analysis of the role of the alkaline and earth alkaline elements when a glassy liquid is

108

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

formed from a random configuration. This second approach is useful to complete the information provided by the direct simulation of a displacement cascade. Indeed, we show that the formation of the basic structural entities are subjected to the same dynamical law in both cases. The general discussion (Section 4) proposes an improved insight on the dynamical phenomena during a displacement cascade in a glass, based on the direct simulations and the complementary studies of the molten glass formation from a random configuration.

2. Part I : Dynamical behavior of a glass submitted to a displacement cascade 2.1. Introduction

The O–U repulsive terms take the following analytic form: ! rij qi qj C þ
6: ð2Þ U2 ðrij Þ ¼ A exp
r qij rij The A parameter value is 1496.9 eV, q is equal to , and C is 65.4 eV A 6 . The qi terms corre0.4 A spond to the formal ion charges. The U-cation potentials are purely Coulombian with no additional repulsive term. Three-body terms are added to account for the directional order imposed by the presence of covalent bonds. The form of the threebody terms is based on work by Stillinger and Weber [18]. The following potential is associated with three atoms i, j and k forming a triplet centered on i:
 cij cik U3 ðrij ; rik ; hjik Þ ¼ k exp þ rij
rc rik
rc 2

 ðcos hjik
cos h0 Þ ; Part I presents the dynamical processes occurring in a simplified nuclear glass submitted to a 16 keV displacement cascade initiated by a heavy particle. This work is the continuation of previous studies dedicated to the analysis of a complex glass behavior under ballistic collisions ([12] and references herein). 2.2. Computational procedure for displacement cascades simulations The glasses considered in this article were simulated by molecular dynamics. A six-element glass was used to simulate the displacement cascades. The size of the cubic glass cell is 18.9 nm and contains 526 847 atoms. The molar composition is the following: 63.79% SiO2 + 17.02% B2 O3 + 13.39% Na2 O + 1.81% ZrO2 + 3.96% Al2 O3 + 0.03% UO2 . The small percentage of UO2 allows to have heavy ions to initiate the displacement cascades. The interaction potentials are Born– Mayer–Huggins pair terms. For two atoms i and j at a distance rij , the analytic formula of the pair term is the following: ! rij qi qj U2 ðrij Þ ¼ A exp
: ð1Þ þ qij rij

ð3Þ

where rij and rik are the pair distances i–j and i–k, and hjik is the angle formed by the triplet. The other terms are adjustable parameters. The adjustable parameters for the pair and three-body potentials are indicated in Tables 1 and 2 (we have included the parameters corresponding to elements used only in the second part i.e. Ca and Cs). The Cs-cation repulsions are purely Coulombian with no additional pair term. These potentials have already been used extensively in studies of glass structure and displacement cascades [19]. Glass structure and alkaline mobilities are correctly reproduced. Unfortunately, the use of formal ionic charges leads to a too rigid structure (bulk modulus [12] and thermal conductivities are too large by a factor probably close to 2. Due to the lack of experimental references, we can only compare with the complete nuclear glass properties). But the simulated structure remains a good model to understand the basic processes occurring during displacement cascades and the relative behaviors of formers and modifiers. To save computational time, the Ewald sum is reduced to its real space component with 
1 and a 20 A  cutoff radius. A recent a ¼ 0:1072 A study showed that this approximation had no appreciable effects on the short range structure [20].

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

109

Table 1  (except q  Adjustable parameter values Aij (eV) for Born–Mayer–Huggins potentials qi
j ¼ 0:29 A O–O ¼ 0:35 A) Species

Si

O

B

Na

Zr

Al

Ca

Cs

Si O B Na Zr Al Ca Cs

834.4

1571.2 352.68

337.7 760.9 121.1

862.0 1396.4 374.5 842.1

2557.4 4805.2 1031.8 2637.9 7822.7

961.3 1734.0 396.4 976.3 2940.7 1100.5

4000.9 687.3 1691.1 3990.8 1223.7 4326.8 1790.0

0 14 740.9 0 0 0 0 0 0

Table 2 ) and h0 for Adjustable parameter values ki (eV), cij ¼ cik (A three-body terms

Table 3 ) between the ZBL potentials, polynomial Coupling distances (A coupling potentials and BMH potentials

Parameter

O–Si–O, O–Al–O

Si–O-Si

O–B–O

Interaction

ZBL potential

Polynomial

BMH potential

ki (eV) ) cij ¼ cik (A h0

149.8 2.6 109.47

6.2 2.0 160.0

11 985.0 2.27 109.47

Other Si–B Si–U; B–B O–Na O–U

r < 0:8 r < 0:7 r < 0:8 r < 0:8 r < 1:1

0:8 < r < 1:1 0:7 < r < 1:1 0:8 < r < 1:5 0:8 < r < 1:4 1:1 < r < 1:3

r > 1:1 r > 1:1 r > 1:5 r > 1:4 r > 1:3

During simulations of displacement cascades, short-range interactions can appear temporarily following ballistic collisions between atoms. These short-range interactions are processed using Ziegler, Biersack and Littmarck (ZBL) potentials [21]. The ZBL potentials and pair potentials are related by fifth-order polynomial terms. Continuity of energy, force and force derivatives is imposed to determine the polynomial parameters. Table 3 indicates the coupling distances between the ZBL and BMH potentials. A cascade is initiated by accelerating an initial projectile: a uranium atom. At the beginning of the cascade the projectile is situated at position (0.1; 0; 0) in reduced coordinates in a cubic simulation cell for which all three dimensions range from 0.0 to 1.0. The direction of the acceleration vector is (1; 0; 0). The outer layer of the simulation cell is maintained at room temperature over a thickness of 3 . The mean temperature of the outer layer is A measured and readjusted to 293 K every 5 time steps by modifying the velocity of the atoms in the outer layer. The time steps are adjusted along a hysteresis loop in which the plateaus are determined by critical atomic velocities. The increasing and de-

Table 4 Velocity thresholds for modifying the time step during a displacement cascade simulation Decreasing velocity threshold (cm s1 ) 2 · 107 5 · 106 2 · 106 5 · 105 3.5 · 105 2 · 105

Time step (fs) 0.002 0.006 0.02 0.06 0.2 0.4 1

Increasing velocity threshold (cm s1 )

Time step (fs)

2.5 · 107 7.5 · 106 2.5 · 106 7.5 · 105 3.75 · 105 3.4 · 105

0.002 0.006 0.02 0.06 0.2 0.4

creasing plateau values are indicated in Table 4. When the largest atomic velocity drops below a decreasing plateau value the time step is increased by the corresponding amount in Table 4. Conversely, when the largest atomic velocity rises above an increasing plateau value the time step is diminished by the corresponding amount in Table 4. This technique ensures that the atom displace per time step during ment never exceeds 0.034 A the cascade simulation, while gradually increasing the time steps as the mean kinetic energy of the atoms decreases.

110

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

Table 5 Cutoff radius used to compute the atom coordination numbers ) Cutoff radius (A

Si–O

B–O

Al–O

U–O

2.1

2.2

2.4

3.4

During the ballistic collision sequence produced by the initial projectile, the degree of polymerization of the glass structure is monitored by computing a depolymerization index at regular intervals: X Npoly ðtÞ ¼ Nj ð0Þ
Nj ðtÞ: ð4Þ j¼Si;O;B;Al;U

For each of the atoms listed beneath the summation sign, we calculate the difference between the initial coordination number Nj ð0Þ and the coordination number at time t, Nj ðtÞ. The subscript Npoly ðtÞ is proportional to the number of formeroxygen bonds broken at time t in the glass structure. The result is positive when the glass network is depolymerized with respect to its initial state. The cutoff radius used to determine the coordination numbers to calculate Npoly are indicated in Table 5. 2.3. Results 2.3.1. Displacement cascades One U ion is accelerated with an energy of 16 keV according to the technique described previously. Fig. 1 shows positions of the atoms displaced by  between the initial instant and the more than 1 A end of the cascade. Inhomogeneities can be observed in the distribution of atom displacements. The displaced network formers and O atoms are found mainly along the projectile path, while the displaced Na atoms form a cloud scattered around the path. Fig. 1 does not distinguish the atoms, but the diffuse displacement cloud filling the entire cell consists essentially of Na atoms. The same phenomenon had already been distinctly observed during previous calculations [12] at lower energies, and were largely amplified as the projectile energy increased. The displaced Na cloud extends to the edges of the simulation cell. The thermal wave associated

Fig. 1. Morphology of a 16 keV cascade. Only atoms displaced  are shown. by more than 1 A

with the projectile and the longer-range Coulombian interactions are probably responsible for this phenomenon. The greater mobility of alkali atoms compared with other atomic species facilitates these displacements. The trail behind the projectile is not homogeneous, but shows alternating zones of higher and lower concentrations of displaced atoms. The coexistence of strongly and weakly perturbed zones has already been observed for lower energy cascades. Some highly modified zones and other zones containing practically no displaced atoms were also observed [12]. Sub-cascades also occur, but they are too limited in size to appear clearly within the overall morphology. The uranium atoms are much heavier than the other atoms in the glass, and the energy transmitted to the atoms initiating the subcascades is too low to create visible branches. Fig. 2 is equivalent to Fig. 1 for O atoms dis between the initial instant placed by more than 4 A and the end of the cascade. Highly damaged zones are delimited in the figure. Studies at lower energies demonstrated that these patterns arise from the structure of the glass itself, with alternating zones of high and low atomic densities. The coexistence of zones enriched by modifiers and formers is consistent with previous simulations

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

111

1000

800

NPOLY(t)

600

400

200

0 0 -200

2

4

6

8

10

12

time (ps)

Fig. 4. Variations of polymerization of the glass network during the displacement cascade. . The Fig. 2. Positions of O atoms displaced by more than 4 A square indicates the limits of the simulation cell.

performed by other authors [13] and with experiment [16]. The probability of collisions is greater in zones of high atomic density, resulting in zones of greater perturbation during displacement cascades. Conversely, high-energy projectiles can pass through zones of low atomic density with less damage. The size of the relatively undamaged zones (about 1 nm) is indicative of the dimensions of the zones forming the glass structure. This order of magnitude is consistent with recent glass structure studies that also reveal heterogeneous zones of comparable size [14–16]. We examined the variation over time of the number of displacements for each species (Fig. 3).

Number of displacements

3500 3000

O

2500 2000

Na

1500

Others

1000 500 0 0

2

4

6

8

10

12

t(ps)

 during the 16 keV Fig. 3. Atom displacements exceeding 1 A cascade.

We discriminated among displacements of alkalis, oxygen and other species, and observed different behavior. Displacements of network formers and oxygen atoms occur mainly at the beginning of a cascade; many alkali atoms are also displaced at the beginning of the cascade, but also with slower dynamics until the end of the cascade. These longer-term displacements are associated with the progressive restoration of the structure after a damage peak. A larger scale view of the displacement curves during the initial instants confirms what we have already observed during previous studies, i.e. the displacement of sodium atoms occurs with a delay of about 0.1 ps relative to the other species [22]. Fig. 4 shows the glass structure depolymerization index (Eq. (4)) versus time. Here again, the behavior is consistent with previously observed phenomena [12]. Strong depolymerization of the structure at the beginning of the cascade is followed by gradual repolymerization. This example shows incomplete structure repolymerization; the polymerization deficit after the cascade corresponds to about 50 atomic bonds. Depolymerization followed by repolymerization is due to the rupture and subsequent reconstruction of entities around network formers, i.e. Si, B and Al atoms. Fig. 5 shows the Si percentage in SiO4 entities. The minimum and the rise from 0 to 2 ps

112

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

placements versus time is probably not an accident, and suggests a correlation between Na displacements and the reconstruction of local entities. The influence of the Na concentration on the formation of local entities is addressed in the second part of the paper.

99.6 99.55 99.5

%SiO4

99.45 99.4

99.35 99.3 99.25 99.2 99.15 0

5 time(ps)

10

Fig. 5. Percentage of 4-coordinated Si versus time during a 16 keV cascade.

correspond to the dissociation and reformation of SiO4 tetrahedrons. Conversions are also observed between BO3 and BO4 entities and between AlO3 and AlO4 . The variations over time of 4-coordinated B and Al are shown in Fig. 6. The mean percentage of 4-coordinated Al returned to its initial level, while the drop in the percentage of 4-coordinated B (replaced by 3-coordinated B) explains the overall depolymerization. The measured internal pressure variation corresponds to an estimated structure swelling of 0.17%. The mean decrease in the boron coordination number has recently been correlated with the overall swelling [23]. The resemblance between the shapes of the curves for tetra-coordinated Si and Al and the shape of the curve for the number of Na dis-

2.3.2. Effect of increasing projectile energy This study of a 16 keV cascade allowed us to extend the energy range investigated to date. Many observations in the (0–7 keV) range remain valid at 16 keV. The overall morphology of the damage zone alternates between highly perturbed and less perturbed regions. After a depolymerization peak, the amplitude of which increases with the energy, the glass structure is restored almost to its initial state. Recent calculations suggest that the depolymerization peak increases with the energy up to 16 keV and begins to saturate beyond this point. The initial and final ring size distributions are identical, confirming the medium-scale reconstruction of the initial order. Sodium atom displacements differ by their extended duration following a rapid initial growth phase. Increasing the energy mainly affects the positions of the displaced Na atoms. The displacement cloud is substantially more dense than observed in the (0–7 keV) range, and extends to the edges of the simulation cell. We attribute the presence of this cloud to the combined effects of the greater Na atom mobility, of the thermal wave created by the projectile, and of long-range Coulombian interactions.

58.3

84.2

58.25

84

%AlO4

84.4

%BO4

58.35

58.2

83.8

58.15

83.6

58.1

83.4 0

5 time(ps)

10

0

5 time(ps)

Fig. 6. Percentage of 4-coordinated B (left) and Al (right) during a 16 keV cascade.

10

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

3. Part II : Analysis of the formation of basic glass entities from a disordered configuration 3.1. Introduction This part is dedicated to the dynamical processes occurring during a molten glass formation from a completely random configuration. We will show in the discussion (Section 4) that the dynamical laws governing the molten glass formation form a random configuration and the reconstruction of a glass after a collision cascade obey the same formula. This study is therefore very useful for a better understanding (by a faster method than the complete simulation of a displacement cascade) of the role of the ionic species, especially the mobile ions, during the glass reconstruction after a local perturbation. 3.2. Computational procedure We performed the following calculations. Beginning with a random configuration in which the only constraint was that all the atoms were initially  apart, we examined the formation rates of 1.2 A SiO4 entities in the following systems: SiO2 –Na2 O, SiO2 –CaO, SiO2 –B2 O3 , SiO2 –Na2 O–B2 O3 , SiO2 – Na2 O–CaO, SiO2 –Na2 O–Cs2 O. The molar compositions are indicated in Table 6. The simulation cells contained 5184 atoms. The interatomic potentials are already exposed in part I (see Section 2.2). The adjustable parameters are given in Tables 1 and 2. The Coulombian interactions are computed using the Ewald formalism [24]. The a parameter and the cutoff radius in real and reciprocal space were selected to ensure a precision of 10
3 on the energy [25].

In each case, the temperature was maintained between 2000 and 3000 K, i.e. whenever the temperature exceeded 3000 K it was reduced to 2000 K by adjusting the atomic velocities. Many corrections were required at the beginning of each simulation, but as the structure neared equilibrium the number of corrections dropped to zero. The simulation cell volume was adjusted each time so the final equilibrium pressure was near zero (between )3 · 108 and 3 · 108 Pa). Several tests were necessary in some cases. Each simulation lasted ten thousand time steps of 10
3 ps each. The cutoff radius used to determine the coordination numbers of network-forming atoms are indicated in Table 5. 3.3. Results 3.3.1. Molten silicate glass formation dynamics The SiO4 percentage is plotted versus the time steps in Fig. 7. Only the final variations between 90% and 100% are shown; the rapid growth phase from 0% to 90% was always identical, and of no particular interest. Each of the curves plotting the SiO4 percentage versus time shows three characteristic elements: rapid growth for a short time, a curvature, and a final plateau around the SiO4 equilibrium percentage. The slope of the rapid growth section gives the rate of formation of SiO4 tetrahedrons and the difference at 100% of the final equilibrium value gives the number of residual defects in the liquid. In compositions containing a single mobile cation the formation rate of SiO4 entities increases as the mobile cation concentration rises (the slope of the increasing SiO4 percentage before the

Table 6 Molar compositions of systems formed from a random configuration SiO2 95SiO2 Æ 5Na2 O 90SiO2 Æ 10Na2 O 85SiO2 Æ 15Na2 O 80SiO2 Æ 20Na2 O 75SiO2 Æ 25Na2 O 85SiO2 Æ 5Na2 O Æ 10B2 O3 85SiO2 Æ 7.5Na2 O Æ 7.5B2 O3 85SiO2 Æ 10Na2 O Æ 5B2 O3

113

95SiO2 Æ 5CaO 90SiO2 Æ 10CaO 85SiO2 Æ 15CaO 80SiO2 Æ 20CaO 75SiO2 Æ 25CaO 85SiO2 Æ 5Na2 O Æ 10CaO 85SiO2 Æ 7.5Na2 O Æ 7.5CaO 85SiO2 Æ 10Na2 O Æ 5CaO

95SiO2 Æ 5B2 O3 90SiO2 Æ 10B2 O3 85SiO2 Æ 15B2 O3 80SiO2 Æ 20B2 O3 75SiO2 Æ 25B2 O3 85SiO2 Æ 5Na2 O Æ 10Cs2 O 85SiO2 Æ 7.5Na2 O Æ 7.5Cs2 O 85SiO2 Æ 10Na2 O Æ 5Cs2 O

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121 100 99 98 97 96 95 94 93 92 91 90

100

25%

25%

10% 5% 0%

10%

98 %SiO4

%SiO4

114

0%

96 94 92 90

0

5000 step

0

10000

5000 step

10000

Fig. 7. Formation rate of SiO4 entities in SiO2 + Na2 O glass (left) and SiO2 + CaO glass (right). The mobile cation percentage is indicated on the curves.

100

2 3

97 95

91

%SiO4

94 92

1 2

98 97

3 4

96 95 94 93 92 91 90 0

2000

4000

6000 step

8000

10000

12000

Fig. 9. Formation rate of SiO4 entities in SiO2 + Na2 O + Cs2 O glass. Subscripts 1 to 4 correspond to concentrations with 15% Na2 O, 10% Na2 O + 5% Cs2 O, 5% Na2 O + 10% Cs2 O and 15% Cs2 O, respectively.

pure SiO2

96

93

99

entities. The slope of the rising SiO4 percentage over time diminishes when Cs2 O replaces Na2 O. No non-linear effect of mixed alkalis is observed on the formation of SiO4 tetrahedrons in the presence of these two oxides.

1

99 98

100

%SiO4

curvature is higher). Saturation levels between 10% and 15% Na2 O or CaO are observed. Moreover, the final number of defects diminishes as the mobile cation percentage increases. Fig. 8 shows the time necessary to form SiO4 entities when two mobile cations (Na and Ca) are present simultaneously. When Na2 O is progressively replaced by CaO, the number of defects increases but the SiO4 formation rates are comparable. In each case, SiO4 entities are formed more readily in the presence of modifier elements. In the presence of CaO, the structure is restored at the same rate, but is of lower quality than in the presence of Na2 O. The continuous replacement of Na by Cs (Fig. 9) diminishes the efficiency of formation of SiO4

90 0

2000

4000

6000 step

8000

10000

12000

Fig. 8. Formation rate of SiO4 entities in SiO2 + Na2 O + CaO glass. Subscripts 1, 2 and 3 correspond to 85% SiO2 + 15% Na2 O, 85% SiO2 + 5% Na2 O + 10% CaO and 85% SiO2 + 15% CaO glasses, respectively.

3.3.2. Specific formation times of SiO4 Modifiers such as Na or Ca atoms or a mixing of both elements clearly enhance the rate of formation of SiO4 tetrahedrons. We have tried to extract from the Figs. 7–9 a specific formation time of the SiO4 tetrahedrons by the following method. The different curves are fitted by the formula (5):

rffiffiffiffi t ; ð5Þ A 1
exp
t0 where A and t0 are constants, and t is the time.

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

99 98

%SiO4

97 96 95 94 93 92 0

2000

4000

6000 step

8000

10000

12000

Fig. 11. Fitting Eq. (5) to the curves for the rate of formation of SiO4 entities in 90SiO2 Æ 10Na2 O glass (top) and pure SiO2 (bottom).

120 115 110 105 100 95 90 85 80 75 70 0

5

10

15 20 %Na2O ou%CaO

25

30

Fig. 12. Specific formation time of SiO4 tetrahedrons versus Na2 O percentage in SiO2 + Na2 O glass (solid symbols) and versus the CaO percentage in SiO2 + CaO glass (open symbols).

ishes to a saturation level above 10% Na2 O or CaO. The shorter the time, the faster the formation of SiO4 entities. The presence of modifiers clearly favors the SiO4 tetrahedron formation rate. The glass is less viscous and, as shown in Fig. 13, the mean square displacements for O increase with the Na2 O concentration. The presence of modifiers enhances oxygen mobility and thus facilitates the formation of the local entities constituting the glass network.

16 14 12 RDF Si-O

100

t0(fs)

The analytic expression of this formula arises from the following reasoning. Let D be a coefficient expressing the mobility of the O atoms, N the number of 4-coordinated Si, and N0 the number of 4-coordinated Si once the final state is reached. Within a time interval dt, the number dN of Si atoms acquiring a coordination number of 4 is assumed proportional to the number of O atoms entering the Si coordination sphere, which in turn is proportional to the distance dr covered on average by the O atoms the time interval dt qffiffiduring ffi pffiffiffiffiffi (i.e. dr ¼ d Dt ¼ 12 Dt dtÞ. dN is also assumed proportional to the number of Si atoms that are not yet tetra-coordinated, i.e. N0
N . Hence N0dN
N is proportional to pdtffit; integrating this expression yields formula (5). We may refer to a coordination sphere for Si even at high temperatures because the radial Si–O distribution functions clearly show a first peak corresponding to the first-neighbor shell followed by a minimum equal to 0 (Fig. 10). t0 in formula (5) is considered to be the specific formation time of the SiO4 entities. Formula (5) was fitted to the curves between 94% and 100% to eliminate the first points corresponding to a state too disordered to have any physical reality. Fig. 11 shows two examples of fitting Eq. (5) to the simulation results. Fig. 12 shows the variation of the specific time t0 versus the Na2 O or CaO composition in the SiO2 + Na2 O and SiO2 + CaO systems. It dimin-

115

10 8 6 4 2 0 0

2

4 r(Ang)

Fig. 10. Radial distribution 0.85SiO2 Æ 0.15Na2 O system.

function

6

for

8

Si–O

in

a

3.3.3. Borosilicate glass When two network formers were present in the structure (SiO2 + B2 O3 ) we compared the SiO4

116

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

1.1E-15 4 3 2 1

r2(cm2)

9E-16 7E-16 5E-16 3E-16 1E-16 0

2000

4000

6000 step

8000

10000

12000

Fig. 13. Mean square displacement of O atoms in ð1
xÞSiO2 + xNa2 O glass. Curves 1, 2, 3 and 4 correspond to x ¼ 0:05, 0.15, 0.2, 0.25, respectively.

100 99 98 97 96 95 94 93 92 91 90

pure SiO2 1

0

3

5000 step

%SiO4

%SiO4

formation rates (Fig. 14). Boron alone has a limiting effect on the formation of tetrahedrons around Si atoms. Non-linear behavior is observed when the B2 O3 concentration is increased. The SiO4 formation rate drops to a minimum at about 15% B2 O3 and then increases again. A mixture of B2 O3 + Na2 O favors the formation of Si tetrahedrons (Fig. 14); the SiO4 formation rate diminishes when the Na2 O concentration decreases. SiO2 –B2 O3 glasses are the only compositions among those simulated here exhibiting non-linear behavior when the B2 O3 concentration increases. Fitting the SiO4 formation rate in SiO2 + B2 O3 glass using Eq. (5) is much less accurate than for SiO2 glass containing modifiers. The simple model

expressed in Eq. (5) cannot account for the shape of the curves. In the absence of mobile cations, more complex movements such as the correlated construction of entities around Si and B could become predominant. This would require a model describing the joint construction of SiO4 , BO3 and BO4 entities. Despite the absence of mobile cations in the pure borosilicate SiO2 –B2 O3 , the B atoms are at coordination number 3 or 4; charge compensation is provided by 3-coordinated O atoms. O atoms are found with three possible coordination numbers: 1, 2 or 3; coordination number 2 largely predominates in each of these systems. Nevertheless, the percentage of 1-coordinated O atoms sharply declines when the B2 O3 concentration increases. It drops from 1.3% to 0.1% between pure SiO2 and a composition containing 25% B2 O3 . At the same time, the concentration of 3-coordinated O atoms rises from 0.23% to 2.3%. Coordination numbers 1 and 3 are O atom defects that increase the degrees of freedom for the formation of the glass structure. The sum of the 1and 3-coordinated oxygen concentrations is plotted versus the B2 O3 percentage in Fig. 15. It is interesting to note that the concentration of these defects drops to a minimum for around 10–15% B2 O3 concentration, which corresponds to the longest SiO4 tetrahedron formation time. In other words, the greater the number of defects, the easier the formation of SiO4 . The non-linearity of the defect concentration accounts for the non-linear behavior of SiO2 –B2 O3 systems.

2

10000

100 99 98 97 96 95 94 93 92 91 90

1 2 pure SiO2

0

5000 step

3

10000

Fig. 14. Formation rate of SiO4 entities in SiO2 + B2 O3 glass (left) and SiO2 + B2 O3 + Na2 O glass (right). Subscripts 1, 2 and 3 on the left-hand curves correspond to glass containing 5%, 15% and 25% B2 O3 , respectively. Subscripts 1 to 3 on the right-hand curves correspond to concentrations with 15% Na2 O, 10% Na2 O + 5% B2 O3 and 15% B2 O3 , respectively.

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

importantly, the rate of disappearance of SiO3 entities in favor of SiO4 entities followed curves of the same shape, which we were able to fit using the following relation:
rffiffiffiffi t þ B: ð6Þ A exp
t0

3 2.5 %O1+%O3

117

2 1.5 1

The origin of this equation is discussed in the preceding section. A, B and t0 are constants, t is the time. The same formula can be used to obtain a correct fit of both types of curves, but with quite different values for the A parameter. The compositions, the initial states and the annealing thermal conditions are different, hence the differences in the fitting parameters. Nevertheless, the fact that the curves can be fitted using the same formula suggests that the dynamics of formation of SiO4 groups from a disordered glass rich in SiO3 groups involves the same processes after the passage of a high-energy projectile as during the emergence of a structure from a random configuration. We verified that the model also satisfactorily describes displacement cascades in simpler glasses. Fig. 17 shows the variation in the SiO3 percentages and the results fitted using Eq. (6) for 1 and 1.2 keV cascades in SiO2 glass and in 90% SiO2 + 10% Na2 O glass. As in Fig. 16, only the decreasing portion of the curve corresponding to the restoration of SiO4 tetrahedrons and the disappearance of SiO3 entities was adjusted. The model in Eq. (6) is thus applicable to cascades in simpler glasses. The reconstruction of glass after a collision cascade and the formation of glass from a random

0.5 0 0

10

20

30

%B2O3

Fig. 15. Percentage of 1- or 3-coordinated O in SiO2 Æ B2 O3 versus the B2 O3 molar concentration.

4. General discussion 4.1. Comparison of the dynamics of SiO4 formation under irradiation and in a random configuration under construction

100

1

80

0.8 %SiO3

%SiO3

Fig. 16 shows the typical variation over time of the SiO3 concentration during glass formation (in this case the glass composition was 90% SiO2 + 10% Na2 O). The SiO3 concentration reaches a peak before diminishing with respect to SiO4 entities. The formation of SiO3 entities is an intermediate step leading to the appearance of SiO4 entities in structures formed from a random configuration. Fig. 16 also shows the rates of formation and disappearance of SiO3 entities during a 16 keV displacement cascade. In both cases the concentration peaked (for different reasons) but, more

60 40 20

0.6 0.4 0.2

0

0 0

0.5

1 time(ps)

1.5

2

0

0.5

1 time(ps)

1.5

2

Fig. 16. Percentage of SiO3 entities versus time in 90% SiO2 + 10% Na2 O glass formed from a random configuration (left) and in glass subjected to a 16 keV cascade (right). The grey curves represent the fit from Eq. (6).

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121 2.4 2.3

1.2

2.2 2.1

0.8

1 %SiO3

%SiO3

118

2 1.9 1.8

0.6 0.4 0.2 0

1.7 0

0.5

1 time (ps)

1.5

2

0

0.5

1 time (ps)

1.5

2

Fig. 17. Percentage of SiO3 entities versus time in SiO2 glass (left) and in 90% SiO2 + 10% Na2 O glass (right) subjected to 1 and 1.2 keV cascades, respectively.

configuration follow the same dynamical trend. It was therefore founded to rest on the second method analysis to understand the dynamical phenomena during a displacement cascade.

second-neighbor shells. And we have calculated the concentration of O atoms which have only bridging oxygens in their first and second-neighbor shells. These oxygens are the center of strongly connected entities. Fig. 18 shows the percentage of oxygen atoms having only bridging oxygen in their first- and second-neighbor shells in SiO2 + Na2 O and SiO2 + CaO systems. This percentage diminishes with the modifier concentration, reaching zero between 10% and 15% Na2 O or CaO, a value near the threshold mentioned above. In pure SiO2 , the percentage of oxygen atoms having only bridging oxygens in the first two neighbor shells is not 100% but about 60%; the small residual percentage of non-bridging oxygens (1.3%) is sufficient to account for this difference. The formation rate of SiO4 entities reaches its maximum when there are no more strongly connected local groups approximately 1 nm in diameter (a central O and two layers of bridging oxygen neighbors) that do not contain a non-bridging oxygen atom. Increasing the concentration of

4.2. SiO4 formation rate threshold We looked for possible correlations between the glass structure in the liquid state and the threshold observed for the SiO4 formation rate at about 10% Na2 O or CaO. One of the consequences of increasing the Na2 O concentration is the formation of non-bridging oxygens by the rupture of Si–O–Si bonds. We analyzed the changes in the environment of the oxygen atoms according to the Na2 O concentration. The O–O radial distribution function (not  shown here) reveals two minima: one at 3.4 A delimiting the first-neighbor shell and a second at  delimiting the second-neighbor shell. 5.85 A So, for each O atom we are able to calculate the number of bridging oxygen inside the first and

70 50 40

% of O

% of O

70

SiO2+Na2O

60

30

40 30

20

20

10

10

0

SiO2+CaO

60 50

0 0

10

20 %Na2O

30

0

10

20

30

%CaO

Fig. 18. Percentage of O with only bridging oxygens in the first two neighbor shells in SiO2 + Na2 O glass (left) and SiO2 + CaO glass (right).

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

4.3. Mixed alkali effect The mixed alkali effect is well known with regard to diffusion in oxide glasses, but was not observed during the structure formation in the systems studied here. The specific time related to the formation rate of SiO4 entities increased continuously as 15% Na2 O was replaced by 15% Cs2 O. In addition, the mean square displacements of O atoms (not shown here) increase continuously during the transition from a system of 85% SiO2 + 15% Na2 O to a system of 85% SiO2 + 15% Cs2 O. Two surprising phenomena are worth mentioning in this regard. The specific time related to the formation rate of SiO4 entities remained relatively constant for all the systems of the type ð100
xÞSiO2 + xCs2 O (where x varied from 0% to 25%). Unlike silicate systems with Na2 O or CaO, the rising molar percentage of Cs2 O did not shorten the formation time of the SiO4 entities.

1.6E-15 1.4E-15

4

1.2E-15 r2(cm2)

mobile cations (Na2 O or CaO) gradually eliminates these highly connected pockets, enhancing the oxygen mobility and thus the glass formation rate. The disappearance of the polymerized zones allows the glass to form at the maximum rate. It is interesting to note that the thresholds in SiO2 + Na2 O and SiO2 + CaO systems are very close. It takes only half as many Ca atom to eliminate the polymerized zones: the +2 charge of the Ca atoms breaks two Si–O–Si bonds while the +1 charge of the Na atoms can break only one. Despite disappearance of the polymerized zones the glass structure remains heterogeneous, as shown by the morphology of the displacement cascades initiated in glass containing more than 13% Na2 O (Fig. 1). Although the Na2 O concentration corresponds to the disappearance threshold of the polymerized zones, the cascade morphology remains heterogeneous with the coexistence of highly and slightly perturbed zones. Structural disparities persist within the glass. A visual examination suggests that these disparities can be attributed to zones of variable atomic density and to nanometric fluctuations in the Na concentration. It is reassuring to note that restoration after the passage of a projectile occurs at near the maximum rate in nuclear glasses.

119

1E-15

3 2

8E-16

1

6E-16 4E-16 2E-16 0 0

2000

4000

6000 8000 Step

10000 12000

Fig. 19. Mean square displacement of O atoms in ð1
xÞSiO2 + xCs2 O glass. Curves 1, 2, 3 and 4 correspond to x ¼ 0:05, 0.15, 0.2, 0.25, respectively.

The second point concerns the mean square displacement of O atoms in SiO2 –Cs2 O systems (Fig. 19). The total mean square displacement of O atoms clearly increases with the Cs2 O percentage. Comparing Figs. 13 and 19 shows that, for the same molar percentage of Na2 O or Cs2 O, the mean square displacement of O atoms glass is higher in SiO2 –Cs2 O glass than in SiO2 –Na2 O glass, while the time of formation of SiO4 entities is longer for SiO2 –Cs2 O glass. These results suggest that the formation of SiO4 entities does not require long-range diffusion of oxygen atoms, but rather more local reorganizations. Analysis shows that the atom displacements during simulated displacement cascades are consistent with this interpretation. We did not observe any long-range migration of alkali atoms, but rather coupled movements among network formers, oxygens and mobile atoms. The peripheral cloud of displaced Na atoms around the highly perturbed zone corresponds to collective displacements of chains of formers and oxygen atoms accompanied by Na atom movements (Fig. 1). Long-range diffusion of mobile cations is unnecessary for the restoration of the glass structure; this accounts for the lack of a mixed alkali effect. The fact that the formation time of SiO4 entities in SiO2 –Cs2 O is longer than in SiO2 –Na2 O glass while the O atom displacements also occur over longer distances can be attributed to a structure effect. The densities of systems containing Cs are

120

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121

lower than those of systems containing Na, and on average the distances between SiO4 tetrahedrons are longer. This simple density effect explains why the construction of SiO4 tetrahedrons from a random configuration takes longer in a glass containing Cs than in a glass containing Na. The fact that the specific formation time of the SiO4 tetrahedrons remained relatively constant irrespective of the Cs2 O molar concentration shows that there are two competing mechanisms. Increasing the mobile cation concentration not only enhances the mobility of the O atoms but also reduces the density and thus lengthens the path covered by the O atoms to form the structure. The first mechanism predominates in glass containing Na or Ca atoms; the two mechanisms cancel each other in glass containing Cs atoms, however, and the SiO4 formation time remains constant.

5. Conclusion The effect of ballistic collisions in a simplified nuclear glass was investigated by molecular dynamics. Systematic results were obtained in the 0– 16 keV energy range. Following a damage peak, reconstruction of the glass structure was observed in terms of the overall degree of polymerization. The reconstruction was facilitated by the presence of mobile cations as shown by additional calculations on simple compositions. The dynamics of restoration of the SiO4 tetrahedrons during displacement cascades and during the formation of a structure from a random configuration can be fitted to curves corresponding to the same analytic formula. This similarity allowed us to examine the influence of mobile cations (alkalis or alkaline earths) on the formation of SiO4 entities in smaller systems formed from random configurations. The presence of alkalis and alkaline earths accelerates the formation pffiffiffiffiffiffiffi of SiO4 tetrahedrons. By fitting an exp
t=t0 curve we defined a specific formation time for SiO4 entities in different glass compositions. The formation rate accelerates with the Na2 O or CaO percentage to reach a maximum rate above about 10 mol% Na2 O or CaO. This

threshold corresponds to the disappearance in the glass structure of zones comprising a central bridging oxygen surrounded by two first-neighbor bridging oxygen shells. At the same time, the oxygen mobility increases in glass containing alkali or alkaline earth elements. No mixed alkali effect was observed in systems containing both Na and Cs because the formation of SiO4 does not require long-range diffusion of mobile cations.

References [1] W.J. Weber, R.C. Ewing, C.A. Angell, G.W. Arnold, A.N. Cormack, J.-M. Delaye, D.L. Griscom, L.W. Hobbs, A. Navrotsky, D.L. Price, A.M. Stoneham, M.C. Weinberg, J. Mater. Res. 12 (1997) 1946. [2] H.J. Matzke, E. Vernaz, J. Nucl. Mater. 201 (1993) 295. [3] Y. Inagaki, H. Furuya, Y. Ono, K. Idemitsu, T. Banba, S. Matsumoto, S. Muruoka, in: Mater. Res. Symp. Proc., XVI Symp. Scientific Basis of Waste Management, Boston, MA, 1992, p. 199. [4] G.W. Arnold, Mater. Res. Soc. Proc. 44 (1985) 617; G.W. Arnold, Radiat. Eff. 98 (1986) 55. [5] R.C. Ewing, W.J. Weber, F.W. Clinard Jr., Prog. Nucl. Energy 29-2 (1995) 63. [6] J.F. DeNatale, D.G. Howitt, Nucl. Instr. Meth. Phys. Res. B 1 (1984) 489; J.F. DeNatale, D.G. Howitt, G.W. Arnold, Radiat. Eff. 98 (1986) 63. [7] B. Boizot, G. Petite, D. Ghaleb, G. Calas, J. Non-Cryst. Solids 283 (2001) 179. [8] G. Battaglin, G.W. Arnold, G. Mattei, P. Mazzoldi, J.-C. Dran, J. Appl. Phys. 85 (1999) 8040; G. Battaglin, G. DellaMea, Nucl. Instr. Meth. Phys. Res. B 19/20 (1987) 948. [9] A. Abbas, Y. Serruys, D. Ghaleb, J.-M. Delaye, B. Boizot, B. Reynard, G. Calas, Nucl. Instr. Meth. Phys. Res. B 166– 167 (2000) 445. [10] E. Lell, J. Kreidl, J.R. Hensler, in: J.E. Burke (Ed.), Progress in Ceramic Science, 4, Pergamon, 1966; A. Wootton, B. Thomas, P. Harrowell, J. Chem. Phys. 115 (2001) 3336. [11] J.-C. Dran, Solid State Phenom. 30&31 (1993) 367. [12] J.-M. Delaye, D. Ghaleb, Phys. Rev. B 61 (2000) 14481. [13] A.N. Cormack, in: K. Wright, R. Catlow (Eds.), Microscopic Properties and Processes in Minerals, Kluwer, Dordrecht, 1999, p. 577. [14] P.H. Gaskell, D.J. Wallis, Phys. Rev. Lett. 76 (1996) 66. [15] L. Cormier, P.H. Gaskell, G. Calas, A.K. Soper, Phys. Rev. B 58 (1998) 11322. [16] C. Chemarin, B. Champagnon, J. Non-Cryst. Solids 243 (1999) 281.

J.-M. Delaye, D. Ghaleb / Journal of Non-Crystalline Solids 330 (2003) 106–121 [17] R.S. Averback, T. Diaz de laRubia, Solid State Phys. 51 (1998) 281. [18] F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262. [19] J.-M. Delaye, V. Louis-Achille, D. Ghaleb, J. Non-Cryst. Solids 210 (1997) 232; J.-M. Delaye, D. Ghaleb, J. Nucl. Mater. 244 (1997) 22; F. Angeli, J.-M. Delaye, T. Charpentier, J.-C. Petit, D. Ghaleb, P. Faucon, J. Non-Cryst. Solids 276 (2000) 132. [20] A. Abbas, J.-M. Delaye, D. Ghaleb, G. Calas, J. NonCryst. Solids 315 (2003) 187.

121

[21] J.F. Ziegler, J.P. Biersack, U. Littmark, The Stopping Range of Ions in Matter, Pergamon, New York, 1985. [22] J.-M. Delaye, D. Ghaleb, Radiat. Eff. Def. Solids 142 (1997) 471. [23] J.-M. Delaye, D. Ghaleb, Nucl. Instr. Meth. Phys. Res. B 191 (2002) 10. [24] P.P. Ewald, Ann. Phys. Z. 19 (1921) 253; S. Deleeuw, J.W. Perram, E.R. Smith, Proc. R. Soc. London A 373 (1980) 27. [25] D. Fincham, Molec. Simul. 13 (1994) 1.