%~Mouliné par MaN_auto v.0.27.3 2023-02-06 14:12:22
\documentclass[CRPHYS, Unicode, XML, screen, Thematic, HideArticleType,published]{cedram}


\newcommand{\clf}{\mathcal{f}}
\newcommand{\clN}{\mathcal{N}}
\newcommand{\FSDP}{\mathrm{FSDP}}

%\sethlcolor{green}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}
\ItHasTeXPublished

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\RubricEF{From everyday glass to disordered solids}{Du verre quotidien aux solides désordonnés}


\title{Topological ordering during flexible to rigid transitions in disordered networks}
\alttitle{Ordonnancement topologique lors des transitions flexible-rigide dans les réseaux désordonnés}

\author{\firstname{Matthieu} \lastname{Micoulaut}}
\address{Sorbonne Université, Physique théorique de la matière condensée, CNRS UMR 7600, 4 place Jussieu 75252 Paris cedex 05 France}
\email[M. Micoulaut]{matthieu.micoulaut@sorbonne-universite.fr}
\thanks{The author acknowlegdes base and recurrent funding from Sorbonne Université and CNRS. He no longer wishes to be supported by research grants from the controversial Agence nationale de la recherche (ANR)}



\keywords{Glasses, structure, network rigidity, elastic phase transitions, molecular dynamics simulations}
\altkeywords{verres, structure, rigidité du réseau, transitions de phase élastiques, simulations de dynamique moléculaire.}

\begin{abstract}
This contribution focuses on the structural origin of flexible to rigid transitions and the possible underlying intermediate phase which have been reported to occur in a variety of network glasses such as chalcogenides or modified oxides. Here, using molecular dynamics simulations of densified glass-forming liquids, 2SiO$_2$-Na$_2$O, which are known to display a numerical reversibility window as a signature of an intermediate phase, we focus on structural functions emphasizing topological ordering using the Bhatia--Thornton formalism. Results not only reveal that densified silicates display topological ordering on lengthscales of about 25~\AA, but also display obvious threshold behaviors close to the isostatic condition when the network undergoes a flexible to rigid transition. The mechanical constraint count of the atomic network structure reveals that a typical lengthscale characterizing the decay of topological correlations emerges for stressed rigid systems at $\simeq$3.5~\AA, whereas small wavevector oscillations are found to be minimal when the isostatic condition is merely satisfied. An additional analysis building on diffusivity and liquid entropy suggests that the {\em locus} of flexible to rigid transitions has also connections with water-like anomalies of densified tetrahedral liquids.
\end{abstract}

\begin{altabstract}
 Cette contribution se concentre sur l'origine structurale des transitions de flexible à rigide et une possible phase topologique intermédiaire sous-jacente qui se produisent dans une variété de verres structuraux tels que les chalcogénures ou les oxydes modifiés. Ici, en utilisant des simulations de dynamique moléculaire de liquides vitreux densifiés, 2SiO$_2$-Na$_2$O, qui sont connus pour présenter une fenêtre de réversibilité en pression lors d'une transition vitreuse numérique, nous nous concentrons sur des corrélations structurales mettant l'accent sur l'ordre topologique en utilisant le formalisme de Bhatia--Thornton. Les résultats révèlent non seulement que les silicates densifiés présentent un ordre topologique sur des échelles de longueur d'environ 25 \AA, mais présentent également des comportements de seuil évidents à proximité de la condition isostatique lorsque le réseau subit une transition de flexible à rigide. Le comptage des contraintes mécaniques de la structure du réseau atomique révèle qu'une échelle de longueur typique caractérisant la décroissance des corrélations topologiques émerge pour les systèmes rigides sous-contraints à $\simeq$3,5 \AA, alors que les petites oscillations du vecteur d'onde sont minimales lorsque la condition isostatique est simplement satisfaite. Une analyse supplémentaire basée sur la diffusivité et l'entropie du liquide suggère que le lieu des transitions de flexible à rigide a également des liens avec les anomalies observés dans les liquides tétraédriques densifiés comme l'eau sous pression.
\end{altabstract}

\makeatletter
\def\cdr@preabstracthook{%
%\vspace*{15pt}
\vskip15pt plus 8pt minus 1pt
}
%\def\cdr@interabstracthook{%
%\vskip2pt plus 3pt minus 1pt
%}
\makeatother
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\dateposted{2023-04-20}
\datepublished{2024-04-26}
\begin{document}
\maketitle

\section{Introduction : from mean-field rigidity to reversibility windows}
The atomic scale description of non-crystalline solids and glasses has remained a challenging field of inquiry in materials science because of the lack of translational periodicity that prevents from using powerful techniques offered by e.g. crystalline band theory such as Brillouin zones, k-space and Bloch functions. In addition, the structural information gathered from light, neutron or X-ray scattering often remains at the level of short order (bond distances, coordination numbers) as long-range order made of the repetition of a well-defined unit cell is notoriously absent. Some information can be gathered from the low wave-vector region of e.g. the structure factor and permits to correlate salient features of scattering functions such as the first sharp diffraction peak (FSDP) to certain features of atomic ordering~\cite{bychkov_prb2005,salmon_nature2005,salmon_jpcm2007}. Only a limited number of glasses where structural groupings are present, essentially ring structures, can be unambiguously analysed and their structure characterized from explicit spectroscopic signatures using Raman or Nuclear Magnetic Resonance (NMR) spectroscopy~\cite{galeener_prl1976,ferlat_prl2008,vignarooban_epl2014}. On the theoretical side, atomistic simulations can be extremely valuable for understanding structural and dynamic properties including aspects of the glass transition and relaxation of the supercooled state~\cite{liu_cras2022}. They also appear to be useful for checking the validity of microscopic models, empirical approximations or theoretical assumptions. The present contribution clearly belongs to this general scope, and attempts to link certain features derived from simulations with results from rigidity theory.


The concept of rigidity in amorphous networks and glasses traces back to the early work of Maxwell on the stability of macroscopic structures such as bridges or roof frameworks~\cite{maxwell_philmag1864}, and to the introduction of mechanical constraints by Lagrange. These ideas and results were then extended to atomic networks by Phillips~\cite{phillips_jncs1979} who highlighted the notion of mechanical isostaticity as promoting glass-forming tendency of covalent alloys. It was recognized that so-called ``good'' glass formers usually form at an optimal network connectivity, or mean coordination number satisfying the Maxwell stability criterion of isostatic structures, i.e. $n_c$ = 3 in 3D, where $n_c$ is the count of atomic constraints per atom arising from relevant neighbor interactions.

