Fragile-to-strong transition in liquid silica

Network-forming liquids, such as H₂O, SiO₂, Si, Ge, Sb, Bi, and Ga,¹ show complex structural and dynamical features. Their specific properties are due to the fact that they form nearest-neighbor bonds that are strongly directional, with specific values for the bond lengths and bond angles. The number of such bonds is smaller than the usual number of neighbors in a liquid; it has for instance the value of four in tetrahedral networks, as formed by water and silica.² Therefore, the density of the liquid decreases as the bond network is formed, and the structure becomes locally more ordered. Since the formation of the network implies a decrease in entropy and in potential energy, the formation of the low-density networks is favored for lower temperatures and pressures, leading to various anomalies. The most famous of these anomalies, the density maximum of water, has been known for a long time.³ 

Water shows additional anomalies, for instance, in the isothermal compressibility⁴ and the heat capacity.⁵ They are often particularly pronounced for the supercooled liquid.^2,6,7 The unusual structural and thermodynamical properties of water are accompanied by anomalous dynamical behavior, such as a diffusion anomaly.^2,7 Moreover, it was proposed that water exhibits a fragile-to-strong (FS) transition in the supercooled regime,⁸ which means that the temperature dependence of the structural (α) relaxation changes from Vogel-Fulcher-Tammann (VFT) to Arrhenius behavior upon cooling, but an experimental observation of this phenomenon is hampered by a high tendency for crystallization. All these features of water can be explained by using simple models that allow for three nearest-neighbor states with different energy and entropy.^9,10

Silica resembles water in several respects. SiO₂ and H₂O have the same stoichiometry with two A atoms and one B atom. Furthermore, the A-B-A bond angle (O-Si-O or H-O-H) is about the tetrahedral angle which, together with 2:1 stoichiometry, allows for a formation of tetrahedral networks. Finally, like water, silica is polar because of substantial partial charges of the silicon and oxygen atoms. Therefore, one may expect that silica shares many properties with water.¹¹ Consistently, a density maximum was also seen in computational^12,13 and experimental^14,15 studies on silica. It occurs around 1800 K at 1 atm.¹⁵ In simulations, also other anomalies were reported,¹⁶ for example in the specific heat.^12,17 Yet, as compared to water, silica has a lower tendency for crystallization and, hence, it is a better glass-former.

Silica is usually considered as the paradigm of a strong glass former.^18,19 Over a wide range, the temperature dependence of self-diffusion coefficients and structural relaxation times is described by an Arrhenius law with activation energies of ca. 5 eV. However, experimental^20–22 and computational^23–29 studies found that the strong behavior changes to a fragile one when silica is heated to temperatures above 3000 K. Similar to the situation for water, it is currently highly debated whether this FS transition is a consequence of a true phase transition or a smooth crossover. Evidence for a transition from a high-density liquid (HDL) to a low-density liquid (LDL) in silica^30–32 supports the idea that the system may exhibit a liquid-liquid phase transition in the viscous regime. In fact, the above-mentioned simple three-state models¹⁰ lead generically to a thermodynamic phase transition between the HDL and LDL silica forms. However, with different parameter choices the phase transition is replaced by a crossover, and therefore the question is yet unsettled. For silica, crystallization does not interfere with studies of the liquid state in the relevant temperature range. Thus, experimental studies should, unlike for water, be able to determine whether a liquid-liquid phase transition occurs for silica. The relatively high temperatures of ca. 3000 K are a serious drawback to the application of many experimental techniques, though. Conversely, a FS crossover occurring at high temperatures and fast dynamics, i.e., for structural relaxation times of the order of 1 ns, provides ideal conditions for molecular dynamics (MD) simulations approaches.^23–29 

In this paper we extend the existing theoretical studies on silica by performing MD simulations for system sizes and temperature ranges that are considerably larger than before, and by combining the investigation of structural and dynamical properties. Our main finding is a clear change from Arrhenius to VFT behavior for the structural relaxation times as temperature is increased. This change occurs approximately at the same temperature as the density turning point and a change from a low-density silica with high local tetrahedral order to a high-density silica with less local order. For very high temperatures, the structural relaxation reverts again to an Arrhenius law. The change from Arrhenius behavior to VFT and back to Arrhenius is accompanied by a transition from a classical Stokes-Einstein relation to a fractional one and back to a classical one as the temperature increases.

