Topological Constraint Theory for Network Glasses and Glass-Forming Liquids: A Rigid Polytope Approach

A variation of the topological constraint theory is proposed where an atomic network is modeled as a collection of rigid polytopes, and which explicitly distinguishes the bond angle constraints as well as rigid bond angles from flexible ones. The proposed theory allows for direct quantitative estimation of the fraction f of zero-frequency or floppy modes of the network. A preliminary model is proposed to connect the theory to the two key experimental observables that characterize glass-forming liquids, i.e., the glass transition temperature T g and fragility m . The predicted values are tested against the literature data available for binary and ternary chalcogenides in the Ge-As-Se system. The T g is related to f in this model by the activation entropy associated with the bond scission-renewal dynamics that is at the heart of transport and relaxation in glass-forming liquids. On the other hand, the large and temperature-dependent conformational entropy contribution of the 1-polytopes, i.e., the selenium chain elements in these chalcogenide glass-forming liquids, plays a key role in controlling the variation of m with f .


I. INTRODUCTION
The deformability, entropy and even temperature-dependent disintegration modes of a random atomic network are largely determined by the network's rigidity. Seminal works by Phillips [1,2] and Thorpe [3,4] have shown the existence of a floppy-to-rigid transition in such networks in the form of rigidity percolation. The prevailing model suggests that the network is floppy when the number of degrees of freedom per atom exceeds the number of interatomic force field constraints, and the transition to rigidity occurs when these two quantities are equal. The network's rigidity therefore depends on the average coordination number ⟨r⟩, and various arguments suggest that the floppy-to-rigid transition occurs at the critical value r p = 2.4. Although Phillips [1,2] originally attempted to relate the glassforming ability of a network to its entropy and deformability, Thorpe [3] placed the idea of network rigidity on more formal footing by considering the number of zero-frequency vibrational modes, i.e., continuous deformations with no energy penalty. The fraction f of such modes is predicted to decrease with increasing ⟨r⟩ to zero at or near r p where rigidity percolates, and increasing the connectivity further only serves to overconstrain the network. The rigidity percolation model suggests that there could be an underlying universal dependence of the physical properties of glasses on composition that is effectively insensitive to their chemical details. The chalcogenide glasses provide an ideal testbed in this regard, being characterized by energetically similar covalent bonds and a wide range of network connectivity. As it is practically impossible to directly obtain f for these complex amorphous networks, typical studies in the literature instead report the variation of physical properties with ⟨r⟩ [5][6][7][8][9][10]. However, it remains unclear what would be the signature of the rigidity percolation in the physical properties of amorphous networks as a function of the connectivity.
In fact, a wide variety of physical properties including elastic moduli, refractive index and the glass transition temperature T g of glasses in the Ge-Se system show continuous and monotonic change through r p [8]. On the other hand, Boolchand and coworkers have suggested that the floppy-to-rigid transition may not be sharp at r p , but rather that the network is optimally constrained and therefore stress-free over a range of ⟨r⟩ around the threshold value of 2.4 [11][12][13][14]. These authors used modulated differential scanning calorimetry (MDSC) to show that Ge-Se glasses with such optimally constrained networks are characterized by vanishing non-reversing enthalpy, and identified this as a property of the "Intermediate Phase" 2 (IP). It was further suggested that these IP glasses would be completely resistant to aging or relaxation below T g . However, subsequent direct aging studies of IP glasses in the Ge-Se and Si-Se systems by other researchers have shown the presence of both structural and enthalpy relaxation [15][16][17], thus questioning the validity of the existence of such IP.
Considering that T g depends strongly on the connectivity of a network, its monotonic variation with ⟨r⟩ raises further questions about the use of the average coordination number as a governing variable for rigidity percolation. For example, the Gibbs-DiMarzio model of the chain polymer glass transition [18,19] predicts that T g (along with the available free volume, chain stiffness and degree of polymerization) is described by the monotonic function T g = T 0 /(1 − κX) where X is the cross-linking density of the chains and κ is a universal constant. Note that while the chain length in chain polymers is not affected by cross-linking, this does result in progressively shorter chains for chalcogenides. Despite this difference, Varshneya et al. [20] found that for chalcogenide networks the compositional variation of T g could be described by the same equation if κX is replaced by β(⟨r⟩−2) where β is a system-dependent parameter. Naumis [21,22] included the effect of floppy modes on the vibrational density of states to provide theoretical justification for this observation by suggesting that the Lindemann criterion for atomic displacements at the melting point could also be applied to the glass transition. Other previous studies [23,24] have shown the existence of a fundamental connection between the temperature dependence of the atomic mean square displacement and viscous flow or shear relaxation in glass-forming liquids.
However, the glass transition is not a true thermodynamic transition and is only significant for the fact that, by definition, the structural relaxation timescale at T g is on the order of ∼100 s. Therefore, the validity of applying the Lindemann criterion to the glass transition remains questionable.
On the other hand, several studies have reported sharp changes in the activation energy of viscous flow near T g , or in the related fragility parameter [25] where η is the viscosity, in the vicinity of T g . That said, the original form of the rigidity percolation model cannot explain temperature-dependent dynamics in supercooled liquids because it does not include any mechanism for the thermal energy to overcome the interatomic constraints. Gupta, Mauro and coworkers attempted to address this by explicitly 3 introducing temperature-dependent network constraints in the rigidity model [26]. This involves assigning each network constraint a switching temperature such that the constraint becomes active only when the system temperature drops below the switching point. Although this approach is used extensively in the literature to model the compositional variation of the thermophysical properties of a wide range of glass-formers, direct experimental verification of the key assumption about the switching behavior has yet to be made.
Naumis [22,27] attempted to connect the rigidity percolation model to statistical mechanics by suggesting that the floppy modes provide channels in the potential energy landscape that serve as pathways for the network to explore many local minima, and hence make a contribution S c to the configurational entropy. This entropy was calculated as where Ω is the number of accessible microstates per atomic degree of freedom and is independent of f . Floppy modes indeed provide an important source of entropy, though their connection with structural relaxation is not obvious from the standpoint of the energy landscape where metabasin hopping and vibrational excitations within the metabasins are expected to be temporally decoupled. A possible connection between the short-and longtimescale processes has been proposed by Dyre and coworkers with their elastic "shoving" model, where the rapid increase with cooling of the activation energy for structural relaxation in a fragile glass-forming liquid is attributed to a corresponding anharmonic increase of the high-frequency shear modulus [28]. Additionally, an exact solvable glass transition model by Naumis and coworkers [24] has related the short-time process to the frequency of probing transition states of the energy landscape, while the long-time relaxation process represents the transition between metastable states.
The rigidity percolation model is undoubtedly useful as means of understanding the behavior of random networks around the glass transition. Nevertheless, the model makes some predictions that are not consistent with experimental observations, and the average atomic coordination number ⟨r⟩ is not always an appropriate measure of the network connectivity (Section II A gives an example), at least not without qualifications described by Thorpe [3].
Moreover, a connection between the model and the statistical mechanics of glass-forming liquids has yet to be conclusively established. That is to say, the topological constraint theory of glasses should not be considered complete and inviolable.
This article proposes a variation of the topological constraint theory where, for the purposes of calculating the available degrees of freedom, the network is considered as being 4 constructed from a collection of rigid polytopes. This view is supported by, e.g., Sidebottom's observation that the mean connectivity of a network's weakest links is more closely related to the fragility than the mean coordination number [29,30], and is related to the concept of rigid unit modes in the context of negative thermal expansion materials [31,32].
The theory is motivated and developed in Section II, along with a preliminary model that relates the available degrees of freedom to the relaxation time of a glass-forming liquid. The resulting equations are tested and evaluated for the family of chalcogenide glasses in Section III, and suggest that this variation of the topological constraint theory could be useful more generally.

