Entropy favors heterogeneous structures of networks near the rigidity threshold

The dynamical properties and mechanical functions of amorphous materials are governed by their microscopic structures, particularly the elasticity of the interaction networks, which is generally complicated by structural heterogeneity. This ubiquitous heterogeneous nature of amorphous materials is intriguingly attributed to a complex role of entropy. Here, we show in disordered networks that the vibrational entropy increases by creating phase-separated structures when the interaction connectivity is close to the onset of network rigidity. The stress energy, which conversely penalizes the heterogeneity, finally dominates a smaller vicinity of the rigidity threshold at the glass transition and creates a homogeneous intermediate phase. This picture of structures changing between homogeneous and heterogeneous phases by varying connectivity provides an interpretation of the transitions observed in chalcogenide glasses. The mechanical and thermodynamic properties of amorphous materials are governed by their disordered network at microscales, but the detail remains elusive. Yan shows that the vibrational entropy induces a floppy-rigid phase separation near the rigidity onset and thus favors heterogeneous structures.

L acking long-range order, amorphous materials are fully governed by their microscopic structures. Increasing evidence indicates that the structural elasticity of such materials correlates with their dynamical properties and mechanical functions, such as the suddenly slowing relaxations of glasses [1][2][3][4] and the allosteric regulation of proteins 5,6 . A crucial factor behind the structural disorder that controls the linear elasticity of a structure is the average number of constraints n of its interaction network and the rigidity transition associated with tuning n 7 . At the Maxwell point n c = d 8 , which is the minimum number of constraints per particle to avoid floppy modes in spatial dimension d, both the elastic moduli and self-stresses vanish, accompanied by a vanishing onset frequency ω * of the soft vibrations on the socalled boson peak [9][10][11] . However, it is questionable whether these results obtained in homogeneous networks apply to heterogeneous network structures, which may be fundamental.
Chalcogenides, for example, are network glasses composed of chemical elements with different covalent valences r, proportional to which the number of covalent constraints n varies. Rather than a point threshold r c = 2.4 12,13 , a range of singular features, named the intermediate phase, bridges the well-connected stressed and poorly coordinated floppy phases, as observed in experiments [14][15][16] and reproduced in molecular dynamics simulations [17][18][19] . Inside the phase, the non-reversible heat, a glasstransition equivalent of the latent heat, vanishes 14 , which is associated with a vanishing stress heterogeneity 15 and a minimal molar volume 16 . All of these measurements are discontinuous when entering the phase from either side 16 . The critical point observed in random networks [20][21][22] (Fig. 1a), which allow fluctuations in local connectivities, fails to capture the nature of the intermediate phase. Emerging in self-organized networks to reduce the energetic costs of self-stressed states 23,24 (Fig. 1b), the rigidity window with distinct onsets of rigidity and self-stress promisingly maps to a critical range like the intermediate phase; however, the stronger heterogeneity inside the critical window actually contradicts the experimental observations, and the window is also sensitive to the appearance of prevailing perturbations such as van de Waals forces 25 . In fact, a rather odd feature is the heterogeneous nature away from the threshold, outside of the intermediate phase. What causes the heterogeneity beyond the local fluctuations?
Recent achievements 26,27 indicate that the entropy, a synonym of 'disorder', leads to order and heterogeneity in many cases, including the gas-crystal phase separation in colloid-polymer mixtures 28,29 and the open lattice structures of patchy particles 30,31 . The key components that allow for this comprehensive role are the high degeneracy of configurations and the separation of degrees of freedom carrying entropy from the ones assembling structures. In amorphous networks, configurations are inherently degenerate. Floppy and soft modes on boson peaks store significant amounts of vibrational entropy 32 , particularly close to n c ; thus, they inevitably shape the network structures.
In this communication, we investigate the role of entropy in regulating network structures and show the appearance of phaseseparated heterogeneous structures ruled by a critical point at the rigidity threshold. We then confirm the appearance of a homogeneous intermediate phase when stress energy dominates at low temperature. Finally, we apply the results to chalcogenides and discuss several experimental evidence of phase separation.