We used the NAMD simulation software package. The MD simulations of SiO₂ in the bulk state were performed with N = 17496 atoms (5832 silicon and 11664 oxygen atoms) in a cubic volume centered about the origin, with periodic boundary conditions. All simulations were performed in the NpT ensemble. To keep the temperature T and the pressure p constant, the Langevin thermostat and the Langevin-Piston-Nosé-Hoover barostat were used. The time step of the integration was set to 1 fs. For temperatures above 5000 K the time step had to be adjusted to 0.5 fs, because otherwise atoms come too close to each other within one integration step. Simulations were performed for up to 135 ns, depending on temperature. The pressure was set to 1 bar. The masses of oxygen and silicon atoms were set to 15.9994 u and 28.0850 u respectively. Prior to data acquisition the system was equilibrated. At all temperatures, we ensured that the equilibration times exceed the structural relaxation times.

We used a modified version of the interaction potential proposed by van Best, Kramer, and van Santen (BKS).³³ The BKS potential proved well suited to reproduce physical properties of silica.^23,24,30 It describes the particle interactions by the combination of a Coulomb term and a Buckingham term V_B. The latter is given by

The use of a modified BKS potential was necessary because the Coulomb and the Buckingham interactions between oxygen and silicon are attractive for small distances. Therefore, we had to ensure that the atoms do not get unphysically close to each other, in particular, at high temperatures where their kinetic energies suffice to cross the energy barrier separating this fully attractive region from the usual minimum region of the Buckingham potential. The modified potential therefore consists of the regular Buckingham potential at larger distances with a continuously differentiable transition to the repulsive part of a Lennard-Jones potential at smaller distances,

Specifically, we switch between both potentials at a distance r_max = 1.09 Å for silicon-oxygen and at r_max = 1.50 Å for oxygen-oxygen interaction, corresponding to the maximum of the Buckingham potential. The parameters used in our simulations are listed in Table I and Table II.

We checked how many atoms penetrated into the region r < r_max. At 10000 K, we found that 0.50 % of the atoms reside in the region where the Buckingham potential was modified. This percentage is relatively small, nonetheless it can change the behavior of the system. However, the percentage decreases fast with lower temperatures, e.g, at 7000 K, this percentage amounts to only 0.07%, suggesting that the modification of the potential does not play a role at T ≤ 7000 K.

We simulated the following temperatures: between 2300 K and 3500 K in 100 K steps and additionally 3800 K, 4100 K, 4400 K, 4700 K, 5000 K, 5500 K, 6000 K, 6500 K, 7000 K, 8000 K, 9000 K, and 10000 K.

The extent to which the atoms order locally can be quantified by the tetrahedral order parameter. The tetrahedral order q_i of the silicon atom i is defined as^35,36

where Θ_ijk denotes the angle between the silicon atom i and two of its four nearest silicon atoms j and k. Possible values of q_i lie between -3 and 1, where the latter value is found for a perfect local tetrahedral structure. The tetrahedral order parameter Q of the whole system is defined as the average over all silicon atoms. If f(q_i) denotes the distribution of the q_i in the system, Q is obtained from the integral

A value of Q = 1 means that all atoms are arranged in a perfect tetrahedral structure, while Q = 0 for an ideal gas. Closely related to the tetrahedral order parameter is the tetrahedral entropy. It is defined as³⁷ 

Here, S₀ is a constant and k_B is the Boltzmann constant.

In order to analyze the structural relaxation of the studied silica melt, we calculate the incoherent intermediate scattering function (ISF). For isotropic systems it can be obtained according to

