Properties of irradiated sodium borosilicate glasses from experiment and atomistic simulations

With a combination of atomistic modeling and experimental techniques, we have investigated the structural and elastic parameters of sodium borosilicate glasses, including irradiation- induced changes. Both approaches show that the Young's modulus depends linearly on the density of material. The simulated glass density and boron speciation match also the estimates by independent, elemental glass composition-based models, indicating that atomistic simulations could be used in validation of theoretical models and experimental results. This allows us to formulate Young's modulus— density relationships for investigated borosilicate glasses and test the existing empirical model for description of Vickers hardness of these materials. The simulation of irradiation reveals a change of B [4] content under irradiation. By applying a simple defects accumulation procedure, we are able to correctly reproduce the measured critical irradiation dose of ~0.1 dpa and provide reasonable information on density change and stored internal energy. With the obtained agreements between the experimental and simulation results, we obtained superior insights into the atomic-scale structural evolution of irradiated borosilicate glasses.

B 2 O 3 by mass. 6 After detailed systematic studies of the effects of different additions and substitutions to a very common soda lime glass, in the late 19th century, Otto Schott and Ernst Abbe found that a borosilicate glass with significant amount of boron has superior chemical durability and improved thermal shock resistance. For instance, it has three times lower thermal expansion coefficient than the soda lime glass. [6][7][8][9] Because of the enhanced performance, borosilicate glasses are technologically applied at extreme conditions, with wide range of applications. These include car headlights, cookwares, such as water bottles, coffee pots or microwave glasses, laboratory apparatus and glassware, and special windows in washing machines or at the international space station. 5,10,11 They are also used as a host of High-Level Nuclear Waste (HLW). 1,3 In the nuclear applications, the borosilicate glasses are widely used as a primary matrix/host for the immobilization/ vitrification of HLW because of good chemical durability and physical parameters, including good thermal stability, resistance to neutral or acid solutions, and the ability to incorporate about 30 different elements or compounds found in spent nuclear fuel. 12,13 They have large waste load (40% HLW by volume) and could be easily formed using low-temperature routes. 1 They are also inexpensive to produce, as boron is abundant worldwide and easily accessible. Besides, the manufacturing technologies are well established. 14,15 Vitrification has been accepted as one of the most suitable methods to solidify HLW. 1,3,12,14 However, due to safety concerns, it is necessary to estimate the parameters of borosilicate glass under residual irradiation. In the past years, many groups have extensively studied the thermodynamic, chemical and mechanics properties of borosilicate glasses, and significant experimental effort has been devoted to evaluation of the properties of irradiated borosilicate glasses. [16][17][18][19][20][21][22] Yet, these experimental results require detailed interpretation, which could be achieved by correct understanding of the atomicscale processes associated with radiation damage and defects accumulation. It could be achieved with the aid of atomicscale simulations using supercomputing resources and simulation methods of computational chemistry and materials science. 23,24 In the last decades, different glasses have been investigated at the atomistic scale by atomistic modeling methods, including molecular dynamics. There are various such studies of nuclear glasses that aim to understand the structural and physical-chemical properties of these materials. 1,4,5,7,[25][26][27][28][29][30][31][32] Atomistic modeling methods have been applied to understand the local structure of the glassy materials, 25,30 speciation in glasses or quenched melts, 26 or elastic parameters of such materials. 28,29 At the same time, various techniques were developed to simulate glassy materials and other material compounds under irradiation. 24,26,33,34 In this contribution we used a combination of experimental measurements and atomistic scale modeling to understand the behavior of irradiated sodium borosilicate glasses and to derive the associated physical parameters. We focus on derivation of the structural, elastic, and thermodynamic properties, such as change in volume, Young's modulus, internal (stored) energy and other parameters that are determined by boron speciation, as well as on understanding the atomic structure of the irradiated glasses and finding structureproperty relationships. The reported results also represent a good validation test for the performance of atomistic modeling methods on the simulation of properties of irradiated, structurally complex and disordered materials.

