The Origin of Deformation Induced Topological Anisotropy in Silica Glass

The atomic interactions were modelled by the first principles based polarizable force field developed by Kermode et al. [27]. Additional simulations using the BKS potential [28, 29] showed no qualitative differences; see supplementary material Fig. S1. The initial, pristine glass sample was prepared by deforming and rescaling an $\alpha-$quartz structure containing 72000 atoms into a cubic simulation box to match with the experimental density of 2.2 g cm^-3 of silica glass [30]. Periodic boundary conditions (PBC) were used along all directions. This sample was molten at 5000 K in the NVE ensemble for 10 ps, kept at this temperature for 10 ps in the NVT ensemble, and 100 ps in the NPT ensemble at 0 MPa while maintaining a cubic simulation box. The sample was cooled to 300 K in the NPT ensemble at a cooling rate of 10¹³ K s^-1 and 0 MPa, resulting in a glass with a density of 2.39 g cm^-3, see also Tab. S1 in the supplementary material. For this quenching rate, the glass transition temperature was determined to $T_{g}=3000$K. Additional glasses were prepared with cooling rates of $2\times 10^{13}$ K s^-1 and $1\times 10^{14}$ K s^-1. Nosé-Hoover thermo- and barostats were used to perform the constant temperature and pressure simulations [31, 32]. A timestep of 0.1 fs was used for the MD simulations at temperatures above 1000 K and 1 fs for temperatures below.

The equilibrated stress-free isotropic glass sample (’pristine sample’) was subjected to uniaxial deformation by tension or compression at 300 K. The deformation was applied by homogeneously scaling the samples along the $x$-direction at a strain rate of 5$\times 10^{8}$ s^-1, while maintaining zero stress in the orthogonal directions. Additional deformations were performed with strain rates of 1$\times 10^{9}$ s^-1 and 5$\times 10^{9}$ s^-1, see supplementary material Fig. S2. At strains of $\epsilon_{x}$ (where $x$ = 2%, 5%, 7% and 17%, during tension and $x$ = 2%, 5%, 7%, 17%, 20%, 25%, 30%, 35% and 40%, during compression) the sample was unloaded back to the stress-free state. Additional simulations to investigate the effect of the cooling rate and deformation temperature on the stress-strain response and anisotropy index are shown in Fig. S3 and Fig. S4 in supplementary materials.

Deformations at higher temperature to mimic the wire-drawing process were performed at $T_{g}$. Here, the melt was first cooled under temperature and pressure control from 5000 K to 3000 K with a cooling rate of 10¹³ K s^-1 while keeping zero stress in all axes. The stress-free sample at 3000 K was subjected to uni-axial tension or compression up to 17% strain along the $x$-axis while keeping zero stress on the other two orthogonal axes. After the deformation, the samples were quenched to 300 K while allowing the box and atoms to relax the accumulated stress using a cooling rate of 10¹⁴ K s^-1. The resulting glasses have densities of 2.35 g cm^-3 for the glasses pre-deformed in tension and 2.43 g cm^-3 for the one pre-deformed in compression at $T_{g}$.

The change of the orientation of a pair of atoms during the loading and unloading is monitored and quantified by the anisotropy index ($\alpha$). The anisotropy index is a scalar value obtained by calculating the eigenvalues ($\lambda_{i}$) of the fabric tensor $F=<n\otimes n>$ [11] as given in eq. 1, and is normalized in a way that it will vary between 0 and 1,

For an isotropic material, $\alpha$ = 0, and in the case of full anisotropy $\alpha$ = 1, higher values mean higher anisotropy in the medium. We computed $\alpha$ for Si–O, Si–Si, and O–O pairs, and it is referred to as $\alpha(Si-O)$, $\alpha(Si-Si)$, and $\alpha(O-O)$, respectively.

To capture the short-range order’s directionality, a suitable expansion of the RDF onto spherical harmonics must be utilized. Particularly, we take the following expansion into spherical harmonics Y${}_{lm}(\theta,\phi)$

where $l$ and $m$ are the degree and order of the spherical harmonic and $g_{lm}$ is the expansion coefficient. During a uniaxial deformation, the deformation axis will be the symmetry axis. Checking the spherical harmonic functions, $Y_{20}$ best captured this behavior. The expansion coefficients of the RDF into the $Y_{20}$ spherical harmonic is given as,

