Topological pruning enables ultra-low Rayleigh scattering in pressure-quenched silica glass

Silica glass is the most indispensable material in optical communication applications due to its superior optical properties. The transmission loss of silica glass has been reduced over the past 30 years by continuous efforts toward decreasing density fluctuations by lowering of fictive temperature, e.g., through improvements in processing or doping. A recent study has shown that shrinkage of structural voids by hot compression is a promising way to further decrease the loss. However, an atomic understanding of the pressure effect is still lacking. Here, using molecular simulations, we connect the void shrinkage to topological pruning of silica network. Two physical models predict that the Rayleigh scattering loss of pressure-quenched silica glass can be reduced by >50% when the glass is quenched at an appropriate pressure (4 GPa in our simulation). Our studies are consistent with available experimental results and demonstrate topologically optimized structure can give desirable properties for optical applications of silica as well as other glasses with similar network structure.


INTRODUCTION
Silica (SiO 2 ) glass is the single most important material for fiber optic communications due to a combination of several favorable properties 1 including high optical transparency, high tensile/ bending strength, good processability, ability to tune refractive index by doping, and low cost. Because of this unique combination of properties, silica glass has revolutionized global telecommunications through the development of high-bandwidth optical communication networks. During the 1970s, research was vigorously conducted to reduce the loss in silica glass by a factor of a hundred [2][3][4] . A 10-20% further decrease in optical attenuation has been achieved, e.g., in fluorine-doped silica glass fiber 4,5 . Interests in low optical loss are also seen in SiO 2 -based waveguides that are used in planar lightwave circuits [6][7][8] . Notwithstanding the extreme importance of SiO 2 glass in these applications, there is a lack of atomistic level understanding of optical attenuation in the glass. In addition, the lack of recent progress in reducing optical loss is due to limitations in the processing methods that have been considered. More than 80% of the optical loss in silica glass fiber is due to Rayleigh scattering 3,9,10 , which can be suppressed by reducing the density fluctuations in glass, e.g., by reducing its fictive temperature. The fictive temperature of silica can be reduced by: (a) reducing the cooling rate 1 ; (b) introducing an annealing process 11 ; or (c) reducing the viscosity of the glass by doping with an alkali metal oxide 12,13 . However, reductions in cooling rate or annealing cause a significant decrease in production yield, and compositional variations can lead to an increase in loss due to chemical concentration fluctuations 14 . Recent experimental reports have shown that controlling the pressure history of silica can affect the density fluctuations or heterogeneity and the Rayleigh/X-ray scattering 15,16 . This method of controlling fictive pressure has provide a new train of thought to achieve even lower propagation loss of silica glass. Although the silica glass structural change due to hot compression has been shown to be very different from cold-compressed silica glass with comparable density as well as the pristine silica glass 17 , a systematic understanding of the pressure effect on the silica network structure in terms of homogeneity of the glass at the atomistic level is not available which is required for the design of ideal glass structure with optimal optical properties.
Hence, in this study, we use molecular dynamics (MD) simulations to examine how density fluctuation of SiO 2 glass can be suppressed under conditions of varying pressure, including conditions beyond those reported in the experiments 15 . The MD simulations provide a direct probe to the SiO 2 network structure which can be compared with experimental results for density and void size within the experimentally achievable applied pressures. The model shows that the optical loss can be improved significantly without using an extremely high pressure. We specifically show that pressure quenching leads to an enhanced homogeneity of the silica glass characterized by a decreasing void size which can be understood using topological pruning picture, bringing a significant reduction in Rayleigh scattering intensity bỹ 56% at a quench pressure of 4 GPa in MD. We examine details of the change in homogeneity of SiO 2 glass due to hot compression between 0 and 12 GPa and find that it is controlled by two different mechanisms depending on the magnitude of the pressure. While the homogeneity at large pressure is compromised due to an increase of network connectivity, at lower pressures, the dominant homogenization mechanism is related to the topological pruning of the network, which is found to be an effective way for tailoring of SiO 2 properties without introducing unwanted defects. These findings provide an atomistic understanding of the pressure-induced network homogenization and also shed light on topological optimization of the silica glass network for optical applications.

