Percolation transitions in compressed SiO2 glasses

Amorphous–amorphous transformations under pressure are generally explained by changes in the local structure from low- to higher-fold coordinated polyhedra1–4. However, as the notion of scale invariance at the critical thresholds has not been addressed, it is still unclear whether these transformations behave similarly to true phase transitions in related crystals and liquids. Here we report ab initio-based calculations of compressed silica (SiO2) glasses, showing that the structural changes from low- to high-density amorphous structures occur through a sequence of percolation transitions. When the pressure is increased to 82 GPa, a series of long-range (‘infinite’) percolating clusters composed of corner- or edge-shared tetrahedra, pentahedra and eventually octahedra emerge at critical pressures and replace the previous ‘phase’ of lower-fold coordinated polyhedra and lower connectivity. This mechanism provides a natural explanation for the well-known mechanical anomaly around 3 GPa, as well as the structural irreversibility beyond 10 GPa, among other features. Some of the amorphous structures that have been discovered mimic those of coesite IV and V crystals reported recently5,6, highlighting the major role of SiO5 pentahedron-based polyamorphs in the densification process of vitreous silica. Our results demonstrate that percolation theory provides a robust framework to understand the nature and pathway of amorphous–amorphous transformations and open a new avenue to predict unravelled amorphous solid states and related liquid phases7,8. Amorphous–amorphous phase transitions in silicon dioxide are shown to proceed through a sequence of percolation transitions, a process that has relevance to a range of important liquid and glassy systems.

Understanding the physical mechanisms controlling the transformation from one amorphous 'phase' to another is an open fundamental issue in materials science 1,2, 8,9 . The short-range structures of amorphous solids such as SiO 2 , GeO 2 , Si, Ge and the chalcogenides are very similar to those of their crystalline counterparts and it is the random nature of the inter-unit connections that produces disorder at long length scales. In such systems, it has been argued that amorphous states can be described in terms of changes in the electronic and structural properties at short-and medium-range order. These properties include electronic bonding 2,4,10 , the coordination number 1,3,4,11 and the ring distribution 11,12 . However, to describe the passage from one amorphous 'phase' to another at the critical threshold, an explicit scale-invariant quantity (that is, an order parameter) needs to be defined. This concept has not been considered for pressure-driven amorphous solid transformations thus far.
When an amorphous solid is pressurized, it is commonly assumed that transformations from low to high density occur gradually, with coexisting low-and high-fold coordinated polyhedra. This structural evolution, sometimes referred to as polyamorphism, differs from the polymorphism of crystals, in which transitions occur from a specific phase to another at a critical pressure. The change from tetrahedral (Si-O coordination number Z = 4) to octahedral (stishovite-like structure, Z = 6) structures in vitreous silica (v-SiO 2 ) is accompanied by a well-known mechanical anomaly at 3 GPa (refs. 13,14 ), a percolation phenomenon at about 7 GPa (ref. 15 ) and a two-step elastic-to-plastic transformation, one around 10 GPa (refs. 12,14,16 ) and a second around 20 GPa (ref. 14 ).
A recent ab initio machine learning method applied in combination with empirical force fields in compressed amorphous silicon highlights a three-step transformation sequence for amorphous silicon under increasing pressure 4 . However, the passage from one 'phase' to another was not addressed. For v-SiO 2 at ambient pressure, ab initio methods complement experimental data 17,18 , but, despite many efforts 12,19,20 , application to the above-mentioned issues has been limited due to the prohibitive calculation times required for a reliable thermodynamic sampling 21,22 .
We were thus motivated to use an ab initio-based approach, the self-consistent-charge density functional-based tight-binding (SCC-DFTB) method 23 . Combined with the molecular dynamics (MD) technique, it is able to reproduce many features of crystalline silica (c-SiO 2 ) and v-SiO 2 with an accuracy similar to those of ab initio methods, but at least three orders of magnitude faster. This allowed us to simulate more than forty SiO 2 glasses and explore the full pressure range up to 82 GPa at room temperature. We also decompressed the samples at 8 GPa and 10 GPa (for details, see Methods). Our samples contain 1,008 atoms (336 SiO 2 units) and are large enough to observe and describe the percolation critical phenomena governing the aforementioned transformations. Figure 1a illustrates the v-SiO 2 transformation from Z = 4 to 6 with increasing pressure p. At ambient temperature and 0 < p < 3 GPa, similar to c-SiO 2 , the Si-O bond of v-SiO 2 is characterized by an sp 3 hybridization that favours the formation of SiO 4 tetrahedra. For 3 < p < 10 GPa, the densification of v-SiO 2 includes changes in the electronic structure, with an increase in the Fermi energy and of the atomic charges (Extended Data Figs. 1 and 2). These changes produce an increase in the Si-O bonding ionicity (Extended Data Fig. 1c), boosting the formation of higher-fold coordinated polyhedra (Fig. 1a, b). As a consequence, beyond 10 GPa, the inter-polyhedra connectivity changes from purely corner-sharing (CS) tetrahedra to a more complex connected network (Extended Data Figs. 3 and 4). The latter involves corner-sharing and edge-sharing (ES) SiO 5 -SiO 6 and SiO 6 -SiO 6 ( Fig. 1a, c) and a lower extent of face-sharing (FS) polyhedra (Extended Data Fig. 5).