II. NETWORKS OF RIGID POLYTOPES
This section revisits the rigidity percolation model. First, the limits of this model are Thorpe's rigidity percolation model [3] begins by classifying atoms in a random network by the number of bonds in which they participate, without distinguishing bonds of different types. If n r is the number of atoms with r bonds, then 3 ∑︁ r n r is the number of degrees of freedom. Specifying the bond lengths and bond angles around an atom with r bonds imposes r/2 and 2r − 3 constraints, respectively, on the degrees of freedom. If all of these constraints are independent, then the fraction of unconstrained degrees of freedom f is given The average coordination number ⟨r⟩ is defined to be ∑︁ r n r r/ ∑︁ r n r (despite the misprint in Ref. [3]), reducing this to the standard equation Since f = 0 when ⟨r⟩ = 2.4, the model is often believed to predict that rigidity should percolate through the network and various properties of the system should be discontinuous around this critical value. Thorpe's original article included several qualifications to this conclusion that are unfortunately not consistently discussed in the literature, but that indicate that the situation is more nuanced than suggested by Eq. 3.
Consider applying this model to the regular network of corner-sharing octahedra in a perovskite material, as shown in Figure 1a. If the formula for this material is ABO 3 , then a single octahedron contains one B atom with six bonds and three O atoms with two bonds each. The average coordination number for the octahedral network is then three, and the model would seem to predict that the network is not only rigid but highly overconstrained.
Nonetheless, it is widely recognized that various patterns of octahedral tilting [33], including the one in Figure 1b, are low-frequency modes of the network. In fact, such soft modes are believed to be responsible for displacive phase transitions and negative thermal expansion 6 in a number of systems [31,32]. Possible resolutions to this apparent contradiction are considered below to help to clarify the limitations of the model as frequently interpreted, and to suggest refinements to be included in the subsequent derivation.
Our first attempt begins with the observation that the symmetry of the ideal system allows only a constant number of octahedral tilting modes, independent of the number of atoms. Hence, as the system size approaches the thermodynamic limit, the fraction of zerofrequency modes goes to zero and the model gives the expected result. This argument is not entirely satisfying though; the number of octahedral tilting modes is constant because they extend throughout the entire network, but it is difficult to exclude the possibility of localized zero-frequency modes whose number scales with the size of a disordered system. Thorpe discussed this possibility, saying that "the number of zero frequency modes is not zero at ⟨r⟩ = r p because floppy inclusions still exist" [3]. Perhaps then the presence of localized zero-frequency modes is less relevant to the system properties than a percolating network of identical ideally rigid bonds. The difficulty with this is that bonds in physical systems are neither identical nor ideally rigid; van der Waals forces allow systems with ⟨r⟩ < r p at room temperature to solidify at low temperatures, and bond breaking allows systems with ⟨r⟩ > r p at room temperature to liquify at high temperatures. The conclusion here is the same as the one arrived at by Thorpe [3], namely, one should neither expect to observe discontinuous behavior in the number of zero-frequency modes, nor in the properties of glass-forming liquids, at any particular value of the average coordination number.
Our second attempt considers the possibility that the number of constraints is overestimated in Eq. 3. Initially observe that if the bond angles of the oxygens at the corners of the octahedra were fixed, then the octahedra would not be able to tilt. Thorpe anticipated this possibility as well: "For example in Si x O 1−x ...it is reasonable not to count the angular force at the oxygen atoms" [3]. Since the bond angle constraints do not explicitly appear in Eq. 3, such a modification needs to be performed using Eq. 2. Unfortunately, excluding the oxygen bond angle constraints is insufficient to make f nonnegative, and only increases the apparent value from −0.5 to −0.25. This idea should instead be taken much further, with only independent constraints included in Eq. 2 [3]. The difficulty with this is that while many constraints are dependent for the octahedral network, distinguishing the dependent from the independent ones is not at all obvious due to the interactions being nonlocal (e.g., a set of bonds could form a ring). While there do exist approaches to identifying the set of 7 independent constrains (e.g., the pebble game [34,35]), these require much more detailed knowledge of the network configuration than is generally available in practice.
Perhaps since all of the quantities in Eq. 2 are defined using only local information, f is more closely related to the properties of an average local environment than the overall network. Indeed, f could be interpreted as an estimate for the average fraction of uncon- Section II B proposes a revised topological constraint theory that is motivated by the first and second points, while the third is the subject of Section II C.