Structure of hot-compressed silica glass from MD
To demonstrate the effect of quench pressure (P Q ) on silica glass, we have conducted pressure-quench simulations of silica using the recently developed SHIK potential 18 (see "Methods" for justification of this potential and Supplementary Fig. 1a). Density obtained by experiment is well reproduced by MD simulations using SHIK potential as shown in Fig. 1 (also see Supplementary  Fig. 1b). The increase of silica glass density with quench pressure is accompanied by the reduction of structural void size in the glass as illustrated in Fig. 2a. Here, isosurfaces with a minimum distance of 0.1 nm (about the size of the largest voids in α-quartz) from the surface of an atom are visualized as "voids". At P Q = 0 GPa, many of the voids are interconnected in 3D space and their sizes are relatively large. From P Q = 0 to 1.0 GPa, these voids become more disconnected and their sizes become significantly smaller. At P Q = 4.0 GPa, only disconnected voids remain, and the number of the voids (with a radius > 0.1 nm) are few. At P Q = 8.0 GPa, the voids even disappear. Such change persists during the whole cooling procedure (see Supplementary Fig. 2). The void size distributions are shown in Fig. 2b. It shows a rapid extinction of large sized voids (radius > 0.2 nm) from ambient pressure up to 2 GPa. The largest void size in the silica glass decreases significantly with increasing P Q as plotted in Fig. 2c. Despite the fact that the experimental data (red points) from positron annihilation lifetime spectroscopy 15 are only available up to P Q = 0.2 GPa, the largest void size results from MD up to this range show a very good agreement in both the absolute value and the trend.
The significant modification of void morphology and distribution in the pressurized silica glass suggests that some network structure change in the silica glass is triggered by the pressure quench. Here, to characterize the structural change in the SiO 2 framework, we make use of the cluster model 19,20 in which we ignore all O atoms and consider each Si atom and its Si coordination sequence CS = {N k } = N 1 , N 2 , …, N k , where N k indicates the number of Si atoms in "sphere" k which are connected to Si atoms in sphere k − 1. A cutoff of 3.4 Å (see the "Methods") is used to find the neighbor Si atoms of each Si atom. Note that here the N 1 , i.e., CN Si-Si is not the same as the Si coordination, CN Si-O (i.e., the coordination number of Si-O), often referred in experiments 17,21,22 , which is very close to 4 in our simulations (CN Si-O ≤ 4.05 at P Q ≤ 8 GPa, see Supplementary Fig. 4) and consistent with other simulations 23 . It has also been shown 23,24 that a slower cooling rate such as the 1 K ps −1 used in current study tends to generate a CN Si-O = 4 network when the SiO 2 glass density is less than 2.7 g cm −3 (or P Q < 8 GPa). An example of CS is shown in the inset of Fig. 3b (left panel) wherein Fig. 1 Densities of the hot-compressed silica glass. The experimental data (red squares) is from Guerette et al. 17 . The compression is done at 1373 K at each corresponding pressure. The density of the MD simulated glass is on par with the pressure-quenched SiO 2 glass in experiments (the open blue circles and the filled blue circles correspond to the simulated glasses before and after the P Q is released, respectively). For MD simulation, quenching is done from 4000 K, and same as the experiment, the quench pressure P Q is released down to 1 atm at room temperature. The standard deviation of the density of MD glasses is around 0.007 g cm −3 . the central Si atom has N 1 = 4 and N 2 = 11. The spatial distribution of N 1-2 (i.e., N 1 + N 2 ) in the silica glass and its statistics are shown in Fig. 3a, b from which a significant increase in the dispersion of N 1-2 appears when P Q increases from 4 to 8 GPa. Likewise, the dispersion of the N 3-4 ( Fig. 3b, right panel) is also seen to increase more significantly at P Q > 4 GPa. The normalized standard deviation (Std (N k ) <N k > −1 ) in N k is shown in Fig. 3c, where a larger increase of variation in the Si coordination sequence is seen in N 1-2 (the 1st and 2nd neighboring Si atoms) than in N 3-4 and so forth in N 5-6 . The change in the 1st (locating at~3 Å away) and the 2nd (locating at~5 Å away) neighboring Si atoms is consistent with the major variation in void size as shown in Fig. 2 particularly in those voids whose diameter is between 2 and 4 Å (note the atomic radius of Si (0.26 Å) and O (1.35 Å) are taken into account). The large variation in N k in Fig. 3b, especially in N 1-2 also suggests that the homogeneity of the silica glass would deteriorate when P Q is too high and different CN Si-Si appears.
Calculation of Rayleigh scattering of densified silica glass As a result of decreasing void size, Rayleigh scattering is expected to be significantly reduced with improvement of the glass homogeneity. The Rayleigh scattering can be quantified via either Mie theory or density fluctuations (see "Methods"). The Rayleigh scattering intensity R sp of all voids can be calculated using the following equation based on the Mie theory 25 which sums contribution from all N voids in the material by treating them as scattering spherical particles, where V is the volume of the material, n is the relative refractive index, D i is the diameter of the ith void and λ is the wavelength the incoming light. R sp calculated from all voids (see "Methods" and Fig. 2b) is plotted against P Q in Fig. 4a by the blue curve. The plot shows a clear suppression of R sp (e.g.,~57% from 0 to 4 GPa or~65% when 8 GPa ≤ P Q ≤ 12 GPa). While it is difficult to directly obtain the Rayleigh scattering coefficient α scat from R sp , we have attempted to convert R sp to α scat at 1.55 µm via two SiO 2 samples from MD whose largest pore sizes are similar as two SiO 2 samples from experiments 15 whose α scat values are known (see "Methods"). Such obtained α scat value for each MD SiO 2 sample is shown using the right axis of Fig. 4a. The inset of Fig. 4a compares the MD results with the available experimental α scat available up to P Q = 0.2 GPa 15 where a suppression in α scat is seen in both cases. The amount of change in α scat from MD seems to be smaller than the experiments at 0.2 GPa, which may result from several factors: (1) the inaccuracy of the conversion from R sp to α scat (only the largest pore is considered), (2) the simple treatment of pore as spherical particles in the Mie theory while their shapes are irregular as shown in Fig. 2a, and/or (3) the much shorter time scale used in MD under high pressure than experiments. Alternatively, we may calculate α scat from the density fluctuations of the silica glass. As noted in "Methods", α scat can be written in terms of density fluctuations hN 2 i À hN 2 i hNi as 10,26 : where λ is the wavelength, n is refractive index, p is the photoelastic coefficient, and n is the number density of particles. N is the number of particles in volume V and the angle brackets indicate an average. The density fluctuation hN 2 i À hNi 2 can be calculated from density/volume variation of multiple samples melt-quenched using, e.g., the NPT ensemble. Given the empirical relation between n and ρ (the density of a glass in g cm −3 ) 17,27 , i.e., n = (ρ − 2.21) × 0.193 + 1.46, the photoelastic constant p for an isotropic solid can be calculated as 2ρ n 3 ∂n ∂ρ 10,28 which, when substituted in Eq. 2, would reduce the dependence on n from the eighth power to the second power. Use of this empirical relation is justified because the relative error of n is <2% based on ref. 27 which would at most cause a relative error of~4% in α scat . α scat as a function of quench pressure using Eq. (2) is shown in Fig.  4b. The Rayleigh scattering is found to be suppressed at P Q = 4 GPa with a reduction of 56% compared with at 0 GPa which is very close to the reduction in R sp (57%) from Fig. 4a. Specifically, a scattering loss~0.081 dB km −1 due to Rayleigh scattering can be achieved at 4 GPa, which is in the same ballpark as the Rayleigh scattering loss of conventional/hot-compressed silica in experiments 11,15 and is extremely low given the much higher cooling rate used in MD (1 K ps −1 ). However, a further increase of P Q to 12 GPa will instead enhance the scattering by more than 10 times, in contrast with Fig. 4a, suggesting that factors other than voids should be considered at higher pressure. Figure 4b also shows the change of the n 8 p 2 term in Eq. (2) with a monotonic increasing trend with P Q (e.g., an increase of 37% is seen from 0 to 4 GPa), confirming that the lowering of α scat comes solely from change of glass homogeneity. We note here that even though it seems that predictions using Mie theory gives smaller α scat than those from density fluctuation, both Fig. 4a, b show that a reduction of Rayleigh scattering by~56% can be achieved using a moderate quench pressure~4 GPa. The calculated value of α scat using voidscattering and density fluctuation models are, respectively, 0.069 and 0.21 dB km −1 under 0.2 GPa. Although the former value matches well with the experimental results (~0.06-0.07 dB km −1 ) as shown by the inset of Fig. 4a, we believe the latter value is more accurate when considering the high cooling rate and the short P Q holding time used in the MD simulations. The fictive temperatures (>2100 K) of the SiO 2 glasses in MD simulations are several hundred degrees higher than the experimental samples, suggesting that an extremely low Rayleigh scattering coefficient may be achieved by tailoring the thermal history of the glass in combination with its pressure history as reported by ref. 15 . The results of α scat calculation using void-scattering and density fluctuation both clearly indicate that, by increasing P Q from 0.0 to 4 GPa, a continuous decrease of α scat is expected. Such improvement in α scat should come from amelioration of homogeneity of silica glass as evidenced in Fig. 2.
Understanding the homogenization of silica glass To understand the homogenization mechanism of the silica network as a result of P Q , we examine the structural differences among those glasses. We evaluate the change of ring size distribution (whose calculation method is described in "Methods" and see Supplementary Fig. 3) for the simulated silica glass. For P Q < 6 GPa, CN Si-O in the silica glass remains close to 4.0 and the Si-Si distance scarcely changes (see Fig. 5c, Supplementary Fig. 4 and Supplementary Methods). The improvement in SiO 2 homogeneity can be understood using the topological pruning mechanism (Fig. 5a) 29,30 which, to our knowledge, has only been used to account for the densification of SiO 2 network. According to refs. 29,30 , ring formation reduces the density of the network by pruning a calculable number of tetrahedrally coordinated atoms Fig. 5 Homogenization mechanism of silica glass. a Configurational change to account for the homogenization effect (via elimination of a three-membered ring 29,30 as denoted in the red circle) and dehomogenization effect (via introduction of different CN Si-Si ) at high P Q in the SiO 2 glass. The yellow area denotes a void whose size decreases after the exemplary three-membered ring is eliminated. All spheres denote the Si atoms (the O atoms are not considered). b Proximity of free space to 3-memberd rings (Si: blue atoms and O: red atoms) in silica glasses whose P Q = 0 and 4 GPa. The yellow regions are isosurfaces with a minimum distance of 0.11 nm from the surface of any atoms in the glass. c Characteristic ring size K * (red) and CN Si-Si or CN Si-O at different densities. CN Si-Si is defined as the average number of first Si neighbor atoms to a central one using a cutoff of 3.4 Å. The dashed red line and black dashed line are guides for eyes. Multiple samples were used to calculate the standard deviation for the error bar. Fig. 4 Rayleigh scattering of pressure-quenched SiO 2 glasses at a wavelength of 1.55 µm. a Calculated R sp (left axis) at different P Q using the Mie theory and the corresponding α scat at λ = 1.55 μm (right axis). Experimental results from ref. 15 are also included for comparison. Note that the experimental data in Fig. 3a contain SiO 2 glasses quenched using the same pressure but different holding times or different glass transition temperature; therefore, multiple values exist at P Q = 0.0 and 0.2 GPa. The error bars are estimated by the difference between the largest R sp and the smallest R sp at each P Q . b Calculated α scat and n 8 p 2 at different quench pressures using density fluctuations. Each α scat is calculated from density variation of eighteen samples (see "Methods"). Additional data from SiO 2 samples with different P Q are added to obtain the trend for α scat . The blue dashed line is an eye guide to show the trend of α scat . The black dashed line corresponds to 0.1 dB km −1 .
(Si atoms) and the relative ability of different sized rings to reduce the network density decreases exponentially with the ring size K. Therefore, a higher density of SiO 2 can be achieved by eliminating smaller rings as illustrated in Fig. 5a: (i)-(ii) and increasing the characteristic ring size K * (see Supplementary Methods for calculation of K * ). At P Q ≤ 4 GPa, small ring ratios (K ≤ 7) are found mostly to decrease with P Q (see Supplementary Fig. 3B, D) and K* approaches 6 as shown in Fig. 5a. Note a linear increase of K * in the simulated glass with pressure at P Q ≤ 4 GPa is observed, which means that the glass network should become more homogeneous in the sense that its K * is closer to that of the quartz (whose K * is between 6 and 8) whose α scat is extremely small (<0.01 dB km −1 using k B Tβ Τ 31 where T = 300 K and β Τ = 2.68 × 10 −11 Pa 32 ). In order to further correlate such change to the void result (Fig. 2), we propose an alternative way to account for the homogeneity improvement via the proximity of free space, i.e., voids, to smaller rings (K ≤ 7) which, partially, can be pruned to densify the glass. An example is given in Fig. 5a: (i)-(ii) in which the free space (yellow region) near the 3-member ring shrinks after the ring is pruned to allow additional Si atoms (light blue) to connect with the group of dark blue Si atom. Similar mechanism can also occur at the proximity of 4 and 5 membered rings. The proximity of the free space to 3-membered rings is evidenced in Fig. 5b (also for K = 4 in Supplementary Fig. 5) at 0 GPa. Unlike the shrinkage of pores due to distortion of the rings under pressure which tends to recover once P Q is released, this achievement in homogenization requires breaking and forming of bonds which shall persist after P Q is released at the room temperature.
On the other hand, large rings (K > 7) increase, and particularly, small rings (K < 6) are found to significantly rebound at P Q > 4 GPa (see Supplementary Fig. 3b). This change in the constituent ratio of the ring size comes from the increase of CN Si-Si in the silica glass as shown in Fig. 5a: (ii)-(iii), leading to an increase in density fluctuation as shown by the different number of Si atoms enclosed by the red and green circles in Fig. 5a (iii) and the local Si neighbor number shown in Fig. 3. We note that the increase in the constituent ratio of small rings is also seen before P Q is released ( Supplementary Fig. 3c), which is consistent with experimental results under high pressure by Hemley et al. 33 and ab initio calculations by Trave et al. 34 and Karki et al. 35 of the silica liquid.
Although there is no consensus on the change of the populations of the 3-and 4-membered rings with pressure in silica glass/liquid from previous MD simulation results 17,29,[34][35][36] probably due to the usage of different simulation settings, e.g., the interatomic potential, the cooling rate, the pressure condition, the ring concepts, etc., our ring statistics seem to be supported experimentally by the change of the two defect peaks, D 1 and D 2 , arising from symmetric stretch of regular 4-membered rings and 3-membered rings 37,38 . Figure 6a shows our Raman measurements from the samples used in ref. 15 . Both D 1 and D 2 peaks are found to decrease in intensity without any frequency shift at higher pressure and at longer holding time, when the silica glasses are quenched from 2073 K with one at ambient pressure and two others held at 200 MPa for 1 and 4 h prior to the quench, suggesting a smaller population of the two types of the rings due to hot compression. The weakening (between 0 and 4 GPa) of the two peaks is also observed by the Raman spectrum results of SiO 2 glass by Guerette et al. 17 (Fig. 6b) quenched from a lower temperature 1373 K. Instead, at large pressures, e.g., 8 GPa, the intensities of the two peaks seem to be enhanced according to Fig. 6b using the right-most point as the reference. By comparison, the intensity of the D 2 peak of silica glass subjected to cold compression increases monotonically with the maximal pressure reached 39 . Here, we have several remarks. First, although experiments show that the CN Si-O of hot-compressed silica remains at close to 4 with P Q up to 8 GPa 17,21 , they do not preclude CN Si-Si from being higher than 4 as shown in Fig. 5c. Second, due to the much higher cooling rate used in MD, it is likely that the increase of CN Si-O and CN Si-Si may be postpone at some P Q value larger than 4 GPa if a slow cooling rate is used. However, heterogeneity shall occur as long as a coexistence of different CN Si-O /CN Si-Si starts to appear which will enhance optical scattering 16 . Lastly, based on the works by Bridgman 40-42 , Vukcevich 43 has shown a maximal compressibility is reached for cold-compressed silica glass under a pressure between 3 and 4 GPa. However, this is different from the observed lowest density fluctuation of the hot-compressed silica glass in our simulations which corresponds to a minimal compressibility at around 4 GPa (see Supplementary Fig. 1d).