\subsection{Mean-field molecular rigidity}
It is useful to review basic results from the Phillips--Thorpe mean-field rigidity theory~\cite{thorpe_jncs1983}. In structural glasses the relevant interactions can be near-neighbor bond-stretching (BS) and next-near-neighbor bond-bending (BB) forces that characterize bonds and angles typical of the network topology. In a mean-field approach, the atomic density $n_c$ of mechanical constraints can be exactly computed, and is given by:
\begin{equation}\label{nc0}
n_c={\frac {\sum\limits_{r\,\geq\,2}n_r\left[{\frac {r}{2}}+2r-3\right]}{\sum\limits_{r\,\geq\,2} n_r}}
\end{equation}
where $n_r$ is the concentration of species being $r$-fold coordinated. Here it should be noted that the temperature is absent. The contribution of the two terms in the numerator is straightforward because each bond is shared by two neighbors, and one has therefore $r$/2 bond-stretching (BS) constraints for a $r$-fold atom. For BB (angular) constraints, one notices that a 2-fold atom involves only one angle, and each additional bond needs the definition of two more angles, leading to the estimate of (2$r$-3). For one-fold terminal atoms, a special count~\cite{boolchand_prb1994} is achieved as no BB constraints are involved, and in certain situations some constraints may be ineffective~\cite{zhang_science1994,micoulaut_amin2008}. By defining the network mean coordination number $\bar r$ of the network by:
\begin{equation}\label{eq1}
\bar r={\frac {\sum\limits_{r\,\geq\,2}rn_r}{\sum\limits_{r\,\geq\,2} n_r}}
\end{equation}
one can reduce~\eqref{nc0} to the simple equation:
\begin{equation}\label{nc_final}
n_c={\frac {\bar r}{2}}+2\bar r-3
\end{equation}
It is easy to verify that the Maxwell isostatic criterion ($n_c$=3) corresponds to a threshold average coordination number~\cite{phillips_jncs1979,thorpe_jncs1983} of $\bar r$=2.40 in 3D, corresponding usually to a non-stoichiometric composition where glass-forming tendency has been found to be optimized experimentally (e.g. GeSe$_4$ instead of the silica analog GeSe$_2$ in the Ge-Se system), including in systems with a strong crystallization tendency~\cite{micoulaut_prb2021}.



\begin{figure}[tbp]
\includegraphics[width=0.6\linewidth]{rect1554.png}
\caption{Fraction of floppy modes $\clf$ as a function of average network coordination number $\bar r$ (adapted from Ref. \cite{jacobs_jncs2000}). The mean-field Maxwell line $f=3-n_c=6-5\bar r/2$ vanishes exactly at $\bar r$=2.4 (solid red curve, eq.~\eqref{nc_final}), whereas percolation of rigidity is obtained from bond-depleted networks at a somewhat lower connectivity (solid black, \cite{thorpe_jncs1983}). At T=0, a self-organized removal of bonds (avoiding the presence of stressed rigidity, solid and broken red curve) leads to two percolative stress and rigidity transitions (green squares), and define a so-called intermediate phase (IP). See also Figure~\ref{thorpi_fig}.}
\label{mike}
\end{figure}

The nature of the underlying phase transition associated with the isostatic criterion has been revealed from a vibrational analysis of bond-depleted random networks~\cite{thorpe_jncs1983,he_prl1985} constrained by harmonic bond-bending and bond-stretching interactions. In such networks, the number of zero frequency (floppy) modes $f$ = 3-$n_c$ (i.e. the eigenmodes of the dynamical matrix) is vanishing for $\bar r$=2.38 when rigidity percolates through the network (Figure~\ref{mike}). The Maxwell condition $n_c$ = 3 therefore defines a mechanical stiffness transition above which redundant constraints produce internally stressed networks, identified with a stressed-rigid phase~\cite{he_prl1985}. Here, $\bar r$ acts as the external control parameter of the transition, whereas $f$ is the order parameter. For n$_c$ $<$3 however, floppy modes can proliferate, and these lead to a flexible phase where local deformations with a low cost in energy (typically 5~meV~\cite{kamitakahara_prb1991}) are possible, their density being given by: $f=3-n_c$. There have been early experimental probes of this peculiar transition from Raman scattering~\cite{feng_prl1997}, stress relaxation~\cite{boehmer_prb1992} and viscosity measurements~\cite{tatsumisago_prl1990,gunasekera_jcp2013,ravindren_jcp2014}, vibrational density of states~\cite{kamitakahara_prb1991}, Brillouin scattering~\cite{vaills_jpcm2005,micoulaut_amin2008}, resistivity~\cite{novita_prl2007,micoulaut_prl2010}, and Kohlrausch exponents~\cite{boehmer_prb1992,micoulaut_prb2009}.



\subsection{Beyond mean-field}

\begin{figure}
\begin{center}
\includegraphics[width=0.55\linewidth]{mike2022.pdf}
\end{center}
\caption{\label{thorpi_fig}
Evidence for a stress-free intermediate phase (IP) from the Pebble Game analysis (adapted from~\cite{jacobs_jncs2000,yan_prl2014}). Fraction of sites on isostatically rigid and stressed rigid percolating cluster in a self-organized network as a function of the network mean coordination number $\bar r$. The intermediate phase which is rigid but unstressed exists in these classes of models between $\bar r_{c1}$=2.375$\leq \bar r\leq$ 2.392=$\bar r_{c2}$, and coalesces in random networks.}
\end{figure}

The construction leading to the mean-field constraint estimate (eq.~\eqref{nc_final}) neglects fluctuations that have been identified numerically at finite temperature~\cite{bauchy_epl2013,bauchy_prl2017} and considers only network averaged quantities, i.e. one assumes that all atoms of a given type have the same constraint density arising from an identical coordination number per chemical species. However, as in ordinary phase transitions fluctuations in constraints or in the order parameter $f$ may be expected close to the critical point at $n_c$ = 3. It has been suggested (for a review, see~\cite{micoulaut_advx2016}) that this non-mean field scenario manifests by the onset of an isostatic or intermediate phase (IP) 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 during the glass transition (red curves, Figure~\ref{mike}).


A certain number of scenarios have been proposed to describe the observed behaviors, and these emphasize either the role of fluctuations in the existence of a double threshold/transition defining an IP between the flexible and the stressed rigid phase~\cite{jacobs_jncs2000,barre_prl2004,micoulaut_prb2003, micoulaut_prb2006,briere_pre2007}, or on a non-trivial elastic coupling between atoms with limited coordination fluctuations~\cite{yan_prl2014,yan_natcom2018}. In the latter scenarios, such approaches lead to the emergence of locally distinct configurations that promote an IP.



Following the path based on coordination fluctuations, several authors have modified the initial mean-field theory~\cite{phillips_jncs1979,thorpe_jncs1983} to account for structural adaptation or self-organization in order to avoid stress from additional cross-linking elements with an energy penalty if stressed rigid elements are considered~\cite{jacobs_jncs2000,micoulaut_prb2003,barre_prl2004,briere_pre2007,moukarzel_pre2013}. This adaptation has been also revealed from molecular simulations showing that rigid angles will soften (adapt) with increasing stress/pressure in order to reduce the global increase of the constraint density upon cooling the melt~\cite{bauchy_natcom2015}. In spirit, these models (built from lattices~\cite{jacobs_jncs2000,moukarzel_pre2013}, clusters\cite{micoulaut_prb2003,micoulaut_prb2006}, spin cavity expansions~\cite{barre_prl2004}), share in common that the addition of bonds in a network with increasing average connectivity will be accepted only if this leads to isostatically rigid regions. Once the network has undergone percolation of rigidity at a first transition at $\bar r_{c1}$, an unstressed (isostatic) structure dominates, and characterizes the IP. The addition of redundant bonds then will contribute to the occurrence of stressed rigid regions clusters that finally percolate at a second transition at $\bar r_{c2}$, identified with a stress transition. Both transitions then define a window in connectivity $\Delta \bar r$=$\bar r_{c2}-\bar r_{c1}$ (Figure~\ref{thorpi_fig}).


