Concepts and applications of rigidity in non-crystalline solids: a review on new developments and directions

Abstract We present the most recent developments and applications of Rigidity Theory in atomic-scale disordered materials such as glasses, liquids or amorphous solids. While the framework of rigidity, elastic phase transitions and topological constraint counting has been extensively used in glass science, its potentially wide applicability has remained difficult because of some inherent limitations of the theory. Here, we review promising methods that have been recently introduced, and which allow to overcome such limitations, and to extend these methods to the case of densified liquids or partially ordered solids, and also to the non-mean-field case. We introduce and describe a several methods based on an explicit account of temperature effects. In addition, the calculation of bond length and bond angle standard deviations is presented, accessed from molecular simulations, which allows to estimate atomic stretching and bending topological constraints as a function of pressure or temperature. These new methods are illustrated through a set of applications ranging from several archetypal glass-forming liquids (SiO2, GexSe1−x, …) to semiconductors (Si-C:H), phase change materials (Ge-Te) and cement, while also emphasizing the role of isostatic composition for the design of stable glasses with reduced ageing. Results furthermore reveal that in certain situations, adaptive isostatic phases appear, and these lead to remarkable physical properties over a finite interval in composition defining an intermediate phase. These features do not seem to be restricted to glasses, and, instead, connect to universal features found in disordered networks.


Introduction
Because of the lack of translational periodicity that prevents from having longrange order for the atomic arrangement and a structure made of a repetition in three dimensions of a well-defined unit cell, the description of non-crystalline ABSTRACT We present the most recent developments and applications of Rigidity Theory in atomic-scale disordered materials such as glasses, liquids or amorphous solids. While the framework of rigidity, elastic phase transitions and topological constraint counting has been extensively used in glass science, its potentially wide applicability has remained difficult because of some inherent limitations of the theory. Here, we review promising methods that have been recently introduced, and which allow to overcome such limitations, and to extend these methods to the case of densified liquids or partially ordered solids, and also to the non-mean-field case. We introduce and describe a several methods based on an explicit account of temperature effects. In addition, the calculation of bond length and bond angle standard deviations is presented, accessed from molecular simulations, which allows to estimate atomic stretching and bending topological constraints as a function of pressure or temperature. These new methods are illustrated through a set of applications ranging from several archetypal glass-forming liquids (SiO 2 , Ge x Se 1−x , …) to semiconductors (Si-C:H), phase change materials (Ge-Te) and cement, while also emphasizing the role of isostatic composition for the design of stable glasses with reduced ageing. Results furthermore reveal that in certain situations, adaptive isostatic phases appear, and these lead to remarkable physical properties over a finite interval in composition defining an intermediate phase. These features do not seem to be restricted to glasses, and, instead, connect to universal features found in disordered networks. Glasses  solids has remained a challenging field of inquiry in materials science. It has been impossible to describe their electronic properties in terms of the powerful techniques offered by crystalline band theory such as Brillouin zones, k-space and Bloch functions. In addition, this absence of a repeating unit cell has often limited to short-range order the detailed structural information gathered from standard experimental techniques. Still, diffraction methods provide a partial characterization of the intermediate range order [1][2][3][4][5][6] with a particular emphasis on the low wave-vector region revealing a first sharp diffraction peak (FSDP) in the structure factor and particular correlations in reciprocal space [7][8][9][10][11] that can be even related to a species dependence using isotopic substitution [12][13][14]. However, there is only a limited number of cases where the structural groupings, essentially ring structures, of the glass structure can be unambiguously analysed and characterized from their explicit spectroscopic signature using Raman or Nuclear Magnetic Resonance (NMR) spectroscopy [15][16][17][18].
Despite these inherent lack of a full structural information, the possibility to tune by appropriate alloying of the composition of non-crystalline solids has led to a certain number of important recent applications among which scratch-resistant glasses with increased hardness [19], amorphous films for electrical and optical data storage [20], ultra-transparent optical fibres designed at the desired wavelength [21], solid-state electrolytes [22], metallic glasses for bio-implants or military applications [23], etc.
From a more fundamental and thermodynamic viewpoint, non-crystalline solids are in most of the situations out of equilibrium systems so that concepts and techniques from statistical mechanics cannot be applied in a straightforward fashion. The breakdown of the relationship between spontaneous fluctuations and dissipation that indicates the absence of ergodicity will, thus, lead in timedependent quantities and ageing phenomena that can be characterized from appropriate theoretical frameworks [24]. This has also some important consequences for the experimental investigation of amorphous solids, and the reproducibility of glass measurements has been regularly questioned because materials properties are found to vary quite substantially with the thermal history of the melt [25], the waiting time before the experiment/measurement [26], the sample homogeneity [27], the preparation method [28], the dryness [29], etc.
Given these intrinsic drawbacks, atomistic simulations have proved to be extremely valuable for understanding structural and dynamic properties including ageing phenomena, as well as for checking the validity of microscopic models, empirical approximations or theoretical assumptions. However, while the available computing power has increased exponentially over the past decades, a large enough computing power for direct molecular dynamics (MD) simulations of glass on a realistic laboratory time scale should not be available in a foreseeable future. Furthermore, in the search for property optimization or new functionalities ruled by changes in composition, a huge amount of time has to be spent on detecting trends in physical properties in the compositional phase diagram, and this also burdens the possibility to use MD for this purpose. Among alternatives, there is another path forward, albeit carrying its own limitations, Rigidity Theory [30], which focuses on the key microscopic physics governing the thermal, mechanical and rheological properties of glass and amorphous materials, while filtering out unnecessary details which ultimately do not affect the overall properties.
The concept of rigidity in disordered solids such as glasses, amorphous networks or sphere packing traces back to the early work of Maxwell on the stability of trusses and macroscopic structures such as bridges, and to the introduction of mechanical constraints by Lagrange. These ideas and results were then extended to atomic networks by Phillips [31] who highlighted the notion of mechanical isostaticity as promoting glass-forming tendency of covalent alloys. Specifically, it was recognized that so-called 'good' glass formers usually form at an optimal network connectivity, or mean coordination number r =r c satisfying the Maxwell stability criterion of isostatic structures, i.e. n c = n d , where n c is the count of atomic constraints per atom, and n d is the dimensionality in which the network is embedded, usually n d = 3 in 3D.
Aspects of rigidity have a large interest in the field of condensed matter physics, elastic phase transitions, glass science and, more generally, in materials science, and we will address here the important theoretical progresses that have been achieved recently such as the extension of constraints to the liquid state using an appropriate temperature dependence or an enumeration from MD simulations. The latter allows establishing the rigidity status of a material under various thermodynamic conditions (composition, pressure, temperature) or excitations (mechanical, irradiation). This has also led to promising applications such as glasses with superior toughness designed for cell phones, and to the discovery of stable isostatic compositions with remarkable properties, the search of thermal stability being a crucial issue in e.g. high-tech glasses or optoelectronic applications leading to e.g. phase change recording.