Polyhedra coordination and connectivity
On increasing the pressure from ambient conditions, c-SiO 2 changes from a tetrahedral local structure to an octahedral one through structural transformations including α-quartz to coesite I (Z = 4) at 3 GPa, then coesite I to stishovite (rutile structure, Z = 6) at 9 GPa, with a higher ionic bonding contribution 10 . Regarding the local environment of the O atoms, CS tetrahedra yield to divalent OSi 2 structures, and ES octahedra to trivalent OSi 3 structures. In the crystalline phase at ambient pressure, each CS tetrahedron is connected to four neighbouring tetrahedra, leading to Z SiSi = 4 and Z OO = 6. In stishovite, each SiO 6 is connected to ten octahedra (Z SiSi = 10), two of which are ES and eight CS, and Z OO = 12. In v-SiO 2 , the coordination numbers exhibit similar trends with pressure (Fig. 1b), as well as jumps that allow us to identify different regimes at intermediate pressures (Fig. 1b, c, top arrows).
For 13 < p < 20 GPa, the compressed glasses contain similar fractions of SiO 4 , SiO 5 and SiO 6 , with a total number of shared edges per polyhedron of n ES ≈ 1.4 and a small contribution of FS polyhedra. This is similar to the crystalline polymorph discovered recently, the so-called coesite IV 6 , which is rarely observed for high-valence and low-coordinated cations in crystals, as it is at odds with the fifth Pauling rule 24 . For 20 < p < 45 GPa, the fraction of SiO 6 octahedra increases with a jump, while the fraction of SiO 4 tetrahedra decreases to almost zero. In parallel, n ES increases to about 1.7 and a small fraction of FS polyhedra persist, which looks similar to coesite V, the second recently discovered c-SiO 2 polymorph 6 . For 45 < p < 70 GPa, the network connectivity is fully dominated by SiO 6 octahedra, either CS or ES, and n ES ≃ 2, which reveals a stishovite-like 'phase'. For p > 70 GPa, n ES > 2 due to the emergence of SiO 7 polyhedra, suggesting the onset of a post-stishovite-like 'phase' 20,25 .
Finally, we calculated the fraction of Si polyhedra and the coordination number for unloadings from 8 GPa and slightly above 10 GPa. In the former case, the decompression is reversible, but it is not in the latter (Fig. 1b, small symbols). The pressure threshold at which irreversibility occurs agrees with experimental observations 14 . When unloading from 10 GPa to ambient pressure, a hysteresis is clearly observed due to persisting residual SiO 5 pentahedra (hence Z ≠ 4, Fig. 1b), confirming previous expectations and pointing out the key role of five-fold coordinated polyhedra in the plastic-to-elastic transition of v-SiO 2 (refs. 12,14,16 ).