from the particle displacement vectors $r i → ( t 0 + t ) − r i → ( t 0 )$ during a time interval t. Here, the angular brackets denote the averages over all atoms of a given species and various time origins t₀. F_s(q, t) probes particle displacements on a length scale determined by the absolute value of the scattering vector q. We use q = 2.0 Å ⁻¹ corresponding to the nearest neighbor distance between silicon atoms, which is roughly 3.1 Å. We define the structural relaxation time τ as the time after which the ISF has decayed from 1 to 1/e.

For strong glass formers, the structural relaxation time τ follows an Arrhenius law

Here, E_A is the activation energy and τ₀ is related to the attempt frequency. In the case of a fragile liquid, the temperature dependence of τ is often well described by a VFT behavior

implying a divergence of the structural relaxation time at the temperature T_∞.

The self part of the van Hove correlation function $G s ( r → , t )$ is well suited to explore the mechanism for the structural relaxation. It is defined as³⁸ 

Thus, for isotropic systems the quantity 4πr²G_s(r, t) measures the probability density that an atom has travelled the distance $r=| r → |$ during the time interval t.

The fractional Stokes-Einstein relation for supercooled liquids is^39,40

Here D is the diffusion coefficient and C a constant. A breakdown of the classical Stokes-Einstein relation occurs if θ ≠ 1.⁴¹ 

The temperature-dependent density of the studied silica melt is shown in Fig. 1. In contrast to most other liquids, where the density decreases monotonously with increasing temperature due to the increasing kinetic energy, BKS silica shows a local density minimum at about 3000 K, followed by a local maximum at approximately 4700 K. The steepest increase in density occurs at T ≈ 3400 K. The density anomaly of silica melt was studied in both simulations^12,13 and experiments.^14,15 An experimental study reported a density maximum at 1800 K for ambient pressure,¹⁵ which is much lower than our value. This difference illustrates the fact that force fields in MD simulations are always optimized for certain quantities. Obviously, the density maximum of silica was not a reference value for the parametrization of the BKS potential. Nevertheless, this model correctly reproduces the qualitative behavior.

The structural changes that cause the density anomaly can be seen in the tetrahedral order parameter q_i of the silicon atoms. Fig. 2(a) shows the distribution f(q_i) at various temperatures. For low temperatures the distribution has a sharp peak at q_i ≈ 0.80. With increasing temperature f(q_i) develops a shoulder at lower q_i values, which evolves into a new peak at q_i ≈ 0.43. At 4500–5000 K this new peak becomes the main maximum of the distribution. These findings indicate that the degree of local order changes considerably with temperature. At low temperatures the structure is basically tetrahedral, similar to supercooled water, which exhibits roughly the same value of the tetrahedral order parameter.⁴² At higher temperatures the local structure is significantly less ordered. These results confirm and improve previous results by Shell et al.¹³ who simulated a smaller system comprised of 450 atoms over a narrower temperature range and worked at constant density rather than at constant pressure.

The structural changes are equally well visible in the tetrahedral entropy S_Q(T). The change in tetrahedral entropy $S Q − S 0 k B$ is shown in Fig. 2(b). We see that the tetrahedral entropy strongly increases in the range 3000–5000 K and has a turning point at T ≈ 3100 K, while the temperature dependence is much weaker at both lower and higher temperatures. The observed growth of $S Q − S 0 k B$ in the intermediate temperature range indicates the breakup of the tetrahedral network upon heating. Interestingly, the temperature range of this loss of tetrahedral order agrees well with that of the density increase, providing strong evidence that these structural changes are closely related to the density anomaly.

Fig. 3 displays F_s(q, t) for a wide range of temperatures (2300 K-10000 K). The ISF is shown for the silicon atoms, but it shows qualitatively similar behavior for the oxygen atoms (not shown), consistent with findings in previous studies on BKS silica.^24,28 We see that F_s(q, t) develops the characteristic two-step signature of glass-forming liquids upon cooling. A short-time decay due to vibrational dynamics is followed by a long-time decay due to structural relaxation, which strongly slows down upon cooling resulting in an extended intermediate plateau regime at sufficiently low temperatures. At the crossover between the vibrational and plateau regimes, i.e, at ca. 10⁻¹ − 10⁰  ps, the ISF shows an oscillatory behavior, which was related to the boson peak by some scientists.⁴³ 

