Microscopic Model of Intermediate Phase in Flexible to Rigid Transition

We introduce a lattice gas model with a modified Hamiltonian considering different energy for cycles of connected atoms. The system can be interpreted as a chalcogenide glass with pollutants forming floppy and rigid structures. We consider an energetic penalization for redundant bonds in the network. This penalization allows us to incorporate the topology constraints of rigidity in the network to study the thermodynamics of the system. We observe, depending on the parameter used for the penalization, that the system exhibits a typical first-order phase transition, or a stepped transition between the low and high density while varying the chemical potential. We also observe a hysteresis loop in the density and energy of the system. We use the area of these loops to calculate the irreversible enthalpy. There are two regimes, one where the enthalpy decreases linearly and the other with almost constant enthalpy. As the enthalpy is almost constant and very low, we interpreted this as the intermediate phase of the chalcogenide glasses.


INTRODUCTION
Lattice gas models are among the simplest thermodynamic models that exhibit a phase transition with an exact solution in 2D. The nearest neighbor interaction and limitation of occupation in volume allows emulating a real gas in potentials such as Lennard-Jones [1]. This model has proven to be useful in different contexts, such as condensation of DNA [2], or the absorption in controlled-pore Glasses [3]. There is also a direct relationship with the Ising model, which was first used to study ferromagnetic materials [4], and then many other materials as spin-crossover materials [5], or spin glasses [6] among others. On the other hand, chalcogenide glasses seem to be some of the most promising materials for future technology, with important applications [7][8][9], ranging from solid state batteries [10,11] to optics and photonics infrared devices [12,13]. Topology and rigidity of the network are characterizing properties of these glasses [14,15].
Experimental modulated differential scanning calorimetry (MDSC) and computational molecular dynamics (MD) studies near the glass transition over chalcogenide materials have found anomalies in the behavior of the macroscopic variables of these materials, giving rise to what is known as the intermediate phase [16][17][18]. Although theoretical explanations regarding the significance and existence of this phase exist, as far as we know, no microscopic model which recovers the thermodynamic macroscopic properties of the system has been constructed to this date. The purpose of the paper is to provide a simplified microscopic model that reproduces the behavior of thermodynamic variables in the intermediate phase.

INTERMEDIATE PHASE AND RIGIDITY IN CHALCOGENIDE GLASSES
The anomalies mentioned previously, along with most of the bulk properties of chalcogenide glasses, have been related to a microscopic property of their covalent network called rigidity [19]. It can be defined as the property of atoms being able to move without deforming current bond angles and lengths. More precisely, a whole mathematical formalism can be developed to study rigidity [20].
Independently of the mathematical theory of rigid networks, simplified models have been extensively used to study the rigidity of glasses, with particular focus on reproducing the transition between rigid and floppy modes [21]. The most relevant of them is the percolation rigidity model, based on constraint counting, which can be exactly solved using mean-field approximations [22].
Besides the obvious complexity arising from the absence of long-range order in amorphous solids, the challenge in building a model that recovers macroscopic properties relies on the difficulty of representing the vibrational entropy accumulated near the boson peak and the inherent degeneracy of most of the configurations [23].
Such difficulty is directly related to the problem of effectively quantifying the rigidity of the network so it can be incorporated into the microscopic model. Mean-field theories are incapable of describing the microscopic scale accurately. In the particular case of two-dimensional networks, the pebble game algorithm [24] is capable of decomposing a network onto its rigid components with sufficient speed, but it does not give a method of relating rigidity to thermodynamic variables.
The algorithm relies on Laman's theorem [25] which characterizes exactly the rigidity of a network embedded in a two-dimensional Euclidean space. The referred theorem hasn't been successfully extended to other dimensions due to the difficulty of the exactly characterizing the rigid components of a network embedded in an arbitrary geometry. For dimensions greater than two, approximations are commonly used [26]. As we are only interested in quantifying the rigidity as a function of thermodynamic variables in a more accurate way than the meanfield theories, we could aim to use a non-precise but a simplified model of the network by using a modified lattice gas model. This model should take into account the results obtained by the pebble game algorithm when simulating the transition to a rigid system.

