Mechanical Failure in Amorphous Solids: Scale Free Spinodal Criticality

TL;DR

This paper proposes that the scale-free criticality observed in the mechanical failure of amorphous solids is due to a spinodal point of a thermodynamic phase transition, characterized by a diverging correlation length and system-spanning instabilities.

Contribution

It introduces a universal theoretical framework linking critical failure in amorphous solids to a spinodal point, with new order parameters and correlation functions applicable across systems.

Findings

Critical exponents for correlation length divergence estimated.

Universal order parameter and correlation functions proposed.

Spinodal criticality observed in athermal systems and analyzed in thermal systems.

Abstract

The mechanical failure of amorphous media is a ubiquitous phenomenon from material engineering to geology. It has been noticed for a long time that the phenomenon is "scale-free", indicating some type of criticality. In spite of attempts to invoke "Self-Organized Criticality", the physical origin of this criticality, and also its universal nature, being quite insensitive to the nature of microscopic interactions, remained elusive. Recently we proposed that the precise nature of this critical behavior is manifested by a spinodal point of a thermodynamic phase transition. Moreover, at the spinodal point there exists a divergent correlation length which is associated with the system-spanning instabilities (known also as shear bands) which are typical to the mechanical yield. Demonstrating this requires the introduction of an "order parameter" that is suitable for distinguishing between disordered amorphous systems, and an associated correlation function, suitable for picking up the growing correlation length. The theory, the order parameter, and the correlation functions used are universal in nature and can be applied to any amorphous solid that undergoes mechanical yield. Critical exponents for the correlation length divergence and the system size dependence are estimated. The phenomenon is seen at its sharpest in athermal systems, as is explained below; in this paper we extend the discussion also to thermal systems, showing that at sufficiently high temperatures the spinodal phenomenon is destroyed by thermal fluctuations.

Equations

$$\Delta U\sim N^0,$$

$$\Delta U\sim N^{\beta },\beta =2/3.$$

$$Q_{12}\equiv \frac{1}{N}\sum _i^N\theta (a-|{\mathbf{r}}_i^{(1)}-{\mathbf{r}}_i^{(2)}|),$$

$$S\equiv \underset{\gamma }{\max }\left(-\frac{d\overline{⟨Q_{12}⟩}}{d\gamma }\right).$$

$$S\approx C{N}^{\theta };\theta =0.41\pm 0.09.$$

$$P({\gamma }_c,N)=\sqrt{N}\tilde{P}\left(({\gamma }_c-{\gamma }^{\ast })\sqrt{N}\right)$$

$${\Delta }_{{\gamma }_c}\sim N^{-1/2},$$

$$Q_{12}(\mathbf{r})\equiv \sum _{i=1}^N\theta (\ell -|{\mathbf{r}}_i^{(1)}-{\mathbf{r}}_i^{(2)}|)\delta (\mathbf{r}-{\mathbf{r}}_i^{(1)}),$$

$$G_L(\mathbf{r})=2G_R(\mathbf{r})-{\Gamma }_2(\mathbf{r}),$$

$$G_R(\mathbf{r})$$

$${\Gamma }_2(\mathbf{r})$$

$${\tilde{\Gamma }}_2(\mathbf{r})=\frac{{\sum }_{i\ne j}(u_i^{(12)}-Q_{12})(u_j^{(12)}-Q_{12})\delta (\mathbf{r}-({\mathbf{r}}_i^{(1)}-{\mathbf{r}}_j^{(1)}))}{{\sum }_{i\ne j}\delta (\mathbf{r}-({\mathbf{r}}_i^{(1)}-{\mathbf{r}}_j^{(1)}))}$$

$${\tilde{G}}_R(\mathbf{r})=\frac{{\sum }_{i\ne j}[u_i^{(12)}{u}_j^{(12)}-2u_i^{(12)}{u}_j^{(13)}+Q_{12}{Q}_{34}]\delta (\mathbf{r}-({\mathbf{r}}_i^{(1)}-{\mathbf{r}}_j^{(1)}))}{{\sum }_{i\ne j}\delta (\mathbf{r}-({\mathbf{r}}_i^{(1)}-{\mathbf{r}}_j^{(1)}))}.$$

$$u_i^{(12)}\equiv \theta (\ell -|{\mathbf{r}}_i^{(1)}-{\mathbf{r}}_i^{(2)}|).$$

$${\chi }_{{}_{G_L}}(\gamma )\equiv \int d^{2}x{G}_L(x,y;\gamma ).$$

$$G_R(x=0,y;\gamma )\approx C+A\exp (-y/\xi ),$$

$$\xi \approx ({\gamma }_c-\gamma {)}^{-\nu },\nu \approx 2.4\pm 0.35.$$

$$X^a\equiv \{{\mathbf{r}}_i^a{\}}_{i=1}^N,X^b\equiv \{{\mathbf{r}}_i^b{\}}_{i=1}^N,$$

$$f[T,T_g]$$

$$Q_{ab}=\frac{1}{N}\sum _{i=1}^N\theta (\ell -|{\mathbf{r}}_i^a-{\mathbf{r}}_i^b|).$$

$$f[T,T_g]=\underset{s\to 0}{\lim }{\partial }_s\Phi [T,T_g;s],$$

$$\Phi [T,T_g;s]=-\frac{1}{\beta N}\log \int d{X}_{0}d{X}_1\cdots d{X}_s{e}^{-{\beta }_{g}V(X_0)}{e}^{-\beta V(X_1)}\delta (q-Q_{01})\cdots e^{-\beta V(X^s)}\delta (q-Q_{0s}),$$

$$\Phi =-\frac{1}{\beta N}\int \mathcal{D}{Q}_{ab}{e}^{-dS(Q_{ab})}.$$

$$\mathcal{D}{Q}_{ab}\equiv \prod _{a<b}^{0,s}d{Q}_{ab},$$

$$F(h=0,T)=\underset{m}{\min }G(m,T)$$

$$\frac{\partial G}{\partial m}=0⟹m=\tanh (\beta m)$$

$$M_{ab;cd}\equiv \frac{{\partial }^{2}S}{\partial Q_{a<b}\partial Q_{c<d}},a,b,c,d\in [1,s]$$

$$\sum _{e\ne f}{M}_{ab;ef}{G}_{ef;cd}=\frac{{\delta }_{13}{\delta }_{bd}+{\delta }_{ad}{\delta }_{bc}}{2}$$