Instead, using models of soft spheres in the context of jamming transitions, different authors~\cite{yan_prl2014,yan_natcom2018} have suggested that the IP could result from the presence/absence of weak noncovalent interactions in addition to harmonic elastic interactions. In a strong force r\'egime, an IP can be obtained that is characterized by an isostatic cluster spanning the whole system, driven by fluctuations in coordination. However, when weak interactions are present, these features vanish below a certain temperature suggesting that the transitions become mean-field at low temperature and coalesce. These results are partially supported by MD simulations~\cite{bauchy_natcom2015} taking into account long-range interactions (Coulomb, Van der Waals) allowing to probe the weak-force r\'egime.



\subsection{Reversibility windows}

\begin{figure}[b]
\includegraphics[width=0.47\linewidth]{review_ip1.pdf}
\includegraphics[width=0.45\linewidth]{review_ip2.pdf}
\caption{Experimental and numerical verification of anomalies in the intermediate phase region (cyan background) for different glassy systems : (a)~Raman threshold pressure (red) and molar volume (right axis) as a function of constraint density in Ge-Se glasses~\cite{wang_prb2005,bhosle_ijags2012}. (b)~Ionic conduction (red) and dielectric permittivity $\varepsilon(0)$ (right axis) as a function of AgI content $x$ in (1-$x$)AgPO$_3$-$x$AgI glasses~\cite{novita_prl2007,micoulaut_prb2009}. (c)~Calculated enthalpy and volume hysteresis (red, right axis) as a function of pressure in MD simulations of densified silicates~\cite{bauchy_natcom2015}. (d)~Calculated enthalpy and internal network stress (red right axis) as a function of the calculated constraint density in (1-$x$)Na$_2$O-$x$SiO$_2$~\cite{song_prb2020}. The cyan areas represent the reversibility window determined experimentally~\cite{bhosle_ijags2012} or numerically~\cite{bauchy_natcom2015,song_prb2020} in the glass transition region.}\label{ip_exp}
\end{figure}

Experimental support to the predictions of non-mean-field rigidity is obtained from various probes (Figure~\ref{ip_exp}), the main one being of calorimetric nature. Modulated differential scanning calorimetry at the glass transition permits, indeed, to measure an enthalpy of relaxation (a non-reversing heat flow $\Delta H_{nr}$) which exhibits a deep minimum over a finite compositional interval, sometimes exhibiting even sharp boundaries and a square-well behavior with composition defining a reversibility window (RW) for which $\Delta H_{nr}\simeq$0. The latter behavior is thought to be a $T\neq$0 manifestation of the intermediate phase characterizing the $T=0$ disordered network representing the model glass~\cite{micoulaut_advx2016}. There is actually a strong experimental support for these findings connecting RW with the isostatic nature of the network glass structure, and a vast literature has been accumulated on this topic during the last decade for a variety of systems, modified oxides or chalcogenides, which exhibit a salient phenomenology of thermal, relaxation, and vibrational anomalies (Figure~\ref{ip_exp}). A certain number of measured anomalies are provided for Ge-Se (Figure~\ref{ip_exp}a) \cite{bhosle_ijags2012} and AgI-AgPO$_3$ (Figure~\ref{ip_exp}b) \cite{micoulaut_prb2009}. The former indicates that for isostatic compositions, a near stress-free character is suggested from Raman pressure experiments~\cite{wang_prb2005} and manifests by a zero threshold pressure in the IP, as in crystals without residual stresses. Corresponding volumes are found to be minimum and indicate space-filling tendencies (Figure~\ref{ip_exp}a, right axis~\cite{bhosle_ijags2012}). In a very different glass (AgI-AgPO$_3$), the boundaries of the RW~\cite{novita_prl2007} coincide with a marked jump in ionic conductivity that grows exponentially once the network has become flexible (Figure~\ref{ip_exp}b, right axis~\cite{micoulaut_prb2009}). It indicates that local deformation (floppy) modes will promote ion hopping. These changes are accompanied by a maximum in dielectric permittivity (left axis). Glasses with various bonding types display RWs (Figure~\ref{review_ip}), from ionic (silicates~\cite{micoulaut_amin2008}) to iono-covalent, covalent (Ge-Se, \cite{gunasekera_jcp2013}), or semi-metallic (Ge-Te-In-Ag~\cite{varma_ssc2014}). In certain of these systems, e.g. for the simple binary network glasses such as Ge$_x$S$_{1-x}$ or Si$_x$Se$_{1-x}$, the experimental boundaries of the RW are found to be all very close~\cite{gunasekera_jcp2013,selvenathan_prb2000,chakraborty_jpcb2014}, i.e. located between 20~\%~$<x<$~25~\%, and aspects of topology fully control the evolution of rigidity with composition, given that there is a weak chemical effect in case of isovalent Ge/Si or S/Se substitution. This compositional interval defining the RW connects to the mean-field estimate of the isostatic criterion (eq.~\eqref{nc_final}) satisfying $n_c=3$ because coordination numbers of Ge/Si and S/Se can be determined from the 8-${\clN}$ (octet) rule to yield an estimate of the constraints $n_c$=2+5$x$ using eq.~\eqref{nc0}. In fact, for these IV-VI glasses, the lower boundary of the RW ($x_c$=20~\%) coincides with the Phillips--Thorpe~\cite{phillips_jncs1979,thorpe_jncs1983} mean-field rigidity transition $n_c$=3 and $\bar r$=2.4.



\begin{figure}
\begin{center}
\includegraphics[width=0.8\linewidth]{rw2022.png}
\end{center}
\caption{\label{review_ip}
Experimental location of reversibility windows thought to be the manifestation of the topological intermediate phase predicted from models : modified oxides (phosphates~\cite{novita_prl2007,mohanty_front2019}, borates~\cite{vignarooban_epl2014}, germanates~\cite{rompicharla_jpcm2008}, silicates~\cite{micoulaut_amin2008}) and chalcogenides (Si-Ge-Te~\cite{das_ssc2012, gunasekera_jap2014}, Ge-Te-In-Ag~\cite{varma_ssc2014}, Ge-Sb-Se~\cite{gunasekera_jpcb2013}, Ge-P-S~\cite{vempati_jpcm2004}, Ge-P-Se~\cite{chakravarty_jpcm2005}, Ge-As-S~\cite{qu_philmag2005}, Ge-As-Se~\cite{chbeir_front2019}, Ge-S-I~\cite{wang_prl2001}, Ge-Se-I~\cite{wang_jpcm2007}, P-S~\cite{boolchand_jncs2009}, P-Se~\cite{georgiev_epl2003}, As-S~\cite{chakravarty_front2019}, As-Se~\cite{ravindren_jcp2014}). Only a limited cases can be represented as a function of constraint density given the obvious atomic coordinations fulfilling the octet rule : Ge-S~\cite{chakraborty_jpcb2014}, Ge-Se~\cite{bhosle_ijags2012}, and Si-Se~\cite{selvenathan_prb2000}.}
\end{figure}

For most of the systems however, a direct constraint counting cannot be performed because uncertainties persist regarding the relevant coordination numbers so that an enumeration of active/inactive constraints must be derived from specific structural models (see, e.g. an example for silicates~\cite{micoulaut_amin2008}). This feature becomes obvious once Group V selenides/sulphides are being considered (Figure~\ref{review_ip}). Here, different RW locations are found for isovalent compounds, e.g. differences emerge between As- and P-bearing chalcogenides, and between sulphides and selenides (e.g. P$_x$S$_{1-x}$ and P$_x$Se$_{1-x}$, \cite{vempati_jpcm2004}) so that differences can only be explained from structural models validated from e.g. spectroscopic studies~\cite{georgiev_epl2003}. The above statements are also valid in tellurides for which the increased electronic delocalization of the Te atoms lead to mixed local geometries that are now composition dependent~\cite{micoulaut_prb2021}, so that a proper constraint count must rely on molecular simulations, in conjunction with dedicated constraint counting algorithms.