Percolation transitions
The above observations only include an evaluation of the local and medium-range order. However, the nature of the structural transformations from one 'phase' to another in compressed v-SiO 2 should also be substantiated by the structure at a long length scale. Indeed, we found that clusters built from SiO n polyhedra (including their mixtures, (SiO n -SiO m ) eventually emerge. One of these becomes dominant and percolates by spanning its structure from one side of the simulation box to the other (Fig. 2a). This effect has important implications for the glass rigidity 15,26 and can be monitored by calculating the percolation probability, P ∞ . The latter is often used as an order parameter for the analysis of a percolation transition. It tends to 1 (or 0) if the largest cluster percolates (or not), and the percolating cluster is representative of the new 'phase' 27 . The notation (SiO n -SiO m ) ∞ used in the following corresponds to a percolating cluster composed of alternating n-and m-fold coordinated polyhedra.  Figure 2b shows the percolation probability, P ∞ , of the largest cluster for connected networks of SiO 4 , SiO 5 and SiO 6 polyhedra, and their mixtures, as a function of pressure in v-SiO 2 (density-dependent P ∞ is shown in Extended Data Fig. 6a). For the (SiO 4 -SiO 4 ) ∞ cluster, P ∞ decreases from 1 to 0 at 13 GPa. In addition, between 3 and 8 GPa, a percolating cluster emerges that is composed of a connected skeleton of alternating tetrahedra and pentahedra, (SiO 4 -SiO 5 ) ∞ . This coincides with the percolation transition observed at 7 GPa (ref. 15 ). With increasing pressure, two more percolating clusters appear at 10 GPa, (SiO 5 -SiO 5 ) ∞ and (SiO 5 -SiO 6 ) ∞ , coexisting with the first one. The percolating cluster structures containing a mixture of SiO 4 and SiO 5 between 8 and 10 GPa, as well SiO 6 between 10 and 13 GPa, recall the pressure-induced post-quartz amorphous states 28 . At 13 GPa, a (SiO 6 -SiO 6 ) ∞ cluster appears and percolates. The coexistence of all these percolating clusters beyond 13 GPa is similar to coesite IV, the crystalline structure of which combines SiO 4 -SiO 5 , SiO 5 -SiO 5 , SiO 5 -SiO 6 and SiO 6 -SiO 6 planes in different crystallographic directions. By analogy, the vitreous state for 13 GPa < p < 20 GPa is labelled v-coesite IV in Fig. 2b.

Article
Slightly above 20 GPa, SiO 4 tetrahedra vanish and, accordingly, the P ∞ of the (SiO 4 -SiO 5 ) ∞ cluster decreases to 0. Simultaneously, a percolating cluster, (SiO 6 -SiO 5,6 ) ∞ , emerges, which is composed of ES SiO 6 -SiO 5 and SiO 6 -SiO 6 polyhedra, together with some FS contribution. In this pressure range, the amorphous state (labelled v-coesite V in Fig. 2b) possesses the same dominant ES polyhedra connectivity as well as alternations of SiO n -SiO m polyhedra equivalent to those in coesite V 6 . Our result correlates with the irreversible structural behavior observed experimentally above about 20 GPa in v-SiO 2 (ref. 14 ), suggesting that ES structures remain stable during decompression.
Analysis of the percolation probability of a cluster composed of pure ES octahedra, similar to stishovite, provides evidence that such a cluster percolates around 40 GPa and, accordingly, a v-stishovite state replaces the v-coesite V one (Fig. 2a, b). The pressure at which this occurs is in good agreement with the value reported for compressed v-SiO 2 when Z becomes equal to 6 (ref. 25 ). Therefore, our results demonstrate that, instead of a single and gradual transition, the mechanism of the structural transformation from Z = 4 to 6 includes a series of percolation transitions between well-defined amorphous states (Fig. 2). Moreover, estimation of P ∞ as a function of the SiO n fraction shows that the critical fractions lie in the range expected from percolation theory, that is, around 2/Z (Methods and Extended Data Fig. 6b). Finally, an analysis from the viewpoint of O atoms also reveals that OSi 2 -OSi 3 and OSi 3 -OSi 3 clusters percolate around 3 GPa and 10 GPa, namely, at pressures similar to those for the SiO n -SiO m structures (Methods, Fig. 2b and Extended Data Fig. 7c).
The correlation length ξ has been estimated in the pressure range at which the different percolating clusters appear. This quantity accounts for the maximum length at which scale invariance exists. At the thermodynamic limit, this quantity is expected to diverge at  the critical pressure and it should reach a maximum for finite size systems, which is indeed the case for the different situations described above (Fig. 2c). Beyond 70 GPa, a large cluster of ES (SiO 6 -SiO 7 ) ∞ polyhedra starts to grow; this, together with the fact that at this pressure the total number of shared edges per Si polyhedron exceeds 2, supports the post-stishovite 'phase' suggested previously. The jump in the percolation probability P ∞ from 0 to 1 provides evidence of the dominant character of a well-defined amorphous state in a given pressure range. However, the non-negligible partial overlap between the ξ curves beyond 3 GPa may indicates a 'phase' coexistence between the dominant percolating cluster and a 'metastable' emerging state (Fig. 2c) when approaching the critical pressures. This is similar to the metastability observed in coesite IV and V polymorphs 5,6 . When unloading from 8 GPa, the percolating (SiO 4 -SiO 5 ) ∞ cluster vanishes at a pressure similar to that at which it percolates during compression (Fig. 2b). This confirms the reversible character suggested above by analysis of the local and medium-range structure. Conversely, when the sample is decompressed from 10 GPa (dashed curve), P ∞ exhibits a hysteresis, that is, the (SiO 4 -SiO 5 ) ∞ cluster continues to percolate at pressures much lower than for the percolating compression route (Fig. 2b). A similar plastic behaviour is observed for the other percolating clusters, that is, (SiO 5 -SiO 5 ) ∞ and (SiO 5 -SiO 6 ) ∞ (data not shown). This demonstrates that the formation of percolating clusters containing SiO 5 pentahedra for p ≳ 10 GPa prevents recovery of the pristine glass structure and probably explains the irreversible behaviour.
To explore the impact of these percolation transitions on glass rigidity, we calculated the bulk modulus K of our pressurized silicas (Fig. 2d). The increasing connectivity at short and long length scales yields a strong increase in K, that is, from approximately 30 GPa at ambient conditions up to approximately 375 GPa at p = 82 GPa. Its evolution parallels those of its crystalline counterparts: SiO 4 polymorphs (except coesite I) below 10 GPa, coesite IV and V at intermediate pressures and stishovite above [35][36][37][38][39][40]29 ). Our calculations also reproduce the well-known minimum around 3 GPa (Fig. 2d, inset) and we associate the increased rigidity starting above 4 GPa to the series of percolating clusters that emerge with increasing pressure. Finally, our estimation of K is very close to high-frequency experimental data except perhaps around 10 GPa, where the simulations exhibit a minimum, when (SiO 5 -SiO 5 ) ∞ and (SiO 5 -SiO 6 ) ∞ clusters emerge. These results suggest that the mechanical properties relate primarily to the peculiar connectivity describing the system at all length scales, rather than to the crystalline or amorphous nature of the network.
We have also analysed the structure of the percolating clusters at the critical pressure thresholds. For this, we expanded these clusters out of the box, taking into account the periodic boundary conditions, and computed their structure factor S(q). As expected from percolation theory, at the critical thresholds the clusters are fractals, that is, ment with the value predicted theoretically (Fig. 3).