Description and Behavior of Chalcogenide Glasses
Chalcogenide glasses are amorphous solids built upon members of the 16 group of the periodic table (S, Se, Te) by doping them with members of another group, most commonly group 15 (As,Sb). To a molecular level they can be completely described by a continuous random network (CRN) [27]: a molecular network where each edge represents a covalent bond. Although van der Waal forces between pairs of free electrons are normally present in the system, those interactions are weak enough to be left out of the CRN model. Raman Spectroscopy allows us to obtain the resonant frequencies of the vibrational modes of the network, which is related to the different molecular structures (components) of the network. This information allows to calculate the entropy by using the formula S x i log(x i ) where x i is the relative fraction of each component [28]. This approach has been particularly useful in numerical studies [23].
The relationship between the topology of the network and its rigidity to the macroscopic properties of the system via the changes in its density of states was first proposed by Phillips [14,15] and further confirmed from mean-field constraint counting approaches to rigidity by Thorpe [19,29]. Experiments have also shown that when examining glasses of the same compound but different stoichiometries, which is equivalent to changing the mean coordination number of the network, macroscopic properties change as a function of the stoichiometries and present a transition when passing for coordination numbers similar to the theoretically predicted by mean-field theory [30,31]. In addition, when performing experimental MDSC calorimetry studies of chalcogenide glasses, we can measure the heatflow during endothermic and exothermic processes of the system. With these measurements, the irreversible enthalpy when passing through the glass transition can be obtained. When analyzing such data as a function of the stoichiometry of the glass, a reversibility window is found, in which the irreversible enthalpy vanishes [16]. Such a window can not be directly associated with a rigid or floppy phase of the CRN; it forms a new phase called intermediate phase [16].
The importance and existence of the intermediate phase are confirmed by the exotic behavior of other macroscopic properties of the glasses in such a window [32], such as ionic conduction [33] and infrared reflectance [34]. The intermediate phase has also been observed by measuring the configurational entropy of the system [23,28,35,36]. Studies on other chalcogenide glasses and oxide glasses also exhibit an intermediate phase with similar anomalies in the macroscopic variables [17], and it can also be observed in molecular dynamics simulations [18].
Such experiments have also measured quantities that are directly related to the average coordination number and the number of floppy modes (see for example [17] and references inside), which can be used to quantify the constraint density [16]. Those quantities have been also developed in analytical and numerical treatment of this materials [14,17].
Despite the experiments and simulations carried out the intermediate phase is still controversial especially due to contradicting experimental studies in which a structural origin of the phase isn't found [37][38][39][40][41].

Microscopic Models
Besides the experimental controversy, some efforts have been made in order to construct a microscopic (structural) model that reproduces the exotic behavior of the macroscopic variables of the system. Outside of the intermediate phase context, several microscopic models have been developed [27,42], more recently putting effort into describing the glass transition and the Arrhenius-like behavior [43,44].
One of the first models [45], uses a random bond network that allowed to control the number of rings of bonds, which are problematic in the rigidity percolation theory. The authors used a computer simulation that allowed constructing the network beginning in a floppy mode and growing it by adding bonds that restrict movement in a non-redundant way (i.e., the bonds reduce the degrees of freedom of the system). After all of the independent bonds have been placed, the bond allocation continues and its performed at random. With this model, the authors were able to observe a transition from a floppy state to a rigid, non stressed (i.e., without redundant bonds) and from there to a stressed network. The key element in the model is in the procedure of the building that allows for self-organization in the network.
This model was further developed by Chubynsky et al. Ref. 46 in order to work with networks in thermodynamic equilibrium, as the network construction process of the original model could lead to highly atypical networks for the system. This model changes the original only in a small way: it grows the network from a floppy state by adding independent bonds, but every time it adds a bond it deletes and creates a bond that doesn't change the stress (i.e., the redundancy of another bond) of the network.
Another model Ref. 47 also uses the self-organization of the network, in this case by explicitly describing the Hamiltonian as the number floppy modes. The network is restricted to a Bethelike lattice (of finite size). The system is then studied using Monte Carlo simulations, switching configurations by rewiring two randomly selected nodes. The intermediate phase is found in terms of changes in the probability of a stressed cluster that exists and percolates through the entire lattice.
It must be acknowledged that the three mentioned models depend on the Pebble Game algorithm in order to describe the independent bonds of the network. Other models have been studied before those mentioned [48]. They do not depend directly on the exact description of the rigid components of the network, but instead approximate those components via loops or cycles of the covalent graph. They are also known as tree-like percolation models. Most of these models do not produce a uniform ensemble with equal probability for all tree-like structures, although they can be treated as if they were in thermodynamic equilibrium. Tree-like percolation models also include self-organization by avoiding the building of loops.
Models that directly attack the thermodynamic properties from the hamiltonian either analytically [49] or numerically via stochastic descriptions [23] have also been developed. Another important aspect of the models in this subsection is the fact that most of them are designed only to describe the system in moments where we can accurately characterize the vibrational entropy of the atoms. Such an assumption implies that the temperature T of our systems is much smaller than the Debye temperature T D .