Results
Network model. To illustrate our main result of an entropyinduced phase-separated connectivity range near the rigidity threshold, we consider a network model on a two-dimensional triangular lattice with periodic boundaries, nevertheless, the result and its derivation depend neither on the dimension nor on the lattice structure. On lattice, a particle at each of N nodes can be wired to at most all of its six neighbors, corresponding to the maximal constraint number n m = 3. Following ref. 20 , we randomly perturb the locations of lattice nodes to avoid straight lines that lead to non-generic singular modes, as shown in Fig. 2a. The key assumption of the model is the separation of energy scales such that we can consider the network of the stronger interactions such as the covalent bonds in chalcogenides and treat the weaker ones such as van der Waals forces as perturbations. In the simplest construction, a network configuration Γ is defined by the allocation of N s ≡ nN linear springs of identical stiffness k on the n m N possible links. Different configurations are probed by relocating one random spring (red solid) to an unoccupied (blue dashed) link at a time, as illustrated in Fig. 2a, such that their number is fixed by a given average number of constraints n, similar to rearranging atoms of different valences in network glasses. The different configurations are sampled with probabilities proportional to the Boltzmann factor expðÀF=TÞ using the Metropolis algorithm, which is documented together with the model parameters in Methods section. Given configuration Γ, its free energy is where vibrational entropy S vib quantifies the volume of thermal vibrations near the mechanical equilibrium of Γ 31-33 , which depends on ω 2 -the eigenvalues of Hessian M and a Γindependent number c. H 0 is the self-stress energy of Γ at equilibrium. We introduce frustrations by imposing that the rest length of the spring γ positioned at the link ij h i, l γ ¼ r i;j h i þ ϵ γ , differs from r i;j h i , the spacing between neighboring nodes i and j, by a mismatch ϵ γ assigned from a Gaussian distribution of zero mean and variance ϵ 2 . In the small frustration limit, where ϵ is much smaller than the lattice constant, we compute M and H 0 in the linear approximation, as derived in the Supplementary Note 1 and refs. 4,25,34 .
We include perturbations of non-specific but weaker interactions by connecting all six second neighbors on the lattice with springs of stiffness k w ( k. At this high connectivity, they act approximately as isotropic potentials of effective stiffness α ¼ Entropy favors phase separation. As shown in Fig. 1c, in the limit of no self-stress penalty ϵ = 0 and thus no energy regulation H 0 ¼ 0, entropy-favored networks present a phase separation into two phases, a highly coordinated stressed cluster (n > n c dark green) and a floppy phase formed by the remainning clusters (n < n c blue), near n c , distinct from the homogeneous structures in Fig. 1a, b, where the percolating rigid cluster would appear indistinguishable from the remainder if the color code and the pivots are removed in Fig. 1. This phase separation is captured by a long-range correlation of the local constraint number and a bimodal cluster size distribution (a system-size stressed cluster plus small ones in the floppy phase) in contrast to a continuous one 35 , as shown in Fig. 3a, b.
Due to the phase separation, the network rigidity arises in a discontinuous fashion as the stressed cluster percolates-growing from an island inside the floppy sea to a continent enclosing floppy lakes. This percolation occurs at a constraint number n * different from n c , which is captured by a discontinuous P ∞ , the probability of springs in the percolating cluster, as shown in Fig. 3c. In Fig. 3d, the bulk modulus K shows a trend to jump at n * , whereas the shear modulus G vanishes.
Phase diagram. Why does entropy alone favor a floppy-rigid phase separation? As the degrees of freedom carrying vibrational entropy (particles) disconnect from the ones coding the configuration (springs), the total entropy increases by creating floppy modes in the floppy subpart of the network by confining springs in the stressed counterpart, particularly when this spring redistribution costs little configurational entropy near the rigidity threshold. When the self-stress energy is not participating, the balance between the vibrational entropic gain and the configurational cost determines the stability of the separation.
Consider a separation into a homogeneous rigid phase and a floppy phase of volume fractions V r and V f controlled by the constraint numbers n r and n f , as illustrated in Fig. 4a. The configurational entropy is the entropy of mixing springs and vacancies summed over the two phases, plus s c,0 , the entropy from the boundary contribution, which vanishes in the thermodynamic limit. As the extra vibrational entropy gains from the floppy modes, let us assume that the vibrational entropy is proportional to the number of floppy modes, changing by Λ per floppy mode. As shown in the Supplementary Note 2, this assumption is approximately valid in the model and per mode entropy gains λ ¼ À the spectrum-average entropy of non-floppy modes. Henceforce, we use the convention of the large Λ as a parameter in the formalism and the small λ as the actual entropic gain in the model. Constrained on the total volume V f + V r = 1 and the average constraint number n f V f + n r V r = n, the total entropy S vib + S conf is optimized with ; ð5Þ Since V r ∈ [0, 1], the heterogeneous phase exists in the selfconsistent range n ∈ [n f , n r ], which is very wide n r Àn f nc $ λÀ  Fig. 3 Features of network structures optimizing the total entropy. a Spatial correlation of connectivity C n ðrÞ = nðrÞnð0Þ h iÀn 2 À Á / nð0Þ 2 D E À n 2 for entropy favored networks, N = 576. b Probability distribution of rigid cluster sizes ρ(S), collapse for a wide range in n < n c . c Probability in the percolating cluster P ∞ , and d Bulk modulus K and shear modulus G vs. constraint number n for various system sizes N, α = 0.0003. The black solid line is theoretical prediction for the thermodynamic limit N → ∞. λ ≈ 3.3, so n f ≈ 0.94, n r ≈ 2.76, and n * ≈ 1.85, fitted by Eqs. (5), (6) and (7) Analogous to the classical spontaneous magnetization and gas-liquid phase separation, the entropy-induced floppy-rigid separation is governed by a critical point at Λ = 0 and n = n c , where the entropy gain Λ plays as the relevant parameter like temperature T c − T and the average constraint number n is the order parameter akin to the mean magnetization M in ferromagnet or the mean density ρ in gas-liquid separation. Close to the critical point Λ = 0 and n = n c , the free energy follows, The counting approximation, Eq. (4), gives α = −1. The order parameter scales as, The mean-field solution, Eqs. (5), (6) and (7) implies β = 1. Both exponents are different from the standard Landau theory. A typical size is thus determined by the critical scaling approaching (0, n c ). In the separation range, the network structure presents the dominant phase (V > 1/2) with droplets of the subdominant one of this characteristic size. Global rigidity arises when the rigid phase becomes dominant at n * = (n f + n r )/2, as indicated by the yellow line in Fig. 4a.

Self-stress penalty and homogeneous intermediate phase.
When creating self-stressed states is prohibited 23,24 , phase separation can still arise for n < n c due to an entropy gain of additional soft modes on the boson peak in isostatic structures. Per degree of freedom in isostatic volume V c , the vibrational entropy increases Λ ′ ∂S vib =N d∂Vc , positive as shown in the Supplementary Note 2. This gain from isostatic structures leads to a separation between an isostatic phase and a floppy phase, as illustrated in Fig. 1j. The corresponding phase boundary follows shown as the white dashed line in Fig. 4a. Because reducing the self-stress energy tends to level the connection distribution 25 , when the energetic cost H 0 competes with the entropic gain, a homogeneous intermediate phase can develop inside the heterogeneous gap at low temperature. In Fig. 1d, we depict the typical network structures equilibrating the total free energy Eq. (1) at the glass transition temperature T g . From left to right, which correspond to below, at, and above n c , the networks are floppy-isostatic heterogeneous, homogeneous, and floppy-stressed heterogeneous, respectively.
At temperature T (in the energy unit kϵ 2 ≡ 1), each selfstressed state contributes an independent direction to store energy 4, 34 . Noticing the duality between self-stressed states and (see the Supplementary Note 3 for the derivation). The selfconsistent condition of floppy-rigid phase separation breaks down when λ F (T) ≤ Λ(n), the phase boundary in Eq. (5). Relying on the insights of the elastic models 37 , we apply a glass transition temperature that is proportional to the shear modulus, T g ∝ G, whose analytical form is derived in the Supplementary Note 4. When n > n c , T g~n − n c 4, 34 , λ F , shown as the blue solid line in Fig. 4b, reenters the homogeneous phase when n decreases close to n c , n r − n c~α , defining the threshold of the homogeneous intermediate phase on the rigid side.
When n < n c , T g $ α ( 1 4, 34 , the self-stress prohibited situation applies. Derived from a flat mode density approximation 4,38 in the Supplementary Note 3, the free energy loss per isostatic volume, shown as the blue dashed line in Fig. 4b, surpasses the heterogeneous boundary Eq. (10) in the dashed line at n c À n f ≳ ffiffiffi α p , giving the transition from the intermediate phase on the floppy side. Altogether, as the connectivity increases, the network structures change from homogeneous floppy to heterogeneous floppy-isostatic to intermediate homogeneous marginal to heterogeneous floppy-stressed and finally to homogeneous stressed, as depicted in Fig. 4b.
Relative entropy. This floppy-rigid phase separation has a general information theory implication. Rewriting the phase boundaries n f (Λ) and n r (Λ) in Eqs. (5), (6) and (7) in terms of relative entropies 39 , DðpjqÞ = p ln p q + ð1 À pÞln 1Àp 1Àq , we find that The connection distributions of the floppy and rigid phases obey the conditions that the relative entropy density from the rigid phase balances the density from the floppy one to the critical network and the entropic gain per unit volume of the floppy phase compensates the relative entropy from the rigid phase to the floppy one. Similarly, when any self-stress structure is forbidden, the phase boundary follows The entropic gain per unit volume of the critical structure compensates the relative entropy from the floppy phase to the critical phase. As derived and numerically verified in the Supplementary Note 5 and Supplementary Fig. 2, these balances, as well as the main results on the phase separation, hold in general for networks of multiple types of interactions, which is the case of real chalcogenides and proteins 40 , as long as the vibrational entropy gain is approximately linear in probability distributions of interactions.
Segregation in network glasses. In network glasses, the degrees of freedom and the covalent constraints, both of which are associated with the atoms, depend differently on different chemical elements. The entropy-induced heterogeneous phase develops by segregating different elements. For illustration purposes, we derive the phase boundaries of compounds A x B 1−x , where x is the number fraction of atoms A, the knob equivalent to the number of constraints n. Both A and B atoms, as isotropic particles, possess d degrees of freedom. The number of constraints, counting both bond stretching and bond bending, satisfies n B < n c = d and n A > n c , so that both floppy and rigid networks can be produced by different compositions. We perturb near the segregation of a stressed rigid phase of volume fraction V r , B concentration ρ B r and A concentration ρ A r and a floppy phase of V f , ρ B f and ρ A f . The vibrational entropy obeys, with Λ the vibrational entropy gain from each floppy mode. The configurational entropy of two segregated regions is, Optimizing entropy with the following constraints, we end up with following phase boundaries, The boundary of the heterogeneous phase when self-stress is prohibited is determined by, As many constraints are associated with a high valence atom, the configurational entropy cost to generate phase separation is lower than in the network model by a factor of 1/n m . So the transition boundary Eq. (21) is at a much lower value than Eq. (14), and the segregation occurs easier. In particular, we plot the phase diagram in Fig. 4c for chalcogenides Ge x Se 1−x , where valences r Se = 2 and r Ge = 4 correspond to the number of covalent constraints n Se = 2 and n Ge = 7 counting both bondstretching and bond-bending contributions 13 . Segregations occur above the critical point (Λ c = 0 x c = 0.2), and five phases with four homogeneous-heterogeneous transitions appear at the glass transition in varying x.

Discussion
This comprehensive structural behavior provides a natural interpretation for the four transitions with discontinuous features, including transitions to the intermediate phase, as observed in chalcogenides when changing the chemical compositions 16 . Out of the intermediate phase, the micron-sized stress bubbles 15 are direct evidence of the heterogeneity. Its consequence on elasticity, the weakened shear modulus, is faithfully recorded in Raman scattering experiments 15 . Distortions of micro-structures shift the Raman peaks proportional to the global elasticity, Δν 2 ∝ G. As shown in Fig. 5, the jump of the Raman shift of the transversal optical branch in the intermediate phase 15 maps to the change of shear moduli between a homogeneous media and a heterogeneous mixture of two components 41 . In addition, high dynamical fragility out of the intermediate phase 16 is consistent with the appearance of very floppy structures 4 , and the Einstein relation breaks down with a floppy-phase-dominated diffusion and a stressed-phase-limited relaxation 19 , which results in a very stretched exponential relaxation 42 .
According to the model, ruling the transitions is predominantly the entropic gain λ, which is negatively correlated with α, the strength of the perturbing interactions relative to that of the strong ones forming the network. The width of the heterogeneous range is Δn / λ~À 1 2 ln α, whereas that of the homogeneous intermediate phase is Δn $ ffiffiffi α p . Thus the larger is the entropic gain, that is, in terms of experimental parameters, the stronger are the covalent bonds or the weaker are the van der Waals forces, the easier is the glass being frozen in a heterogeneous structure and the narrower is the intermediate phase.
This rule provides a general reference to the componentdependent widths of the intermediate phase 19 . Stabilizing the floppy parts as the weak interactions 43 , the pressure should be another experimentally approachable knob. Starting from a heterogeneous structure, increasing pressure effectively increases α and leads to a transition to the homogeneous phase 18 . However, further pressure that distorts the strong interactions, α~1, breaks our premise on the separation of energy scales and thus ends up in new physics 19 .
In conclusion, we have shown that the entropy favors heterogeneous structures in the vicinity of the rigidity threshold of networks. Based on the counting approximation 8 [14][15][16] . The counting approximation simplifies the entropic gain as a single parameter independent of the configurations. To go further, it is necessary to treat the entropic gain more carefully and study the global minimum and the dynamics toward it in a rougher free energy landscape induced by the complex entropic consequences of structures such as long chains. Meanwhile, it is important to test the separation in molecular dynamics simulations 17 for various temperatures and non-specific weak forces. Finally, it is useful to apply the role of entropy in protein foldings and self-assembly, where flexible units appear vital for elastic functions 5,45 .

Methods
Metropolis algorithm and chosen parameters. We equilibrate network structures Γ using the Metropolis algorithm. From an initial configuration Γ, a new configuration is proposed by the random relocation of a spring, as illustrated in Fig. 2. By comparing the free energy Eq. (1) between the current and the new configurations, we sample and reset to the new configuration with probability min 1; exp À Data availability. The authors declare that the data supporting the findings of this study are available within the article and its Supplementary Files, or are available from the authors upon reasonable request. The numerical data of the network model were generated by a home-written code on MATLAB interface. This code is available upon request. Theory for homogeneous Heterogeneous ref. [41] Data from ref. [15] G (a.u.) -0.5 0 0.5 n − n c Fig. 5 Experimental evidence of structural heterogeneity. Shear modulus G in arbitrary unit predicted for homogeneous networks (blue) and for heterogenous networks composed of two phases at ends (red). Black circles are Raman shift data of the transverse optical (TO) branch ν 2 TO =ν 2 0 À 1 in ref. 15