DISCUSSION
The change of density fluctuations in SiO 2 glass with P Q can now be clearly understood in terms of the two topological modification of SiO 2 network. At P Q < 4 GPa, decrease of population of smaller rings and increase of larger ring population drive the network to be closer to a perfect Bethe lattice (Fig. 5a) whose density fluctuations is extremely low (all lattice points have equal chemical sequence). Such tendency of SiO 2 structure to approach the Bethe lattice can lead to a decrease in void size (Fig. 5a). At higher P Q , the equality of the Bethe lattice points is undermined by appearance of Si atoms with different CN Si-Si , leading to an increased density fluctuation. Thus, the highest homogeneity of silica glass is achieved at some intermediate value of P Q (~4 GPa) in our MD simulations. Lastly, the high cooling rate used in MD may shift the critical pressure to a smaller value above which CN Si-O and CN Si-Si start to increase rapidly. This suggests that further improvement in Rayleigh loss can be achieved by using a slower cooling rate at a higher P Q .
Previous MD results by Kaiki et al. 35 show that the Grüneisen parameter and heat capacity of the liquid SiO 2 have distinct values in this low-pressure regime (P ≤ 4 GPa). While a complete analysis of these properties is beyond the scope of this study, the SiO 2 structure between 0 and 4 GPa changes solely in the ring structure without introduction of defects, e.g., under-or over-coordinated Si/O atoms, with adverse effect on optical properties. Our study, therefore, shows a topological way of improving the optical properties of silica glass which may be applicable in other network glasses such as GeO 2 and B 2 O 3 . Lastly, while it is tempting to connect the characteristic ring size K * to the void size in SiO 2 , according to our results (Figs. 2b, 5b), they vary with P Q in opposite directions at P Q < 4 GPa, which suggests that the ring pruning should play a critical role in the homogenization.
To sum up, we have shown that, using the pressure-quench method and two different methods to calculate Rayleigh loss, the homogeneity of silica glass can be significantly improved, achieving a reduction in Rayleigh scattering by more than 50%. Given the different simulation conditions compared to experiments, the improvement in loss is expected to be even greater if lower cooling rates are used. Such improvement is due to the homogenization in the pressurized silica network. Specifically, at P Q ≤ 4 GPa, the reduction in density fluctuation results from a change in intermediate-range order that is supported by Raman spectrum and ring structure analysis in our MD simulations and can be understood as an approaching of the network's characteristic ring size to a perfect Bethe lattice. Higher values of P Q , however, are found to potentially bring about coexistence of different network coordinations which will undermine the homogeneity of the network. Our findings have shown, through an atomistic understanding, how the homogeneity of SiO 2 is improved by pressure; topological pruning from decrease of small rings can drive the network to be more like a Bethe lattice and at the same time decrease those voids close to the pruned rings without introducing unwanted defects. This provides not only the detailed understanding of the structural change induced by hot compression, but also an example of ideal SiO 2 glass structure of high homogeneity which is desirable for optical applications. Such ideal glass structure is useful not only for silica glass, but also for other types of inorganic glass with similar network structure such as GeO 2 . Although a high P Q at high temperature is still technically challenging, the ideal glass structure might be approachable using other techniques such as high-pressure chemical vapor deposition 44,45 , which might be applied to improve the optical transparency of light waveguide.

