Revealing the role of molecular rigidity on the fragility evolution of glass-forming liquids

If quenched fast enough, a liquid is able to avoid crystallization and will remain in a metastable supercooled state down to the glass transition, with an important increase in viscosity upon further cooling. There are important differences in the way liquids relax as they approach the glass transition, rapid or slow variation in dynamic quantities under moderate temperature changes, and a simple means to quantify such variations is provided by the concept of fragility. Here, we report molecular dynamics simulations of a typical network-forming glass, Ge–Se, and find that the relaxation behaviour of the supercooled liquid is strongly correlated to the variation of rigidity with temperature and the spatial distribution of the corresponding topological constraints, which ultimately connect to the fragility minima. This permits extending the fragility concept to aspects of topology/rigidity, and to the degree of homogeneity of the atomic-scale interactions for a variety of structural glasses.

O nce the melting point has been bypassed, both viscosity (Z) and structural relaxation time (t a ) of a supercooled liquid increase enormously as the temperature is lowered. Fragility (M) of glass forming melts has been introduced to quantify such temperature dependences, and a convenient way of obtaining M 1 is to estimate the slope of log 10 (Z) or log 10 (t a ) with temperature, near the glass transition temperature T g where, by definition, one has Z ¼ 10 12 Pa s, and t a B100 s.
This dimensionless slope M, examined for a broad range of glass forming liquids, has been found to vary between a value of 214 (ref. 2) to 14.8 (ref. 3), the latter low value being associated with so-called 'strong' glass formers, whereas liquids with a high M are termed as 'fragile'. When represented on an appropriate plot (that is, log 10 (t a ) versus T g /T), strong glass formers (for example, silica) follow an Arrhenius behaviour of the form exp[E A /k B T], while fragile ones will exhibit an important nonlinear behaviour and a dramatic increase of Z (or t a ) over a limited temperature interval once the supercooled liquid approaches T g . A certain number of correlations between the fragility index M and materials properties have been emphasized, and relationships with T g 2 , structural features 4,5 and elastic properties 6 , have been established. These contributions simply reveal that various aspects of physical properties, structure and interactions influence the glassy relaxation and the fragility.
Here we show by using molecular dynamics (MD) simulations that a typical network-forming system with changing composition displays a certain number of relaxation dynamics anomalies that can be correlated with the strong or fragile nature of the corresponding supercooled liquids. The relaxation properties in network glasses has been quite extensively studied, as in for example, silica 7 . We adopt here, however, a slightly different approach, albeit in line such previous work 7 . Off-stoichiometric glass formers, for example, Ge x Se 100 À x with GeSe 2 (x ¼ 33%) being isochemical to silica, afford, indeed, a unique venue for exploring the dynamics, and understanding ultimately the physical origin of fragility and relaxation in network glasses. First, such studies permit to tune the composition in a continuous fashion, and to compare compositional trends in physical properties, leading to the identification of new correlations. Second, by focusing on a combination of several methods, molecular simulations and rigidity theory 8,9 , one can accumulate direct probes and statistical time-dependent quantities that matter during the glass transition. Third, the focus on key physical behaviour, that is, the small temperature evolution of network rigidity for certain compositions is seen to cause the strong character of corresponding glass-forming melts and is driven by the spatially homogeneous distribution of constraints/ interactions at select compositions. This offers a new interpretation of the physics driving fragility of network glasses.
It is detected that the compositions displaying an anomalous dynamics merely satisfy the well-known Maxwell stability (isostatic) criterion 10,11 , while also weakly depending on temperature. Given that such isostatic compositions also lead to a homogeneous distribution of mechanical constraints at the atomic scale, we argue that fragility deeply connects to the rigidity properties of the underlying network structure, and that fragile glasses are characterized by a spatial distribution of heterogeneous interactions. These conclusions are successfully compared with available experimental data that support the present findings, while also emphasizing the general character of the results for a variety of structural glasses.