$$G_{ab;cd}=\overline{⟨(Q_{ab}-⟨Q_{ab}⟩)(Q_{cd}-⟨Q_{cd}⟩)⟩},$$

$$M_{ab;cd}$$

$${\lambda }_R$$

Full text

I Introduction

Mechanical failure of amorphous solids is an unwanted and often catastrophic event, occurring when enough strain and stress accumulate due to external loading. The phenomenon is ubiquitous in nature in the form of earthquakes due to tectonic activity and in material engineering due to shear or tensile strains. The phenomenon is known to be “scale-free” in the sense that the statistics of energy release upon failure appears to have no typical scale, a characteristic that is exemplified by the Gutenberg-Richter law Gutenberg (2013) in the geophysical context. Many authors commented that this scale-free nature indicates that material failure should be a critical phenomenon with power-law scaling, but until recently the precise origin and the actual character of this criticality remained unknown. Precisely three decades ago P. Bak and coworkers Bak et al. (1987) offered the idea of “Self Organized Criticality” to explain the ubiquity of such scale-free statistics, but the correspondence to the microscopic structure of amorphous solids and the particle-scale mechanisms that are responsible for the phenomenon remained mysterious. Recently Jaiswal et al. (2016); Parisi et al. (2017) the source of the criticality was revealed in the form of a spinodal criticality which appears to be quite universal in athermal conditions independently from the detailed microscopic interactions between the particles forming the amorphous solid. This criticality is not at all ‘self-organized’, rather it is forced on the system by the external loading. The aim of this paper is to review the pertinent features of this phenomenon and extend its exploration from athermal systems to amorphous solids at finite temperatures. Among other issues discussed below it will be shown that when the temperature becomes high enough the spinodal characteristics are destroyed by thermal fluctuations.

Solids are states of matter capable to respond elastically to a small externally applied shear deformation Landau and Lifshitz (1959). However when the external strain grows the response of all solids becomes mixed with plastic deformations, and eventually they suffer a mechanical yield. In crystalline solids plasticity and yield involve defects and dislocations. In amorphous materials such as molecular and colloidal glasses, foams, and granular matter there is no long range order with respect to which defects can be defined. Thus the mechanisms of plasticity and yield in amorphous materials need be understood along different lines from those of crystalline matter. The physics near the yielding point of this vast class of materials, as reported in a host of strain-controlled simulations Varnik et al. (2004); Maloney and Lemaître (2004); Demkowicz and Argon (2005); Tanguy et al. (2006); Lemaître and Maloney (2006); Lerner and Procaccia (2009); Rodney et al. (2011) and experiments Subhash et al. (2006); Kara et al. (2013); Noradila et al. (2013) shows a high degree of universality despite the different nature of the systems involved. Importantly, one finds that at the onset of flow at yielding, there appear typical system spanning excitations referred to as shear-bands Berthier and Biroli (2011); Dasgupta et al. (2012). We refer to a plastic event as a shear band when previously homogeneous shear strongly localizes, leaving the rest of the material less perturbed. This phenomenon is of capital importance for engineering applications as it is responsible for the brittleness typical of glassy materials, in particular metallic glasses Ashby and Greer (2006), whose potential for practical use is stymied by their tendency to shear-band and fracture Dasgupta et al. (2012, 2013a, 2013b). Measurements of plastic events occurring after yield reveal scale-free energy or stress drops, typically characterized by power-law statistics Karmakar et al. (2010a); Franz and Spigler (2017). The aim of this paper is to present the current understanding of this scale-free behavior which, as said above, is suspected to be related to some criticality.

This paper is organized as follows: in Sect. II we discuss the universal features of mechanical yield, explaining that any appropriate theory must use generic order parameters which are equally applicable to a large variety of amorphous solids. This is crucial. After introducing the “overlap” order parameter, we turn to using it to investigate the physics of yield in athermal conditions. The key result will be that yield is tantamount to a spinodal point in the emerging phase transition that is associated with the phenomenon. In Sect. III we follow up on the identification of the precise criticality that is implied by the spinodal point, and we study the correlation functions that are expected to exhibit a divergent correlation length. We then show in Sect IV that the correlation length associated to these correlators diverges as a power law in the distance from the spinodal point, cf. Eq. (17) below. Section V explores the modifications caused by having a finite temperature. Not too surprisingly, we will discover that at sufficiently higher temperatures, fluctuations destroy the spinodal characteristics, forcing a cross over to different statistics of the energy drops. In Sect. VI we offer a summary of the paper and thoughts about the road ahead. We also comment on the notion of ‘self organized criticality’ which is a very vague notion and explain how it is related to the results of this paper.

II Mechanical Yield, Universality and Order Parameter

II.1 Universality of Mechanical Yield

To introduce the main issue consider Fig. 1 showing a typical stress vs. strain curve obtained using standard numerical simulations in a strain-controlled athermal quasistatic (AQS) shearing protocol. This particular figure pertains to a Kob-Andersen 65-35% Lennard Jones Binary Mixture Kob and Andersen (1994) of 4000 particles in $2d$. Similar curves were computed and measured in a large variety of simulations and experiments. The universal features that need to be observed are the following: (i) For very small strain values the stress increases linearly according to the laws of linear elasticity. One should note that the region of purely elastic behavior is expected to reduce with the system size, shrinking to nonexistence in the thermodynamic limit Hentschel et al. (2011). Nevertheless, before a value of the strain known as the “yield strain” $\gamma_{{}_{\rm Y}}$, the plastic events are “small”, in the precise sense that the energy drop $\Delta U$ associated with them is system size independent

[TABLE]

where $N$ is the number of particles in the systems. The nature of these plastic events is identified as quadrupolar displacements, known also as Eshelby Eshelby (1957) events, which can release stress locally in regions that are particularly susceptible to the type of loading employed. The important point is that, whether elastically or punctuated by plastic events, the stress $\sigma$ continues to increase with the strain $\gamma$ until the latter exceeds $\gamma_{{}_{\rm Y}}$ which in Fig. 1 is about $\gamma_{{}_{\rm Y}}\approx 0.07$. After that point, in strain controlled protocols the strain increases without increasing the average stress - the material “flows” keeping an average “flow stress”. In stress controlled experiments, exceeding the average flow stress results in a mechanical collapse of the material. In athermal conditions it was found that the transition around $\gamma=\gamma_{{}_{\rm Y}}$ is associated with a change in the plastic response which is no longer localized, but rather exhibits system spanning events, known also as micro shear bands, in which the energy release becomes sub-extensive Karmakar et al. (2010a),