where $\alpha$ and $\beta$ are the atom types, $N_{\alpha}$ and N_β are their corresponding number of atoms, $V$ is the volume of the simulation box, $x$, $y$, and $z$ are the atomic position in Cartesian coordinates. This function is positive when more bonds between atom types $\alpha$ and $\beta$ of length $\langle|r_{i}-r_{j}|\rangle$ are oriented in parallel to the loading axis than orthogonal to it, and negative when there are more bonds of a certain length orthogonal than parallel to the loading axis.

To the best of our knowledge, the stereographic projection is used for the first time to characterize the structure of oxide glasses in this work. The steps that we used to adapt this method to our needs are detailed in the supplementary materials. First, the tetrahedral face normal calculation was made by finding the nearest neighbors of each Si atom. Then, the four-fold coordinated Si atoms were identified, each Si atom and its four nearest neighboring O atoms were separated, these five atoms together were referred to as tetrahedron, each tetrahedron has four faces, and each face has three atoms. Next, the normal vector of each face of the tetrahedra was computed using the position vectors of the three O atoms. Lastly, these vectors were normalized and then stereographically projected on the required plane as described by the procedure detailed in the supplementary materials.

For the ring normals, a similar approach was used where all rings of a given sample were identified using the RINGS code [33]. Each ring’s Si and O atoms were separated; see Fig. S5(d) for an example with a six-membered ring. The center of mass of each ring was computed based on the position of the atoms, see Fig. S5(e). Then, each ring was discretized into triangles with 2 neighboring Si atoms and the center of mass of the corresponding ring. After that, the normal vector of each triangle was computed, see Fig. S5(f). The ring normal is determined as the normalized average of all its triangle normals, see Fig. S5(g). Lastly, all the normals of the same-sized rings were projected on the required plane using stereographic projection by following the procedure above (see Fig. S5(h - j)).

Bond statistics following [34] was used to determine the number of atoms with new, switched, or broken bonds. This was achieved by comparing the deformed configuration with the pristine glass configuration and finding the differences in bonding topology between them. Each atom has a unique ID. This allows to determine atoms with broken or added bonds by comparing the number of their nearest neighbors with respect to the reference configuration. Atoms with switched bonds have the same number of nearest neighbors as in the reference configuration but have at least one nearest neighbor with a different ID than in the reference configuration.

The usual structural characterization of the simulated silica glass at RT are shown in the supplementary material, Fig. S6, and these are consistent with the data reported in the literature from both experiments [5] and simulations [27]. The stress-strain curves of pristine silica glass subjected to tension and compression at RT are shown in Fig. 1 and show brittle fracture in the case of tension and plastic flow in compression. Unloading simulations to zero stress were performed starting from different values of the maximum strain as depicted in Fig. 1(a and b). Unloading from strains beyond 12% tensile strain resulted in a permanent plastic strain. For compressed samples, a permanent plastic strain was observed for unloading from strains beyond 10%. The densities after unloading from different loads are provided in the supplementary material, Tab. S1.

During straining, the radial distribution function and the bond angle distribution are clearly changed from their pristine state, see Fig. S6 in the supplementary material. However, after unloading, only negligible differences are observed, and these isotropic measures cannot be used to characterize anisotropy. Therefore, following earlier studies [11], the anisotropy was first quantified by the anisotropy index $\alpha$ (see method section for details), which is shown in Fig. S7 in the supplementary material. The increase of $\alpha$(Si-Si) and $\alpha$(Si-O) during the loading up to the yield stress shows that the deformation induces some alignment in the Si–Si and Si–O bonds. This transient anisotropy during loading was present irrespective of the deformation mode. Significant persistent anisotropy was only observed in compression, as shown in Fig. S7(d). Simulations performed at different strain rates and on samples produced with different cooling rates showed that the observed trends regarding transient and persistent anisotropy do not depend on those parameters. However, the degree of anisotropy and the remaining plastic strain depends on the cooling and strain rates; see Figures S3 and S2 in the supplementary material. Deformation at higher temperatures showed that anisotropy is affected by the deformation temperature (see supplementary material Fig. S4).