Fig. 4 shows the correlation times τ of the silicon and oxygen atoms as a function of inverse temperature. For both types of atoms, we observe a crossover in the temperature dependence at T_cross = 3200 − 3300 K. While the data are well described by a VFT fit above this crossover, they follow an Arrhenius law at lower temperatures. Table III shows the fit parameter. The VFT fits suggest a divergence of the correlation times near T_∞ ≈ 2180 K for silicon and T_∞ ≈ 2340 K for oxygen atoms. Alternatively, a power-law divergence predicted by the mode-coupling theory was used to describe the slowdown of the dynamics in the high-temperature regime.^24,28 In any case MD simulations consistently show that BKS silica cannot be considered as a strong glass-former at T > T_cross. Our Arrhenius fits at T < T_cross yield activation energies E_A of 5.3 and 5.1 eV for the silicon and oxygen atoms, respectively. These values compare reasonably well with the activation energies from previous experimental^44,45 and computational^23,24,28,29 approaches. Interestingly, at temperatures T ≥ 5000 K correlation times obey again an Arrhenius law (see dashed line in Fig. 4). The activation energies E_A in this regime are 1.2 eV, which is far lower than in the supercooled regime.

The observed crossover from VFT to Arrhenius behavior at T_cross = 3200 − 3300 K was also reported in previous MD work on BKS silica.^23–29 It yields evidence for the existence of a FS transition of this model. Previous studies rationalized this effect by a crossover between two liquid phases of silica²⁵ or an existence of a cutoff of the potential energy landscape.²⁹ Our simulations cover for the first time extended dynamic and temperature ranges that include clear Arrhenius behavior at lower temperatures and clear VFT behavior at higher temperatures.

To investigate the mechanism for the structural relaxation, we analyze the self part of the van Hove correlation function. Fig. 5 shows the time evolution of 4πr²G_s(r, t) for oxygen atoms at 2400 K and 3500 K. For short time intervals t, the distribution is dominated by a sharp peak at around 1 Å, which is due to vibrational motion within local cages formed by neighbouring atoms. When extending the time interval this peak broadens. While a single-peak signature is retained at 3500 K, a secondary maximum grows at the expense of the primary maximum at 2400 K. This secondary maximum is located at r ≈ 3 Å, corresponding to the oxygen-oxygen distance. These observations indicate that a hopping motion of the oxygen atoms sets in when the temperature is decreased through the crossover region at T_cross = 3200 − 3300 K. In contrast we do not observe hopping motion for silicon atoms. These results agree with the literature.²³ The fact that only oxygen atoms show hopping motion can be understood by realizing that oxygen atoms have less mass, are twice as numerous, and interact only with half as many neighbors compared to silicon atoms in the tetrahedral structure.

While the van Hove correlation function gives information about the type of motion within the liquid, the average displacement helps to identify the regimes of ballistic and diffusive motion, see Fig. 6, where the mean square displacement (MSD) for silicon atoms is shown. At very short time scales below 0.1 ps, the curves have a slope of 2 due to ballistic motion. At longer times diffusive motion sets in. Between these two regimes a plateau can be seen for temperatures below ca. 3200 K. The plateau value can be interpreted as the maximum squared distance an atom can travel within its cage. Motion over larger distances requires the rearrangement of several surrounding atoms.

Based on the data displayed in Fig. 6, we calculated the diffusion coefficient by fitting the expression 6Dt to the MSD at sufficiently long times. The resulting diffusion coefficient D is plotted as a function of the correlation time τ, divided by the temperature, in Fig. 7. Over a wide range of correlation times the data are reasonably well described by $D∝ τ T − θ$ with θ = 0.89 rather than θ = 1. The corresponding temperature range is 2300 K–10000 K. We note that this breakdown of the classical Stokes-Einstein relation in BKS silica was not observed in previous studies covering a narrower T range.⁴¹ We did not observe this breakdown either when we used a NVT instead of a NpT ensemble. Closer inspection of Fig. 7 reveals, however, that the Stokes-Einstein relation is recovered at low as well as high temperatures, where the solid lines with a slope of θ = 1 fit the data well. This means that the breakdown of the Stokes-Einstein relation occurs in the same temperature regime as the non-Arrhenius behavior.