Results
Dynamic and relaxation anomalies. To establish our conclusions, we performed First Principles Molecular Dynamics (MD) of Ge x Se 100 À x liquids and glasses at various Ge content (10rxr33%). The modelling scheme 12 reproduces accurately structural data accessed from neutron diffraction in glasses and liquids (Supplementary Figs 1,2 and Supplementary Note 1). The structural changes taking place in Ge-Se are well documented [13][14][15] , and generic to the family of Group IV chalcogenides, that is, glassy systems made of two-fold Se chains are progressively cross-linked by the addition of fourfold germanium atoms that transforms the basic chain structure into a highly connected network (Fig. 1a,b). The addition of such higher coordinated species leads to a global stiffening of the network structure, and by using the notion of Maxwell rigidity 10,11 and an enumeration of constraints n c (refs 8,9), a flexible-to-rigid transition has been identified at the composition x c ¼ 20%, confirmed experimentally at room temperature 15 , and identified with the mean-field Maxwell isostatic stability criterion n c ¼ 3.
The key findings of our study are shown in Fig. 1c,d. Different quantities encoding the dynamics and relaxation behaviour of the Ge-Se melts (at 1,050 K allowing for a validation of the structural models) have been calculated. From the mean-square displacement or i 2 (t)4, one obtains the diffusion constants in the long time limit using the Einstein relation D ¼ or i 2 (t)4/6t. A global decrease of the Ge and Se diffusivity is obtained with increasing Ge content, in line with the stiffening of the network structure that progressively hinders atomic motion. However, a striking feature is the presence of a narrow interval in composition (typically 18-22% Ge) for which D Ge displays an even more marked decrease, about a factor 2-3 less with respect to what would be expected from a smooth decrease of D Ge (broken red line). Further characterization of the dynamics and detecting the anomalous behaviour of the Ge diffusivity is provided by the self part of the Van Hove correlation function 4pr 2 G(r,t), which represents the probability density of finding an atom at time t knowing that this atom was at the origin (r ¼ 0) at time t ¼ 0. In fact, a minimum for the corresponding jump probability is found (inset of Fig. 1c) revealing that the atomic motion will be substantially reduced for this particular interval in composition.
Topological constraints. A central result is the detection of the crucial role played by rigidity in the anomalous relaxation dynamics of Fig. 1, and the indication of a fragility minimum. Within temperature-dependent rigidity theory 16 , the assumption of an Adam-Gibbs form for the relaxation time lnt a p1/TS c , and the contribution of the number of topological degrees of freedom f(x,T) ¼ 3 À n c (x,T) to the configurational entropy 17 S c (x,T) leads to the prediction of the fragility index M: which depends only on the scaling of f(x,T) with temperature. The validity and predictive power of equation (2) has been checked on a certain number of glasses by building structural models, establishing their constraint count, and comparing with a composition dependence of fragility measurements such as for borates 18 or phosphates 19 . MD-based constraint counting algorithms (Supplementary Note 2) are used 20 to yield n c (x,T) from Ge/Se bond-stretching (BS) and bond-bending (BB) motions, and they show that (i) for the 300 K glasses, n c follows exactly the mean-field estimate n c ¼ 2 þ 5x ( Fig. 2a) 8,9 , leading to an isostatic condition (n c ¼ 3) for B20% Ge that agrees with the experimentally measured threshold value 15 , and, importantly, (ii) the calculated n c (x,T) displays a minimum change with temperature for 20-22% Ge. The obtained minimum change is robust (see Supplementary Note 3). Following equation (2), this leads to a minimum for the liquid fragility with composition given the weaker variation of the term dn c /dlnT at xB22%, and coincides with n c ¼ 3.
Viscocity anomaly. Independent support for this conclusion is given by a Stokes-Einstein estimate of the liquid viscosity.
Given the system size (250 atoms), and First Principles MD timescale (ps), a direct calculation of the Green-Kubo viscosities is out of reach. However, the validity of the Stokes-Einstein Z ¼ k B T/6prD where 2r ¼ r À 1/3 involves the liquid densities r, can been established from calculated diffusivities for the supercooled state, similarly to a previous application 21 , and this provides a validation of the dynamic behaviour of the models (see Supplementary Note 4). Next, one can derive from the calculated diffusivities D Ge (x,T) (Fig. 1a) corresponding viscosities using the Stokes-Einstein relationship. Results using such calculated viscosities are shown in Fig. 3, and a direct comparison in an Angell plot represents calculated and measured log 10 (Z) as a function of T g /T. Note that we have checked that the Stokes-Einstein remains valid for the considered range of temperatures. The agreement is found to be very good for the two compositions for which a direct comparison is possible (Ge 20 Se 80 and Ge 33 Se 67 ). Having validated the approach, we then use it to extract the viscosity Z(x,1,050 K) along the 1,050 K isotherm (Fig. 1d).
A maximum for Z is obtained that correlates close to the isostatic composition (Fig. 1d), with the maximum of the structural relaxation time t a calculated from the a-relaxation regime of the intermediate scattering function (Supplementary Note 5 and Supplementary Fig. 5). This provides an additional indication that strong glass-forming liquids with a minimum fragility are obtained for the 22% composition.