MODEL
Similar to tree-like percolation models, we can argue that in twodimensions the redundant constraints put by adding a bond to an atom with zero degrees of freedom is equivalent to whether the connectedness of the graph will depend on such bond. In a lattice gas model, an independent bond will be a bridge edge of the graph, and rigid components become equivalent to components isomorphic to cycles. To obtain the macroscopic variables from the microscopic model, we need its hamiltonian. We will base our model in the lattice gas model, which has the following hamiltonian: where μ is the chemical potential, n i is 1 or 0 depending if the site i is or not occupied, and J is the energy of the bound between two nearest neighbors represented by 〈i, j〉. This hamiltonian will be modified to take into account the energy cost of stress (redundant constrains). In our model these redundant constrains are the bonds forming cycles. Then, adding a cost energy of the bonds belonging to a cycle to Eq. 1, we obtain: (2) here L is the set of all the nearest neighbors that form clusters without bridge edges. This is equal to the set of all the edges between vertices in cycles. C > 0 represents a penalization for forming rigid clusters that delay the normal phase transition of the system between a low-density and a high-density state. This allows our system to find new non-rigid configurations and stay near them for large simulation times. Figure 1 shows occupied sites as circles, while the bridge edges are in black, and cycles are plotted in purple.
Our model differs from the original tree-like percolation [50] due to the fact that, although it is also an Ising-like model, in the tree-like model the connectedness problem arises as an interpretation of the system, meanwhile that in our model we are explicitly modifying the hamiltonian with a new term.
In the limit T → 0 we can approximate how the system in a square lattice will behave in order to orient the range of parameters of the Hamiltonian for which we will see certain transitions.
For C ∼ μ ∼ J, if we have N atoms with a fraction l of them forming rigid clusters, we can suppose than adding a new atom adds three rigid bonds and change two normal bonds to rigid bonds on average (see Figure 2), which becomes equally likely when which gives a critical value μ c1 for the chemical potential when the inequality becomes an equality: and suggests a change between medium and high density. Physically, this medium density corresponds to the case where adding particles in rigid, non stressed components is favored. This can be related to the constrain density n c of the network, so that μ > μ c1 will correspond to high constrain density n c (fully connected system). For μ, J ≫ C or C ≈ 0 we have the same behavior than a normal lattice gas, and the critical value for μ is expected at μ c2 −2J with a jump between low and high density.
For sufficiently high values of C ≫ μ, J energy minimization will be achieved by having as less cycles as possible. For N atoms in the system not belonging to a cycle, it will be equally probable to add or remove atoms if: which happens for μ c3 −J and predicts a transition between high and low density states. μ < μ c3 will correspond to a floppy system with low constrain density n c ∼ 0. For μ c1 aμ c3 , which is obtained for C c a2J/5 we can expect only a double transition in the system. Summarizing, μ c1 the critical value for the transition from a medium to high density (rigid to stressed), μ c2 is the typical critical value from zero density to one in the lattice gas model, and μ c3 correspond to the value where appear a transition from zero to non zero (floppy to rigid) density.