We performed a comprehensive MD study of the structural and dynamical features of liquid silica and found coherent evidence of a transition between two qualitatively different liquids in the vicinity of 3300 K. The structural features show an anomalous behavior of the density, which increases with temperature between 3000 K and 5000 K, with the maximum slope near 3400 K. Around this temperature, the tetrahedral order parameter changes from a high value (Q = 0.8) to a considerably lower value (Q = 0.4). The temperature at which the tetrahedral entropy changes fastest is around 3100 K. These structural changes are accompanied by dynamical changes. The correlation time shows a qualitative change from an Arrhenius behavior, i.e., from a strong glass at lower temperatures (below 3300 K) to a Vogel-Fulcher-Tammann behavior, i.e., a fragile glass at higher temperatures. At very high temperatures (T ≥ 5000 K), dynamics shows again Arrhenius behavior, which means that silica becomes a simple liquid at sufficiently high temperatures. In the temperature range where the Arrhenius behavior breaks down, the Stokes-Einstein relation breaks also down, and the liquid shows anomalous diffusion with an exponent θ ≃ 0.9 for the mean-square displacement.

All our results were obtained with a pressure of 1 bar. We also performed simulations at higher pressures and found that the transition from strong to fragile behavior moves to lower temperatures. For a pressure of 100 kbar the liquid is fragile for T ≥ 2700 K. This change is expected since the phase with higher density is thermodynamically favored at higher pressures.

Taken together, our results extend considerably earlier evidence for a fragile-to-strong (FS) transition in silica. Such a transition was observed in experiments^20–22 and simulations.^25,26,28 Saika-Voivod et al. found a FS transition in the potential energy surface²⁶ and in the energy landscape.²⁵ Vogel et al. discovered a FS transition in the transport coefficient.²⁸ Horbach et al.²⁴ evaluated the correlation time, but could only see the onset of the transition from VFT to Arrhenius behavior, due to the much narrower range of temperature and time scales. The breakdown of the Stokes-Einstein relation in silica has not been reported before. The transition back to simple dynamical behavior at very high temperatures has not been reported either.

The finding of a FS transition in simulations of small systems can be interpreted in several ways. The different possible thermodynamic scenarios compatible with the obeserved anomalies were discussed in Ref. 46 for liquid water. Heckmann et al. showed that these scenarios are expected for any liquid with the same structure as water,^9,10 and which of the possible scenarios is realized depends on the values of the parameters that characterize the interactions between the molecules, and on the entropy associated with the different possible nearest-neighbor configurations. First, there can be a true thermodynamic phase transition between two liquids with different density and different degree of order. In a small system, this phase transition would be visible as a smooth crossover between the two phases, and the transition temperature can be expected to be in the vicinity of the temperature where the changes are fastest, i.e. around 3300 K. Second, the system may have a true phase transition, but the investigated range of temperature and pressure values is beyond the critical end point of the transition line. The line where changes are fastest is the so-called Widom line, and it is the extension of the phase transition line beyond the critical point. Third, the critical point may be at zero temperature, so that the system shows no true phase transition, but it shows the anomalies caused by this critical point. Our data are compatible with any of these three situations. Future simulation studies or experiments will hopefully be capable of settling this interesting question whether silica or other network-forming liquids show a true thermodynamics phase transition in the liquid phase. If the observed FS transition is a true phase transition in the thermodynamic limit, this transition should become sharper with increasing system size. In contrast to water, silica has the advantage that it can be studied in the supercooled phase in experiments because it does not crystallize that easily. The reason is that the ion bonds in silica have a higher energy than the hydrogen bonds in water, making the barriers higher that have to be overcome when forming the crystal. For the experimental studies this means in particular that no nanoconfinement is required in order to prevent crystallization.

This work was performed within the DFG research unit FOR 1583 and supported by grant numbers Dr300/11-2 and Vo-905/9-2.