Molecular dynamics simulations
Silica glass samples were generated by randomly distributing 3000 atoms in a cubic box and melted in the liquid state at 4000 K at different pressures for 100 ps, followed by a quenching process to 300 K at a rate of 1 K ps −1 . It has also been shown 23,24 that a slower cooling rate such as the 1 K ps −1 used in current study tends to generate a CN Si-O = 4 SiO 2 network when the SiO 2 glass density is less than 2.7 g cm −3 (P Q < 8 GPa according to Fig. 1). The quench pressure P Q = 0-12 GPa is maintained using a Nosé-Hoover thermostat/barostat [46][47][48] in the isothermal-isobaric ensemble (NPT). Larger samples of 24,000 atoms were used to get the cooling curves (e.g., Fig. 1) and measure the density fluctuation. After the meltquench, the pressure is reduced to 1 atm in 100 ps followed by an additional run for 100 ps at 1 atm. The SHIK potential 18 used for the atomic interaction has been found to accurately reproduce radial distribution function (RDF) of silica liquid 18 , vibrational density of states 18 , and the density of glass at 300 K at various pressures. The Si speciation in liquid SiO 2 predicted by this potential also compares well with ab initio calculation and TS potential considering O polarizability (see Supplementary Fig. 1c). For comparison, BKS potential 49,50 broadly applied for SiO 2 simulations is also used to generate silica glass by pressure quench of the liquid from 7000 to 300 K using the same cooling rate (see Supplementary Fig. 1b). Note a higher cooling rate is used in ref. 24 , making the resulted glass densities different from us. All MD simulations and structure visualizations are carried out using LAMMPS 51 and OVITO 52 , respectively.