Conclusions and outlook
The structural transformation from low-to high-density amorphous SiO 2 solid occurs through a sequence of percolation transitions. Each emerging amorphous phase (polyamorph) actually corresponds to a state of a megabasin of the configurational energy landscape 30 , and each megabasin is associated with a specific crystalline polymorph with the same local structures and connectivities as those of the percolating cluster. When the pressure is increased and the energy barrier separating two megabasins is overcome, a transformation is induced to an amorphous state with higher connectivity. Other pressure-or temperature-driven amorphous transformations in oxide glasses, chalcogenides and metallic glasses could be explained in these terms. A percolation phenomenon has been proposed to explain the transition between low-and high-density liquid phases of water 31 . The same has been argued for the glass transition in model systems [32][33][34] . Changes in the local connectivity similar to those observed here have been demonstrated in pressurized SiO 2 , GeO 2 and H 2 O liquids, but for different pressure windows and interpolyhedra distances [35][36][37][38] . The question whether amorphous states are related to underlying and unidentified liquid phases can now be addressed by considering a scale-invariant quantity as the percolation probability. This will help to catalogue amorphous solids and unveil their affinities with liquid and crystalline phases.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-03918-0.

Glass preparation
The initial glass configuration of 1,008 atoms (336 SiO 2 units) in a cubic simulation box was prepared by carrying out classical MD simulations using periodic boundary conditions. We used the melt-and-quench approach and the constant volume-constant temperature (NVT) ensemble, and the atomic interactions were modelled using the Van Beest-Kramer-VanSanten (BKS) pair potential 39 . We then carried out MD simulations within the framework of the SCC-DFTB method for the electronic structure calculations 23 . These calculations were also performed using the NVT ensemble. An Andersen thermostat was used to maintain room temperature and the samples were relaxed for more than 40 ps, with a time step of 2 fs. Within the SCC-DFTB approximation, the electronic integrals corresponding to the Hamiltonian matrix elements were replaced by analytical functions. The coefficients of these functions were obtained by fitting DFT calculations corresponding to systems obtained at different physico-chemical conditions. We used the Periodic Boundary Conditions (PBC)-0-3 set of parameters for SiO 2 (ref. 40 ). This resulted in an electronic structure calculation method at least three orders of magnitude faster than DFT, which is currently considered the best-performing ab initio method in materials science.
The orbitals were occupied according to a Fermi-Dirac distribution with a temperature of 300 K. We included Broyden charge mixing in the SCC cycle with a mixing parameter of 0.1. The Γ-point approximation was used for the electronic band structure energy and the projected electronic density of states (DOS) calculations. For the latter, a Gaussian broadening of 0.4 eV was considered.