[TABLE]

The mechanism for the creation of these micro shear bands was elucidated in Dasgupta et al. (2012), and it has to do with the preferred appearance of concatenated series of Eshelby quadrupoles (lines in $2d$ or embedded in a plane in $3d$) that organize the displacement field to localize the shear on narrow lines or planes respectively. The interested reader is referred to Ref. Dasgupta et al. (2013a) where detailed energy estimates were offered to explain the energetic preference of single Eshelby quadrupoles at low strains vs. the appearance of a density of such objects at higher values of the strain.

The key observation is that after yield strain $\gamma_{{}_{\rm Y}}$ the stress cannot grow on the average, no matter how much the strain is increased. What remained obscure for a long time is what is the difference in the material before and after the yield point; why the stress could continue growing with the strain before yield, but it cannot do that after yield. Since the phenomenon is ubiquitous, the universality of this basic phenomenology of yielding begs an explanation in terms of a universal theory, in the sense that such a theory should rely on a statistical-mechanical framework and be independent of details such as chemical composition and the production process of the material.

II.2 Order Parameter and Transition

In Ref. Jaiswal et al. (2016) it was made clear that the difficulty in making a distinction between the pre- and post-yield configurations lies in the fact that there is really no distinction. The crux of the matter is not in the nature of configurations but in their number. The yield takes place because of a sudden opening up of a vast number of marginally stable configurations that are not at the system’s disposal before yield. To demonstrate this one needs to employ an order parameter that is designed Franz and Parisi (1995); Cammarota et al. (2010); Yeo and Moore (2012); Berthier (2013); Berthier and Coslovich (2014); Biroli et al. (2014); Berthier and Jack (2015); Ninarello et al. (2015) to compare two different glassy configurations $\{{\bm{r}}^{(1)}_{i}\}_{i=1}^{N}$ and $\{{\bm{r}}^{(2)}_{i}\}_{i=1}^{N}$

[TABLE]

wherein $\theta(x)$ is the Heaviside step function. The parameter $a$ is of the order of the microscopic interaction length, and is determined by trial and error. The quantity $Q_{12}$ is called an “overlap” since it has a value that goes from [math] (completely decorrelated configurations) to $1$ (identical particle coordinates within the tolerance of $a$). Its purpose is to measure the degree of similarity between configurations.

Let us now consider a glass, made by quenching a super-cooled liquid with $N$ particles down to a certain temperature $T\geq 0$ at a suitable rate. A glass is an amorphous solid wherein particles vibrate around an amorphous structure. So, if we take two configurations $\{{\bm{r}}_{i}^{(1)}\}_{i=1}^{N}$ and $\{{\bm{r}}_{i}^{(2)}\}_{i=1}^{N}$ from this glass at two different times, they will be most likely close to each other with $Q_{12}$ of the order of unity. If one is able to obtain a good sampling of the typical configurations visited by the particles in the glass, one can measure the probability distribution of the overlap $P(Q_{12})$, which will be strongly peaked around an average value $\langle Q_{12}\rangle$ close to unity. The configurations visited by the particles will then form a small connected “patch” in the configuration space of the system, selected by the amorphous structure provided by the last configuration that was realized by the liquid glass former before it fell out of equilibrium while forming a glass.

Things will change once we begin straining this glass. While the stress increases, there appear plastic events that are associated with irreversible displacements in the particle positions. The average order parameter $\langle Q_{12}\rangle(\gamma)$ responds to these displacements, reducing from $O(1)$ to lower values. An important point to understand is that before reaching the yield strain $\langle Q_{12}\rangle(\gamma)$ tends to remain around unity, but as the mechanical yield takes place a sharp phase transition occurs, whereupon sub-extensive plastic events Lerner and Procaccia (2009); Dasgupta et al. (2012, 2013b) begin to take place. These are sufficiently large, cf. Eq. (2), to cause substantial displacements, allowing different regions of the configuration space to affect the order parameter. In such a situation, the distribution $P_{\gamma}(Q_{12})$ may develop two peaks: one at high $Q_{12}$ corresponding to configurations in the same patch and one for a smaller value of $Q_{12}$ corresponding to configurations that were “ergodized” by the mixing of the sub-extensive plastic events.

To demonstrate this fundamental idea we can use any model glass, since this order parameter description is expected to be universal. Here we review molecular dynamics simulations of a Kob-Andersen 65-35% Lennard Jones (LJ) Binary Mixture in $2d$, using five system sizes, $N=500,1000,2000,4000$ and $N=10000$. We chose $Q_{12}$ with $a=0.3$ in LJ units, but verified that changes in $a$ leave the emerging picture invariant. As a first step, we prepared a glass by equilibrating the system at $T=0.4$, and then quenching it (the rate is $10^{-6}$ in LJ units) down to $T=1\cdot 10^{-6}$ into a glassy configuration. The sample is then heated up again to $T=0.2$, and a starting configuration of particle positions is chosen at this temperature. Note that while at $T=0.4$ equilibration is sufficiently fast, at $T=0.2$ the computation time is much shorter than the relaxation time. The configuration is then assigned a set of velocities randomly drawn from the Maxwell distribution at $T=0.2$, and these different samples are then quenched down to $T=0$ at a rate of $0.1$. This procedure can be repeated any number of times (say $n$ times), and it allows us to get a sampling of the configurations inside one single “patch”. We verify that the typical overlap of the ensemble of configurations so obtained in one patch is close to $\langle Q_{12}\rangle=1$, signaling that indeed the ensemble is completely located in a single patch.

Having generated one such patch, we repeat the procedure starting from another equilibrated configuration of the liquid to create another patch. The process is then repeated to generate as many patches (say $m$ patches) as needed to obtain good statistics, depending on the system size.

We then apply to each configuration in a given patch an AQS protocol as described above. This will create for each value of $\gamma$ a strained ensemble of configurations in the patch. The order parameter Eq. (3) is computed by using all the $n(n-1)/2$ unique pairs of configurations generated in the strained ensemble at a given $\gamma$. We stress that we do not compare configurations at a given value of $\gamma$ to the reference configuration at $\gamma=0$, but rather the overlap between pairs of configurations at the same value of $\gamma$. Having computed the $\gamma$ dependence of an average $\langle Q_{12}\rangle$ from the $n(n-1)/2$ of configurations in one patch, we average the results over $m$ patches to obtain the average order parameter denoted as $\overline{\langle Q_{12}\rangle}$, wherein the acute brackets denote the average over a single patch and the overline denotes the average over all patches. We present the results for $N=10000$ in Fig. 2. Note that the initial ensemble for $\gamma=0$ shows a value of the averaged order parameter $\overline{\langle Q_{12}\rangle}=1$, signifying that our initial ensemble is indeed composed of close-by configurations. As the ensemble is strained, the value of the order parameter gets lower, dropping towards zero when the strain is increased beyond the yield strain. Below we will show that the sharpness of the transition depends on the system size $N$, getting sharper and sharper when $N$ increases, as expected.

