⇦Chinese Physics Letters, 2019, Vol. 36, No. 11, Article code 116102 Molecular Dynamics Study of the Structural Modification of Porous Silica from Low-Energy Recoils * Ji-Hua Zhang (张继花)1,2, Ye Tian (田野)2, Wei Han (韩伟)2, Fang Wang (王方)2, Fu-Quan Li (李富全)2, Xiao-Dong Yuan (袁晓东)2, Xia Xiang (向霞)1** Affiliations 1School of Physics, University of Electronic Science and Technology of China, Chengdu 610054 2Research Center of Laser Fusion, China Academy of Engineering Physics, Mianyang 621900 Received 29 July 2019, online 21 October 2019 *Supported by the National Natural Science Foundation of China under Grant No U1830204.
**Corresponding author. Email: xiaxiang@uestc.edu.cn Citation Text: Zhang J H, Tian Y, Han W, Wang F and Li F Q et al 2019 Chin. Phys. Lett. 36 116102    Abstract Molecular dynamics simulations are performed to investigate the effects of low-energy recoils on the microscopic structure of porous silica. Exhibiting a logistic growth with the recoil energy, the displacement probability of Si is shown to be smaller than that of O at the same primary knock-on level. Computations of pair distribution functions and bond angle distributions reveal that this material upon irradiation with energies around the displacement thresholds mainly undergoes structural changes in the medium-range order. In the porous network, while the formation of nonbridging oxygen defects tends to induce shorter Si–O bonds than those formed by bridging oxygen atoms, a remarkable increase of inter-tetrahedral bond angles created by multiple recoils can be observed and associated with the rearrangement of ring statistics. DOI:10.1088/0256-307X/36/11/116102 PACS:61.80.Hg, 61.43.Er, 61.43.Gt, 02.70.Ns © 2019 Chinese Physics Society Article Text Sol-gel-derived porous silicas are widely used as anti-reflection films in inertial confinement fusion (ICF) facilities due to their excellent damage resistance in the ultraviolet regime.[1] Typically, such a system with $\sim $50% porosity has been shown by off-situ experiments to be damage-free under 6–8 J/cm$^{2}$ irradiation at 355 nm.[2] However, the optical lifetime of porous silica films can possibly be reduced by the prompt neutron radiation from the fusion target, which tends to induce irreversible modifications in the silicate structure. Although the mechanism of laser damage has been studied for decades, this radiation effect on porous silica—which is related to the subsequent response of radiated materials to laser exposure—remains poorly understood. Details of atom displacements around their thresholds in porous silica are thus of fundamental importance to the assessment of the primary damage production, which is the first step towards exploring collision cascades in the radiation process. Determining the radiation effects in porous silica with respect to the structural modification is challenging on account of the complex coordination environment, where the randomly distributed pores could provide channels for an atom to escape from its neighbors. Given the success in atomic-scale modeling, molecular dynamics (MD) simulation is known as an effective tool to investigate the structural and thermomechanical properties of materials. Fitted to experimental measurements or ab initio simulation results, various force fields have been demonstrated to feature atom interactions in silicate systems, such as the BKS, TTAM, and Teter potentials, which makes this numerical method a convenient approach to model neutron radiation damage in nuclear materials, especially for low-energy recoil events wherein electronic stopping is not significant.[3] Up to now, there have already been several MD studies of neutron-induced recoils in vitreous silica, which is a commonly used material that also approximately serves as the backbone of porous silica. For instance, Mota et al. have confirmed the presence of various coordination defects in the radiated silica structure, wherein permanent densification of vitreous silica has been revealed by Kubota et al.[4,5] Following previous studies, MD simulations have been performed in this work with the large-scale atomic/molecular massively parallel simulator (LAMMPS) code to gain insights into low-energy displacements and the resultant structural modifications in porous silica. Generated in MD steps of 1.0 fs, a 6000-atom system of porous silica was created by the charge-scaling method with the Teter potential.[6] To begin 2000 Si$^{+0.48}$ and 4000 O$^{-0.24}$ atoms were randomly inserted into a $56.57\times56.57\times56.57$ Å$^{3}$ periodic supercell, in which way the initial density was chosen as 1.1 g/cm$^{3}$ to eventually achieve a porosity of around 50%, close to that of the sol-gel films applied on the ICF final optics. Thermalized at 300 K in the NVT ensemble, the ionic charges were successively increased to 2.4$e$ and $-1.2e$ for Si and O respectively to allow for atomic movements with limited columbic interactions. After further equilibration at 300 K and 0 GPa, the system was stabilized at 1.4 g/cm$^{3}$, corresponding to 40% porosity since the density of vitreous silica predicted by the Teter potential is 2.3 g/cm$^{3}$. Given that the structural differences for porous silica of porosities 40%–60% have been shown to be small, the slight overestimation of the simulated density in comparison with the experimental films is thus not problematic.[7] Figure 1 depicts the final structure of the simulated porous silica, which consists of an amorphous Si–O network and several channeled pores. Classified with a cutoff of 2.0 Å, the coordination number of Si and O atoms (marked as red and blue spheres respectively in Fig. 1) are 4.0 and 2.0, respectively, suggesting the presence of Si–O tetrahedrons that are typically formed in various glassy silicates.
cpl-36-11-116102-fig1.png
Fig. 1. An atomic view of the porous structure.
To scale the recoil simulations properly, the threshold displacement energy of porous silica was determined in advance by the binary-search algorithm in the 0–128 eV range. Despite the sufficiency of Teter in reproducing the static silica structure, such a potential is not optimally suited for recoil situations that involve bond breakage and reformation as it assumes an unphysical potential well at ultrashort internuclear distances. Therefore, the bond-ordered reactive force field (ReaxFF) previously developed for Si–O systems was used for the radiation simulations instead, given that recent advances in the researches of silica-based chemical reactions have validated the advantage of the ReaxFF in detailing bond state transitions and defect formations in ab initio qualities.[3] Considering the amorphous nature of porous silica, individual displacements on multiple atoms (100 Si and 100 O) were performed to obtain a probability distribution with the recoil energy, in which process the displacements were conveniently indicated from bond-state transitions. Although a Langevin thermostat was applied to the 5-Å-thick vicinity of the boundaries to maintain a constant system temperature of 300 K, the primary knock-on atom (PKA) was randomly chosen from the inner shell to guarantee recoils within the unthermostatted region, where the atoms moved freely in the NVE ensemble. Beginning from 64 eV, the kinetic energy of the PKA was scaled to be the median of the search range, which would be chopped depending on the subsequent bonding changes in the recoil event. Specifically, if the atoms pre-bonding with the PKA differed from those after recoils, then the upper half in which the threshold cannot lie would be eliminated by decreasing the upper limit to the test value. Otherwise, the lower limit would be increased instead. The search procedure of the displacement threshold then repeated by scaling the recoil energy on the remaining half to narrow the search scope until a precision of 1.0 eV was reached. Note that all of the recoils were simulated for at least 5 ps with a timestep as small as 0.01 fs to achieve converged results. The displacement probability functions for both Si and O are plotted versus the recoil energy in Fig. 2. In general, the probability-based curve for Si exhibits a wider distribution range than that of O atoms. Although no displacement of O occurs below 5 eV, a small number of Si can be pushed away from their respective balance positions by recoil energies as low as 2 eV. Moreover, the displacement probabilities for both species tend to increase rapidly with energy until they start to approach unity. The saturation point for Si, where almost all the atoms can be displaced, is about 15 eV larger than that of O. For most recoil energies, the displacement probability for Si is significantly higher than that of O. Intuitively, this can be further quantified by the energy that corresponds to a 50% displacement probability, which is 12 eV and 20 eV for O and Si atoms, respectively. This observation is unsurprising due to the fact that a Si atom is tightly constrained by 4 Si–O bonds in the center of a (SiO$_4$) tetrahedron, whereas one O atom only bridges two Si atoms. Hence, when a Si–O bond is broken by a specific energy, the corresponding O atom tends to rotate with the remnant bond while the excess bonds with Si can still keep the Si atom steady. Furthermore, the larger atomic mass of Si also contributes to the lower displacement probability.
cpl-36-11-116102-fig2.png
Fig. 2. Displacement probability plotted versus the PKA energy of Si and O atoms.
Once the displacement probability was determined, low-energy recoils were induced by initiating various numbers of atoms simultaneously in the unthermostatted region of the system. The recoil energies for this process were chosen to be around 40 eV, which is the smallest value to produce permanent displacements of each atom. Depending on the kinetic energy for excitation, the atoms were expected to interact dynamically with their surroundings, resulting in breakage and reformation of the amorphous network. Each recoil procedure which typically lasted for several picoseconds was followed by extensive NVT relaxations of about 100 ps at 300 K to stabilize the system. With trajectories traced and averaged over 600 configurations in the NVE ensemble, systemic structural analysis in both the short and medium ranges was performed for all the disturbed systems to gain information on the radiation effects. Figure 3 shows the primary peaks of pair distribution functions for the porous silica with accumulated recoils of various Si or O atoms to reveal changes in bond lengths or nearest pair distances. For further comparison of the recoil effects, the numbers of atoms used as PKA were chosen as 200 and 400 for each species, which turn out to display similar characteristics indeed. For the initial structure, the Si–O curve peaks at 1.6 Å, which is consistent with the bond lengths reported previously in neutron diffraction experiments.[8] After the 200-atom or 400-atom recoils, there is a conspicuous protrusion occurring on the left toe of the Si–O peak, suggesting the presence of shorter Si–O bonds. To account for this phenomenon, direct measurements of the Si–O bond lengths were performed. It is found that this small peak arises from the short bonds formed by Si and non-bridged O that induced due to the recoils, its length is only 1.38 Å in average. For the Si–Si and O–O pairs, on the other hand, the relevant peak positions remain unchanged after the recoils.
cpl-36-11-116102-fig3.png
Fig. 3. Comparison of pair distribution functions of porous silica.
In addition to the pair distribution functions plotted for bonds or nearest neighboring atoms, the bond angle distributions of the radiated porous silica are shown in Fig. 4(a) for Si–O–Si and in Fig. 4(b) for O–Si–O, which describe the connectivity and regularity of SiO$_4$ tetrahedrons, respectively. It is found that the average O–Si–O angle at 103$^{\circ}$ slightly increases by about 2$^{\circ}$ for all the structures, whereas the average Si–O–Si angle increases more significantly from 142$^{\circ}$ to $\sim $$146^{\circ}$. Moreover, the 400-atom recoils appear to generate a larger average Si–O–Si angle than the 200-atom case due to the stronger structural modification. In our previous work, such an angular increase has been shown to originate from both the increasing degree of plane reorientation and bond swapping.[9] Since atom displacements were invoked in the recoil process, it is thus natural to expect an overall change of ring structures in porous silica. Computed by searching the shortest close path that formed by Si–O bonds, the ring size distributions are shown in Fig. 4 for each sample. The initial structure has a Gaussian-type distribution centering around 6, which is similar to the ring distribution of bulk vitreous silica.[10] When recoils occur, the rings tend to shift towards larger sizes. Specifically, there is a prominent decrease of 5–7-membered rings and increase of those larger rings. Hence, the recoils mainly contribute to a repolymerization process that eliminates the small rings. On the other hand, the numbers of 3- and 4-membered rings change only slightly, differing from the previous studies on vitreous silica that radiation could increase these rings instead.[5]
cpl-36-11-116102-fig4.png
Fig. 4. Comparison of bond angle distributions of porous silica.
cpl-36-11-116102-fig5.png
Fig. 5. Comparison of ring size distributions of porous silica.
In summary, low-energy recoils in porous silica have been investigated with molecular dynamics simulations to explore the effect of neutron radiation on the structure of this material. Indicated by the breakage and reformation of Si–O bonds, atom displacements mainly occur with recoil energies of around 12 eV and 20 eV for O and Si, respectively. As the number of permanent displacements increases, the oxygen defects tend to generate shorter Si–O bonds that correspond to the small peak at 1.38 Å of the Si–O pair distribution function, whereas comparisons of the pair distribution between the same species reveal no significant changes in the short range. As evidenced by the increase of Si–O–Si bond angles and large Si–O rings, the recoils are found to mainly influence the medium-range structure instead. This study could offer an atomistic level picture of structural evolutions in porous silica when neutron radiation is induced in optical films during fusion experiments. References National Ignition Facility Laser System PerformanceUV laser conditioning of the sol-gel SiO2 film coated fused silica opticsThermal conductivity of vitreous silica from molecular dynamics simulations: The effects of force field, heat flux and system sizeMolecular dynamics study of structure transformation and H effects in irradiated silicaThreats to ICF reactor materials: computational simulations of radiation damage induced topological changes in fused silicaMolecular dynamics simulation of nanoporous silicaHydrogen gas Mixture Separation by CVD Silica MembraneNeutron scattering from vitreous silica IV. Time-of-flight diffractionDensification effects on porous silica: A molecular dynamics studyThe medium range structure of sodium silicate glasses: a molecular dynamics simulation
[1] Spaeth M L, Manes K R, Bowers M et al 2016 Fusion Sci. Technol. 69 366
[2] Jiang Y, Chen J, Yang K et al 2017 Optik 139 178
[3] Tian Y, Du J, Han W et al 2017 J. Chem. Phys. 146 054504
[4] Mota F, Caturla M J, Perlado et al 2009 J. Nucl. Mater. 386 75
[5] Kubota A, Caturla M J, Stolken J et al 2003 Nucl. Instrum. Methods Phys. Res. Sect. B 202 88
[6] Beckers J V L and de Leeuw S W 2000 J. Non-Cryst. Solids 261 87
[7] Rimsza J M, Du J C and Chen L Q 2014 J. Am. Ceram. Soc. 97 772
[8] Grimley D I, Wright A C and Sinclair R N 1990 J. Non-Cryst. Solids 119 49
[9] Tian Y, Du J et al 2018 Scr. Mater. 149 58
[10] Du J and Cormack A N 2004 J. Non-Cryst. Solids 349 66