SIMULATIONS
We studied the model performing Markov Chain Monte Carlo simulations, based on the Metropolis-Hastings algorithm. The stability and convergence of chains were analyzed to determine the number of Monte Carlo Sweeps (algorithm steps per lattice site) needed to achieve convergence. For τ 100 Monte Carlo Sweeps all the systems are thermalized. This is also confirmed after seeing a drastic reduction of the standard deviation of samples between the whole system and the last τ/2 states.
We simulated for 41×31 equally spaced values of (μ, C) in the [−3.5, 0.0] × [0.5, 1.2] interval and fixed J 2, K B T 0.5 over four independent square lattices with periodic boundary conditions of side L 32 and 40. Macroscopic variables were calculated as the average of the thermalized values as shown in Figure 3.
The simulation was coded in the Julia Language [51], and can be found online in a public repository at: https://github.com/ sayeg84/latticeModels

Macroscopic Variables
In a canonical ensemble, the probability to be in a particular state is proportional to exp(−H/KT) exp(−βH). This give us a relation between temperature and the parameters of the hamiltonian. Developing this expression we obtain exp(−βH) exp(βμ n i )exp(−β(H int + H rig )).
Given so, varying the parameter μ is roughly equivalent to varying β 1/KT, by this we mean both variables can be used to induce the transitions in a similar way, but the corresponding critical exponents are different [52]. Varying the other two parameters is not equivalent to varying β, since these are not multiplied not by the sum of n i , but by the sum of the product (for many configurations, the product βC n i n j 0). However, these parameters can be useful to extend the model to study the effects of pressure or to perform simulations varying T instead of μ.
Changing pressure on the system is equivalent to deforming the lattice (change the volume V) which translates into a change of the potential energy because J J(V) and C C(V). The probability to be in a particular state in the NPT isothermalisobaric ensemble is proportional to exp − β(H + PV − Nlog(V)/β) [53]. This probability looks like the probability for where A 0 is the harmonic contribution to the elastic interaction energy between two neighbors [54]. However, we still don't have a good model for C(V). This ensemble is also probably more suitable to variate the temperature than the NVT-ensemble, since the experiments are performed at constant pressure, and if the volume is fixed, when the temperature increases the probability to change a node tends to 1, which means that the density tends to 1/2, so the whole transition is not reached, the reason why we used μ instead. We would now like to give an interpretation of C. Average coordination number 〈r〉 and rigidity are correlated and are considered as the defining parameters for the study of the intermediate phase. Redundant bonds have an energy cost that we related with C in Eq. 2. In this case, the total energy cost is the product of all the redundant bonds times C. We can variate then the energy cost by varying the number of redundant bonds (so, by varying 〈r〉) or by varying C. So, we can interpret the variation of C as the variation of the percentage of pollutants if the reached configuration is the same. Then, small values of C correspond to flexible configuration and high values of C to rigid or stressed configurations.
The susceptibility χ and the specific heat C v show two large jumps near values μ c3 ≈ − 2.5 and μ c1 > μ c3 . We observe that μ > μ c1 shows a high density state, while μ c1 > μ > μ c3 gives a medium density state and μ c3 > μ returns to low density.
Because a medium-density state could be unexpected for a model of this kind, we checked that its existence is independent of system size by performing simulations for smaller sizes. It occurs for every size and even becomes more stable as N grows. We define μ c1 as the value of μ where there is a transition from medium to high density. In Figure 3 we can see that this also corresponds to the value of μ that maximizes the derivatives of the variables, that in this case result in the density susceptibility χ and the specific heat C v . To better visualize these plots, in Figure 4 we show a cut of each macroscopic variable when $C 0.92$. We define μ c1 (C v ) argmax(C v (μ)) and μ c1 (χ) argmax(χ(μ)). In Figure 5 is displayed the relationship between μ c1 (C v ), μ c1 (χ) and C.
The linear fit of the data displays parameters close to the theoretical analysis done previously (section 3, and equation refeq: mu1) with a high correlation coefficient even if the analysis was very rough. An analysis of hexagonal and triangular 2D lattices revels that the macroscopic variables exhibit the same behavior for those values of μ but for a different interval of C. This is expected from a theoretical point of view due to the critical point dependence of the average rigid bonds added when going from the spanning tree configuration to a rigid configuration.