| Experimental methods
We used the same sample preparation method and multi-energy ion irradiation experiment procedure as Guan et al. 19 The irradiation with multi-energy ions has been applied to obtain uniform damage in the samples, in order to simulate the α-decay-based radiation environment of the vitrified nuclear waste. The irradiation experiment was performed at the Institute of Modern Physics in Lanzhou. The highly energetic Xe ion beams of different energy (1.6, 3.2, and 5 MeV) were induced/selected from an electron cyclotron resonance (ECR) ion source and then redirected to the surfaces of the glass. With ions of multi-energy, uniform damage was created in the sample at room temperature and at the pressure of 7 × 10 −5 Pa. The penetration range of incident ions was within 2 μm. The used irradiation setup guarantees dominance of nuclear collision processes over electronic effects, 35 which is important for the comparison of the experimental results with the ballistic events-based simulation approach. The dose in displacements per atom (dpa) was derived using standard procedure outlined in detail in our previous studies of irradiation effects on borosilicate glasses. 18 After irradiation experiments, the nanoindentation measurements were carried out at the Suzhou Institute of Nanotechnology and Nano-bionics with the MTS G200 nanoindenter with a three-sided pyramid tip (Berkovich diamond tip). The maximum load on the indenter was ~500 μm and to avoid the influence of glass surface, the indentation at depth larger than 0.3 μm was investigated. These measurements were done in continuous stiffness mode. In Yuan et al. 30 and Guan et al., 19 we reported the experimental results of the properties of irradiated NBS1 and NBS2 (sodium (Na)-Boron-Silicate) glasses. In addition to those measurements, here we report results of irradiation experiments on one more glass: NBS4 a . These glass compositions were selected so all three have constant molar ratio of [

| Computational methods
The atomistic simulations were performed with the aid of force field-based molecular dynamics simulations using the LAMMPS 37 and GULP 38 codes, with the later used to derive the elastic properties of the simulated glass structures. Below we provide detailed description of the applied simulation procedures.

| Interaction potentials
A key factor for successful simulation-based studies is the application of a reliable description of interactions between atoms constituting the simulated medium. Because the radiation damage simulations require simulation of a large number of atoms (usually thousands), ab initio simulations are prohibited and classical molecular dynamics with a simple force field has to be applied. In our studies, we used the Buckingham-type interatomic, pair interaction potential of the form: where r is the distance between the interacting atoms and A ij , ρ ij , and c ij parameters for different pairs of interacting i and j atoms are provided in Table 2. These parameters come from the studies of Guillot and Sator 39 and were fitted to reproduce the thermodynamic, structural and transport properties of natural silicate melts. For the interaction between B and O atoms, we used the force field developed recently and tested on borosilicate glasses by Wang et al. 7 The parameters for the Zn-O interaction were fitted using GULP code, 38 so the interaction potential reproduces the measured structural parameters of ZnO oxide. We note that the used force field parameterization differs from the widely used Kieu et al. 40 scheme, but according to Wang et al. 7 it substantially improves the description of density and boron speciation of borosilicate glasses.

| Simulations of glass structure
In all our simulations, the glasses were represented by models containing ~3000 atoms, with correct number depending on the glass composition (see Table 1 for details). This size of the simulated systems is consistent with the number of atoms used in simulation studies by Wang et al. 7 The glass structures have been virtually generated following the procedure applied by Wang et al. 7 The procedure is illustrated in Figure 1. We started with a randomly distributed set of atomic positions and such a system has been initially equilibrated for 100 ps at a temperature of 3000 K using NPT (constant pressure-temperature) ensemble. Then, we slowly cooled the system down to 300 K by applying the cooling rate of 1 K/ps. Such simulations ran for 2700 ps. The cooled glasses were further equilibrated for 100 ps at a temperature of 300 K using NPT ensemble. For some alkali metal-rich NBS glasses, we also used NVT (constant volume-temperature) ensemble for cooling, because of problems with obtaining correct glass densities with NPT approach. This procedure follows recommendation of Wang et al., 7 who noticed that for some glass compositions, NPT ensemble results in uncontrollable expansion in volume. We found that this happens for glasses rich in Na due to formation of a Na-oxide phase, which at high temperatures tends to separate and drastically expand in volume. This could be avoided by lowering the equilibration temperature to below 2000 K. Such a low temperature, however, does not guarantee a well equilibrated structure. The simulations were performed using LAMMPS simulation package. 37

| Simulations of elastic parameters
The elastic Young's modulus as well as bulk (B) and shear (G) moduli have been simulated on a set of 10 snapshots uniformly selected along the equilibrated trajectory. We used the Voigt-Reuss-Hill approach 41 and the simulations were performed using GULP code. 38 The hardness has been estimated as Vickers hardness, H v , from the following empirical formula derived by Chen et al. 42 for polycrystalline materials and bulk metallic glasses: In the above equation, k = G/B. We notice that it is not intuitive to expect such a simple relationship between hardness, which describes resistance of a material to a localized plastic deformation, and elastic moduli that describe the material's non-permanent elastic resistance to the applied stress. Here we only test if such a simple empirical relationship holds for the measured data on the considered pristine and irradiated borosilicate glasses.

| Simulations of accumulation of radiation damage
The radiation damage has been simulated using the defects accumulation procedure applied in previous atomistic simulation studies. 43,44 We performed the damage accumulation molecular dynamics runs in an iterative loop. Each single cation defect formation run was 2 ps long. It consisted of NPT equilibration simulation run (assuming ambient condition), which was followed by a random direction and distance displacement of a randomly selected cation in such a way that the displacement distance was at least 6 • A. This is essential to assure formation of a permanent defect, to which the oxygen sublattice effectively readjusts. Such a simple procedure of gradual defect accumulations has been applied before, for instance, in simulation of radiation damage effects in ceramic materials. 43,45 As a part of our studies, we intended to check if this approach could correctly predict the critical irradiation dose, stored internal energy, or change in Young's modulus Step 1 NVT Step 2 Step 4 Step 5

NPT NPT NVT
Step 3 NPT of the measured glasses. We note, however, that in the simulations that are based on simple interaction potentials, we account only for the ballistic effects. Although used experimental setup guarantees dominance of collision events, 35 in reality, electronic structure effects may also impact the response of a material to irradiation. Due to very short simulation times, the simulated dose rate is also orders of magnitude larger than the one realized in experiments 44 and the overall irradiation effects may depend on type and energy of the applied radiation.

| Simulated system
We simulated three different series of borosilicate glasses that reflect the measured compositions of: (G1) Wang et al., 7 (G2) Zhao et al., 32 and (G3) Guan et al. 19 The exact compositions are given in Table 1 7 for silicon-rich glasses, and both simulations overestimate the glass density for these compositions. The offsets between the three datasets could result from slightly different boron speciation in both simulations and experiments. As indicated in Figure  3, in our simulations, there is less B [4] species (by up to 10 % of total B [4] content) for boron-rich cases. We carefully discuss this discrepancy with the aid of the existing theoretical models.
The formation of B [4] species occurs through the presence of alkaline metals, Ca, and Na atoms in our case. These atoms act as charge-compensating network modifiers. 36 Thus, one should expect that the number of B [4] species will be limited by the amount of available alkali elements. For the case without silicon, we could estimate this limit as: For the G1 75B glass, this results in B [4] content of 0.3333. The value of 0.4 reported by Wang et al. 7 exceeds this limit (Figure 3).
(3) [B [4] [4] for Si-rich compositions. This is understandable as fraction of the alkali metals should participate in saturation of Si-O chains. 36 According to Smedskjaer et al., 36 there exist three models that describe the B [4] content as a function of composition in a borosilicate glass, similarly to Equation (3). The ideal counting (IC) model 36 assumes the equal chargecompensating effect of modifiers (Na and Ca 2+ ions) on B [4] speciation, and B [4] occurs in corner-sharing pairs that are fully bounded (connected) to B [3] or SiO 4 tetrahedra units. In that model, the excess of alkali atoms create non-bridging oxygen (NBO) on these B [3] or SiO 4 units. As discussed by Smedskjaer et al., 36 [4] content but still overestimates the amount of B [4] in G1 glasses at high Si content. Both models predict experimentally unseen saturation (as limit of B [4] = 1 has not been reached experimentally). The experimental results indicate competition between the formation of B [4] species and NBO for high Si content/low B content. 36 It is accounted for in a two-state model of Smedskjaer et al., 36 where the enthalpy difference (ΔH) between the formation of B [4] species and NBO on SiO 4 is taken into account. The results of all three models for G1 glasses are reported in Table 3 and Figure  3, which are equivalent to Figure 6 of Smedskjaer et al. 36 Interestingly, the simulated data could be well fitted with the two-state model assuming ΔH = 0.072 eV (as proposed by Smedskjaer et al. 36 ). The fact that all the models reproduce well the simulated data at low alkali content validates our simulations and points toward potential overestimation of B [4] content in data of Smedskjaer et al., 36 as also indicated in Figure 6 of Smedskjaer et al. 36 In the next step, we performed simulations of the series of SNBS G2 glasses investigated by Zhao et al. 32 Zhao et al. 32 measured and simulated (with the similar force field-based T A B L E 3 Model-predicted, simulated, and experimental fraction of B [4] and density (in g/cm 3  method) the Young's modulus and hardness of these materials. However, in these simulation studies, the potential of Kieu et al. 40 was used, which according to Wang et al. 7 cannot correctly capture the density and boron speciation of borosilicate glass. The results of simulations of density and boron speciation with the Wang et al. 7 potential are reported in Figures 4 and 5. Our simulations predict a density of ~2.40 -2.47 g/cm 3 for this class of glasses. This is by 10% larger than the experimental densities reported by Zhao et al. 32 However, our result is very close to other experimental and simulation studies by Barlet et al. 49 of glass of similar composition, and in general within a density range expected for borosilicate glasses. 49 In order to shed light on this discrepancy, we applied the models of glass densities as a  Table 3 and in Figure 4. We note that the simulated densities are consistent with the model of Inoue et al. 52 On the other hand, model of Budhwani and Feller 51 is consistent with the measured densities, but as discussed by Barlet et al., 49 it, in general, underestimates the densities. With these results, we thus suspect that the experimental densities reported by Zhao et al. 32 are somehow underestimated or that the real parameters of the investigated glasses were slightly different. The predictions of the models for the fraction of B [4] and density, together with the experimental data are reported in Table 3.

) of G1-G3 glasses
The results obtained for the two series of glasses (G1 and G2) indicate that the applied simulation procedure results in production of glass structures with reasonable density and boron speciation, which is important in terms of analysis of the elastic properties of glasses and their performance under irradiation.

| Elastic properties
The simulated and measured Young's modulus for the series of G2 SNBS glasses of Zhao et al. 32 are reported in Figure  6. The simulations show the Young's modulus of ~65-70 GPa, with no clear dependence on the glass composition (B content). This is in contradiction to the experimental results, which show much lower modulus for B-rich compositions (~45 GPa), but similar for Si-rich glasses, ~70 GPa. Interestingly, for Si-rich compositions, we observe the best agreement between the measured and simulated densities (Figure 4), and our simulated values are well consistent with the simulation results of Zhao et al. 32 In Figure 7,  where ρ is the density given in g/cm 3 . The simulated data show much less variation with density, because of nearly identical densities ( Figure 6; Table 3). The difference between the simulated and measured Young's modulus must reflect the differences in the measured and simulated glass structures, namely B [4] speciation. The simulated boron speciation is given in Figure 5 and Table 3, and it is well consistent with the models discussed in Section 3.1.1. The hardness of G2 SNBS glasses estimated as Vickers hardness from the computed bulk and shear moduli is provided in Table 4. We note that Zhao et al. 32 performed direct molecular dynamics simulations of hardness of these glasses and obtained values that are one order of magnitude larger than the measured values. Both datasets are provided in Table  4. On the other hand, our estimate matches reasonably well the measured values. This indicates that a simple, empirical relationship between hardness and elastic moduli derived by Chen et al. 42 for polycrystalline materials and bulk metallic glasses can be successfully applied to estimate the hardness of glassy materials. On the other hand, our analysis indicates some problems with the data on hardness simulated explicitly by Zhao et al. 32 We note, however, that such direct simulations as performed by Zhao et al. 32 are not trivial and several factors can impact the quality of the simulated values.

| Irradiated glasses
Here we intended to perform simulations of irradiation of series of G3 glasses, in order to interpret experimental data published recently 19,20,21,30 and produced within the scope of this paper (NBS4 glass). The computational procedure is tested on data on irradiated PNL 76-68 glass. 17,46 The presented comparison between the experimental and simulation results should be considered having in mind all the limitations of the simulation method and difference in simulated conditions from those realized in experimental irradiation studies, as mentioned in Section 2.2.4.

| Properties of PNL 76-68 glass
In order to validate the simulation procedure for simulation of irradiated glasses, we simulated the PNL 76-68 glass for which the irradiation-induced (stored) internal energy and volume change have been experimentally measured. 17,46 We selected this glass composition as this is the least compositionally complex of the measured glasses and the required potential parameters are known for all the elements, except for Zn (which we fitted here). However, this glass is a good test case because experimentally it shows uniquely negligible volume (density) change and significant stored internal energy (of ~90 J/g; Figure 8). We thus tested how our computational setup reproduces these data. The results are given in Figure 8. Indeed, our simulation shows negligible volume/ density change (within 1%) and the simulated stored energy of ~80 J/g is also very consistent with the experimental value,  glasses. The results on Young's modulus are given in Figure  9. Interestingly, the data show saturated accumulation of radiation damage at critical irradiation dose of ~0.1 dpa. In Figure 9, we also present the experimental data for one more glass compositions (NBS4). All of the considered compositions have similar K parameter, which allows for investigation of the impact of silicon content on the property of the irradiated glass. Regarding the elastic and mechanical properties, on the experimental side, the irradiation causes the decrease in the Young's modulus and hardness. For instance, the Young's modulus of NBS2 glass decreases from 85 GPa (unirradiated case) to 70 GPa. In view of our analysis of behavior of borosilicate glasses provided in previous sections, this indicates decrease in glass density and B [4] content. In order to understand these effects, we performed molecular dynamics simulation of irradiation process.

| Properties of G3 glasses: Simulations
The simulated properties of irradiated G3 NBS1, NBS2 and NBS4 glasses are reported in Figures 9-12 and Tables 5 and  6. First, we report the simulated change in Young's modulus.
Interestingly, the Young's modulus of the pristine glasses is substantially smaller than the measured value (~50 GPa vs. ~90 GPa), and the simulations show increase in the modulus during irradiation (Figure 9). On the other hand, experiments and simulations give consistent value of ~70 GPa for irradiated glass after the critical dose of 0.1 dpa. The difference seen between the experimental and simulated cases could be explained by the difference in measured and computed glass density and boron speciation. As indicated in Figure 11a, the density of the "virtually" irradiated glasses increases by ~10%. As a result, the simulated Young's modulus increases by ~50%, from ~50 to 75 GPa. This is consistent with the model (Equation 4). The increase in density is accompanied by increase in the B [4] content (see Table 5). The increase in the Young's modulus in simulated irradiated glass is thus clearly correlated with the increase in B [4] content of the irradiated glass. However, we noticed that the simulated densities of G3 glasses are inconsistent with the best theoretical prediction (ρ in in Table 3), and the simulated glasses have substantially lower densities. The discrepancy is more pronounced for the glass compositions with large sodium content. In fact, we noticed some convergence problems with simulations of NBS1 glass (highest sodium content) and suspect that high abundance of sodium affects the simulated cooling procedure (see Section 2.2.2). This results in not well equilibrated final glass structures. In order to correct for the density effect, we thus performed the simulations by applying NVT ensemble for cooling (glass preparation), assuming glass density provided by the model of Inoue et al. 52 The results with so "virtually" prepared glass are provided in Figures 10 and 11b, and Table  6. Now, we get much better agreement with the experiment regarding Young's modulus of pristine glasses. The simulated properties of irradiated glasses, however, are not that sensitive to the simulation setup (Figures 9-12 In Figure 12, we provide information on dependence of the Young's modulus on the glass density during the irradiation simulation runs. We observe strong correlations that resemble these seen experimentally ( Figure 6 and Equation  4). Both experiment and simulation thus show linear-like increase of Young's modulus with increase of density. The Vickers hardness estimated from the simulated elastic moduli is also well consistent with the measured values (Tables  5 and 6).
In Figure 11, we present the simulated stored internal energy of G3 glasses. We notice that, although the two applied glass preparation schemes result in slightly different densities and structure of pristine glasses (Tables 5 and 6) and thus resulting stored internal energy, the density of irradiated glasses and critical irradiation doses are nearly identical. The obtained values are also qualitatively consistent with the experimental and simulated results for PNL 76-68 glass ( Figure 8).

| CONCLUSIONS
We performed extensive simulations and experimental measurements of the structural and elastic parameters of series of borosilicate nuclear glasses under irradiation. With the larger composition-dependent set of simulated and measured glass parameters, we performed in depth analysis of the structure-property relationships in these materials, focusing on glass density, boron speciation, elastic Young's modulus and hardness. We found that the properties of borosilicate glasses are mainly determined by the glass elemental composition and the density. The simulated glass density follows well the prediction of the model by Inoue et al., 52 and the boron speciation, namely the content of B [4] , is well described by the two-state model of Smedskjaer et al. 36 We found a clear linear relationship between density and Young's modulus. The hardness of the measured glasses is also reasonably well estimated from the knowledge of the bulk and shear moduli by a simple empirical relationship proposed by Chen et al. 42 This is an interesting result, as it is rather "non-intuitive" that hardness could be estimated from parameters that describe only the resistance of a material to elastic deformations.
The simulations of glasses under irradiation show that a simple defect accumulation procedure can capture well the irradiation process, correctly predicting the critical irradiation dose (~0.1 dpa), change in density and elastic/ mechanical properties, and give reasonable estimate for the stored internal energy. The simulations results and the existing models for glass density and B [4] content indicate potential problems with characterization of the measured samples. This shows the power and importance of atomistic simulations of glassy materials for correct interpretation of the experimental data.