Molecular dynamics simulations have been able to substantiate the characteristics of RWs~\cite{bauchy_natcom2015,mantisi_prb2015,song_prb2020,trachenko_prl2004}. Using a similar methodology as the one used in experimental calorimetry, cooling-heating cycles have been performed in pressurized silicates (Figure~\ref{ip_exp}c~\cite{bauchy_natcom2015}) and sodium silicates ((Figure~\ref{ip_exp}d). Concerning the latter, the hysteresis resulting from enthalpic relaxation during a numerical cooling-heating cycle has been found to be minimized for 12~\%$\leq x\leq$ 20\% Na$_2$O (Figure~\ref{ip_exp}d~\cite{song_prb2020}), which echoes with the experimental observation of the RW (18~\%$\leq x\leq$ 22\% Na$_2$O~\cite{vaills_jpcm2005}). For both investigated silicates, the location of the RW corresponds to an isostatic character of the network that is independently computed. In addition, a certain number of other anomalies have been found from simulations, such as minimal activation energies in viscosity and diffusivity~\cite{micoulaut_pssb2013}, minimal relaxation time~\cite{bauchy_natcom2015}, growth of internal stress (Figure~\ref{ip_exp}d, right axis~\cite{song_prb2020}) with a marked jump at the stress transition.

\subsection{The structural signature of the intermediate phase}

In the present work, we examine structural features that could be associated with the onset of stressed rigidity or the IP. The question continues to be actively debated in the literature given the rather elusive structural signatures of rigidity transitions in spite of thermal and spectroscopic evidence. Direct structural signatures have been searched by different authors. For instance, X-ray diffraction in Ge-Se glasses~\cite{wang_jncs2004,sharma_physicab2005} have revealed that the area and the inverse position of the first sharp diffraction peak of the structure factor $S(k)$ exhibited a plateau found in the same compositional window as the reported IP, i.e. suggesting some effect of intermediate range order for isostatic compositions. These conclusions were not confirmed by a similar study~\cite{sartbaeva_prb2008} and the reported x-ray absorption fine structure analysis furthermore found no correlation between the short-range order and the IP, as also evidenced from various NMR studies~\cite{lucas_prb2009,gjersing_jpc2010}. Another x-ray absorption experiment on Ge-Se (near-edge structure) on Ge-Se glasses led to an opposite conclusion~\cite{chen_apl2010} as spectra revealed compositional plateaus coinciding with the location of the RW, providing structural evidence for the IP. More recently, the possibility of a structural signature of the IP in Ge-Se has been examined from high-resolution neutron diffraction to measure different quantities associated with topology~\cite{zeidler_front2017}, and results did not point to an obvious structural origin of the IP. For instance, they did not confirmed the possible deviation of atomic coordinations from the 8-N rule reported from first-principles molecular dynamics simulations~\cite{inam_jpcm2007}.


Molecular simulations have, indeed, attempted to characterize the origin of the IP. On the archetypal Ge-Se system, first principles simulations have shown that isostatic glasses originate from a competition between amorphous GeSe$_2$ and amorphous Se clusters to produce measurable signatures in electronic properties~\cite{inam_jpcm2007,chen_apl2010}. On the same system, it was shown that while no compelling evidence for signatures were found in structure functions (pair correlation $g(r)$, structure factor $S(k)$, in contrast with simulated As-Se glasses~\cite{bauchy_prl2013b,bauchy_jcp2014}), constraint-counting algorithms showed that broken bond-bending constraints are associated with the stressed-rigid phase, a result that was also recovered for As-Se~\cite{bauchy_jcp2014} and Ge-S glasses~\cite{chakraborty_prb2017}. The softening of bending constraints under increased applied stress \emph{via} cross-linking or applied pressure turns out to be a specific feature in stressed rigid systems and this feature has been also found in different densified tetrahedral liquids~\cite{pacaud_jcp2015,bauchy_prl2013a}.



\begin{figure}[t]
\begin{center}
\includegraphics[width=0.9\linewidth]{FSDP.pdf}
\end{center}
\caption{\label{fsdp}
(a) Position k$_{\FSDP}$ of the FSDP (in \AA$^{-1}$) of different partial structure factors as a function of the calculated constraint density $n_c$. (b) Width $\Delta$k$_{\FSDP}$ of the FSDP (in \AA$^{-1}$). Adapted from Ref. \cite{micoulaut_pssb2013}}.
\vspace*{-15pt}
\end{figure}

The most noticeable structural anomaly related with onset of rigidity has been determined for selected systems only. While principal peak characteristics (position, width) of total scattering functions do not display any specific trend related to the underlying phases~\cite{zeidler_front2017,bychkov_prb2005}, partial correlations in Fourier space reveal that some change in structural ordering takes place as networks become rigid. This is exemplified in densified 2SiO$_2$-Na$_2$O for which the FSDP found at ambient pressure (at wavevector position $k_{FSDP}$=1.53~\AA$^{-1}$) evolves monotonically up to 1.9~\AA$^{-1}$ with pressure/rigidity (black curve, Figure~\ref{fsdp}a). The behavior contrasts with the one calculated for corresponding network partials (Si-O, O-O~\cite{micoulaut_pssb2013}) which peak at $n_c\simeq$3. Corresponding trends are obtained for the width of the FSDP (Figure~\ref{fsdp}b). As the position of the FSDP reflects some repetitive characteristic distance between structural units, the result indicate that an intermediate range order lengthscale of distance 7.7/$k_{FSDP}$ emerges, and then decreases. Here the factor 7.7 relates to the first maximum of the spherical Bessel function $j_0(kr)$. Similarly, the broadening of the FSDP is indicative of a correlation length, following the well-known Scherrer equation for microcrystals which connects the width $\Delta k_{FSDP}$ of a Bragg peak with the average size of the microcrystals. Similar features have been obtained for partials in the As-Se system~\cite{bauchy_prl2013b,bauchy_jcp2014}. While, the total structure factor characteristics did not showed any specific trend across the rigidity transitions as in experiment~\cite{bychkov_prb2005}, it was found that $k_{FSDP}$ of the As-Se partial structure factor displayed a threshold behavior at the Se-rich boundary of the IP, whereas $\Delta k_{FSDP}$ displayed a minimum value in the IP. These findings indicate that structural features might be material specific, and the question of structural correlations in relationship with flexible to rigid transitions continues to be questioned.


Here, inspired by the work of Wilson~\cite{wilson_pccp2012} and using large scale molecular dynamics we are able to probe such structural correlations up to 20-25~\AA\ in a network (sodium silicate) that is known~\cite{bauchy_natcom2015} to undergo flexible-intermediate-stress transitions/thresholds, to display a RW (Figure~\ref{ip_exp}c) and a series of anomalies in dynamic heterogeneities~\cite{bauchy_epl2013,bauchy_prl2017}. By focusing on structure functions describing topological ordering at long-range (Bhatia--Thornton formalism), we show that a typical lengthscale of about 3.5~\AA\ characterizing the decay of correlations emerges once the system becomes rigid. Similarly, the wavevector characterizing the long-range oscillations of topological ordering is found to be minimal for a calculated $n_c\simeq$ 3 which signals that periodic ordering takes place on larger lengthscales in isostatic networks. Since we are focusing on structural correlations, we also calculate the pair correlation entropy which measures excess contributions with respect to a non-structured reference state, and maximizes for $n_c\simeq$3. Taken together, these results indicate that isostatic systems bear striking similarities with densified tetrahedral liquids displaying water-like anomalies, i.e. correlated extrema in transport coefficient (diffusivity), thermodynamics, excess entropy and structural features. These commonalities are discussed and contrasted.

\section{Simulation methods}

\subsection{System characteristics}

The system under investigation is a densified sodium silicate supercooled liquid (2SiO$_2$-Na$_2$O, or NS2) that has been investigated from classical MD simulations (integration time using the Verlet algorithm with a 1~fs time step). The atoms interact with a rigid ion two-body potential parametrized by Teter~\cite{cormack_jncs2003} of the form :
\begin{equation}\label{eq1}
\Phi_{ij}(r)=\Phi_{ij}^{SR}(r)-{\frac {C_{ij}}{r^6}}+{\frac {q_iq_j}{r}}
\end{equation}
where $\Phi_{ij}^{SR}(r)=A_{ij}\exp{-\rho_{ij}/r}$ is a short-range term arising from electronic repulsion between atoms. The accuracy of the potentials has been demonstrated in various studies on structural~\cite{du_jncs2004,krishnan_actam2019,bauchy_jcp2014b,du_jncs2006} and dynamic properties. Noteworthy is the fact that calculated diffusion and viscosity are comparable to experimental data~\cite{bauchy_chemgeol2013} so that the efficiency is substantially improved with respect to alternative potentials~\cite{kramer_jacs1991}. The cut-off for both Buckingham and Coulomb part of the Teter potential has been set to 18~\AA.

\begin{figure}[!htbp]
\includegraphics[width=0.7\linewidth]{snap.png}
\vspace*{-14pt}
\caption{A snapshot of the simulated system: 16767 atoms representing the NS2 liquid (red: silicon, yellow: oxygen, blue: sodium). At zero pressure, the cell length is 65~\AA.}\label{snap}
\end{figure}

In order to investigate long-range information of the partial pair correlations functions, the initial configuration has been obtained from a rather large crystalline disilicate sample consisting of 16767 atoms (1863 Na$_2$Si$_2$O$_5$ units). The system has been melted at an equilibration stage at 4000 K over 1 ns, prior to sequential quenches at 3500~K, 3000~K, and finally the target temperature of 2000~K corresponding to the numerical glass transition region where the reversibility window has been obtained numerically~\cite{bauchy_natcom2015}. The resulting structure (Figure~\ref{snap}) has a cell length at zero pressure of about $L$=65~\AA\ which permits to safely consider distances up to 20~\AA$<L$/2 from simulations up to the most elevated pressures. Densification has been achieved by increments up to 20.12~GPa in NPT Ensemble ($L$=55.39~\AA).



\subsection{Molecular dynamics based rigidity}

In order to connect the structure to features of rigidity, and particularly the constraint density $n_c$ with pressure and temperature, we enumerate mechanical constraints from the atomic scale time dependent trajectories obtained at a given thermodynamic condition (2000~K, $P$). Previous applications of this methodology have been performed for the particular NS2 system~\cite{bauchy_natcom2015,bauchy_jncs2011}, but also on a variety of different chalcogenide network glasses~\cite{chakraborty_prb2017,micoulaut_prb2021,micoulaut_prb2013} and liquids~\cite{yildirim_natcom2016}.


To determine the number of BS constraints, one focuses on neighbour distribution functions 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_c^{BS}$=$r_i$/2. This number usually correspond to small radial motion between two neighboring atoms.


To determine BB constraints, one uses partial bond angle distributions (PBADs) $P({\theta_{ij}})$ which split the usual bond angle distribution into partial contributions defined by a central atom $0$ and the N first neighbours which define N(N-1)/2 possible triplets or angles $i0j$ ($i$=1..N-1, $j$=2..N), i.e. 102, 103, 203, etc. The standard deviation $\sigma^{ij}_\theta$ of each distribution P($\theta_{ij}$) gives a quantitative estimate of the angular excursion around a mean angular value, and provides an indication of the bond-bending strength. Small values for $\sigma^{ij}_{\theta}$ correspond to an intact bond-bending constraint which maintains a rigid angle at a fixed value, whereas large $\sigma_{\theta}^{ij}$ correspond to a bond-bending weakness giving rise to an ineffective or broken constraint~\cite{bauchy_jncs2011}.


The use of such standard deviations permits to derive the constraint density under various thermodynamic conditions $n_c(x,T,P)$ (x-axis of Figure~\ref{ip_exp}d) without any need of dedicated phenomenological structure models. Model potentials (such as eq.~\eqref{eq1}) are validated from a successful comparison of structure functions ($g(r)$, $S(k)$) obtained from scattering experiments.

\section{Results}
\subsection{Bhatia--Thornton formalism and the decay of correlation functions} In order to concentrate on long-range topological ordering, we use the Bhatia--Thornton (BT) formalism~\cite{bhatia_prb1970} which focuses on the number-number correlations (via a corresponding structure factor $g_{NN}^{BT}(r)$) whose definition for binary mixtures can be easily extended to the present ternary (Si,O,Na) system (see also Ref. \cite{grosdidier_jncs2021}):
\begin{equation}
g_{NN}^{BT}(r)=\sum\limits_{i,j}x_ix_jg_{ij}(r)
\end{equation}
where $g_{ij}(r)$ represent the different calculated partial pair correlation functions (not shown). In scattering experiments, this function (or its counterpart $S_{NN}^{BT}(k)$ in Fourier space) probes correlations that are independent of the chemical nature of the scattering centers and, therefore, provides a measure on topological ordering at intermediate and extended lengthscales\cite{salmon_nature2005,salmon_jpcm2007} depending on position/width of the FSDP and principal peak, respectively. Other partials can be derived within this formalism and these focus on e.g. how chemical species are distributed over the scattering centers and might provide information on chemical ordering, i.e. the preference to form unlike bonds (as Ge-Se bonds in Ge-Se glasses~\cite{salmon_nature2005}). For a ternary alloy, concentration-concentration correlations of the silica network can e.g. be characterized via:
\begin{multline}\label{eqbt}
x_{Se}^{-1}g_{CC(SiO)}^{BT}(r)=x_{Si}+x_{O}+x_{Na}\biggl[x_{Si}^2g_{SiSi}+x_{O}^2g_{OO} +2x_{Si}x_{O}g_{SiO}-2x_{O}(x_{Si}+x_{O})g_{ONa}\\
-2x_{Si}(x_{Si}+x_{O})g_{SiNa}+(x_{Si}+x_{O})^2g_{NaNa}\biggr]
\end{multline}
with $x_{Na}$=$x_{Si}$=0.22 and $x_O$=0.56 given by the stoichiometry of NS2. Note that for a three component system there are three different concentration-concentration pair correlation functions and $g_{CC(SiNa)}^{BT}(r)$ and $g_{CC(ONa)}^{BT}(r)$ are obtained by a cyclic operation on (Si,O,Na). In the high wavevector region $r$ all concentration-concentration structure factors converge to the limit $x_i(x_j+x_k)$ ($i\neq j,k$) which is a direct consequence of the definition of $g_{CC(SiO)}^{BT}(r)$, $g_{CC(SiNa)}^{BT}(r)$ and $g_{CC(ONa)}^{BT}(r)$ (eq.~\eqref{eqbt})). We won't comment more on additional correlations (number-concentration, see Ref. \cite{grosdidier_jncs2021}).