Structural analysis
A cutoff of 2.2 Å is used for Si-O pair, and 3.4 Å for Si-Si pair for bond angle and coordination analysis. These values are the position of the first minimum in the corresponding RDFs. The void and its distribution was obtained using Zeo++ 53,54 . Shannon radii 55 for Si, 0.26 Å and O, 1.35 Å are used to define the surface of an atom and for the void analysis, which are very close to the radii determined in ref. 22 . The void radius distribution from Zeo++ does not give the number of voids. To calculate the total amount of voids, we use the surface area function in Zeo++ to obtain the number of voids whose radii are larger than a critical value so that they comprise of the largest~20% voids using Fig. 2b. The total number of voids is then calculated by dividing this number with this percent. The R.I.N.G.S. code 56 was used to determine the primitive ring distribution.

Raman spectroscopy
Raman spectra were collected with a Raman microscope (Nicolet Almega, Thermo Fisher Scientific Inc., Yokohama, Japan) in backscattering geometry. The excitation wavelength was 532 nm. Polished plates with dimensions of 15 × 15 × 1 mm were used for the measurements. 100× microscope and 1800 l mm −1 grating were used. Data were taken with the exposure time set as 60 s, with accumulation of 10 times. All measurements were conducted at room temperature (around 300 K).

Calculation of Rayleigh scattering (R sp and α scat )
The R sp can be directly calculated using Eq. (1) and the total number of voids and the void distribution obtained from Zeo++ 53,54 . Two methods were used to calculate the Rayleigh attenuation α scat . In method I (see Eq. (3)), the following equation is assumed to derive α scat from R sp for the MD simulated silica glasses, α scat ¼ A log 10 where A, B, and C are constants and the second equation is valid when the term inside the logarithmic function is close to 1. Two SiO 2 samples pressure-quenched at 0 and 0.2 GPa in the current MD simulation whose largest pore sizes (0.2435 and 0.2415 nm) and densities (2.210 and 2.215 g cm −3 ) are similar to two of the experimental samples in ref. 15 (0.2436 nm and 0.2414 nm, and 2.230 g cm −3 and 2.256 g cm −3 ) whose α scat are known. C are found to be 4.14 × 10 −3 dB km −1 µm −4 . In method II, we directly calculate the density fluctuation (both configurational and vibrational ones are included) which causes the Rayleigh scattering of pure SiO 2 glass. The density fluctuation is described by 57,58 , where N is the number of particles in volume V, n is the number density, k B is the Boltzmann's constant, T is the temperature, β Τ is the isothermal compressibility, and angle brackets indicate an average. μVT and NPT represent the grand canonical ensemble and the isothermal-isobaric ensemble, respectively. Because all MD simulations were carried out using the NPT ensemble, we have used volume variation of multiple glass samples instead of variation of the number of particles to calculate the density fluctuation. Data reported are from 18 samples prepared under similar conditions but different initial configurations.

DATA AVAILABILITY
Details of the simulations are available within the article. All data supporting the findings of this study are available from the authors upon reasonable request.