To characterize in detail preferred orientations in short-range order, we use the spherical harmonic expansion of the pair distribution function [35]. In the case of an axial symmetry, such as in the current study, the radial distribution function is projected on the spherical harmonic $Y_{20}$; see the methods section for more detail. Figure 1(c and d) shows the spherical harmonic coefficient $g_{20}(r)$, which describes the intensity of spherical harmonic term $Y_{20}$ for the pristine glass, loaded to 17% in tension and compression and unloaded from 17% strain to zero stress. In the following, we will focus on these samples; loading up to and unloading from higher strains resulted in qualitatively similar behavior. For the pristine glass, $g_{20}(r)$ shows that the Si–O bonds are distributed randomly without preferential orientation and indicates that pristine silica glass has an isotropic short-range structure. During tension, Fig. 1(c), more Si–O bonds that are longer than the average bond length in a pristine sample ($\langle l_{\text{Si--O}}\rangle=1.6$Å  see supplementary information, Fig. S6) are aligned in the loading direction than in directions orthogonal to it, while shorter bonds tend to be oriented orthogonal to the loading axis. The isotropic structure is mostly recovered when unloaded back to the zero-stress state, as indicated by the absence of noticeable peaks after the unloading. In contrast, the behavior during compression is different from the one in tension, and only one negative peak centered around $\langle l_{\text{Si--O}}\rangle$ is observed, indicating that there are more Si–O bonds of average length oriented perpendicular to the loading axis than parallel to it. Interestingly, the unloaded sample after compression shows qualitatively the same $g_{20}(r)$ signature as the glass loaded during tension: shorter than average bonds oriented mostly orthogonal to the loading axis, and more elongated bonds along the loading axis (Fig. 1(d)).

The glass structure can be further characterized by determining the orientation of tetrahedra face normals, see the methods section and the supplementary material for more details on the method. Figure 1(e) shows the normalized stereographic projections of tetrahedra face normals on the $yz$-plane in both tension and compression for the pristine, loaded, and unloaded samples. As shown for the pristine silica glass, the normalized stereographic projection of tetrahedra face normals shows that the tetrahedra are randomly oriented in all possible directions without having any preferential orientation. In the loaded sample in tension, the normalized stereographic projections indicate that most tetrahedra faces are oriented perpendicular to the tensile axis. At the same level of strain in compression, the normalized stereographic projection shows that most of the tetrahedra faces are oriented perpendicular to the compression axis. In the unloaded samples, the normalized stereographic projections show a uniform orientation of the tetrahedra.

The ring statistic presented in Fig. S6 shows that the number of 6-membered rings, which represent the majority of the rings in our glass, is only slightly affected by the deformation. Figure 1(f) shows the normalized stereographic projection of 6-membered ring normals for the pristine silica, loaded, and unloaded samples (see the methods section and supplementary materials Fig. S5 for the computation of ring normals). Although a uniform orientation of ring normals is expected in pristine silica, insufficient averaging due to the small sample leads to small deviations of the projected 6-membered ring normals from a uniform density, see Fig. 1(f). However, a clear orientation of the 6-membered ring plane normals away from the loading axis is visible for the sample loaded in tension. Under compression, the loaded samples show a strong preferential orientation of 6-membered ring-plane normals along the compression axis. After unloading to zero stress, a uniform orientation of 6-membered ring normals is observed for both tension and compression.

Figure 2(a), shows the evolution of new, switched, and broken bonds during tensile loading to 17% strain and unloading from it. The methodology of computing the bond statistics follows ref. [34] and is described in the method section. The percentage of atoms which switched neighbors before yielding at around 12% tensile strain is less than 0.5%. With yielding, this fraction increases to 1.5% at 17% strain. During unloading, some more bond switching occurs, and then the percentage of atoms with switched bonds stays roughly constant, indicating that the new bonding topology is maintained with unloading. As expected for tensile loading, the percentage of atoms with broken bonds increases, but less strongly than the number of switched bonds, and only a negligible number of atoms adds bonds. While the number of newly formed bonds stays constant, some broken bonds are healed during unloading. During compression, Fig. 2(b), only very few bonds are broken, and more bonds are formed than switched until yield when the percentage of bond switching significantly increases from about 10% strain on. However, even in compression, the percentage of atoms with changed bonding topology remain below 3% at 17% strain. After unloading from 17% strain to zero stress, the percentage of atoms with broken bonds decreases slightly, and the percentage of atoms with new bonds reduces from 0.9% to 0.3%, while the number of switched bonds further increases during the initial stages of unloading and then only slightly decreases.