B. Revised topological constraint theory
The concept of a metabasin requires separating two types of relaxations in glass-forming liquids, namely, local rearrangements of particles and substantial structural relaxations [36].
If the available kinetic energy is such that the system is confined to a single metabasin, then the system is expected to behave as a solid. Conversely, a system that often experiences the substantial structural relaxations associated with transitions between metabasins explores more of the configuration space and is expected to behave as a liquid. This suggests that characterizing the metabasins is essential to understanding the glassy state, and much effort has been expended in this direction [37,38].
Suppose that the local rearrangements in this picture can be identified with operations involving a small number of covalent bonds, and structural relaxations with changes to the network conformation or connectivity. Constraints associated with individual bonds would then be less relevant to understanding the glass transition than constraints imposed by the network connectivity on larger structural units. This serves as motivation for the topological constraint theory developed in this section, where the network is considered as being composed of rigid polytopes rather than individual atoms, and any relaxations involving bonds within a rigid polytope are explicitly not considered. While others have used rigid polytopes for this purpose before [39][40][41], and even general networks composed of multiple types of rigid polyhedra with edge and face sharing [42], those authors specifically considered the restricted question of whether topologically-disordered networks could exist.
Moreover, those theories consider the number of degrees of freedom per vertex (rather than the fraction of unconstrained degrees of freedom) and require considerably more detailed knowledge of the constituent polyhedra than the theory developed here.
A polytope is a generalization of polygons and polyhedra to any nonnegative dimension.
More precisely, an i-polytope is a finite region of i-dimensional space that is bounded by a finite set of (i − 1)-dimensional hyperplanes. The boundary of a polytope is itself comprised of polytopes, and a j-polytope of this type is called a j-facet [43]. For example, a polyhedron .
It is tempting to say that ϕ(d) is the number of degrees of freedom made available to the network by a direct d-polytope, but this is not necessarily the case; Figure 4 shows several situations where a direct d-polytope contributes fewer than ϕ(d) degrees of freedom. For example, the direct 3-polytope on the right constrains the distance between the connected intersection 0-polytopes, but can rotate about the relevant edge without affecting the configurations available to the network (direct polytopes not connected by intersection polytopes do not interact). That is, this direct 3-polytope effecively functions as a direct 1-polytope.
Along with the other examples in Figure 4, this suggests that a direct d-polytope actually contributes only ϕ(δ) degrees of freedom, where δ ≤ d is the dimension of the convex hull of the connected intersection polytopes and is called the reduced dimension.
. f . This will imply functional dependencies of the glass transition temperature T g and the fragility m on f that will be compared with experiments in Section III.
Given a single primary relaxation mechanism, a relaxation time τ can be modeled as having an Arrhenius dependence on temperature, at least over a limited temperature range near T g . The relaxation time τ can then be simply written as The free energy of activation ∆G can be written in terms of an activation enthalpy and entropy. If only the the activation entropy ∆S depends explicitly on f , the expression for τ can be re-written as This leaves only the construction of a model for ∆S(f ) to be able to relate f to an experimental observable. where A, µ and γ are constants [44,45]. Using Boltzmann's formula S = k B ln Z, the conformational entropy of a single chain depends on n as where k B is Boltzmann's constant. Since entropy is an extensive quantity, the change in configurational entropy when breaking a bond to convert one chain of length n into two chains of length ∼ n/2, as on the right of Fig. 5, is This is interpreted as the activation entropy for the bond scission-renewal dynamics. While many distinct microscopic events are expected to contribute to the effective bond-breaking mechanism in practice, the reasoning above is intended only to motivate the form of the dependence of ∆S on n.
The dependence of n on f for the system on the left of gives for the change in configurational entropy, where a, b and c are constants. The value of a is expected to be around 0.157 from the literature [44,45], and b is expected to be in the interval 0.4 ≤ b ≤ 0.5 by the above reasoning. Finally, substituting the equation for ∆S(f ) into Eq. 5 gives for the relaxation time, where τ ′ 0 = τ 0 exp(−c). If there is a characteristic relaxation time τ g ≈ 100 s that is a universal constant at the glass transition temperature T g , then Eq. 8 can be inverted to find T g as a function of f : The fragility m, on the other hand, is related by the Adam-Gibbs model of relaxation [46] to the temperature dependence of the configurational entropy S c : Sidebottom [29] has recently suggested that the temperature dependence of bonds [47]. The analysis of S-rich glasses is further complicated by the presence of S 8 rings that do not participate in the network [48]. The chalcogen-deficient glass-forming compositions are not considered either, since such networks are over-constrained and the estimation of f for such networks is not physically sensible. Moreover, the structure of chalcogendeficient Ge-As-Se networks is known to be complicated by the formation of molecular elements of the type As 4 Se 3 , As 4 Se 4 and As 4 , as well as the appearance of Ge-Ge bonds in ethane-like Se 3/2 -Ge-Ge-Se 3/2 units and Ge-As bonds [7,49,50]. The quantitative estimation of the relative concentrations of these structural units and their effects on the network rigidity are not straightforward and will be considered in the future. Finally, it is well known that the GeSe 4/2 tetrahedral units in the structural network of Ge x Se 1−x glasses over the composition range 0 ≤ x ≤ 0.33 are predominantly corner-sharing, but a small fraction (15-35%) that increases monotonically with x form edge-sharing tetrahedral units [15,51].
Sidebottom [29] suggested that the edge-sharing tetrahedra provide more degrees of freedom to the network compared to their corner-sharing counterparts since the shared edge does not participate in the connectivity of the rest of the network. However, this does not account for the fact that edge-sharing between tetrahedra or other rigid polytopes can severely over-constrain the local degrees of freedom of the polytope itself [17,41]. Here we assume that these effects approximately cancel, and that all of the tetrahedra are corner-sharing at the concentrations considered below.
Now consider the fractional degrees of freedom f of a chalcogenide network with composition Ge x As y Se 1−x−y that contains n a atoms and n d direct polytopes. Let α be the fraction of direct polytopes that are GeSe 4/2 tetrahedra (3-polytopes), β be the fraction that are AsSe 3/2 triangles (reduced 2-polytopes), and (1 − α − β) be the fraction that are Se 2/2 line segments (1-polytopes). The first step to calculate f is to find expressions for the direct polytope fractions α and β in terms of the atomic fractions x and y. Equating the number of Ge, As, and Se atoms with those that appear in the direct polytopes gives the following system of equations: Solving for the number of direct polytopes n d as a function of x and y and subtituting the result into the first and second equations gives for the fractions of direct 3-polytopes and direct 2-polytopes. The second step to calculate f is to find the types and numbers of all intersection polytopes. Given the assumption that all direct polytopes are corner sharing, the network contains only intersection 0-polytopes, each with two connections to direct polytopes and no associated angular constraints. Since there is precisely one intersection 0-polytope for every Se atom, the number n i of intersection 0-polytopes is Evaluating Eq. 4 for a network containing αn d direct 3-polytopes, βn d direct 2-polytopes, (1 − α − β)n d direct 1-polytopes, and n i intersection 0-polytopes gives for the fraction of unconstrained degrees of freedom.
The value of f is compared with the average coordination number ⟨r⟩ = 2x + y + 2 in Fig. 6 for the binary and ternary Ge-As-Se compositions considered in this study. That f should nearly be a function of ⟨r⟩ is surprising, given that both f and ⟨r⟩ are functions of two independent variables relating to the number of tetrahedral and pyramidal cross-linking elements. It is also immediately apparent that the variation of f with ⟨r⟩ is quite different from that obtained by Thorpe using Eq. 3, since that version of f decreases linearly with ⟨r⟩ < r p until the sharp transition at r p = 2.4. The f obtained using Eq. 11 is instead nearly a nonlinear function of ⟨r⟩ that goes to zero at r p = 2.67. As a result, the model developed in Section II B predicts that As 2 Se 3 and GeSe 4 with ⟨r⟩ = 2.4 are underconstrained networks, whereas Eq. 3 predicts that they are isostatically rigid; the source of this difference is in the application of the angular constraints around Se atoms. This prediction of our model calls into question the existence of the IP as defined elsewhere [11][12][13][14]. Furthermore, one should expect a rapid increase in the elastic moduli beyond r p = 2.67 as the network is increasingly over constrained; such behavior is indeed observed for a wide variety of binary and ternary Ge-As/Sb-S/Se chalcogenide glasses [5,8].
The glass transition temperature T g for glass-forming liquids represents an isoviscous temperature where the viscosity is ∼ 10 12 Pa · s and the relaxation time is ∼ 100 s, and is expected to be a monotonic function of the average connectivity of a network. Since  FIG. 7. The glass transition temperature T g for binary Ge-Se, As-Se, Si-Se and ternary Ge-As-Se glasses is nearly a linear function of the fractional degrees of freedom f . The black line is a fit using Eq. 9 for the parameter values reported in the text. T g data are from [7,9,[52][53][54][55].
the bond strengths of the homopolar Se-Se and heteropolar Ge-Se and As-Se bonds are comparable, the T g of Ge-As-Se glasses should monotonically decrease with f . This is borne out in Fig. 7, which shows a nearly universal and linear decrease in T g with increasing f in these glasses. The black line in Fig. 7 is given by fitting Eq. 9 to the data using τ g = 100 s, τ ′ 0 = 10 −13 s, b = 0.45, and a standard conjugate gradient minimization algorithm for the remaining parameters. This procedure gives ∆H = 3.552 ± 0.001 eV and a = 36.32 ± 0.02, with parameter uncertanties estimated by calculating the Hessian matrix at the minimum.
An activation enthalpy ∆H of around twice the bond energy is reasonable enough for a scission event given the simplicity of the model developed in Section II C, and is roughly consistent with the value of ∼ 2.7 eV reported for the Ge x Se 1−x system as reported in the literature [57]. The deviation of a from the expected value of 0.157 is more dramatric, with the activation entropy being ∼ 200 times larger than expected. This could be related to the number of configurations in SAW models generally being calculated for chains on lattices rather than in continuous space, or to the neglect of the effect of a scission event on the conformations available to the surrounding network components.
The dependence of the fragility m on f for the same Ge-As-Se liquids is shown in Fig.   8. A nearly universal functional dependence between m and f is again observed for glassforming liquids in this system, where m is nearly constant ∼ 30 for 0 ≤ f ≤ 0.3 and rapidly increases to ∼ 80 for 0.3 ≤ f ≤ 0.4. Recent studies have suggested that this behavior is related to the disappearance of the conformational entropy of the selenium chain segments as their average length is reduced to ∼ 3 [22,29], and is consistent with the prediction made in Section II C for the corresponding threshold value of 0.30 ≤ f ≤ 0.33. By comparison, the onset of the sharp rise in m coincides with ⟨r⟩ ≈ 2.3, below where the conventional topological constraint theory predicts a transition.

IV. CONCLUSIONS
A rigid polytope model of glass structure is developed and used to calculate the relative fraction of unconstrained degrees of freedom f , the result having a fundamentally different dependence on ⟨r⟩ than that originally obtained by Thorpe [3] using a mean-field approximation. The variation of T g and m of binary and ternary chalcogenide glass-forming liquids in the Ge-As-Se system shows nearly universal dependence on f over a wide range of com-