Phillips-Thorpe formulation
We start by briefly reviewing the basic notions of Rigidity Theory following the Phillips-Thorpe formulation [31,32], i.e. how a disordered network can be viewed as an assembly of mechanical constraints that change with composition (Figure 1), and how an isostatic condition can be achieved that coincides with a flexible to rigid transition [33] having as control parameter for the phase transition the network mean coordination number r, and as order parameter the fraction of zero-frequency (floppy) modes defined by :  = 3 − n c (see below).
In covalent amorphous networks, the most important interactions are usually near-neighbour bond-stretching (BS) and next-near-neighbour bond-bending (BB) forces which constrain the structure at the atomic level [34] and lead to e.g. well-defined bond distances and bond angles between species. Under this assumption, the number of mechanical constraints n c per atom can be exactly computed in a mean-field way at zero temperature on a fully connected network, and is given by: Figure 1. A graphical summary of the rigidity approach: An atomic structure (a) (here a 3000 atom bond depleted amorphous silicon system whose atoms interact via harmonic Bs and bending (BB) interactions, and which can be viewed as a bar network constrained by bonds and angles. An analysis of the dynamic matrix in such simple bond-depleted random networks leads to the fraction of floppy modes  (b, red triangles) that follows essentially a Maxwell constraint count (solid line, equation (3)) which vanishes at the mean-field critical mean coordination number where n r is the concentration of species being r-fold coordinated. The contribution of the two terms in the numerator is obvious because each bond is shared by two neighbours, and one has r/2 BS constraints for a r-fold atom. For the BB (angular) constraint count, it has to be remarked that a twofold atom involves only one angle, and each additional bond will need the definition of two more angles in 3D. For particular species (e.g. onefold terminal atoms or broken constraints [35]), a special count can be achieved. By defining the network mean coordination number r of the network by: Equation (1) can be reduced to the simple equation: n c =̄r 2 + 2r − 3 which leads to the density of floppy modes: Applying the Maxwell stability criterion, isostatic glasses (n c = 3) are expected to be found at the average coordination number [32] of r c = 2.40 in 3D ( Figure  1(b)). For instance, in Group IV chalcogenides (Ge x S 1 − x , Si x Se 1 − x ) with r = 2 + 2x, such a condition is met for e.g. GeS 4 (x = 20%) which is different from the stoichiometric composition (GeS 2 ).
The problem can be posed an alternative way, and the physical origin of this stability criterion can be established [33] from the vibrational analysis of bonddepleted random networks constrained by BB and BS interactions using a harmonic Keating potential (Figure 1(a)). For such systems, the number of zero-frequency (floppy) modes  (i.e. the eigenmodes of the dynamical matrix) is vanishing when rigidity percolates in the network, and this happens for r = 2.38 (Figure 1(b)). The Maxwell condition n c = n d therefore defines an elastic phase transition, above which redundant constraints produce internally stressed networks, identified with a stressed-rigid phase [32,33] which lead, among other features, to a power-law variation for the elastic constants. For n c < n d however, floppy modes do exist, and their presence defines a flexible phase in which local deformations with a low cost in energy (typically 5 meV, [36]) are possible.
There has been a large body of experimental investigations related to these results and conclusions in the literature, and various signatures for rigidity percolation have been detected from Raman scattering [37], calorimetry [37], viscosity measurements [38], Brillouin scattering [39], resistivity [22,40], stress relaxation [41] and many other properties [30].

Limitations
There are, however, a certain number of obvious limitations for the theory, and especially, for Equation (1). First, the approach is, by essence, a topology-based theory which relies essentially on coordination numbers (r in Equation (1)) and on the connectivity of the network. In this respect, there is no distinction made between isovalent systems (e.g. Group IV chalcogenides) or between systems belonging to the same family of modified glasses (e.g. alkali silicates), although there are obvious differences among such systems when experimental data are carefully put in perspective. One major drawback is, indeed, the assumption that r is determined from the 8- rule where  is the number of s and p (outer shell) electrons. This physical situation is fulfilled in covalent glasses such as chalcogenides for which, in addition, the dominant BS and BB interactions are strong as compared to other more weaker interactions (e.g. Van der Waals). These BS and BB interactions involve typical frequencies of 19-36 meV [36], as detected from the vibrational density of states. In this case, equation (1) holds and the prediction  = 0 is recovered experimentally. However, this assumption does not apply when non-directional (e.g. Coulombic) forces are present such as in alkali or alkaline earth containing glasses [39], or when the system involves partially metallic bonding (Ge-Te) that is characterized by a breakdown of the 8- rule and the sp 3 hybridized geometry of germanium [42]. This prevents from estimating n c and r in a simple fashion. Secondly, the enumeration of bonding constraints is performed over a fully connected network, i.e. at T = 0 K when neither bonds nor constraints can be broken by thermal activation. This physical situation may be valid as long as the viscosity η is very large at T < T g and when the bond lifetime is large, but Equation (1) can certainly not be used in a high temperature liquid. Similarly, in Equation (1), it is also assumed that the average constraint count is performed over all the atoms of the network which implicitly neglects the possibility of phase separation or a structural decoupling from the backbone, a situation encountered in Group V chalcogenides [43] or in sulphur-rich glasses [44].
Finally, it must be emphasized that the construction leading to the constraint estimate of Equation (1) is of mean-field type because ensemble averaged quantities are used and fluctuations are neglected, i.e. one assumes that all atoms of a given type have the same number of constraints and the same coordination number. However, as in ordinary phase transitions, large fluctuations in constraints or in the order parameter () may be expected close to the critical point at r =r c . It has been suggested [45] that this non-mean field scenario manifests by the onset of an isostatic phase (see below) in which fluctuations in coordination may be important, as well as network adaptation in order to lower the stress induced by the increasing cross-link density.

Temperature dependent rigidity theory
Mauro and collaborators have proposed a method [46] to account for temperature effects in Rigidity Theory using a two-state thermodynamic function q(T). This function quantifies the number of rigid constraints as a function of temperature and subsequently modifies Equation (1) which becomes: where q r (T) and q r (T) are step functions associated with BS and BB interactions of a r-coordinated atom so that n c now explicitly depends on temperature. This function q(T) has two obvious limits because all relevant constraints can be either intact at low temperature (q(T) = 1, as in the initial Phillips-Thorpe formulation) or entirely broken at high temperature (q(∞) = 0). At a finite temperature, however, only a fraction of these constraints can become rigid once their associated energy is less than k B T. Different forms can be proposed for q(T) based either on an energy landscape approach [47] or on an equilibrium reaction' broken ↔ intact between broken and intact constraints which leads to: where ΔG * is the Gibbs energy characterizing the equilibrium constant K e (T), and this leads to a functional form for q(T): which is found to fit directly estimated q(T) functions from MD simulations (see below, [48,49]).
A certain number of thermal and relaxation properties of network glass-forming liquids can be determined within this new framework, and analytical expressions for fragility and glass transition temperature [46,47,[50][51][52], glass hardness [53], and heat capacity [54] can be obtained under the assumption of a simple step-like (Heaviside) function for q(T) = 1 − Y(T − T α ) with onset temperatures T α for BS and BB constraints which act as parameters for the theory.
For instance, a prediction of a T g variation is obtained by assuming that the Adam-Gibbs model [55] for viscosity, ln ∝ 1∕TS c , holds and that the configurational entropy S c is proportional to the fraction  of floppy modes [56]. By furthermore stating that T g is a reference temperature at which η(T g (x), x)=10 12 Pa s for any composition, one can write: and T g (x) can be determined for a given composition x from a reference compound having a composition x R and a glass transition T g (x R ), knowing the number of temperature-dependent flexible modes  for compositions x and x c from equation (4), and the behaviour of the step functions q(T).
Using the expression for S c (T g (x), x) in Equation (7), and the definition of the fragility index: which describes the degree of change in transport properties close to the glass transition, one can furthermore extract an expression for  as a function of composition: Typical applications ( Figure 2) for the prediction of the glass transition temperature and fragility concern simple chalcogenides [46], borates or borosilicates [54], phosphates [52] or borophosphate glasses [50]. Equation (9) usually leads to a good reproduction of fragility data with composition ( Figure 2(a)), but requires a certain number of onset temperatures T α that can be estimated from model assumptions.

MD simulations and standard deviations
Alternatives routes can be derived for the rigidity properties of liquids, and especially, for liquids under pressure. From accurate structural MD models reproducing e.g. the total structure factor S(k), the pair correlation function g(r) and other experimental observables (viscosity, diffusion [57]), one enumerates mechanical constraints from the atomic-scale time-dependent trajectories. The global strategy can be sketched in Figure 3(a).
To determine the number of BS interactions, one focuses on neighbour distribution functions [58] around a given atom i, the global sum of all such functions yielding an i-centred pair correlation function g i (r) whose integration up to the first minimum gives the coordination numbers r i , and hence the corresponding number of BS constraints n BS c = r i ∕2. Figure 4(a) shows such an application [59] for the As 2 Se 3 glass showing that three neighbours contribute to the first peak of g As (r), very well separated from the second shell of neighbours, and n BS c = 1.5 for the As atom in this case.
To determine BB constraints, one uses partial bond angle distributions (PBADs) P(θ ij ) which split the usual bond angle distribution into partial contributions defined by a central atom 0 and the N first neighbours which define individual N(N − 1)/2 possible triplets or angles i 0 j (i = 1 … N − 1, j = 2 … N), i.e. 102, 103, 203, etc. and these are followed with time ( Figure 3(b)) over the MD trajectory. The standard deviation ij of each distribution individual P(θ ij ) gives a quantitative estimate of the angular excursion around a mean angular value, and provides an indication of the BB strength [48,49,58]. Small values for ij correspond to an intact BB constraint which maintains a rigid angle at a fixed value (Figure 3(a)), whereas large ij correspond to a BB weakness giving rise to an ineffective or broken constraint.
To determine a meaningful limit between intact and broken BB constraints a cut-off value min is needed. The analysis over the whole system of the distribution of individual P( ij )s with their standard deviations leads to a distribution of standard deviations ij (Figure 4(c)). Interestingly, it is seen that at low temperature all constraints can be considered as intact [46] because f ( ) is limited to a sharp peak found at low ij ≃ 7°. At high temperature, all constraints can be considered as broken because of thermal activation and f ( ) displays, indeed, a broad pattern centred at a large value (≃ 25°) indicative of a full breakdown of angular constraints. For intermediate temperatures, a bimodal distribution sets on which reveals a population of broken and intact constraints whose statistics follows a function q(T) (inset of Figure 3(c)) such as the one of Equation (6) or an alternative one proposed by Gupta and Mauro [46], and having the obvious limits q(0) = 1 and q(∞) = 0. A natural choice for min corresponds to the minimum of f ( ) (red curve in Figure 4(c)) that is weakly temperature dependent, and in the amorphous phase one can safely choose min ≃ 12-15°.
For e.g. oxide systems, the PBADs relative to the Group IV (Si, Ge) atom have a low standard deviation σ θ , of the order of 10-20° when the first four neighbours are considered. One finds e.g. ≃ 7° for the PBAD 102 of GeO 2 , which is substantially smaller as compared to those computed from other distributions (105, 106, etc.) which have ≃ 40°. This indicates that n BB c = 5 for Ge/Si atoms, i.e. there are five angular constraints associated with the tetrahedral angles (109°). The count from standard deviations (Figure 4(b)) leads to six constrained angles but only five of them are independent. The red curve defines a boundary σ min between broken and intact constraints, estimated to be about ≃ 15° at low temperature, and serves as a cut-off for the determination of constraints. The inset shows the fraction q(T) of intact BO constraints as a function of temperature. The solid curve is a fit using equation (6). Alternative functionals for q(T) can be used [46].

Pressure-temperature effect
A first central outcome of the coupled MD and constraint counting method is the possibility to extend Rigidity Theory to densified liquids. The introduction of the step functions has led to Equation (4) offering an increased applicability to temperature effects in liquids. The use of standard deviations permits to push the approach to the next level because a pressure-dependent function n c (x, T, P) can now be derived, and this is out of reach from simplified structural models [46,47] given the non-linear behaviour of coordination changes with pressure [66]. Furthermore, this extension of topological constraint count to densified liquids permits to investigate certain thermodynamic (T, P) conditions which lead to anomalous behaviours in structural, dynamical or thermodynamic properties [67,68]. Such anomalies furthermore display similarities with those observed for isostatic compositions at low temperature and ambient pressure [37,44,69].

Adaptive liquids
The first outcome of such applications is the detection that liquids adapt under the pressure increase. In fact, in a densified silicate liquid, the increase in connectivity arising from the tetrahedral (r Si = 4) to octahedral (r Si = 6) coordination change with pressure [57] leads to an increase in the number of BS constraints ( Figure  5(a)). In order to accommodate this additional stress, some of the BB constraints soften and n BB c decreases at low pressure ( Figure 5(a)). The angular adaptation reduces the network stress but, due to the continuously increasing pressure and coordination change, this holds only up to a certain point (≃ 12 GPa) beyond which angles will stiffen and BB constraints will be restored.
Such trends are correlated with a certain number of transport anomalies [57,70], and in the region where angular adaptation occurs, viscosity and diffusivities are found to be minimum (D = D min ) and maximum (D max ), respectively, and corresponding activation energies E A are also found to minimize. This indicates the central role played by the softening of the bending constraints which lead to an increased ease of relaxation given the minimum in E A .

Water-like anomalies
These identified anomalies detected in diffusivities and viscosity [70] have actually a close connection with water-like anomalies detected in densified tetrahedral systems (SiO 2 , GeO 2 , H 2 O, BeF 2 ). For such liquids, it has been found [71][72][73] that the loci D max and D min of the diffusivity maxima and minima are correlated with thermal (temperature of maximum density line), pair correlation entropy and structural anomalies. For the latter, the (T,V or ρ) evolution of local structural order parameters [74] quantifying both translational order (the tendency of pairs of molecules to be separated by a preferential distance) and orientational order (the tendency of a molecule and its nearest neighbours to adopt preferential orientations) shows a direct correlation with the other anomalies.
In liquid densified GeO 2 , an analysis [75] of the constraints with temperature and system volume V shows that the locus of the diffusivity anomaly D max takes place in a thermodynamic region (T, V) where angles soften, reduce n BB c and lead to a global decrease (minimum) in the total n c (V, T) ( Figure 5(b)). Here, the mechanism is identical to the one found for other densified liquids [70,76], i.e. in order to accommodate the stress induced by the density-driven increase in the oxygen and germanium coordination numbers which increase n BS c , larger bond angle excursions are experienced and reduce the number of BB interactions. This provides evidence that such an angular adaptation drives the observed anomalies in diffusivities by reducing the liquid rigidity.

Quantifying the intermediate phase
The angular adapation detected in the liquid phase ( Figure 5) leads to particular relaxation phenomena which induce anomalous glass transitions, i.e. weak thermal changes at the glass transition are obtained during a cooling/heating cycle (Figure 6(a)), and this is also detected in a variety of experiments on oxide and chalcogenide glasses [43,44,[77][78][79][80]. When MD numerical cycles are performed across the glass transition from the high temperature liquid, one finds, indeed, a hysteresis between the cooling and heating curve in accordance with the salient experimental phenomenology. This behaviour simply reflects the off-equilibrium nature of glasses which slowly relax at T < T g , and induce a decrease in volume or enthalpy. However, for selected thermodynamic conditions (i.e. pressure), and at a fixed cooling/heating rate the hysteresis curves become minuscule (Figure 6(a)), and the cooling/heating curves nearly overlap. This defines a reversibility window (RW) that can be quantified by the hysteresis area A H of the enthalpy variation or the volume hysteresis (A V , inset of Figure 6(a)).
Such thermal anomalies are directly linked to the isostatic nature of the glass-forming liquid, as detected from the constraint count using radial and angular standard deviations ( Figure 5(a). At zero pressure, the investigated liquid (sodium silicate) is flexible (n c < 3, Figure 6(b)), in agreement with experimental results [39]. But as the pressure is increased, the total number of constraints reaches for P = P c(1) ≃ 3 GPa a plateau value of n c = 3 that is characteristic of isostatic systems, and which defines an intermediate phase (IP) between the flexible (n c < 3) and the stressed rigid phase (n c > 3). Upon further compression, this nearly constant value of n c which is driven by angular adaptation ( Figure 5) holds up to a threshold pressure of P c(2) = 12 GPa beyond which an important growth in n c (P) takes place, and the locus P c(2) is identified with a stress transition. When the trend in n c (P) is compared to the one of A V and A H , a direct correlation emerges and one arrives to the obvious conclusion that isostaticity must drive the reversible tendency of the glass transition given that both floppy (low frequency) mode and stress relaxation are minimized for an optimally constrained system.
One furthermore notices that a constraint count at low temperature (300 K, [81]) leads to a single threshold n c = 3 close to the centroid of the pressure window (P ≃ 8 GPa), and at zero pressure one has n c = 2.68 (2.56 from a mean-field count [39]). This behaviour is consistent with recent results on the rigidity and RWs of spring networks [82]. For such systems, temperature and weak non-covalent (van der Waals) interactions drive the onset of a RW, the main outcome being that the temperature considerably affects the way an amorphous network becomes rigid under a coordination number increase, and the existence of an isostatic RW sets on only for T ≥ α where the parameter α characterizes the relative strength of the weak forces. This might indicate that the rigidity transitions become mean field at low temperature as depicted in Figure 6(b). However, RWs result mostly from the particular relaxation phenomena at play during the glass transition interval and may, therefore, be sensitive to the timescale of the MD simulations. In such simulations and in contrast with the spring network results, coordination fluctuations are found to be small [57], given the weak abundance of fivefold Si atoms (10-20%), and fluctuations are essentially found in angular constraints [83] which show a non-random distribution (inset of Figure 6(b)) and signal a breakdown of the mean-field estimate of Equation (1).

Insight from alternative theoretical approaches
Alternative scenarii regarding the IP have been reported, and these are not all consistent with the findings of e.g. Figure 6 because different systems, toy model ingredients and/or thermodynamic conditions are used to assess particular conclusions. Broadly speaking, these models can be put into two main categories. Some studies emphasize the central role of fluctuations [45,[84][85][86] in the emergence of a double threshold/transition that defines an IP between the flexible and the stressed rigid phase, and these features bear similarities with ordinary phase transitions given that fluctuations in coordination number will induce corresponding fluctuations in the order parameter (). Alternatively, mean-field aspects of jamming have been considered, and here fluctuations in coordinations are thought to be limited, but atoms are coupled spatially via elasticity [87,88] and can organise locally into distinct configurations that may promote an IP. Physical ingredients (structure, dynamics, interactions, …) leading to the IP are still being debated in the literature, and its experimental characterization can be achieved from a variety of probes (thermal, calorimetric, optical, structural). In this respect, realistic MD simulations on glasses and supercooled liquids may provide an additional and insightful information [89][90][91], and it is interesting to note that thermal effects initially absent in theories of the intermediate phase [45,84,85,88], are now taken into account explicitly [86,87] which should permit to better understand how fluctuations, constraints, elasticity behave with temperature, and how such quantities are coupled to glassy relaxation.

The intermediate phase: an experimental controversy
Although the signature of a double threshold has been established from a variety of models -we here disregard those having a much weaker theoretical basis, the experimental observation and characterization of the IP continues to be a subject of strong controversy due, in part, to its difficult reproducibility.
One of the most obvious signatures of the IP arises from a modulated differential scanning calorimetry measurement [37] which permits one to measure a non-reversing heat enthalpy, ΔH nr , embedding a large part of the kinetic events associated with the slowing down of the relaxation close to the glass transition. In the IP, it has been found that ΔH nr nearly vanishes (Figure 7(a)) over a finite compositional interval, and, under certain sample/relaxation preparation conditions, the boundaries defining this interval can become very sharp. A certain number of experimental studies have challenged the existence of the IP [92][93][94], suggesting that the vanishing of ΔH nr might be the result of an experimental artefact. However, in such studies, alternative signatures of the IP and related anomalous behaviours of physical properties are not discussed, and usually overlooked. Figure 7 displays a survey of some of such anomalies occurring in the IP for different families of glasses : (1 − x)TeO 2 − xV 2 O 5 [95], Ge x Se 1 − x [77,96], (1 − x) AgPO 3 − xAgI [22] and (1 − x)GeO 2 − xNa 2 O [80]. In all systems, the softening of the rigid network is achieved by a change in composition. When the atomic sizes are comparable, it has been realized that because of the absence of stress [85], networks will pack together more easily in the IP as manifested by a minimum in the molar volume (Figures 7(a)-(c)). The stress-free nature of the IP has been detected from pressure experiments [96] showing the vanishing of the threshold pressure (Figure 7(b)) prior to a pressure-induced Raman peak shift which is an indication of residual stresses as in crystals. In ionic conductors (Figure 7(d)), the exponential increase in ionic conduction appears only in the flexible phase when the network can be more easily deformed at a local level because of the presence of floppy modes [22]. However, in the IP, an intermediate conductive regime sets in, which also shows anomalous behaviour for a typical jump distance [97], and dc permittivity ɛ(0) (Figure 7(d)), right axis) [98].
The search of a possible structural signature or structural origin of the IP also continues to be the subject of an intense debate. High-energy X-ray synchrotron radiation and X-ray absorption fine structure measurements have indicated that typical structural parameters such as peak heights in atomic pair distribution functions g(r) or the FSDP in the structure factor S(k) display a smooth trend with composition as glasses evolve from the flexible to the stressed rigid phase [99]. These results seem to be consistent with conclusions drawn from a compilation of diffraction data [11] but are contradicted by similar measurements [90,91] coupled to MD simulations. Similarly, specific Raman mode characteristics do exhibit changes and/or thresholds across the two boundaries of the IP, particularly Raman mode frequencies [30,37,44,79,80,96] which are sensitive to optical elasticity. However, changes in local structural order and topology probed by NMR do not seem to lead to particular anomalies across the IP, and the band assignment in relation to structural models is also subject to some controversy [92,100,101]. Simulations indicate that only certain correlations in real space (e.g. homopolar bondings) and reciprocal space (Faber-Ziman) may be sensitive to the isostatic nature, as revealed from First Principles MD [102]. This means that e.g. the total structure factor accessed from neutron or X-ray scattering may not represent a structural probe for the IP [11,99,103]. In this respect, neutron diffraction studies using the isotopic substitution method to provide the partial contributions to both structure factors and pair correlation functions [104] across the required range of compositions are highly desirable. Preliminary work in this direction has been achieved from anomalous X-ray scattering coupled to Reverse Monte Carlo simulations [105,106], and results seem, indeed, to indicate that the elastic nature (flexible, intermediate, stressed) of the atomic structure leads to particular signatures in partial correlations only. Conclusions are, however, strongly material dependent, and one should keep in mind that there are only a very limited number of systems for which the mean-field constraint counting (Equation (1)) can be achieved with confidence. It should, therefore, be only considered as a general and initial guide prior to experimental and numerical investigations.

New directions
The discovery and the experimental characterization of the IP have also led to some more unexpected developments in the field, due, in part, to the difficulty in reproducing the measurements of ΔH nr which need to be corrected [107] in order to avoid spurious effects arising from the frequency dependence of the specific heat [108]. However, once corrected, it has been furthermore discovered that this quantity is highly sensitive (Figure 8) to impurities, inhomogeneities and the relaxation status of the glass [77]. Careful sample preparation (e.g. free of water) has therefore become an important prerequisite prior to the detection of the IP [109,110] and its characterization.
The second interesting outcome is the recognition that melt homogeneity greatly affects ΔH nr , and, more generally, measured properties of network glasses [27,44,109]. This issue is often overlooked in the literature but also appears to be crucial as percolative elastic phase transition [84] is to be considered. It has been found, indeed, that the dynamics of melt homogenization during glass melting is not only strongly composition dependent but also limited by the slowest diffusing intermediate phase species [111,112], i.e. for compositions displaying the lowest possible melt fragility.

Rigidity phase diagram
One interesting perspective of MD-based constraint counting is the possibility to identify rigidity transitions for materials in which the 8- and Equation (1) cannot be used straightforwardedly. This is the case for telluride materials used for phase change recording (e.g. Ge-Sb-Te [113]) which display semi-metallic bonding, and lead to more complicated local structures, as highlighted both from experiment [114] and simulations [20,28,42,115]. Using a precise enumeration for constraints arising from BS and BB interactions via First Principles MD trajectories [58], amorphous phase-change materials can be investigated in a similar fashion to network glasses (e.g. Ge-Se). Results reveal a phenomenology that is close to lighter chalcogenides, albeit mean-field rigidity percolation thresholds of Ge-Sb-Te are found to be shifted in composition with respect to Ge-Sb-Se [116] due to the increased complexity of the chemical bonding. The phase diagram of Ge-Sb-Te system can, indeed, be separated into a flexible chalcogen-rich phase (here Te) and a (Sb,Ge)-rich phase that is stressed rigid. A rigidity transition (isostatic) line is expected close to the GeTe 4 -SbTe 4 join, and the most commonly compositions used Ge-Sb-Te phase-change materials belong to the stressed rigid phase in which glass-forming tendency is usually rather weak [30]. , low temperature annealed (filled circles) and artificially water contamined (red curve). The increase in ΔH nr outside of the window (arrows) is due to floppy mode relaxation in flexible phase and to stress relaxation in the stressed rigid phase. (b) evolution of ΔH nr in Ge-se glasses for different melt homogenization conditions at 950° prior to a quench at low temperature: two days [37], four days [142] and seven days [77,109,110], using for the latter a Raman profiling probe to check for homogeneity at the micron scale.

Quantifying tetrahedra
The neat estimate of constraints from MD simulations [117,118] also permits to address an important issue in the field of phase change materials, i.e. the detection of Ge tetrahedral units and the measurement of their population with e.g. composition. Early studies of tellurides have, indeed, emphasized the presence of different structural geometries [114], Ge tetrahedra or octahedra, the latter being found in crystalline Ge 2 Sb 2 Te 5 whose structure is made of a rhombohedrally distorted rock-salt lattice containing two embedded Te-and Ge/Sb sublattices with vacancies [113]. Tetrahedra have been found in the amorphous phase [119,120] and have been proposed to contribute to the phase switching phenomena because of the participation of different electronic orbitals (as compared to octahedra), the sp 3 hybridized tetrahedra typical of covalent semi-conductors, involving only s electrons, whereas p-bonded octahedral bonding may facilitate conduction.
In recent numerical models accessed from First Principles MD, the fraction of tetrahedra can be estimated qualitatively by different ways, i.e. from a selection of the typical bond angles [115] but there must be an obvious overlap between the two populations (centred at 90° and 109°) in the total Ge-centred bond angle distribution. Alternatively, a bond length argument can be used [121] which reveals that the fourth neighbour distribution of a central Ge atom is bimodal, and can be associated with either three-and tetrahedrally four-coordinated germanium. A local order orientational order parameter distribution is introduced [122] for this purpose but the fraction of tetrahedra will dependent on the integration boundaries of the distribution.
The detection and quantification of tetrahedral germanium can be made from angular standard deviations which are associated with BB interactions [123,124]. The method offers, indeed, the advantage to focus on angular excursions rather than working directly on angles [115]. Angles are followed individually during the simulation from the N(N − 1)/2 possible triplets i 0 j defined by a set of N first neighbours around a Group IV atoms. If the number of low standard deviations around such atoms is six, a tetrahedron is identified because this geometry is defined by six rigid angles that give rise to corresponding low standard deviations (Figure 3(c)), and also leads to associated angles that are close to 109°. Averages over the entire system then leads to a precise fraction of Ge or Si tetrahedra in amorphous phase change materials such as Ge-Te and Ge-Si-Te, which can be compared to an experimental measure from Mössbauer spectroscopy of 119 Sn substituted (1% wt) tellurides ( Table 1). The latter technique allows to probe the local geometry, tetrahedral vs. octahedral, and also permits determining the fraction of Sn tetrahedra, but not directly the fraction of Ge tetrahedra if one does not assume a perfect chemical substitution between Sn and Ge. It yields a higher correlation of structure, not just bonds as accessed from e.g. neutron diffraction, but how these bonds are locally configured because the nuclear hyperfine signals are directly transferable [125] to those in reference compounds which contain either a purely tetrahedral (crystalline Si) or a purely octahedral structure (crystalline SnTe).
In binary Ge-Te and Si-Te glasses [123], the chalcogen-rich region (14-15%) is dominated by tetrahedral-rich structures (Table 1) but dramatic changes are acknowledged in the former system once the Ge content reaches 18%, and Ge atoms are then predominantly found in an octahedral geometry (e.g. 58.4% in Ge 18 Te 82 ). This situation contrasts with the Si-Te binary which contains a large fraction of SiTe 4/2 tetrahedra reminiscent of the local structure of the crystalline Si 2 Te 3 polymorph [126] which is 100% tetrahedral. In the ternary Ge-Si-Te system [124], the calculated tetrahedral fraction is found to be somewhat lower (33.7%) when compared to its isovalent compound (Ge 20 Te 80 , 54.6%) but is very close to the one determined experimentally from Mössbauer spectroscopy (40.8%). The detail of the numerical analysis permits to determine the different contributions, and it is found that Si sites dominate the population of tetrahedra (42.74% for Si against 24.6% for Ge).
Using the same analysis tools, the remaining non-tetrahedral Ge geometry can also be characterized in Ge-Te systems, and in e.g. Ge 20 Te 80 only three angles are found to have small angular excursions. These involve the first three neighbours of Ge at distances of about 2.69 Å [118,123], and with angles found at (≃ 98°). This defines a pyramid with a triangular basis having the Te-Te bonds as edges, and a Ge at the remaining vertex, similar to the pyramidal geometry found in As 2 Se 3 for which three rigid angles are also obtained [59]. However, in contrast with the selenide compound, the first coordination shell contains a fourth neighbour that is found at a slightly larger distance (2.96 Å) but associated angles have their BB constraints broken which indicates, on the overall, a much softer geometry as compared to the one encountered in As-Se or Ge-Se.

Dielectric materials
Rigidity Theory also finds recent applications in the field of insulating materials having lower dielectric constants (i.e. low-k). Such materials are used in order to minimize capacitive related power losses in integrated circuits.
A study on amorphous SiC:H with changing composition [127] shows that a variety of properties are deeply linked to the change in network connectivity accessed from NMR. The determination of the average mean coordination number r) from NMR differs from a numerical estimate [128] which assumes the 8- rule and only sp 3 geometries for both C and Si atoms. However, the same tendency is recovered in the experimental results [127], that is, r steadily decreases as the Table 1. experimentally measured (Mössbauer) tetrahedral fraction in tin substituted amorphous tellurides, and constraint-based calculated fraction (in %) of tetrahedral sites in amorphous Ge-Te and si-Te [123] and Ge-si-Te [124] systems with changing Group iv atom content. hydrogen content increases. Such a controlled incorporation of terminal hydrogen leads to large changes in mechanical (Young's modulus, hardness, intrinsic stress), transport (thermal conductivity, resistivity), optical (refractive index) and materials properties (mass density and porosity), and also point out that the isostatic condition recovered at the critical connectivity r c exists in this material system and places limitations (Figure 9, left) on the range of properties that can be achieved for future low-k amorphous SiC:H materials. For instance, porosity of such systems is found to be deeply affected (i.e. divided by a factor 3) once the network becomes stressed rigid. This might be due to the larger bond density and the vanishing of flexible modes which promote mass transport.

Cement
The nanoscale structure of cement can be also investigated from a combination of MD simulations and Rigidity Theory. The structure of a Calcium-silicate-hydrate (C-S-H), which is the binding phase of cement, also appears to be a challenging material for the extension of the constraint counting using Equation (1) given its partially order showing both crystalline and amorphous features and a quasi-layered structure that is not homogeneous. There is the possibility to have broken constraints [35] and coordination numbers that do not necessarily follow the 8- rule, while some undissociated water molecules are not part of the network. From realistic models of cement [129] that accurately reproduce measured structure factors and pair correlation functions [130], a precise enumeration of topological constraints [65] shows that the rigidity of the C-S-H network can be increased by decreasing the Ca/Si molar ratio (Figure 9, right). The detail of the  [65]. A flexible to rigid threshold is expected for ca/si ≃ 1.5. The inset shows the respective contribution due to Bs and BB constraints. estimate of BS and BB constraints shows similarities with other glass-forming materials, i.e. silicon atoms have n BS c = 2 and n BB c = 5 as in silica [35,61] or silicates [39], and calcium atoms are found to have only stretching constraints (n BB c = 4.8 on average) as in the soda-lime silicate system [60] but with two populations of cations depending on their presence between (n c = 4) or inside (n c = 6) the layers. The analysis finally reveals that the average composition of ordinary Portland cement [131] (C-S-H with Ca/Si ratio of 1.71) has a flexible network, and this floppy character is mostly due to the hydration of the network which disrupts the connected rigid (Si, O) based layered structure. Within the same analysis, an isostatic network can be expected for an approximate Ca/Si ratio of 1.5 which may lead to the achievement of optimally constrained cements displaying the salient properties of IP glass compositions such as maximal fracture toughness [132] and weak ageing phenomena [133].
Having the enumeration of constraints in hand, the hardness H v [53] of cement can be predicted from temperature-dependent Rigidity Theory [46], as a function of the C-S-H composition. It is found [64] that H v does not depend on the total number of constraints n c , but only on the weaker ones which are related to the bending motion. This contrasts the role played by BS and BB constraints which contribute differently to hardness, but also indicates that the assumption of BS and BB contributing in a similar way to the network stiffening cannot resolve mechanical properties, as already revealed from their different associated activation energies [36].

Glasses with reduced ageing
Isostatic network glasses are found to display a significantly reduced tendency to ageing and this has been detected for a certain number of systems. In calorimetric studies of As-Se glasses [134], there is a weak evolution of the enthalpy with waiting time t w for glasses belonging to the RW ( Figure 10). Although this result has been challenged [93,135], additional analysis [136] have shown that such contrasting conclusions might result from a nanoscale-phase separation induced by light exposure.
One should note that the experimental protocol used in [103,134,136] does not follow the one utilized for ageing studies in e.g. spin glasses for which the system is maintained in the glassy state at a fixed T/T g < 1 [137]. In fact, in these ageing experiments on chalcogenides, the ageing temperature is usually found to vary because T g varies with composition. However, one still realized that flexible glasses (e.g. As 15 Se 85 ) which are below the RW age significantly over a three month waiting period (Figure 10), while stressed rigid glasses (As 40 Se 60 , above the RW) age somewhat slower, over ≃ 5 months, an observation that directly results from the slower ageing kinetics connected with higher glass transition temperatures. For RW compositions (As 30 Se 70 ), there is minimal tendency of ageing, even after a t w = 3000 h waiting period.
Calorimetric studies of Ge-P-Se glasses [133] have actually confirmed the weak tendency to ageing in the isostatic phase [134], and a deep and wide RW in these chalcogenides is found to sharpen as glass compositions outside the RW age at 300 K over different periods (Figure 9, inset).

Perspectives and conclusions
Recent theoretical efforts attempting to reconsider Rigidity Theory in a substantially revised version, represent attractive venues for an improved quantitative description of glasses and non-crystalline solids. A promising way to extend such approaches to the liquid state is provided by MD-based constraint counting which permits to calculate various properties from ensemble averages, and to connect to Rigidity Theory via an explicit account of the topological constraints. Given the general use of such simulations in the description of glassy materials, this recent extension now offers the possibility to rationalize the design of new families of other materials using as input the rigidity state of the underlying atomic network, as recently demonstrated for the case of cement or in the experimental studies of semiconducting materials. An alternative path is proposed, which introduces an analytical temperature dependence in the mean-field Phillips-Thorpe formulation. It permits to calculate thermal and mechanical properties as a function of composition, among which, quantities associated with glassy relaxation such as the fragility index or heat capacity. Figure 10. Ageing effect on the non-reversing heat enthalpy ΔH nr in As-se glasses [134]. For RW compositions (As 30 se 70 ), the evolution with waiting time is substantially reduced. The inset shows the same quantity as a function of composition in Ge x P x se 100 -2x [133] for two different waiting times t w .
Relationships between network rigidity/elasticity and the thermodynamics and relaxation of supercooled liquids appear to have an even more general ground as recently emphasized [82], and such concepts can actually also be extended to other glass-forming liquids including fragile ones [138]. Glass elasticity and the presence of soft elastic modes have been found, indeed, to drive many aspects of glassy relaxation, and also to relate to thermodynamic changes across the glass transition [87]. In fact, an abundance of such soft modes permits exploration of the phase space without large changes in energy [56], and it ultimately leads to small changes in the specific heat. This, of course, connects back to the notion of floppy modes [33] that are present in weakly connected (flexible) network glasses. The presence of floppy modes has been also seen as important in the mechanisms of dissolution [139] and subcritical crack growth [140], as recently revealed from a coupled MD/constraint analysis.
In this respect, the ongoing theory and experimental investigation of the particular isostatic compositions and the exploration of the nature of the intermediate phase would be of great interest for applications given its remarkable stability and the resulting anomalous behaviour for a variety of physical properties. A combination of molecular simulations and theoretical models of rigidity transitions has shown that adaptation under stress seems to be the key feature for the onset of an IP, and it is quite remarkable to acknowledge that the most recent results from molecular simulations [89] have a one-to-one correspondence with experimental measurements from calorimetry showing a deep minimum in the relaxation enthalpy once the glass is optimally constrained. In fact, without stress (present at high network connectivity) that induces bond/stress relaxation and without floppy modes (present at low connectivity in chalcogen rich glasses) that lead to low-energy relaxation phenomena, an adaptive isostatic glass has a significantly different relaxation regime that leads to a minimum in the enthalpic overshoot and the hysteresis appearing in glass transition cycles.
Other anomalies in properties both in the liquid and glassy phase are found for isostatic glasses. Among other signatures, it is worth to recall at this stage that structural anomalies might also occur, and these manifest by a maximum in typical correlation lengths, and coherence length as revealed from the non-monotonic evolution with density of the position k FSDP and the width Δk FSDP of the FSDP of partial structure factors [141]. This also suggests that there are structural signatures for the IP, as also demonstrated recently from a numerical study of the archetypal system As-Se [102]. Given the ongoing debate on the structural signature of the IP, undoubtedly, this should motivate further studies on this topic.

Disclosure statement
No potential conflict of interest was reported by the author.