\begin{figure}[!htbp]
\includegraphics[width=0.5\linewidth]{Bhatia.pdf}
\caption{Plot of the Bhatia--Thornton Number-Number pair correlation function $rh_{NN}(r)$ (a) for three different pressures in NS2 : total contribution (black curve) and network contribution (SiO$_2$, red curve) for $P$=0~GPa, $P$=1.96~GPa and $P$=18.80~GPa. (b) The representation $\vert rh_{NN}(r)\vert$~\cite{evans_jcp1994} permits to highlight the exponential decay of the partial pair correlation function (same pressures and same color code for the contributions).}\label{BT}
\vspace*{-8pt}
\end{figure}

Figure~\ref{BT}a now represents $G_{NN}(r)$=$rh_{NN}(r)$=$r[g_{NN}(r)-1]$ for three different pressures corresponding to the flexible phase ($P$=0, $n_c$=2.92), to an isostatic condition ($P$=1.96~GPa, $n_c$=3.01) and to the stressed rigid phase ($P$=18.80~GPa, $n_c$= 3.59, see Figure~\ref{ip_exp}c). Usual features of such functions are obtained and manifest by a very intense principal peak ($r$=1.65~\AA) arising from the Si-O correlations of the base silica network, together with a secondary peak found at $r$=2.81~\AA\ which is mainly due to O-O and Si-Si correlations, the alkali related correlations leading to much broader distributions with principal peaks located at 2.44~\AA\ and 3.66~\AA\ for Na-O and Na-Si correlations, respectively. In the region around $r\simeq$5.0~\AA, the pair distribution function $rh_{NN}(r)$ is dominated by a double peak structure that is due to Si-O secondary correlations and a superposition of different partials (O-O and Si-Si), respectively. Upon increasing pressure, the principal peak shifts to larger distances (1.71~\AA\ for 18.80~GPa) consistently with the usual tetrahedral to octahedral conversion of Si which involves a lengthening of the Si-O bonds~\cite{bauchy_chemgeol2013}. Noteworthy is the merging of the double peak at $r\simeq$5~\AA\ which is produced by a near equivalent contribution of Si-O, Na-O and Si-Si partials at high pressure.