The orientation of the bonds that will be broken at 17% strain and the bonds that will be newly formed at 17% strain are shown in Fig. 2(c and d). Although in tension only relatively few bonds are broken before reaching the maximal tensile stress, the bonds that break were already in the pristine sample clearly oriented parallel to the loading axis. In compression, the very few bonds that break were preferably oriented at about 60° to the compression axis before deformation. Most of the bonds that were newly formed during tension make at 17% strain an angle between 60 and 90° with the tensile axis. The new bonds formed during compression make at 17% compression preferentially an angle between 0 and 30° with the compression axis.

In the elastic regime in compression, the strain is accumulated by reducing Si–O bond lengths in the direction of compression and/or by changing Si–O–Si angles, leading to transient anisotropy. The major portion of densification in this regime is associated with rearranged network topology upon reduced Si–O–Si bond angles, which is in consistent with Trease et al.  [36]. The reduced Si–O bond length is in consistent with Benmore et al.  [37]. When the sample is unloaded from this regime, the stored elastic strain is released, and the silica glass returns to its initial state of isotropy. Unloading from the plastic regime resulted in a finite plastic strain at zero stress, and the unloaded glass showed structural anisotropy. The transition from elastic to plastic regime is correlated with a strong increase in the number of atoms with switched bonds.

Persistent anisotropy can thus be understood in terms of changes in bond topology through plastic deformation and bond switching, as opposed to simple ring deformation reactions. As shown in Fig. 2(c and d), this bond switching is not random. Due to the Poisson expansion orthogonal to the compressive axis, bonds in this direction are elongated and, therefore, more prone to bond switching or breaking. New bonds are then preferentially formed along the compression axis, where atoms are brought closer together. Upon unloading, the new bonding topology is maintained (only few bond switching events during unloading, Fig. 2(b)). However, the newly formed bonds along the loading axis become now elongated as the atoms move further away. The initial bonding partner, however, can just relax, leading to shorter bonds orthogonal to the loading axis. This explains the $g_{20}(r)$ signature of the persistent anisotropy introduced by pre-compression, Fig. 1(d). The same rationale applies to glasses pre-deformed in tension, however, there the significantly reduced tendency for bond-switching compared to compression leads to only weak anisotropy.

Compression of silica glass is well known to result in densification through (1) decreasing the Si-O-Si bond angles, (2) by enlarging the Si-O bond lengths, and (3) by increasing the number of over-coordinated Si [38, 37]. It is, however, important to note that a) densification is usually observed as persistent densification, i.e., after load removal; and that b) densification is commonly examined after samples were subjected to hydrostatic strain states or similar cases, in which the flow of material is restricted to a certain volume like in the case of indentation. Our observations show that under purely uniaxial compression, i.e. by allowing Poisson expansion, the densification is much less than under comparable hydrostatic strains, see Table S1 in supplementary material. The densification in our simulations of uniaxial compression is due to the reduction of Si-O-Si bond angles, see Fig. S6d in the supplementary material and not by bond shortening (Fig. S6b) or by overcoordinated Si atoms (Fig. S8). Changes in bonding topology, bond lengths and bond angles after compression, can, however, due to different boundary conditions, very well be different for uniaxially compressed and indented silica glass. This is important when comparing our results to complementary work, e.g., by He et al. [39, 40]. Furthermore, scalar measures like bond-length distributions can not provide significant information on densification by uniaxial compression, as changes in different directions can cancel each other out. Finally, the notion of persistent densification implies also the possibility of transient densification which, as our work on anisotropy has shown, might have different topological signatures compared to persistent densification.

