⇦Chinese Physics Letters, 2020, Vol. 37, No. 10, Article code 106101 Molecular Dynamics Simulations of the Interface between Porous and Fused Silica Ye Tian (田野)*, Xiaodong Yuan (袁晓东), Dongxia Hu (胡东霞), Wanguo Zheng (郑万国), and Wei Han (韩伟) Affiliations Laser Fusion Research Center, China Academy of Engineering Physics, Mianyang 621900, China Received 25 May 2020; accepted 11 August 2020; published online 29 September 2020 *Corresponding author. Email: tianye8911@outlook.com Citation Text: Tian Y, Yuan X D, Hu D X, Zheng W G and Han W et al. 2020 Chin. Phys. Lett. 37 106101    Abstract Molecular dynamics simulations are performed to gain insights into the structural and vibrational properties of interface between porous and fused silica. The Si–O bonds formed in the interface exhibit the same lengths as the bulk material, whereas the coordination defects in the interface are at an intermediate level as compared with the dense and porous structures. Clustered bonds are identified from the interface, which are associated with the reorganization of the silica surface. The bond angle distributions show that the O–Si–O bond angles keep the average value of 109$^{\circ}$, whereas the Si–O–Si angles of the interface present in a similar manner to those in porous silica. Despite the slight structural differences, similarities in the vibrations are observed, which could further demonstrate the stability of porous silica films coated on the fused silica. DOI:10.1088/0256-307X/37/10/106101 PACS:61.43.Gt, 61.43.Bn, 68.35.-p, 61.43.Fs © 2020 Chinese Physics Society Article Text Due to its large surface-to-volume ratio, light weight and unique physical properties, sol-gel porous silica is widely used in various scientific and industrial applications such as biomedical and chemical engineering. In particular, the excellent laser damage resistance in the ultraviolet regime makes this material the only choice as antireflection coatings used for the final optical elements of inertial confinement fusion (ICF) facilities.[1] In the complex radiation environment with laser and particle beams, the stability of interface between porous silica and the base material remains an ongoing topic as it is associated with the output capability of laser systems.[2,3] Recently, the films coated on fused silica have drawn broad research interest in optics because they turn out to be more stable than those on third-harmonic-generation crystals upon laser exposure, which inspires studies on the relevant interfaces consequently. While there have been extensive experimental and numerical research on the damage mechanism of glasses and films, fundamental knowledge of their interfacial structure is still limited on account of the amorphous nature of these materials. As an effective tool to investigate the structure and thermomechanical properties of materials, molecular dynamics (MD) simulations have been applied to the studies of silicates for decades.[4–6] A considerable number of force fields, such as the BKS, Teter and ReaxFF, have been successfully demonstrated to feature atomic interactions in the amorphous network with acceptable computation costs.[7,8] Notwithstanding the slight differences in depicting the disordered structures, reasonable agreements with experiments have been achieved to validate the modeling capability of these MD potentials. Based on the simulations of vitreous silica, procedures to mimic porous silica have also been established.[9–12] For example, Kieffer et al. and Nakano et al. have reproduced the porous structures by successive or one-time volume expansion of vitreous silica, whereas channeled pores as expected in sol-gel processing have been generated by Beckers et al. by scaling the ionic charges.[11,12] Up to date, numerical models of coating porous silica on the surfaces of vitreous silica have not been universally studied yet. In this regard, MD simulations are performed in this work to explore the interface established by porous and vitreous silica. Consisting of 2250 atoms, the simulated system was generated in steps of 1 fs by the Teter force field, a pairwise potential with fixed partial charges of $+2.4e$ and $-1.2e$ assigned for Si and O, respectively.[7] To begin a 1500-atom vitreous silica bulk was constructed by the conventional melt and quench procedure.[5] With 500 Si$^{+2.4}$ and 1000 O$^{-1.2}$ atoms randomly inserted into a $28.3 \times 28.3 \times 20$ Å$^{3}$ periodic cube, the 2.2 g/cm$^{3}$ supercell was thermostated by the Nose–Hoover algorithm to slowly scaling the temperature from 5000 K to 300 K in the NVT ensemble. The cooling rate was set as 5 K/ps that is slow enough to reduce the concentration of coordination defects. Then, the system was relaxed at 0 GPa in the NPT ensemble, resulting in a final density of 2.3 g/cm$^{3}$ that is close to the silica density under ambient conditions. Based on this bulk structure, a vacuum layer of equal thickness was induced to one surface along the $z$-direction, whereas the atoms in the center slab of the bulk were excluded from the following integrations to build a 10 Å thick frozen layer. To construct steady surfaces, the other atoms on both sides of the frozen layer were cooled from 1000 K to 300 K, followed by relaxations in the NVE ensemble to reach equilibrium. Due to the periodic boundary conditions, the vacuum layer was also imposed on the opposite surface to which it is added, leading to formation of two surfaces instead. Upon creation of the silica slab (denoted as dense silica), the porous film was simulated by the charge-scaling method, which has been considered as a crude analog of sol-gel processing.[12] To represent the sol-gel films applied to ICF final optics, the vacuum layer was filled with 250 Si and 500 O atoms to generate a 1.15 g/cm$^{3}$ homogeneous liquid phase so that the porosity was scaled to 50%. The initial ionic charges were set to be $0.48e$ for Si and $-0.24e$ for O, which were successively increased to $2.4e$ and $-1.2e$ in the NVT ensemble to allow for atomic clustering with enhanced columbic interactions at 300 K. In this process, the appended atoms indeed experienced interactions with the other atoms in both the porous liquid and the dense silica surface. With the effects of pressure on the porous silica and interface excluded by extensive relaxations in the NPT ensemble, trajectories of the system were recorded in the NVE ensemble for further analysis.
cpl-37-10-106101-fig1.png
Fig. 1. An atomic view of the porous structure on silica surface (pink Si and blue O for dense silica, purple Si and green O for porous silica).
Figure 1 depicts a snapshot of the final structures for the simulated system. As can be identified from the amorphous network, the interfaces are sandwiched by the porous and dense phase. A few atoms from the porous silica, however, are inserted in the bulk of dense silica. Due to the presence of voids, the number density of atoms belonging to the porous part shows uneven characteristics, given the observation that there are more atoms appended to the silica surface than in the porous bulk. Despite the different densities, the porous and dense silica both consist of SiO$_{4}$ tetrahedrons as fundamental structural units, which expose no Si atom to the surfaces as a result. Separated by the porous phase, the two interfaces manifest themselves as Si–O pairs with atoms provided by different species. In detail, a bonding pair in the interface can be an O atom from the porous silica group linked to a Si atom in the dense silica group, or vice versa. As will discussed below, these Si–O bonds can be classified by a cutoff distance from the pair distribution functions. The pair distribution functions for the dense and porous silica bulk in the mixed structure are shown in Fig. 2 to reveal changes in the first coordination shell at the mean separation of atoms. In the 0–4 Å range, the primary peak positions for porous silica are identical to those for dense silica, indicating the same short-range order for both systems. Representing the bond lengths or nearest distances, the first peaks for the Si–O, O–O and Si–Si pairs in the mixed structure are at 1.6 Å, 2.6 Å and 3.2 Å, which are all consistent with the previous experimental values for silica.[13] Cutting off at the first minimum of the primary peak for each Si–O pair distribution, the maximum length to determine a Si–O bond was chosen to be 2.2 Å. According to this criterion, it is found that there are a total number of 249 atoms forming the interfaces, and the number may vary depending on the simulation procedures. On the other hand, the average length of Si–O bonds formed by the interfacial atoms is calculated to be 1.6 Å, showing no difference from those bonds in the silica bulk. Therefore, it turns out from the simulations that the bond lengths remain unchanged during formation of interfaces.
cpl-37-10-106101-fig2.png
Fig. 2. Pair distributions for porous and dense silica.
Despite the resemblance of the peak positions, the cross-sectional area of the primary peaks for porous silica appears to differ from that for dense silica, as can be inferred from the different amplitudes and widths shown in Fig. 2. In this regard, growth of porous silica on the surfaces is supposed to generate a defective structure, given that integrations of these peaks over the covered distances are known to depend on the average coordination number. To probe the concentration of coordination defects, the bonding state of each atom is calculated. As is expected, in the dense silica bulk the Teter potential reproduced a realistic coordination environment for Si and O, with 99.3% of the Si being 4-fold coordinated and 97.9% of the O being 2-fold coordinated. For porous silica, there is a larger proportion of over- and under-coordinated atoms in the bulk, where 89.7% of the Si and 90% of the O are normally coordinated. These defects corresponding to oxygen-excess or oxygen-deficiency environments mainly include 3-fold and 5-fold coordinated Si as well as 1-fold and 3-fold coordinated O. By contrast, the atoms forming the interfaces tend to connect with each other in an intermediate manner, given that 96.3% of the Si and 92.1% of the O are identified to be normally coordinated. Except for a slight decrease of Si–O bond length for the nonbridging oxygens that has been reported before, these oxygen-related coordination defects show no significant structural differences as compared with the normal-coordinated atoms.[3] If viewed at the electronic level by first-principles studies, however the under-coordinated or over-coordinated atoms will tend to form defects that have different electronic structures and can bring about enhanced absorption to further degrade the damage resistance of films.[14] Another important feature of the silica interfaces that is related to the structural stability is the clustering of bonds, which can be considered as the number of atoms in one species linked to more than one atom in the other species. From a geometric point of view, clustering corresponds to the intersection of Si–O bonds at the interface. Among all the interfacial atoms, it is found that there are 35 Si and 64 O atoms exhibiting clustering. While the porous silica contributes to the clustered Si atoms by 42.8%, 68.8% of the clustered O atoms come from the dense bulk. Since there are no 2-fold coordinated Si in the original silica slab, the presence of Si atoms that participate in the clusters indicates that in the sol-gel process the SiO$_{4}$ tetrahedrons on the silica surface can be unfolded to accommodate atoms from the film, which will result in a stable connection by the covalent bonds. Similarly, while an undercoordinated O defect can naturally be bonded to the appended Si atom, the O atoms in the dense silica in which each link two Si atoms from the porous bulk in turn validate the reorganization of the silica surface. Upon the formation of the interface, both the lengths and angles of such clustered Si–O bonds have been calculated, which indeed show no difference from the other interfacial Si–O bonds. In addition, for the few atoms from the sol-gel that dive below the silica surface, further computations of the bulk structure containing such atoms show that these bulk clusters just form in the same way as typical SiO$_{4}$ structural units. Assuming the same cutoff of 2.2 Å, the medium-range order is expressed in terms of Si–O–Si and O–Si–O bond angles. The bond angle distributions are shown in Fig. 3 for Si–O–Si and O–Si–O, which are used to describe the connectivity and regularity of tetrahedrons, respectively. While the average O–Si–O bond angle keeps its value at 109$^{\circ}$, pronounced differences in the distribution of Si–O–Si angles can be observed. For the dense silica bulk, the Si–O–Si distribution peaks at around 150$^{\circ}$, which is in reasonable agreement with the experimental measurements.[15] For porous silica, however, the primary band are split into two neighbored peaks at 136$^{\circ}$ and 156$^{\circ}$, which are accompanied by the small peaks around 90$^{\circ}$ that is indicative of 2-membered rings. Based on this distribution, the average value of Si–O–Si angles is calculated to be 145.3$^{\circ}$ for porous silica, smaller than that computed for dense silica. For the interfacial atoms, the average Si–O–Si angle is calculated to be 146.5$^{\circ}$, which thus suggests that the bond angles of the interface are organized in a similar manner to porous silica.
cpl-37-10-106101-fig3.png
Fig. 3. Bond angle distributions for porous and dense silica.
In addition to the static structures, the dynamic property of the amorphous system is explored with the partial vibrational density of states that shows the distribution of phonon modes. Typically, the experimental spectrum measured by neutron scattering reproduces three peaks at 400 cm$^{-1}$, 800 cm$^{-1}$, and 1100 cm$^{-1}$, which represent the rocking, bending, and stretching modes, respectively.[16] In the simulations, the vibrations are analyzed by diagonalization of the dynamical matrix.[17] To compute the Hessian matrix, each atom was displaced from its equilibrium position with a small magnitude of 0.01 Å along the positive and negative Cartesian directions to obtain the relevant force constant elements. Given the large number of atoms, only the gamma point is considered in the simulation. The eigenvectors of the dynamic matrix are broadened with Gaussian functions so that a continuous curve can be achieved. The spectra generated for the dense silica, porous silica and their interfaces are shown in Fig. 4 for comparison, each of which includes a broad band around 100–700 cm$^{-1}$ and a peak split near 1100 cm$^{-1}$. Although the high-frequency states at 1100 cm$^{-1}$ are not pronounced for the interface, the three structures share the same features in the positions that indicate the band edges and peaks. Moreover, it is found that the O atoms dominate most vibrational modes except those near 700 cm$^{-1}$ that are governed by Si instead. The resemblance in the atomic vibrations indicates that vibrations in silica can be smoothly transferred to the porous bulk without bringing about significant surface effects that hinder stability of the interface, thereby suggesting a smooth structural transition.
cpl-37-10-106101-fig4.png
Fig. 4. Partial vibrational density of states for porous and dense silica.
In summary, porous silica coated on the surface of vitreous silica has been investigated with molecular dynamics simulations in this work. The Si–O covalent bonds bridged by the atoms from different species show no difference in bond length, whereas the interfacial Si–O–Si bond angles are found to resemble those in porous silica. The content of coordination defects in the interface lies at a level between those in the dense and porous bulk material, whereas surface reorganization can be inferred from the clustered bonds formed by the interfacial Si and O atoms. As compared with the silica bulks, similarities in the vibrations have been observed for the interface, as evidenced by the vibrational density of states. This study could offer an atomistic level picture of structural stability of porous silica films on fused silica optical elements. References Large Optics for the National Ignition FacilityMolecular Dynamics Study of the Structural Modification of Porous Silica from Low-Energy RecoilsDensification effects on porous silica: A molecular dynamics studyStructural features of sodium silicate glasses from reactive force field‐based molecular dynamics simulationsAn atomic-level look at the structure-property relationship of cerium-doped glasses using classical molecular dynamicsRevisiting silica with ReaxFF: Towards improved predictions of glass structure and properties via reactive molecular dynamicsThe medium range structure of sodium silicate glasses: a molecular dynamics simulationStructural and Mechanical Properties of Nanoporous SilicaFractal dimensions of silica gels generated using reactive molecular dynamics simulationsGrowth of Pore Interfaces and Roughness of Fracture Surfaces in Porous Silica: Million Particle Molecular-Dynamics SimulationsMolecular dynamics simulation of nanoporous silicaNeutron scattering from vitreous silica. V. The structure of vitreous silica: What have we learned from 60 years of diffraction studies?Stable structure and optical properties of fused silica with NBOHC- E ′ defectNMR determinations of Si O Si bond angle distributions in silicaComparison of the neutron, Raman, and infrared vibrational spectra of vitreous Si O 2 , Ge O 2 , and Be F 2 First principles phonon calculations in materials science
[1] Baisden P A, Atherton L J, Hawley R A et al. 2016 Fusion Sci. Technol. 69 295
[2]Ruello P, Vaudel G, Brotons G et al. 2017 Proc. SPIE 10447 1044717
[3] Zhang J H, Tian Y, Han W et al. 2019 Chin. Phys. Lett. 36 116102
[4] Tian Y, Du J, Hu D et al. 2018 Scr. Mater. 149 58
[5] Deng L, Urata S, Takimoto Y et al. 2020 J. Am. Ceram. Soc. 103 1600
[6] Pedone A, Tavanti F, Malavasi G et al. 2018 J. Non-Cryst. Solids 498 331
[7] Yu Y, Wang B, Wang M et al. 2016 J. Non-Cryst. Solids 443 148
[8] Du J and Cormack A N 2004 J. Non-Cryst. Solids 349 66
[9] Rimsza J M and Du J 2014 J. Am. Ceram. Soc. 97 772
[10] Bhattacharya S and Kieffer J 2005 J. Chem. Phys. 122 094715
[11] Nakano A, Kalia R K and Vashishta P 1994 Phys. Rev. Lett. 73 2336
[12] Beckers J V L and de Leeuw S W 2000 J. Non-Cryst. Solids 261 87
[13] Wright A C 1994 J. Non-Cryst. Solids 179 84
[14] Lu P F, Wu L Y, Yang Y et al. 2016 Chin. Phys. B 25 086801
[15] Pettifer R F, Dupree R, Farnan I et al. 1988 J. Non-Cryst. Solids 106 408
[16] Galeener F L, Leadbetter A J and Stringfellow M W 1983 Phys. Rev. B 27 1052
[17] Togo A and Tanaka I 2015 Scr. Mater. 108 1