In order to highlight the asymptotic decay of $rh_{NN}(r)$, we represent on the right side of Figure~\ref{BT} the quantity $\ln\vert rh_{NN}(r)\vert$. While correlations become noisy at long-range distance due to the limited sampling of the structure at long distances ($r\geq$20~\AA), a clear signature of an exponential decay $\exp[-r/\xi]$ is acknowledged (blue curve, bottom panel in Figure~\ref{BT}b). This feature is present for the whole system but also if one focuses only on the underlying Si-O subnetwork red curves). The asymptotic decay of total correlation functions has been interpreted by Ma et al.~\cite{ma_natmat2009} as the result of a self-similar packing of atomic clusters which leads to a medium range order having the characteristics of a fractal network of dimension~$D_f$. From a variety of metallic glasses, the asymptotic decay ($r>$6~\AA) could be fitted using the functional form :
\begin{equation}
G(r)-1=Ar^{D_f-D}\exp[-r/\xi]cos\left(q_1r+\phi\right)
\end{equation}
with $D_f-D$ = $D_f-3$ = 0.69 in 3D, and $\xi$ a cutoff length characterizing the finite cluster size and possible cluster-entanglement effects~\cite{freltoft_prb1986}. Here, the sinusoidal function $sin(q_1r+\phi)$ was introduced to describe the oscillatory correlation of $G(r)$ with $q_1$ the FSDP position. A similar form inspired from the Ornstein--Zernike (OZ) theory of liquids~\cite{hansen} has been proposed for network glasses by Salmon and co-workers~\cite{salmon_jpcm2007,zeidler_pccp2013} without invoking the need of a fractal dimension, and has been successfully applied to network glasses and liquids. The approach builds on a simple rigid-ion pair potential model of the form given in eq.~\eqref{eq1} and an exact solution indicates that at long distance $rh_{NN}(r)\rightarrow r^{-5}$ in absence of the dispersion term ($C_{ij}$=0 in eq.~\eqref{eq1}). A pole analysis~\cite{carvalho_molphys1994,evans_jcp1994} of the $k$ space solutions of the OZ equations then leads for the particular case of dense liquids to an exponential decay of the form:
\begin{equation}\label{oz}
G_{NN}=rh_{NN}(r)=2A_{NN}\exp[-r/\xi]cos\left(q_1r+\phi_{NN}\right)
\end{equation}
with parameters $\xi$ and $q_1$ that are specific to the fitted functions ($h_{NN}$, $h_{CC}$ and $h_{NC}$) in several glasses and liquids (GeO$_2$, GeSe$_2$, ZnCl$_2$~\cite{salmon_jpcm2006}). Importantly, it was found that in such tetrahedral systems $q_1$ is not related to the FSDP position. In the present study, while dispersion terms are present in O-O, O-Si and O-Na interactions~\cite{cormack_jncs2003}, the functional form proposed in eq.~\eqref{oz} still provides an interesting guide for examining structural correlations (see however, the effect of dispersion forces in Ref. \cite{kjellander_jpa2005}). Secondly, although the functional forms have been proposed for binary mixtures~\cite{grodon_molphys2005} and since the Bhatia--Thornton pair correlation function $rh_{NN}(r)$ also displays an obvious asymptotic decay together with an oscillatory behavior, it is tempting to use the proposed functional forms to extract an information on topological ordering with respect to the rigidity status of the considered systems.

\begin{figure}[!h]
\includegraphics[width=0.48\linewidth]{fit.pdf}
\vspace*{-5pt}
\caption{Calculated Bhatia--Thornton NN pair correlation function $G_{NN}(r)$=$r[g_{NN}(r)\linebreak-1]$=$rh_{NN}(r)$ for $P$=4.0~GPa (circles), together with a fit from eq.~\eqref{oz} (red curve). The inset shows the entire function with the fit, and the gray zone defines the $r-$range over which the fit has not been performed.}\label{fitbt}
%\vspace*{-15pt}
\end{figure}

An inspection of Figure~\ref{BT}b shows repeated maxima obeying a straight line $ln\vert rh_{NN}(r)\vert\linebreak\simeq -r/\xi+\textit{constant}$ and clearly indicates an exponential decay of the oscillations. It, therefore, represents a convenient starting point for discussing extended-range ordering, as also performed experimentally~\cite{zeidler_pccp2013}. The fit using eq.~\eqref{oz} on the calculated functions $h_{NN}(r)$ (Figure~\ref{fitbt}) now suggests evidence that at large distance (here $r\geq$7~\AA) the damped oscillations of topological correlations can be reproduced using the OZ form. For the represented system at $P$=4.0~GPa, we find $\xi$=2.86~\AA\ and $q_1$=2.61~\AA$^{-1}$. These values already indicate that $q_1$ cannot be related to the FSDP position $k_{FSDP}$ as suggested by Ma et al.~\cite{ma_natmat2009} given our calculated $k_{FSDP}$=1.95~\AA$^{-1}$~\cite{micoulaut_pssb2013}. Instead, $q_1$ is rather linked to the principal peak position found at $k_{PP}$=2.98~\AA$^{-1}$, and this suggested correlation appears to hold for the experimentally investigated glasses (e.g. $q_1$=2.10~\AA$^{-1}$ against $k_{PP}$=2.04~\AA$^{-1}$ and $k_{FSDP}$=1.14~\AA$^{-1}$ in glassy GeSe$_2$~\cite{salmon_jpcm2006,salmon_jpcm2007}). Similarly, the obtained $\xi$ values for the present silicates appear to be close to those obtained experimentally (e.g. $\xi$=3.00~\AA\ in GeO$_2$~\cite{salmon_jpcm2006}, represented in Figure~\ref{xi_vs_nc}b).
\vspace*{15pt}