Hysteresis
The simulations presented in the previous subsection were all performed by initializing the system with a random low-density configuration and making the simulations over it in a μ-increasing direction. When performing the same procedure but for μ-decreasing and beginning with a high-density configuration, the macroscopic thermodynamic variables show a different path. This difference in the path is called hysteresis.
Hysteresis is usually related to loss of internal energy and the work that the system produces [55,56]. The area of the hysteresis loop in the density of the lattice gas is directly proportional to the work ΔW, while the area of the hysteresis loop in the internal energy corresponds to the loss of energy ΔU in the process. Using the first law of thermodynamics, we can obtain the heat ΔQ ΔW + ΔU, which we associate with the irreversible enthalpy. The area enclosed by the loop can be positive or negative depending on how the loop is walk by varying μ. To correctly define the sign, we will say that the area is positive if the loop is walked counterclockwise, and negative otherwise.
The hysteresis loops and their areas, corresponding to the work produce by the system ΔW and the change of the internal energy ΔU are shown in Figure 6. Figure 7 shows the heat ΔQ ΔW + ΔU as a function of C. The integrals used to calculate the area of the loops were calculated using a trapezoidal rule.
Comparing Figures 3, 6 we can see that the difference in internal energy ΔU is negative when there is a middle step in the transition when varying μ. The interval where this occurs is from C ∈ [0.6, 1.14], and the minimum is reached for C ∼ 1.0. For ΔW, the function decreases in the interval C ∈ [0.5, 0.65], for C ∈ [0.65, 85], a local minimum is reached, then the function becomes increasing until a maximum around C ∼ 1.1 is found, and then it decreases again until C ∼ 1.14, where ΔW ∼ 0.
As a result ΔQ ΔW + ΔU is a linearly decreasing function in the interval C ∈ [0.5, 0.78], while in the interval C ∈ [0.78, 1.14], ΔQ ∼ c remains constant, with c ∼ 0.1. We associate this interval C ∈ [0.78, 1.14] (shaded area) with the intermediate phase reported in chalcogenide glasses when a pollutant is added. This region was also obtained as the region where there is a step on the plot of ρ, passing through ρ 0.5. However for values C > 1.14, we did not observed an increase of ΔQ as reported in literature [17], instead we observed ΔU ∼ ΔW ∼ ΔQ ∼ 0. We associate this with the limitations of the model, where a transition is not fully achieved when C is too high. Recalling the association of the parameter C with 〈r〉 we see that our results have qualitatively a similar behavior to the experimental studies in literature [17]. Comparing with the other numerical results of microscopic models, our results display phase separation in terms of the macroscopical variables related to the hamiltonian. Even if a more concise comparison with experimental results is not able with our work, the model is sufficiently generic to be expanded in order to achieve more precise simulations that could reproduce experimental results.

CONCLUSION
We present a simple model that exhibits a change in enthalpy behavior when varying the rigidity of the system, a behavior similar to that reported in chalcogenide glasses when increasing the amount of added pollutants. The fact that the model is simple results in efficient simulations, that allows us to make large enough systems and study how they change by varying different parameters. In addition, this model allows us to obtain density, internal energy, density susceptibility, and specific heat.
For certain parameters, the model exhibits a step transition in the density, while for others the transition is with no medium values of density. From the values of C v and χ we observe that the first transition is not clear while the second seems to be a first-order transition. We interpreted the two possible states as "solid" and "fluid".
We were able to analytically approximate the parameter range where the transitions would appear. Our numerical results are in close agreement with such rudimentary approximations. Furthermore, our model is the first to connect the microscopic properties with the macroscopic thermodynamic variables.
Studying the hysteresis loops we were able to observe a change in behavior of the enthalpy, which is related to the change in density observed when the chemical potential (or temperature) varies when the transition becomes stepped width a medium density.
Despite success in qualitatively describing the first transition (flexible-rigid) to the intermediate phase, we were not able to reproduce the transition to stressed systems. We speculate that to see such a transition we would need to vary the temperature instead of the chemical potential.   Figure 6, the shaded area corresponds to the region where there is a step on the plot of ρ, passing through ρ 0.5.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AP performed the simulations and both authors contributed to designing the simulations, discussing the results, and writing and revising the manuscript prior to submission.

FUNDING
This work was partially funded by Universidad Nacional Autónoma de México via PAPIIT project IA106618.