Glass compression
At ambient pressure, the glass sample prepared using the above procedure had a density of 2.3 g cm −3 . To mimic quasi-static compression, the box length was reduced in steps of 1% (and atom positions rescaled) every 40 ps, until a pressure of approximately 27 GPa was reached. A second series of glasses was obtained by instantaneous volume reduction of the sample obtained at 15 GPa, thus yielding seven samples at target pressures of up to 82 GPa. Two successive volume reductions of 0.125% followed by a 40-ps relaxation were performed on six samples from this second series. The largest compression rate used here is almost equivalent to 0.02 GPa ps −1 . This corresponds to a compression velocity at least one order of magnitude slower than that used in a previous DFT study of compressed SiO 2 at T = 300 K (ref. 12 ), while the number of atoms in our calculation box is almost seven times larger. To distinguish between the elastic and plastic regime, we also decompressed the samples from 8 and 10 GPa. This was done by increasing the box length and scaling the atom positions by 1% every 40 ps, until ambient pressure was reached. The reported structural quantities were averaged over all configurations generated in the last 10 ps of the sample relaxation time.  Fig. 1a). The effect induced by sample densification in the total DOS of v-SiO 2 also agrees 10,19 (Extended Data Fig. 1c, d). Similarly, the Si and O Mulliken atomic charges and the Si-O bond ionicity change are consistent with previous experimental and DFT findings 10,12 (Extended Data Fig. 2). This favours the formation of pentahedra and higher-fold coordinated Si polyhedra 42 (Fig. 1). The maxima reached in the pressure range of 10-20 GPa where the SiO 5 fraction is dominant (Fig. 1) and the slight decreases beyond 20 GPa highlight the large Si-O bonding ionic character of SiO 5 pentahedra compared to octahedra and tetrahedra.