\subsection{Threshold behavior in topological ordering}

The resulting parameters $\xi$ and $q_1$ are now represented in Figures~\ref{xi_vs_nc} and~\ref{q1_vs_nc} as a function of applied pressure (panels a) or calculated constraint density $n_c$ (panels b).

\begin{figure}[!htbp]
\includegraphics[width=0.5\linewidth]{parameters_fit1.pdf}
\caption{Distance $\xi$ obtained from the fit of the Bhatia--Thornton pair correlation function $rh_{NN}(r)$ (Figure~\ref{fitbt}): (a) $\xi(P)$, (b) $\xi(n_c)$. The gray zone represents the reported numerical reversibility window~\cite{bauchy_natcom2015}. Green symbols represent experimental data for GeO$_2$ and GeSe$_2$~\cite{salmon_jpcm2006}.}\label{xi_vs_nc}
\end{figure}

\begin{figure}[!htbp]
\includegraphics[width=0.5\linewidth]{parameters_fit.pdf}
\caption{Wave-vector $q_1$ obtained from the fit of the Bhatia--Thornton pair correlation function $rh_{NN}(r)$ (Fig. \ref{fitbt}): (a) $q_1(P)$, (b) $q_1(n_c)$. The gray zone represents the reported reversibility window~\cite{bauchy_natcom2015}. Green symbol represent experimental data for GeO$_2$~\cite{salmon_jpcm2006}. The data point for GeSe$_2$ is out of range (2.10~\AA$^{-1}$).}\label{q1_vs_nc}
\end{figure}

When represented as a function of pressure, $\xi(P)$ leads to a smooth increase (Figure~\ref{xi_vs_nc}a), and a near constant behavior if only the SiO$_2$ subnetwork is considered. This typical distance converges to a limiting value of about $\simeq$3.5~\AA\ (for the total network) that is also close to the one determined for the stressed rigid system GeSe$_2$~\cite{salmon_jpcm2006} for which a mean-field estimate leads to $n_c$=3.67~\cite{micoulaut_advx2016}. On the other hand, the other experimental data point (GeO$_2$, $n_c\simeq$3~\cite{bauchy_prb2011}) is compatible with the obtained trend in $\xi(n_c)$. Conversely, a minimum occurs for the parameter $q_1$ that characterizes the periodicity of topological ordering (Figure~\ref{q1_vs_nc}a). This minimum (less pronounced for the underlying silica network) is obtained at the flexible-intermediate boundary of the IP previously characterized~\cite{bauchy_natcom2015}, whereas for intermediate and stressed rigid system, a linear increase of $\xi(P)$ is acknowledged, suggesting a reduction of the periodicity of topological ordering upon further compression.

The obtained behavior obviously suggests a structural signature for the IP as also previously emphasized for peak characteristics in Fourier space at low wavector~\cite{micoulaut_pssb2013}. These signatures are essentially obtained once all Faber--Ziman or BT partial correlations are being considered (Figure~\ref{fsdp}) and are, therefore, barely visible from total correlations accessed from X-ray and neutron scattering studies without isotopic substitution~\cite{zeidler_front2017,sartbaeva_prb2008}. Parameters of the FSDP region display, indeed, an anomalous behavior in some of the partial structure factors. The position $k_{FSDP}$ maximizes for a select number of partials only (Si-O, O-O Figure~\ref{fsdp}a), whereas the width also show anomalies (minimum) close to the isostatic system.
\vspace*{50pt}

\break

\section{Discussion}

We relate the anomalies in topological ordering in isostatic systems with a measure of structural disorder, and with diffusivity extrema which underscore water-like anomalies found in densified tetrahedral liquids.

\subsection{Entropy}

An additional and interesting insight is, indeed, provided by the pair correlation entropy which measures from the atom-atom radial distribution functions $g_{ij}(r)$ representing an excess entropy with respect to the entropy of an ideal gas. The latter consists of a non-interacting mixture of particles having the property $g_{ij}(r)$=1 over all distances
\begin{equation}\label{entro}
S_2/k_BN=-2\pi\rho_0\sum\limits_{i,j}x_ix_j\int_0^\infty\biggl[g_{ij}(r)\ln g_{ij}(r) - \left[g_{ij}(r)-1\right]\biggr]r^2dr
\end{equation}
where $x_i$ is the mole fraction of component $i$ and $\rho_0$ is the number density of the system. The latter quantity permits to capture effects that departs from form random liquid networks and show interesting correlations with calculated thermodynamic and transport extrema occurring in densified tetrahedral liquids, known as ``water-like'' anomalies. These simply signal deviations from regular liquid behaviour that the excess entropy is able to capture, and are characterized by sets of state points ($T$,$\rho$) or ($T$,$P$) for which the diffusivity increases with density/pressure.


Figure~\ref{entropy}a represents the calculated entropy for the present system, as a function of the constraint density $n_c$. Unlike the monotonic decrease in $S_2$ with density $\rho$ as in simple liquids, the present system has a well-defined anomalous regime where the excess entropy rises with increasing density (Figure~\ref{entropy}a) or pressure (not shown). The increase of $S_2$ is large in flexible networks with increasing connectedness, correlated with the increase in density (right axis), and reaches a maximum close at the isostatic condition, while for $n_c\geq$3 entropy decreases moderately. These are linked with other anomalies noticed in transport coefficients when tracked as a function of $P$ or $\rho$ or $n_c$~\cite{bauchy_chemgeol2013,mantisi_prb2015}. Here oxygen diffusivity (and silicon, not represented) displays a minimum $D_m$ at $n_c$=2.85 and a maximum $D_M$ at $n_c\simeq$3, whereas Na is merely constant in the flexible phase (Figure~\ref{entropy}b), exhibits a dramatic jump in the identified intermediate phase~\cite{mantisi_prb2015} and then decreases in a similar fashion to oxygen in the stressed rigid phase corresponding to the high pressure NS2 system.