While the deformation-induced anisotropy was studied here only for silica glass, the results can be used to discuss what can happen when modifiers (e.g. Na₂O) are present in the glassy matrix. Silica glass has a high structural homogeneity, which renders it difficult to induce structural anisotropy. A higher level of structural anisotropy is expected in glasses that contain structural inhomogeneities due to the presence of highly directional covalent bonds, i.e., Si-O and modifiers-O bonds with less to no bond directionality. The presence of modifiers within the glass leads to depolymerization of the structure by transforming the bridging oxygens (BO) to non-bridging oxygen (NBO). For example, meta-silicate and meta-phosphate glasses have a chain-like structure made of glass former tetrahedra that are linked by BOs. These chains are separated by modifier-NBO bonds, which have low bond strength compared to Si-BO or P-BO bonds. Breaking these chains under load can be harder than breaking the bonds associated with the non-bridging oxygen. As reported by Yang et al. [9] and Endo et al. [8], a higher level of anisotropy can be induced by freezing these chains into a preferential orientation.

The mechanisms and structural changes causing anisotropy can not only be inferred from analysis averaging over the entire sample, Figures 1 and 2, but can also be directly observed in the simulations. Figure 3 and movie M1 show the process in detail, where during compression to 17% strain, the central Si atom of the shown tetrahedron switches a bond to a nearby O atom. This leads to a bond loss in the orthogonal direction to the compression axis and a new bond along the loading axis. During unloading, the newly formed bond does not switch back to the original configuration, but the tetrahedron rotates to adjust for the increase in sample length during unloading. The middle and lower row of Fig. 3 show in addition the change of orientation of the ring and tetrahedra normals as well as the relative bond length of differently oriented Si–O bonds of the tetrahedron. This provides direct evidence for our interpretation of Figures 1(d-e) and 2(c, d).

In many industrial processes like fiber- or sheet-drawing, deformation takes place at temperatures above $T_{g}$. In order to mimic such experimental conditions, we have prepared anisotropic samples by applying the deformation at high-temperatures (HT), i.e., near $T_{g}$, similar to the procedure in Ref. [26], see the methods section for more details. The stress-strain curves, the transient and persistent anisotropy for different deformation temperatures are shown in Fig. S4 in the supplementary material. While the transient anisotropy is like the stress response systematically lower with increasing temperature, the magnitude of persistent anisotropy after quenching from 3000 K to 300 K and relaxing can – for tension – be similar to the one at 300K. However, for compressive loading, the persistent anisotropy does not reach the same high values as for 300 and 600 K. The anisotropy of these samples as characterized by $g_{20}(r)$ is shown in Fig. 4(a and b) together with the samples pre-deformed at RT. This plot clearly shows that the persistent anisotropy of glasses prepared through different routes has different structural origins; i.e., while the situation is more complex for the high-T deformed glasses than the RT deformed ones, the majority of short Si-O bonds are oriented after tensile predeformation at 3000 K orthogonal to the deformation axis, and parallel for compression. The situation is exactly the opposite for RT-predeformed glasses.

The measured transient anisotropy of these samples are significantly lower than for the glasses deformed at lower temperatures. This is attributed to the viscous nature of the high-temperature deformation that involves much more bond breaking than during the mostly plastic deformation at RT, see Fig. S9 in the supplementary material. This prevents the formation of larger, anisotropic structures. This is also evident from Fig. S10 of the supplementary material, which shows in contrast to the situation at 300 K, at 3000 K no preferential orientation of Si-O bonds during loading is discernible. It is therefore not surprising that the quenched and relaxed samples have also a generally low value of the persistent anisotropy factor $\alpha_{\text{p}}$, see Fig. S4 in the supplementary material. However, for tensile deformation, the HT-quench method seems to yield slightly higher anisotropy than in compression and also than the tensile plastic deformation at RT, see Fig. S4 in the supplementary material. This method – akin to wire-drawing – thus yielded the highest persistent anisotropy in tension. This method of inducing anisotropy is expected to be even more interesting in glasses with partially-polymerized structures like in meta-silicate and meta-phosphate glasses as discussed in section 3.3. The higher persistent anisotropy after compression at RT is due to the anomalous compressibility of pure silica glass, which allows for high plastic strains. Introducing anisotropy by compression is therefore best performed at RT.