Structure factor
To compare the structural properties of our samples to experimental data, we computed the neutron and X-ray total static structure factors. We first computed the partial structure factors S αβ (q) using the definition 43 Here, α, β = Si, O. Here, f αβ = 1 for α = β and f αβ = 1/2 otherwise, N α is the number of particles of species α and N is the total number of atoms. The total structure factors are combinations of partial structure factors. For the neutron structure factor, we used the relation 43 with the neutron scattering length b α given by b Si = 4.1491 fm and b O = 5.803 fm, respectively 44 . The X-ray total structure factor S X (q) is given by 45 Here f α (s) is the scattering-factor function (also called the form factor), computed as a linear combination of five Gaussians using parameters derived elsewhere 46 .
The calculated X-ray S X (q) and neutron S N (q) structure factors are compared to experimental data 3,11 . The calculated S X (q) reproduces the decrease and shift of the first sharp diffraction peak (FSDP) due to the collapse of the open structure of v-SiO 2 at intermediate length scales (Extended Data Fig. 3).

Local connectivity
We computed the pair distribution function g αβ (r) for all compressed v-SiO 2 samples using the definition 43 where ⋅ represents the thermal average, V is the volume of the simulation box and δ αβ is the Kronecker delta function. For each of the three pairs, the first peak of the corresponding pair distribution function corresponds to the distribution of the first-neighbour shell distances for Si-O, O-O and Si-Si, respectively. For the Si-O pair, a minimum after the first peak located at 2.3 Å defines well the upper limit of this distribution at all pressures. We used this value as a cutoff distance to estimate the Si-O coordination number and the fraction of SiO n polyhedra for the considered pressures. Similarly, a cutoff of 3.5 Å was used to identify the first Si-Si neighbours and to estimate the evolution of the polyhedra coordination number with pressure (Fig. 1b). Identification of the Si-O and Si-Si first neighbours also allowed us to estimate the number of O atoms shared by two Si polyhedra neighbours, that is, one, two, three or more, corresponding to CS, ES and FS polyhedral connectivity, respectively ( Fig. 1c and Extended Data  Fig. 4b, c. The above-mentioned upshift of the FSDP translates into a fast reduction of the angles up to the first plastic transformation around 10 GPa when CS SiO 5 -and SiO 6 -based percolating clusters emerge. With increasing pressure above 10 GPa, and in addition to the primary peak at around 130°, we notice the occurrence of two other peaks pointing approximately at 98° and 68°, respectively, in the bond-angle distribution of Si-O-Si. The positions of the three peaks show a weak pressure dependence, and the two peaks located at smaller angles become prominent in the pressure range of v-coesite IV and V for the first one and v-stishovite for the second. The Si-O and Si-Si bond lengths exhibit anomalies at pressures of approximately 10 GPa and 20-25 GPa, corresponding to the two stages of plastic transformation, and the O-O bond length presents a maximum at a pressure of about 3 GPa, corresponding to the minimum of the bulk modulus. To provide more insights into the medium range of the silica network with increasing pressure, we computed the ring distributions. No changes occur for pressures below 4 GPa. For pressures below 8 GPa, we observe a mild increase in small ring sizes (2-, 3-and 4-membered ones) accompanied by a decrease in the proportion of large ring sizes (7-and 8-membered ones). When approaching the threshold of the plastic regime at around 10 GPa, the same pressure dependences are observed, but with stronger trends. This increase in small rings at the expense of large ones at high pressure confirms previous DFT calculations 12 .

Cluster analysis
To identify the different clusters with specified polyhedra coordination and connectivity, as well the largest cluster, we used a strategy inspired by the friends-of-friends algorithm, which is widely used to characterize dark-matter halos from N-body simulations 47 . In our implementation, the first nearest neighbours correspond to the specific Si-O coordination and polyhedra connectivity we want to analyse. To determine whether the largest cluster percolates or not, we first extract all clusters from the box, except the largest one. We then replicate the box containing this cluster along the three spatial directions. The friends-of-friends algorithm is applied again and, if the resulting largest cluster is larger than the box size in at least two directions, we assume that the cluster percolates. For a given pressure, the above procedure is implemented for each possible SiO n -SiO n+1 connectivity, that is, CS, ES and FS polyhedra.
The percolation probability P ∞ is defined as the number of times the resulting largest cluster percolates during the last 10 ps of the relaxation time, normalized by the total number of explored configurations. For different cluster connectivities, Fig. 2b and Extended Data Fig. 6a, b show P ∞ as a function of pressure, sample density and SiO n fraction, respectively. In Extended Data Fig. 6b we plot only the data for the connected polyhedra that have the same coordination, otherwise the fraction of mixed coordination cannot be properly defined. The critical fractions where CS (SiO n -SiO n ) ∞ (n = Z = 5 and 6) percolating clusters emerge and SiO 4 -SiO 4 ones disappear are consistent with the expected critical occupation probability in bond percolation (P c ≈ 2/Z). For example, in square lattices and diamond networks (Z = 4), P c is equal to 0.5 and approximately 0.39, respectively, and in triangular lattices and simple cubic networks (Z = 6) is equal to approximately 0.35 and 0.25, respectively. For ES (SiO 6 -SiO 6 ) ∞ , the critical fraction of SiO 6 is larger than for CS cases. This is because only a small number of SiO 6 octahedra are ES-connected. Note that ES connectivity is similar to a site percolation mechanism presenting a larger P c than the one resulting from bond percolation.
The where the sums run over pairs of Si atoms belonging to each cluster of size s. The correlation length ξ is averaged over the last 10 ps of the relaxation time.

Mechanical properties
We estimated the bulk modulus K of the compressed glasses by using the relation K = −Vdp/dV, where V is the sample volume and p is the pressure averaged over the last 10 ps of the relaxation time. To obtain K, we used the central difference scheme, where its error bar is calculated by error propagation and the pressure error is estimated from the standard deviation. To minimize errors in the finite difference, when two subsequent pressures lie within the estimated error bars, the second is excluded from the estimation of K. We recall that volumetric measurements give direct access to the static compressibility 13 , whereas high-frequency experiments measure the volumetric variations at the timescale of the relaxational processes probed by the instrument 14 . In compressed v-SiO 2 , the two values are similar in the elastic regime (p < 10 GPa), but strong variations have been observed in the plastic case for 10 < p < 55 GPa (refs. 14,48,49 ). Interestingly, above 12 GPa, our calculated K values are much closer to the high-frequency data than to the static ones, a behaviour that possibly reflects the typical time relaxation (approximately 10 ps, corresponding to gigahertz frequencies) and the small volumetric variation (<0.4 nm 3 ) used for the estimation of K in our MD simulations. For p > 55 GPa, the experimental value (static and high-frequency) is around 420 GPa (refs. 48,49 ), in good agreement with our calculations.

Affinities with c-SiO 2 polymorphs
To identify the SiO 2 polymorph most akin to amorphous states, we built supercells containing 582 atoms for tridymite, cristobalites, coesites and stishovite, and 432 and 486 atoms for β-and α-quartz, respectively. The unit cells of the polymorphs were obtained from crystallographic databases under AMCSD codes 0000789 for α-quartz, 0018071 for β-quartz, 0001629 for α-cristobalite, 0017646 for β-cristobalite, 0000531 for tridymite and 0001306 for stishovite 50 . These supercells were relaxed using the Limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) minimization method implemented in the DFTB+ package 23 . Only the lengths of the lattice vectors were allowed to vary during structure optimization (with the exception of coesites IV and V) to preserve the same sample densities reported in ref. 6 . The criterion used to define the v-coesite IV 'phase' relies on the coexistence of (SiO 4 -SiO 5 ) ∞ , (SiO 5 -SiO 5 ) ∞ , (SiO 5 -SiO 6 ) ∞ and (SiO 6 -SiO 6 ) ∞ percolating clusters. We found that almost all the SiO 5 -SiO 6 and SiO 6 -SiO 6 are ES, with a weak contribution of FS polyhedra, as in coesites IV and V. Indeed, in coesite IV we have validated the coexistence of different crystalline planes that have the above-mentioned connectivities.
The ES SiO 5 -SiO 6 and SiO 6 -SiO 6 polyhedra become dominant at a pressure of 20 GPa. They are in a similar proportion (around 40%) and control the network connectivity as in coesite V. Beyond 40 GPa, the fraction of SiO 5 starts to decrease and the SiO 6 -SiO 6 connectivity then becomes dominant, which allows us to identify the stishovite-like 'phase'.
To validate the above observations, the total DOS, Fermi energy, Mulliken atomic charges and bond ionicity were computed and compared to the results obtained for c-SiO 2 polymorphs (Extended Data Figs. 1 and 2). The electronic properties of v-SiO 2 appear to be closer to α-quartz than to other SiO 4 -based polymorphs at ambient pressure, as suggested elsewhere 10 . With increasing pressure, the electronic properties of the vitreous polyamorphs-v-coesite IV, v-coesite V and v-stishovite-are also compatible with those of coesite IV and V and stishovite crystals, respectively. The ionic character of the Si-O bond is also larger for the five-and six-fold Si coordinated polyhedra 51 . A similar conclusion holds when comparing the polyhedra coordination and connectivity with the corresponding v-SiO 2 percolating clusters, as well as for the calculated bulk moduli. The latter agree very well with those of SiO 4 -tetrahedra-based polymorphs at low pressures and to those of coesite IV and V and stishovite (Fig. 2d) when the pressure is increased beyond 13 GPa.
Finally, Extended Data Table 1 presents the densities of the crystalline polymorphs and the corresponding pressure ranges 6,28,52,53 . A comparison with the corresponding polyamorph states of v-SiO 2 is striking. One notices, for example, that a small pressure of around 3 GPa is enough to compact the open structure of vitreous silica up to a density very close to that of α-quartz. At high pressure, the density ranges of the polyamorphs follow those of their crystalline counterparts.