Discussion
When the details of the constraint contributions are investigated, and their spatial distribution analysed (Fig. 2b), additional important features can be linked to the weak variation of n c with temperature. First, for the composition at which the weakest change in n c is found (22%), the BB constraints are nearly homogeneously distributed in the structure given that n c BB maximizes to its nearly low temperature value (n c BB ¼ 5), which reduce the possible fluctuations (Fig. 2c). This contrasts significantly with all the other compositions (for example, 14.3 and 33%, Fig. 2b), which display a heterogeneous distribution of constraints broken by thermal activation (Fig. 2d). The variance of n c BB is also a measure of the spread in the spatial distribution of atomic interactions leading to rigid constraints. For the compositions in the region 20-23%, one now acknowledges that homogeneity of stress (that is, interactions) will not only induce a cascade of dynamic anomalies as detected from the diffusivity minima and viscosity/relaxation time maxima located in the same range of compositions ( Fig. 1) but also represents a more stable state at the nanoscale. As a consequence, thermal changes must lead to minor changes of n c with temperature, and ultimately will minimize fragility (equation (2)). In this particular region of anomalous behaviour, one furthermore detects weak effects on structure with chemical composition given that both the pair correlation function g(r) in real space and the structure factor in reciprocal space exhibit small differences with Ge content (see Supplementary Note 1 and Supplementary Figs 1,2). For instance, the structural properties of the 22% composition that appears to be so particular in terms of dynamic properties (Fig. 1), are found to be nearly identical to the 20% one, although their relaxation time differs by a factor 2 (Fig. 1b) and their rigidity also varies substantially within a tiny compositional interval. This, ultimately, emphasizes the predominance of rigidity and its spatial distribution on the liquid relaxation properties, over aspects of structure. This physical picture revealing the central role played by network rigidity is consistent with conclusions drawn from a compilation of fragility data on network glasses including chalcogenides and oxides for which the location of the isostatic composition is determined from calorimetric measurements 3,27-30 . These liquid fragilities are represented in Fig. 4, and show that once the composition is rescaled with respect to the centroid x c of the isostatic composition interval (for example, x c ¼ 22.5% in Ge x Se 100 À x ; ref. 3), all the data exhibit a minimum in M at xBx c . This not only emphasizes that rigidity affects the fragility evolution of melts, a qualitative correlation that has been reported in the literature for quite some time (ref. 3) and reference therein), but from Fig. 2b, one also realizes that the spatial distribution of atomic scale interactions/constraints is a key feature for the understanding of transport properties during the glass transition. This conclusion can be only drawn from a detailed inspection of constraints, accessed from MD on individual atoms. Although further detailed analysis on the constraint behaviour of, for example, organic glass formers is necessary, these results may indicate that liquids with nondirectional bonding and, therefore, with a more probable heterogeneous distribution of interactions will lead to a much more fragile behaviour.    22 Senapati et al. 23 and Guegen et al. 24 Note that only Stølen et al. have investigated the high temperature region where numerical data can be directly compared. The solid lines correspond to fits using the Mauro-Yue-Ellison-Gupta-Allen (MYEGA) equation 25 that describes with an increased accuracy the high temperature viscosities as compared with alternative fitting functionals 26 (for example, Vogel-Fulcher-Tamman, VFT).

Methods
First-principles molecular dynamics simulations. Ge-Se liquids and glasses have been investigated using Car-Parrinello molecular dynamics simulations. The system contained 250 atoms, and up to nine compositions Ge x Se 100 À x have been simulated in NVT ensemble with cubic cells of sizes allowing to recover the experimental density 31,32 of corresponding liquids and glasses (that is, 19.97 Å for GeSe 9 ). Density functional theory has been used to describe the electronic structure that evolved self-consistently with time. We have adopted a generalized gradient approach using the exchange energy obtained by Becke 33 and the correlation energy according to Lee, Yang and Parr (LYP) 34 . The BLYP approach was used due to its account on valence electron localization effects, and an increased agreement with experiments on structure as revealed by a previous study 12,35 . Valence electrons have been treated explicitly, in conjunction with norm conserving pseudopotentials of the Trouiller-Martins type to account for core-valence interactions. The wave functions have been expanded at the G point of the supercell on a plane-wave basis set having an energy cutoff E c ¼ 20 Ry which is a standard value for the investigation of chalcogenides 36 . A fictitious electron mass of 200 a.u. was used in the first-principles molecular dynamics approach. The time step for integrating the equations of motion was set to Dt ¼ 0.10 fs. The temperature control was achieved for both ionic and electronic degrees of freedom using Nosé-Hoover thermostats. The initial coordinates of 250 atoms have been constructed using the atomic positions in a GeSe crystal. To achieve the correct compositions, Ge atoms were then replaced with Se atoms in an appropriate fashion, depending on the target composition.
We carried out simulations at T ¼ 2,000 K for a period of 22 ps to lose the memory of initial configuration, and then investigated a certain number of isotherms (for 25-30 ps each): 2,000 K, 1,600 K, 1,373 K, 1,200 K, 1,050 K, 800 K and 300 K. The first 2 ps at all the temperatures have been discarded. The 300 K trajectories result from three independent quenches with starting temperature 1,050 K.
Certain compositions have been particularly checked, for example, the 22 and 23%, and a series of five trajectories at 1,050 K obtained from five independent high-temperature liquid starting configurations at 2,000 K have been analysed. The statistical average calculated quantities (for example, n c , n c BB , D Se , and so on) have been found to be consistent (see Supplementary Note 3). Similarly, effects of size (480 atoms and the present 250 atoms) and excursion in density have been tested and have confirmed that the constraint count does not depend on these different simulation conditions. For instance, for the 23% composition, it has been found that n c ¼ 2.84 for N ¼ 480, a value that is very close to the value found for N ¼ 250 (n c ¼ 2.79).