To determine the yield strain $\gamma_{{}_{Y}}$ accurately, one should construct the probability distribution function (pdf) $P_{\gamma}(Q_{12})$ by hystogramming the values of $Q_{12}$ within a patch of $n$ configurations obtained as explained above, and then average the result over the $m$ available patches. The result is denoted $\overline{P_{\gamma}(Q_{12})}$. We ask at which value of $\gamma$ this averaged pdf has two equally high peaks, see Fig. 3. The resulting $\overline{P_{\gamma}(Q_{12})}$ determines a value of $\gamma_{{}_{Y}}\simeq 0.088$. Note that this criterion implies a sharp definition of “yield” which seems absent in the current literature. If accepted, it indicates that the mechanical yield occurs beyond the stress overshoot in correspondence with the mean-field results of Ref. Rainone et al. (2015). We should also state here the yield point and the spinodal point (denoted below as $\gamma_{S}$) are not identical for finite $N$, although they become closer when $N$ increases, and see below for details.

Once we identify the phase transition point, we can demonstrate the transition itself. In Fig. 4 we display the change in $\overline{P_{\gamma}(Q_{12})}$ in the vicinity of the critical point $\gamma_{{}_{Y}}$ as a function of $\gamma$.

Within a very narrow range of $\gamma$, of the order of $\Delta\gamma\approx 0.017$, we observe a first-order like transition from a pdf with dominant peak at high values of $Q_{12}$ to a dominant peak at low values of $Q_{12}$. We capture a very unambiguous and qualitative change in behavior as the yielding point is reached.

To sharpen the understanding of what is happening in the vicinity of the yield point we examine next how many of our realizations loose the tight overlap and where the loss of overlap is taking place. To this aim we consider, as an example for the system of 4000 particles, all the 50,000 realizations that we have from 100 patches each containing 500 configurations. These are obtained by 100 choices of liquid realizations, each of which is velocity randomized 500 times (chosen with Boltzmann probabilities). When the strain $\gamma$ is increased in our AQS algorithm, we keep computing the order parameter $Q_{12}$ where the first configuration $\{{\bm{r}}^{(1)}_{i}\}_{i=1}^{N}$ in Eq. (3) is chosen randomly from all the available configurations at that value of $\gamma$, and the second is any one of the other available configurations at the same value of $\gamma$. We confirmed that changing the randomly chosen $\{{\bm{r}}^{(1)}_{i}\}_{i=1}^{N}$ does not affect the results. Next, choosing $Q_{12}=0.8$ as a threshold value, we now count how many of our observed configurations cross this threshold and exhibit $Q_{12}\leq 0.8$. The number of configurations that do so as a function of the strain (superimposed on the stress vs. strain curve) is shown in Fig 5.

The conclusion of this test is that in the vicinity of the yield point $\gamma_{{}_{Y}}$ all the configurations lose their overlap with the initial configuration, but not before. The mechanical yield is tantamount to the opening up of a vast number of possible configurations, whereas before yield the system is still constrained to reside in the initial meta-basin of the free energy landscape.

The upshot of these results is that we are able to focus on the essential feature that is responsible for the mechanical yield: a very constrained set of configurations available to the system before yield is replaced upon yield with a vastly larger set of available configurations. This much larger set is generic; we would like to refer to the phenomenon as “stressed ergodization”. The initially prepared close-by configurations are now scattered, but all of them are stressed with stress value close to the yield stress. They are all marginally stable in the sense that they would yield plastically with any increase of strain Karmakar et al. (2010b); Gendelman et al. (2015). We propose this as a universal mechanism for the ubiquitous prevalence of stress vs. strain curves that look so similar in a huge variety of glassy systems.

II.3 System size dependence

In the context of first-order phase transitions one expects that the transition should become sharper as a function of system size. To this aim we consider the dependence of $\overline{\langle Q_{12}\rangle}$ on $\gamma$ for a series of system sizes, see Fig. 6.

Indeed, the sharpening of the transition is obvious to the bare eye. To quantify it we evaluated the derivative of this function, see the upper panel in Fig. 7 for $N=4000$, and computed the maximum of this derivative function, denoted as

[TABLE]

Finally, the value of $S(N)$ is plotted in a log-log plot vs. the system size $N$ as shown in the lower panel of Fig. 7.

This log-log plot indicates the existence of a power law of the form

[TABLE]