\begin{figure}[tbp]
\includegraphics[width=0.7\linewidth]{entropy.pdf}
\caption{(a) Calculated entropy in densified NS2 (eq.~\eqref{entro} as a function of calculated constraint density $n_c$. The calculation is performed either on the total network (black) or on the underlying silica network (red). The right axis correspond to the calculated system density (blue symbols). (b) Calculated oxygen and sodium diffusivity as a function of the calculated $n_c$. Transport properties are obviously affected by the elastic nature of the liquid. Adapted from Bauchy et al.~\cite{bauchy_chemgeol2013} and Mantisi et al.~\cite{mantisi_prb2015}. The location of the reversibility window (cyan) has been determined from NVT cooling-heating cycles~\cite{mantisi_prb2015}.}\label{entropy}
\end{figure}

A closer inspection of the different Figures~\ref{q1_vs_nc} and~\ref{entropy} indicates that the structural properties not all track the {\em locus} of the IP. Parameters characterizing the topological ordering ($q_1$ and $\xi$) appear to display a change in behavior as the system is exactly isostatically rigid ($n_c$=3), i.e. at the flexible to intermediate boundary, whereas structural partial correlations involved in the first sharp diffraction peak merely follow the location of the IP region. On the other hand, previous investigations have shown that the change in dynamic properties (diffusivity) takes place throughout the IP. These features signal that the signatures of the two IP boundaries are essentially detected from calculated dynamic quantities (diffusivity~\cite{bauchy_natcom2015,mantisi_prb2015}, viscosity~\cite{bauchy_chemgeol2013}, enthalpy of relaxation~\cite{song_prb2020}) or from experimental thermal/calorimetric measurements~\cite{micoulaut_advx2016}.

\subsection{Water-like anomalies}

The present behaviors obtained for pair correlation entropy and diffusivity correspond to salient features of water-like anomalies which consist in a series of dynamic/thermodynamic anomalies that manifest by diffusivity maxima and minima in densified tetrahedral liquids such as water~\cite{errington_nature2001}, SiO$_2$~\cite{shell_pre2002} or GeO$_2$~\cite{pacaud_jcp2015,jabes_jcp2010}. Such anomalies have been found in NS2 as well, and extrema in diffusivities with density or pressure for the network forming species (Figure~\ref{entropy}b and~\cite{bauchy_prl2013b,mantisi_prb2015}) or viscosity~\cite{bauchy_chemgeol2013} have been reported. In addition, when both diffusivity and viscosity are studied as a function of temperature in an Arrhenius plot, it has been found that corresponding activation energies for viscosity or diffusion also display a minimum for the same range of pressures/densities where the diffusivity anomalies are found~\cite{bauchy_natcom2015}. The latter behavior indicates the increased ease of relaxation in isostatic network-forming liquids that contributes to heterogeneous dynamics~\cite{bauchy_prl2017} and to the nearly reversing character of the glass transition (Figure~\ref{ip_exp}c).


In studies on structural and dynamic anomalies of densified tetrahedral liquids, it has been stressed\cite{errington_nature2001,shell_pre2002,jabes_jcp2010} that the definition of local structural order parameters could help in understanding the relationships between such diffusivity anomalies, structural and thermodynamic anomalies under temperature and density change. One order parameter (translational) measures the tendency of pairs of molecules to be separated by a preferential distance, while a second order parameter (orientational) measures the tendency of a molecule and its nearest neighbours to adopt preferential orientations.


Here, one now recognizes that the isostatic nature of the network leads to structural correlations which manifest by a certain degree of topological ordering characterized by a decay length $\xi$ and a period $q_1$ that is linked with the presence of an entropy maximum for $n_c\simeq$3. The diffusivity anomalies are driven by constraint $n_c$ softening which is the dominant feature controling the evolution of transport coefficient under density and temperature change~\cite{bauchy_prl2013a}. We, thus, view these transport anomalies as a consequences of structural rearrangements driven by stress adaptation in a connectivity window located between roughly 2.85$\leq n_c\leq$3.0, and the typical features of diffusivity minima $D_{m}$ and maxima $D_{M}$ can be related to the boundaries of the RW (Figure~\ref{entropy}b) \cite{mantisi_prb2015,bauchy_natcom2015}. In light of these correlations, we now interpret the location of $D_{m}$ as the boundary for the onset of a rigid but stress-free network-forming liquid, whereas the location of $D_{M}$ is related to the upper boundary of the IP. The trend with $n_c$ (Figure~\ref{entropy}b) also reveals that local deformation modes which are present in the flexible phase will promote transport. Once the system is becoming more rigid, atomic motion and diffusivity reduce. In the IP, the adaptive nature of the structure (constraint softening~\cite{bauchy_natcom2015,bauchy_prl2013a}) will facilitate transport that is also induced by a reduction of stress~\cite{song_prb2020} (no redundant bonds/constraints), prior to an important decrease once the system has become stressed rigid and the network is locked by an important bond density.

\section{Conclusion}

Here we have reviewed some features regarding flexible to rigid transitions in glasses and disordered networks. The Phillips--Thorpe mean-field theoretical framework inspired by the Maxwell treatment on trusses has been used for decades in glass science. More recent non-mean-field scenarios indicate that measured reversibility windows at the glass transition might be the signature of the ability of networks to adapt under a stress imposed by composition or pressure, that lead to a so-called intermediate phase. The picture of flexible to rigid transition appears, thus, to be more complex than previously believed, and connects to certain anomalous features typical of densified liquids. Molecular dynamics simulations provide an interesting added value, able to substantiate these key results derived from phenomenological models. Simulated sodium silicates appear as attractive systems in this context given that they undergo a rigidity transition with sodium content, and display with pressure a variety of anomalies that connect to the presence of a reversibility window, and an intermediate phase. Here we have focused essentially on possible structural signatures of the IP, and the analysis of the decay of Bhatia--Thornton pair correlation functions with pressure indicates that typical lengthscales for topological ordering display anomalies (thresholds, minima) across the IP or its related boundaries. These results obtained from number-number correlations might suggest that structural signatures of the IP are only visible from partial correlations, as also previously found in parent systems~\cite{bauchy_prl2013b,micoulaut_pssb2013}.


The NS2 liquid displays striking similarities with densified tetrahedral liquids~\cite{errington_nature2001,shell_pre2002,jabes_jcp2010} which display similar anomalous thermodynamic and kinetic properties, the most obvious one being the existence of a regime of anomalous density behaviour where diffusivity and entropy maximizes. In contrast to those liquids which do not contain ions (i.e. Na), no negative isothermal expansion coefficient is found that give rise to temperatures of maximum density in appropriate density-temperature or temperature-pressure planes~\cite{jabes_jcp2010}. Instead, one acknowledges an obvious threshold behavior once $\rho_0$ is tracked with the constraint density $n_c$ (Figure~\ref{entropy}b) which only underscores the enhanced compaction tendency of the flexible phase and the rapid evolution of $\rho_0$ with increasing rigidity.


These anomalies not only correspond to a set of thermodynamic state points for which e.g. the diffusivity increases with density but also underscores the role played by network connectivity and structure. A unified explanation~\cite{sharma_jcp2005,agarwal_jcp2007} for such thermodynamic and kinetic anomalies has been proposed on the basis of the existence of an excess entropy anomaly involving a rise in excess entropy on isothermal compression. The collection of numerical data on NS2, once represented as a function of constraint density $n_c$ highlights now the central role of network rigidity and provides some indication that stretching and bending interaction softening~\cite{bauchy_natcom2015,bauchy_prl2013a} will control the various obtained anomalies including the one detected in pair correlation entropy, while also driving the onset of a near reversible glass transition where the enthalpy of relaxation is minimal. This is the dominant feature of the intermediate phase. Ultimately, these features give rise to some structural ordering that manifests once long-range correlations of the network topology are investigated as a function of the order parameter of the flexible to rigid transition, i.e. the fraction of floppy modes 3-$n_c$.

\section*{Conflicts of interest}
The author has no conflict of interest to declare.


\bibliographystyle{crunsrt}
\bibliography{crphys20220570}
\end{document}