Percolating clusters at the critical threshold
We analysed the structures of the percolating clusters around the critical pressures. For this, we expanded and extracted the percolating clusters out of the box by taking into account the periodic boundary conditions. At the critical threshold, these clusters are usually larger than the box size because of their branched structure. We computed the structure factor S(q) of the isolated clusters by using the corresponding expression for a diluted system 54  pressed v-SiO 2 near the critical pressures, we have found that D f = 2.5, in agreement with the value predicted by percolation theory (Fig. 3).

Connectivity of OSi n structures
With increasing pressure, the evolution of the coordination number of the O atoms, Z′, parallels that of the Si atoms; that is, Z′ passes from 2 to 3 when Si octahedra (Z = 6) replace tetrahedra (Z = 4) (Extended Data Fig. 7a, b). The anomalous steep jumps defining the different amorphous states, as described in the main text, are also reproduced. From a strictly structural point of view, these observations suggest that the two approaches are similar. The abrupt increase in OSi 3 structures around 10 GPa leads to the percolation of (OSi 3 -OSi 3 ) ∞ clusters, in line with the percolation of (SiO n -SiO m ) ∞ clusters (n,m = 5,6) between 10 and 13 GPa ( Fig. 2b and Extended Data Fig. 7c). Interestingly, the pressures at which OSi 2 -OSi 3 and OSi 3 -OSi 3 structures percolate, approximately 3 GPa and 10 GPa, are similar to the characteristic pressures corresponding to the minimum of the bulk modulus and the elastic-to-plastic transition, respectively. The fact that these mechanical properties correlate with the percolation of (OSi n -OSi m ) ∞ clusters probably underlines the role of the Si-Si non-bonded interactions, as stated by O'Keeffe and Hyde in their model 55,56 , based on the observation that the distance between non-bonded first-neighbour cations in many non-molecular crystals is nearly independent of the bridging atom (anion). When applied to silicates, this implies that the structures are constrained by the problem of fitting large Si atoms around small O atoms. However, the model deserves to be reformulated in view of our results. Indeed, it was developed for stable crystalline minerals, whereas with pressurized v-SiO 2 we are facing unrelaxed structures with connectivities that mimic those of metastable coesites IV and V. These structures are at odds with the third and fifth Pauling rules 24 and show non-monotonic behaviour in their electronic and structural properties around 10-20 GPa (Extended Data Figs. 2 and 4a, respectively). These effects can hardly be accounted for by the model in its current formulation. With a view to future developments, we have noticed that the ratio between the glass density and the O-O distance is almost constant over the entire pressure range explored, suggesting that the non-bonded O-O distance could be an interesting marker to follow.

Code availability
The DFTB+ code is publicly available at https://dftbplus.org/. Additional information may be found there. The percolation code is available freely for non-commercial research at Zenodo