The error bars measured here suggest that the exact value of the exponent $\theta$ is $\theta=1/2$. Such an exponent indicates that the width of the transition is not determined by the thermal fluctuations in the parent fluid from which our glassy patches were quenched, as it would be in the case of an ordinary first-order transition. Rather, it is dominated by the disorder fluctuations (i.e. sample to sample fluctuations, due to the fact that each glass is randomly “selected“ at quenching time by a parent configuration in the high temperature liquid. This causes $\gamma_{{}_{\rm Y}}$ to vary from sample to sample. To test this hypothesis we return to our numerical data and compute, for each patch, a yield point $\gamma_{c}$ which we identify as the first value of the strain for which $\langle Q_{12}\rangle\leq 0.5$. Having done so we can evaluate the probability distribution function $P(\gamma_{c},N)$. These functions obviously depend on the systems size as shown in the upper panel of Fig. 8.

To examine the scaling of the width of these distributions we rescale the data according to the ansatz

[TABLE]

where $\gamma^{*}$ is the peak value of each pdf. The data collapse means that indeed the disorder leads to a spread $\Delta_{\gamma_{c}}$ in the values of $\gamma_{{}_{\rm Y}}$ that scales like

[TABLE]

which will end up as the scaling law Eq. (5) with $\theta=1/2$. If we just had a thermal origin to the measured width we could expect rather a scaling law with $\theta=1$, as typical of first-order transitions Binder and Landau (1984). This finding highlights the pivotal role played by the fluctuations over the disorder in the finite-size scaling of the yielding transition.

II.4 Concluding this section

The upshot of this section is that the yield is associated with a first order phase transition such that before yielding the amorphous system is limited to a small patch in the configuration space, very far from any kind of ergodicity. The yielding transition is an opening of a much larger available configuration space, whereupon the system is ergodized subject to the constraint of constant mean stress. The generic configurations that are created by the mixing caused by micro shear-bands include many marginally stable states which yield easily upon the increase of strain. This is why the stress cannot increase further on the average.

This realization does not explain yet where is the criticality. In general first order phase transitions are not characterized by diverging correlation lengths, while critical points associated with second order phase transitions do. The point to understand, as sharpened in the next section, is that first order phase transitions are bordered by spinodal points which do exhibit criticality. To see this pictorially examine again Fig. 4 and focus on the pdf associated with $\gamma=0.097$. At that point the maximum of high values of $\overline{\langle Q_{12}\rangle}$ has been reduced to a saddle. This is a spinodal point that we denote as $\gamma_{S}$ where the slope of the curve vanishes as well as the second derivative. This is where a correlation length is expected to diverge as we are going to explain in the next section. The reader should also take into account that when $N\to\infty$ also $\gamma_{{}_{Y}}\to\gamma_{c}\to\gamma_{S}$.

III Theory of spinodal criticality

The aim of this section is to clarify the identification of the yielding transition as a spinodal point Urbani and Zamponi (2016). This is the point where the metastable, high overlapped glassy patch of configurations, becomes unstable with respect to a new phase with low $Q_{12}$, associated with a stressed ergodized system in the presence of disorder Nandi et al. (2016). A previously known example of such a spinodal is the Mode Coupling crossover Berthier and Biroli (2011), characterized by dynamical slowing down and heterogeneities, whose behavior is characterized by a dynamical lengthscale which can be extracted from suitable multi-point correlators Berthier and Biroli (2011). This kind of critical behavior should also be found at the yielding transition, conditional that one is able to derive the expression of the right correlator to measure. It is important to stress here that the reason that a spinodal point can be exposed and measured is that the glassy time scales and the athermal conditions stabilize the metastable system until the spinodal point is crossed and the system becomes unstable against constrained ergodization. We will see below how thermal fluctuations may destroy the spinodal characteristics.

In statistical mechanics with a suitable Gibbs free energy $G[\phi]$, $\phi$ being the order parameter of choice, stable phases are identified with its points of minimum in $\phi$. Of particular interest are instances for which the curvature of these minima goes to zero, inducing a critical behavior which manifests diverging susceptibilities-fluctuations, critical slowing down of the dynamics, and growing correlation lengths Zinn-Justin (2002). At a spinodal point, for example, one such minimum becomes unstable and transforms into a saddle. In the case of the order parameter $Q_{12}$ the general form of the free energy $s[Q_{12}]$ had been already derived and studied (see De Dominicis et al. (1998) for a review) in the context of the theory of replicas originally developed for the study of spin-glasses, and its properties, at least at mean-field level, are well known (we refer to Rainone et al. (2015); Rainone and Urbani (2016) for the derivation of $s[Q_{12}]$ in the specific case of mean-field hard spheres); the matrix of second derivatives (or, using a more field-theoretic terminology, the mass matrix) is not diagonal in the basis of $Q_{12}$, and after diagonalization is found to have only three distinct modes, or masses De Dominicis et al. (1998). Of these, the most relevant ones are the so called replicon mode $\lambda_{R}$, which for example goes to zero at the newly proposed Gardner transition Charbonneau et al. (2014), and the longitudinal mode $\lambda_{L}$ which is instead related to spinodal points Rainone and Urbani (2016); Urbani and Zamponi (2016) such as our yielding transition. In Appendix A of this paper we review briefly the background theory that is at the basis of the present approach.

III.1 Correlation functions

Based on the introductory discussion, we now derive an expression for the correlator associated with the longitudinal mode, from whence one can extract the diverging correlation length associated with the onset of criticality at the spinodal point, and define an associated susceptibility which will shoot up as the spinodal point is approached. The first step is to “localize” the overlap function and define the ${\bm{r}}$-dependent quantity

[TABLE]

Next, as mentioned above, the expression for the longitudinal correlator in terms of four-replica correlation functions can be found by diagonalization of the correlation matrix $G_{ab;cd}$, (see Appendix A) which is defined as the inverse of the mass matrix $M_{ab;cd}$ of the replicated field theory of the overlap order parameter $Q_{12}$ De Dominicis et al. (1998). The derivation is a matter of standard diagonalization algebra, so we shall not report it here and refer to Appendix A for the details. The expression, employed for example in Baity-Jesi et al. (2014a, b) in the case of a model with spins on a lattice, reads for athermal systems

[TABLE]

with the definitions

[TABLE]

We reiterate that angular brackets denote a patch average and $\overline{(\bullet)}$ indicates an average over different patches. The quantity $G_{R}(\bm{r})$ is the correlation function of the replicon mode De Dominicis et al. (1998) and $\Gamma_{2}(\bm{r})$ is just the garden-variety four-point correlator.

Using these definitions and taking Eq. (8) into account, the quantities we compute in numerical simulation, before taking the ensemble average, are (see Appendix A and Ref. Berthier et al. (2016)):

[TABLE]

and

[TABLE]

with

[TABLE]

These four-replica objects can be computed for any quadruplet of distinct replicas. The ensemble averaged correlation functions are simply obtained as $\Gamma_{2}\equiv\overline{\tilde{\Gamma}^{(12)}_{2}}$ and $G_{R}\equiv\overline{\tilde{G}^{(12)}_{R}}$, and cf. Appendix A for a proof. We stress that one must keep the full space dependence of the correlators in the definitions above, as the shear strain breaks the rotational symmetry of the glass samples and so the correlators are not just functions of a distance $r$.

IV Numerical Results

The three correlation functions discussed in the previous section were computed within the same numerical framework as discussed in Subsect. II.2. Typical results are shown in Fig. 9, here at $\gamma_{S}=0.09405$. One notices the obvious fact that the correlation functions reflect the breaking of isotropy that is caused by the strain. In fact the spatial structure of the correlation function is quadrupolar, precisely indicating where shear bands are bound to appear in simple strain: the $x$ and the $y$ axes are in $45^{\circ}$ to the principal axis of the stress Ashwin et al. (2013).

To demonstrate the strain dependence of the correlators we consider first the susceptibilities $\chi_{{}_{G_{L}}},\chi_{{}_{G_{R}}}$ and $\chi_{{}_{\Gamma_{2}}}$ that can be obtained from the correlators. For example

[TABLE]

In Fig. 10 upper panel we show the susceptibility $\chi_{{}_{\Gamma_{2}}}$ as a function of $\gamma$ for the three system sizes at our disposal. Superimposed are the stress vs. strain curves obtained by averaging the individual curves over all the available configurations and glass samples. One sees very clearly the singularity that develops near the spinodal point as a function of the system size. In the middle and lower panel of the same figure we show the other two susceptibilities $\chi_{{}_{G_{R}}}$ and $\chi_{{}_{G_{L}}}$ as a function of the strain $\gamma$, again with the stress-strain curve superimposed for comparison. As we expected, the susceptibilities show a distinct peak at the spinodal point $\gamma_{S}$ where criticality is reached.

The scaling of the peak of the susceptibility $\chi_{{}_{\Gamma_{2}}}$ with the system size is expected to mirror the scaling of the response as written in Eq. (4), at least if standard fluctuation-dissipation theorems should apply to the present problem. Indeed, plotting the maximal values of $-\chi_{{}_{\Gamma_{2}}}$ as a function of $N$ in a log-log plot, cf. Fig. 11 we find that the maxima $\chi^{\rm max}_{{}_{\Gamma_{2}}}$ scale like $\sqrt{N}$ as expected.

More detailed information is provided by the full dependence of the correlators on their arguments. To see most clearly the change in the correlators as the spinodal point is approached, we consider for example the one-dimensional function $G_{R}(x=0,y;\gamma)$, shown for $N=4000$ in Fig. 12.

We note that the correlator changes both in amplitude and in extent when we approach the critical point. To quantify these changes we fit a 3-parameter function to $G_{R}(x=0,y)$ in the form

[TABLE]

where all the fitting coefficients are functions of $\gamma$. In Fig. 13 we present the $\gamma$ dependence of the amplitude $A(\gamma)$, the constant $C$ and the correlation length $\xi(\gamma)$.

An interesting observation concerns the constant $C$ used in the fit Eq. (16). This constant is also sensitive to the approach of the criticality, cf. the lower panel in Fig. 13. One could worry that integrating this constant over $y$ could contribute to the divergence of the susceptibilities. In fact the rise in $C$ near the spinodal point goes down with the system size and its contribution to the integral is reduced as well, as can be seen in Fig. 14 which presents the integral $\int dyG_{R}(x=0,y)$ from which $C\times L$ is subtracted. The conclusion is that indeed the contribution of $C$ goes down also when integrated over the system size, showing that the main contribution to the divergence of the susceptibility is from the divergence of the correlation length. It is interesting to notice that the constant $C$ decreases with the system size, presumably becoming irrelevant in the thermodynamic limit. The amplitude $A$ is still increasing with the system size, and it is difficult to assert whether it converges or not. On the other hand we can safely conclude that the data present a strong evidence for the increase in the correlation length. This conclusion is substantiated below using the correlation function $\Gamma_{2}(x,y)$. Before doing so we need to discuss the fitting procedure for the correlation function $G_{R}(x=0,y)$. In Fig. 15 we show the full results for this correlation function for all the available values of $\gamma$ and for two larger systems sizes at our disposal. One sees that the exponential decay that is used for the fit is only reliable up to the minima of the functions. The reason for the upward trend is the periodic boundary condition that reflects the correlations. To eliminate this spurious effect we presented in Fig. 12 the fit up to the minimum in the function. One should note however that the distance to the minimum increases with the system size, presumably diverging in the thermodynamic limit. Thus the fit up to the minimum allows a faithful estimate of the correlation length $\xi$.

The dependence of $\xi$ on the distance from criticality and on the system size is not easy to read from Fig. 13. In fact a smoother dependence is available from the correlation function $\Gamma_{2}(x=0,y)$ and $\Gamma_{2}(y=0,x)$. An exponential fit similar to Eq. (16) was applied to these two projections of $\Gamma_{2}(x,y)$ and the correlation length $\xi$ was determined as shown in the upper panel of Fig. 16.

The scaling exhibited in the lower panel of Fig. 16 is not perfect, but a least square fit to all the three curves leads to a scaling law in the form

[TABLE]

The estimated value of $\nu$ is unusually high. The error bars are significant, and it is quite likely that this result indicates that $\nu=2$, although at the present time we cannot offer a theoretical basis for this number.

The result Eq. (17) may have important experimental consequences, predicting the length of micro-shear bands in materials as a function of the distance from criticality. We propose that such measurements should be carried out, providing a possible direct test of the present ideas.

V The effects of finite tempeartures

The mechanical yield in athermal conditions is an excellent conceptual laboratory for clarifying the essence of the yield mechanism, but in reality many yielding amorphous solids operate under thermal conditions, effected by thermal fluctuations. It is therefore interesting and important to assess the effects of temperature on the findings described above.

To assess the effects of temperature we repeat precisely the same protocol described above to create a patch of $n$ replica at $T=0$, including the creation of $m$ such patches. The difference is that presently we warm up all the replica in a given patch to a target temperature. Results will be reported for target temperatures $T=0.1,0.2$. While keeping the strain at $\gamma=0$ each configuration was thermalized by molecular dynamics. Afterwards, each configuration was strained by increasing the strain in steps of $\delta\gamma=2\times 10^{-4}$, allowing the energy to stabilize after each such step before straining again. Typical averaged strain vs stress curves (with averages computed firstly over a patch and secondly over all the patches) for a system with $N=10000$ are shown in the upper panel of Fig. 17.

We see that at the lower temperature $T=0.1$ there is still a stress peak before the yield, but at the higher temperature $T=0.2$ the stress peak no longer exists and the stress reaches the flow steady state stress quite monotonically. At both temperatures the steady state is attained at lower values of the strain than at $T=0$. Computing the average overlap order parameter $\overline{\langle Q_{12}\rangle}$ (cf. lower panel of Fig. 17) we observe a corresponding behavior. At $T=0.1$ a remnant phase transition is still observable, with the order parameter still falling somewhat sharply after $\gamma\approx 0.04$. At $T=0.2$ there is no longer a sharp decrease, but rather a smooth decline of $\overline{\langle Q_{12}\rangle}$ as a function of $\gamma$. It is obvious that temperature fluctuations at $T=0.2$ are sufficient to destroy the spinodal characteristics of our phase transition.

The same conclusion is drawn from examining the pdf of the order parameter. In Fig. 18 we show the function $\overline{P_{\gamma}(Q_{12},T)}$ for three temperatures $T=0,0.1,0.2$ for a system with $N=10000$.

While the phase transition is observed nice and clear at $T=0$, and it still remains observable at $T=0.1$, it changes to a smooth migration of the single peak of $\overline{P_{\gamma}(Q_{12},T)}$ from high to low values of $Q_{12}$ when $\gamma$ is increased. We lose completely the double hump structure which underlies the spinodal criticality.

It is important to stress that with the loss of the spinodal criticality we also lose the qualitative distinction between the pre-yield and post-yield statistics of the energy drops as shown in Eqs. (1) and (2). Observing a stress vs. strain curve for a single realization one finds the same statistics of energy and stress fluctuations before and after the yield, since it is dominated now by temperature fluctuations rather than by mechanical instabilities. The sharp appearance of system spanning events with the yield phenomenon is caused by the spinodal criticality as explained in this paper. Once this gets destabilized by temperature fluctuations there is no increase in correlation length and we remain only with standard temperature fluctuations.

VI Summary and Conclusions

In summary, we have presented evidence that the scale free yielding transition in amorphous solids is governed by a spinodal point with disorder. The associated correlation length is exhibited by suitable four-point correlators whose expression can be obtained from replica theory. The full implications of the theory pertain to an athermal setting, and the full fledged criticality is destroyed by thermal fluctuations Shrivastav et al. (2016). In athermal conditions the transition becomes ever sharper with increasing the system size. We have found that the range of strain values over which the transition takes place goes to zero like $1/\sqrt{N}$. The correlation length $\xi$ appears to diverge following a scaling law, cf. Eq. (17). We have commented above that this prediction may be tested experimentally by examining the lengths of micro shear bands as a function of the strain or the stress while approaching mechanical collapse.

For sufficiently high temperatures the system will generally be able to escape through thermal activation from the high-$Q_{12}$ minimum before this has a chance to flatten and the relative susceptibility to diverge. However, since the nucleation time is expected to be fairly long, one should anyway be able to observe transient shear-bands/heterogeneities, as long as the temperature is low enough that nucleation does not take place until the system is close to the spinodal, which, interestingly, is precisely the behavior of transient shear bands as reported in Shrivastav et al. (2016).

Finally we need to touch upon the notion of ‘self organized criticality’. Not being quite sure what it means, we propose that it refers to the fact that after yield the system remains critical in the sense that the yielded configurations are still maintained at the yield stress, and are therefore marginal; any increase in strain will cause a repeat of the phenomena discussed above. To see this we need to recreate a patch of configurations that are closed to a given yielded configuration and examine what happens upon further straining of such a patch. To create a patch of configurations having almost same value of stress $\sigma\approx\sigma_{{}_{\rm Y}}$, we take a single post-yield configuration, and apply random displacements of randomly selected particles (keeping the displacements infinitesimal, so that the overlap function $Q_{12}$ remains close to unity). We then perform conjugate-gradient minimization to return each configuration to an athermal mechanical equilibrium (T=0). The obtained configurations are close to the selected post-yield configuration and have almost the same stress. This method is repeated to generate as many configurations as we need that belong to an approximately iso-stressed patch. Finally we strain al the configurations in the patch and compute the strain dependent order parameter as explained above. The result of this analysis is displayed in Fig. 19.

We see that the phase transition is now eliminated, the order parameter decreases smoothly and without a sharp decline at any specific value of $\gamma$. The reason is of course that each patch contains many marginal configurations that yield again and again in the strain controlled experiment. We reiterate that in stress controlled experiments the system collapses anyway when the yield-stress is exceeded. The conclusion is that the criticality is not at all self-organized, it is caused due to the mechanical straining by the external agent; the system is driven to a marginal state and is maintained there by continuing to strain the system.

In the future one needs to examine further the universality of the proposed scenario and of the scaling laws found in this paper, examining different amorphous systems and different space dimensionalities. Another interesting future path is the study of the mechanical yield in frictional aggregates. It is not known at the present point in time whether these systems fall in the same universality class or whether they might exhibit totally different behavior.

VII Acknowledgements

We thank Giorgio Parisi for inspiring discussions. We acknowledge interesting exchanges with Giulio Biroli regarding the effect of disorder on the width of the transition. This work has been supported in part by the Minerva foundation with funding from the Federal German Ministry for Education and Research, and by the Israel Science Foundation (Israel Singapore Program).

Appendix A The longitudinal correlation function

Let us start from the expression of the free energy of a glass state, prepared by equilibrating a generic glass former down to a glass transition temperature $T_{g}$ where it can still be equilibrated, and then quenching it out of equilibrium to a given temperature $T<T_{g}$. Such a free energy was first defined in Franz and Parisi (1995) in the context of spin-glass physics. Its definition in the case of structural glasses, and its computation in the particular case of hard spheres were first discussed in Rainone et al. (2015). The definition, in the case of a generic glass former made of $N$ particles is based on comparing two configurations $X^{a}$ and $X^{b}$ of the same glass. Here

[TABLE]

where the labeling ${\bm{r}}_{i}$ refers to the position of the same particle $i$ in the two different configurations. For a generic interaction potential $V(X)$ the definition of the free energy is

[TABLE]

where $\beta_{g}=1/(k_{B}T_{g})$, $\beta=1/(k_{B}T)$ and $q^{*}_{r}$ is the value of $q_{r}\neq 0$ whereupon the free energy attains a local minimum Rainone et al. (2015). The overlap function $Q_{01}$ for any two configuration, say $a$ and $b$ is Jaiswal et al. (2016)

[TABLE]

Here $\ell$ is a coarse graining parameter (in Jaiswal et al. (2016), $\ell\simeq 0.3$ in Lennard-Jones units). The idea is to consider the free energy at temperature $T$ of the glass former, which is constrained to stay close to an amorphous configuration $X_{0}$ which is selected from the equilibrium ensemble, using the canonical distribution when the glass is still at equilibrium at $T_{g}$.

The properties and computation of the free energy (19) are discussed extensively in Rainone et al. (2015); Rainone and Urbani (2016), so we refer the interested reader to those works. The explicit analytic computation is accomplished in the mean field approximation. In our paper we use the results far from the mean field limit, but we ascertain that the relevant correlation functions that are fleshed out in the mean field calculation are the relevant ones also in the general case. Of course, critical exponents can differ. In the sequel we sketch how from this mean-field theory in terms of an overlap order parameter $Q_{ab}$ one can extract the definitions of the correlation functions that are expected to show critical behavior.

The outermost integral in the Eq. (19) can be computed with the replica trick,

[TABLE]

where $\Phi$ is defined as

[TABLE]

so we are considering $s$ replicas of the $X$ configuration. In infinite dimensions for the case of hard spheres it was shown Kurchan et al. (2012) that the functional defined above can be written as

[TABLE]

Here $\mathcal{D}Q_{ab}$ denotes an integration measure over all the distinct $Q_{ab}$s,

[TABLE]

and $d$ is the number of spatial dimensions. The functional $S(Q_{ab})$ is referred to as the “replica action”. In the mean-field limit $d\to\infty$, the integral above can be computed via the saddle point method Bender and Orszag (1999), which means that one must consider the optimum points in $Q_{ab}$ of the replica action $S(Q_{ab})$. This means that $S(Q_{ab})$ plays the role of a Gibbs free energy, i.e. the free energy for fixed order parameter. An illustrative example is the case of a Curie-Weiss model (mean-field ferromagnet) wherein, for the Helmholtz free energy $F$ in zero magnetic field, one has Rainone (2014)

[TABLE]

where $G(m,T)$ is indeed the Gibbs free energy for fixed magnetization $m$. The minimization equation for $G$ is then the celebrated equation for the spontaneous magnetization

[TABLE]

and the ferromagnetic phase transition takes place when the paramagnetic, $m=0$ minimum of $G$ flattens and splits in two degenerate minima with $m\neq 0$, which implies that at the critical temperature $\frac{\partial^{2}G}{\partial m^{2}}=0$. The derivation of the $S(Q_{ab})$ action in the case of mean-field hard spheres can be found in Kurchan et al. (2012).

In the present case the $f[T,T_{g}]$ plays the role of the Helmholtz free energy $F$ and the $S(Q_{ab})$ of the Gibbs free energy $G$. With this analogy, one can understand how the critical properties of glass states are related to the matrix of second derivatives of the replica action $S(Q_{ab})$,

[TABLE]

in the limit $s\to 0$ (we stress that $X_{0}$ is not involved in this definition). The inverse $G_{ab;cd}$ of the tensor $M$, defined as

[TABLE]

is then the covariance matrix of the mean field theory

[TABLE]

wherein the angled brackets denote the thermal average restricted to a single glass sample at temperature $T$ (that is over the canonical distribution of the $X_{1}$ configuration in the (19)), and the overbar denotes the average over all possible glass samples selected at $T_{g}$ (that is over the canonical distribution of the $X_{0}$ configuration in the (19)). This covariance tensor encodes the critical fluctuations of the system near the critical points whereupon the tensor $M_{ab;cd}$ develops a zero mode.

Let us now assume that the glass state under study is a single minimum of the free-energy landscape of the system wherein all replicas from $1$ to $s$ can move ergodically, this means that the replicas are all equivalent and the matrix $Q_{ab}$ must then be invariant by any replica permutation, an hypothesis referred so as replica-symmetric (RS).

In Rainone et al. (2015) it is discussed how this is not true in all cases, i.e. there exist a regime wherein the glass basin undergoes an ergodicity breaking and fractures into sub-basins. Nevertheless, here we stick to the simple RS ansatz. In this case, since the action $S(Q_{ab})$ must in turn be invariant for any replica permutations, the most general form that the Hessian $M$ can take is

[TABLE]

and the same goes for the covariance matrix $G_{ab;cd}$. This form is completely general as it only pertains to the RS symmetry; then the only model-dependence is in the parameters $M_{1}$, $M_{2}$ and $M_{3}$, which must be computed case by case and are generally dependent on the external parameters like temperature or magnetic field.

The diagonalization of the tensor $M_{ab;cd}$ is an exercise of standard linear algebra and has been already carried out many times, see for example Crisanti and Sommers (1992); Bray and Moore (1979); De Dominicis et al. (1998) and Zamponi (2010) where it is proposed as an exercise. It is found that the tensor $M$ has only three distinct eigenvalues

[TABLE]

and the same goes for the tensor $G$. Those three eigenvalues (or modes) are called the replicon, longitudinal, and anomalous, respectively Zamponi (2010). We are interested in the longitudinal mode (which in the limit $s\to 0$ is degenerate with the anomalous one), which becomes soft at the yielding transition Rainone and Urbani (2016); Urbani and Zamponi (2016). Let us consider the $G$ tensor. Because of replica symmetry, there are only three distinct correlators that one can define, namely

[TABLE]

and in the limit $s\to 0$ we know that

[TABLE]

It is then immediate to check that

[TABLE]

which then implies

[TABLE]

with the definitions

[TABLE]

as in the main text. We have used $\Gamma_{2}=G_{12;12}-G_{12;34}$ which derives from replica symmetry, as $\left<Q_{12}\right>\left<Q_{12}\right>=\left<Q_{12}Q_{34}\right>$ in the replica-symmetric phase.

Let us now detail how to transform these definitions into quantities that can be measured in simulation. We start by ”localizing” the definition of the $Q_{ab}$ overlap in the following way

[TABLE]

In a thermal simulation the $a$ and $b$ configurations would depend on the time $t$, and so would the $Q_{ab}(r)$, so one would need to perform the in-state thermal average $\left<\bullet\right>$ by considering the equilibrium value of these quantities. In the present paper we focus un athermal solids under quasi-static shear, so we do not have dynamics and the $a$ and $b$ configurations will simply be two distinct minima of the inter-particle potential obtained through the protocol described in the main text, and the thermal average will be the average over this ensemble of configurations which make up a glassy patch.

We now apply the definition (43) in the (41), (42) to construct the correlators. For illustrative purposes we use the $\Gamma_{2}(\bm{r})$. We get, omitting the overline to lighten the notation,

[TABLE]

with

[TABLE]

as in the main text, and we used that $\left<Q_{ab}(x)\right>=Q_{ab}$. Because of translational invariance, the correlator is actually independent of $\bm{x}$. We can get rid of $\bm{x}$ by performing an integration over this variable, which, using the $\delta$-functions, gives as a result

[TABLE]

then, following Berthier et al. (2016), we omit the terms with $i=j$ (which are anyway relevant only for $\bm{r}=0$) and we normalize the correlator with the pair distribution function of the glass; we finally obtain

[TABLE]

as in the main text. The derivation for the $\tilde{G}_{R}(\bm{x})$ is then an obvious generalization.