Equilibrium avalanches in spin glasses

# Equilibrium avalanches in spin glasses

Pierre Le Doussal, Markus Müller, and Kay Jörg Wiese CNRS-Laboratoire de Physique Théorique de l’Ecole Normale Supérieure, 24 rue Lhomond,75005, Paris, France The Abdus Salam International Centre for Theoretical Physics, P. O. Box 586, 34151 Trieste, Italy
July 5, 2019
###### Abstract

We study the distribution of equilibrium avalanches (shocks) in Ising spin glasses which occur at zero temperature upon small changes in the magnetic field. For the infinite-range Sherrington-Kirkpatrick model we present a detailed derivation of the density of the magnetization jumps . It is obtained by introducing a multi-component generalization of the Parisi-Duplantier equation, which allows us to compute all cumulants of the magnetization. We find that with an avalanche exponent for the SK model, originating from the marginal stability (criticality) of the model. It holds for jumps of size being provoked by changes of the external field by where is the total number of spins. Our general formula also suggests that the density of overlap between initial and final state in an avalanche is . These results show interesting similarities with numerical simulations for the out-of-equilibrium dynamics of the SK model. For finite-range models, using droplet arguments, we obtain the prediction where and are the fractal dimension, magnetization exponent and energy exponent of a droplet, respectively. This formula is expected to apply to other glassy disordered systems, such as the random-field model and pinned interfaces. We make suggestions for further numerical investigations, as well as experimental studies of the Barkhausen noise in spin glasses.

## I Introduction

The low-temperature response of disordered systems often proceeds in jumps and avalanches. crackling (); granular (); earthquakes (); fracture (); LeDoussalWieseMoulinetRolley2009 (); vortex (); mwalls (); barkhausen (); eglass () These processes are beyond standard thermodynamic calculations and are therefore relatively difficult to access and describe analytically soc (); AvalanchesRFIM (); avalanchesFRG (); LeDoussalMiddletonWiese2008 (); RossoLeDoussalWiese2009a (); vives (). In a recent article us (), we succeeded in calculating the statistics of equilibrium avalanches (also called shocks) in a variety of disordered systems described by mean-field theory based on Parisi replica symmetry breaking. This encompasses in particular the canonical Sherrington Kirkpatrick (SK) model for the Ising spin glass SK (); MPV () and elastic manifolds in the limit of a large number of transverse dimensions. Although it has been known for a while that the equilibrium magnetization curve of the SK model undergoes a sequence of small jumps as is increased NonAvOfsuscept (), their statistics had not been obtained previously. The aim of the present article is to provide a detailed derivation of the distribution of avalanche sizes for the SK model. We introduce replica techniques that significantly extend the formalism developped in Ref. BMP, to study velocity correlations in high-dimensional Burger’s turbulence. It also generalizes previous studies of the variance of equilibrium jumps to their full distributions FRGRSB (); FRGlargeN (); RizzoYoshino (). We expect this technique to be useful in several other contexts as well. In particular, it should be helpful to describe the response of complex systems to a small change of parameters, a problem that arises in a variety of fields ranging from condensed-matter physics of complex systems, optimization problems to econophysics vertexcover (); coloring (); 1stepstabilitykSAT (); jpb ().

The main result of our calculation is that the distribution of jumps takes a scale-free form, described by a power law of the jump size. This is intimately tied to the criticality of the spin-glass phase of the models analyzed criticalityMF (), and we conjecture that such a criticality is a feature of a large variety of frustrated glassy systems. Palassini () The exact result obtained in the SK model finds a natural interpretation which allows for an extension to finite dimensions via droplet scaling arguments. Those relate the equilibrium-avalanche exponent to critical properties of droplet excitations.

Our results complement previous numerical simulations by Pazmandi et al. SKnumerics () on out-of-equilibrium hysteresis at in the SK model, which exhibit surprising similarities, as we will discuss. Understanding the relations between these results requires further numerical investigations of dynamic and static avalanches, both in mean-field and finite-dimensional spin glasses. Our results suggest to look for power-law distributed Barkhausen-type noise in spin and electron glass experiments, as will be discussed.

This paper is organized as follows: In Sec. II we revisit the Parisi saddle-point equations in the presence of a small varying external magnetic field. In Sec. III, we generalize the Paris-Duplantier equations to compute the moments of the magnetizations in different fields. From that calculation we extract the distribution of equilibrium jumps in Sec. IV. In Sec. V we consider the case of finite-dimensional spin glasses, and using droplet arguments we obtain a power-law distribution of equilibrium avalanches. In Sec. VI we discuss the connection with previous numerical studies on spin and electron glasses, and propose experimental and numerical investigations.

## Ii Model and method

### ii.1 Model and aim of the calculation

We study the SK spin-glass model of energy

 H=−N∑i,j=1Jijσiσj−HextN∑i=1σi, (1)

where the are i.i.d. centered Gaussian random variables of variance , that couple all Ising spins, and is the external field.

Our aim is to follow the equilibrium state as a function of the applied field at low temperature . We consider small variations of the applied field around a reference value , .

We are interested in the total magnetization in a given sample,

 M(Hext)=∑i⟨σi⟩Hext=−∂HextF , (2)

where is the free energy. Since upon variation of of order one we expect jumps of the total magnetization of order we define:

 mh=1√NM(H+h√N)=−∂hF(h) (3)

where from now on we denote the free energy in the external field . Note that is the sum of a constant part of order , , plus a fluctuating part of order unity.

To characterize the statistics of these order-one jumps in we need to compute the following cumulants in different physical fields , :

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpc=∂h1…∂hpS(p)(h1,h2,…,hp) . (4)

It is obtained from the cumulants of the sample-to-sample fluctuations of the free energy,

 S(p)(h1,h2,…,hp)=(−1)p¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯F(h1)…F(hp)J,c, (5)

where we denote by the average over disorder and its connected average.

These can be obtained from the generating function of replica submitted to different fields ,

 exp(W[h]):=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯exp[−βn∑a=1F(ha)]J . (6)

Note that fields with replica index are denoted with upper index to distinguish it from the physical field with lower index. Hence

 W[h] = ∞∑q=0βqq!∑a1,…aqS(q)(ha1,…,haq) (7)

We now derive a formula for from the saddle-point equations in the large- limit.

One has:

 eW[h] (8) = ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∑{σia}exp[β∑ijσiaJijσja+β∑i(H+ha√N)σia]J = ∑{σia}∫∏a≠bdQab∏iexp(nNβ2J22) ×exp[∑aβ(H+ha√N)σia] ×exp[β2J2∑a≠b(−N2Q2ab+Qabσiaσib)].

Note that on spins , we put the replica-index at the bottom, and the site index at the top. Now we define the local partition sum

 eA(Q,h) (9) :=∑{σa}exp⎡⎣β2J2∑a≠bQabσaσb+β∑a(H+ha√N)σa⎤⎦,

in terms of which we can write

 eW[h] = ∫∏a≠bdQabexp[nNβ2J22−N2β2J2∑a≠bQ2ab (10) +NA(Q,h)].

In the limit of we can perform a saddle-point evaluation. For this is the usual SK saddle-point equation in presence of a field . In the low-temperature phase considered here, it has a set of solutions, denoted . They are obtained from the standard Parisi solution by applying a permutation of the indices. Each saddle point of the path integral over satisfies the self-consistent equation for :

 ⟨σaσb⟩A(q,0)=qab , (11)

where refers to an average with action from Eq. (9). Since changes in the external fields are of size , they alter the solution of the saddle-point equation from to . Hence we can compute the contribution to of each saddle point in perturbation theory. For a given saddle point, each contribution to is of the form , with

 V[q,h]:=nNβ2J22−N2β2J2∑a≠bq2ab+NA[q,h] . (12)

 ∂qabV[qh,h]=0 . (13)

Using this equation, we obtain the following expansion in replica fields :

 V[qh,h] = V[q,0]+∑aha(∂haV)[q,0] +12∑abhahb(∂2hahbV)[q,0]+O(1√N) = V[q,0]+β∑aha¯¯¯¯¯¯¯m0 +12∑abhahbβ2[δab+(1−δab)qab]+O(1√N)

In the first line we used the condition (13), its total derivative w.r.t. , and that to eliminate the cross term ; in the second line we used Eq. (11).

The final expression for is obtained by performing the sum over all saddle points ,

 eW[h]−W[0]−β∑aha¯¯¯¯¯¯m0 =′∑πeβ22∑a(ha)2(1−qaa)+β22∑abqπabhahb . (14)

The prime on the permutation sum indicates that the sum is normalized by . For convenience, we have introduced to be defined later.

Let us define the “non-trivial” part of as

 ~W[h] := W[h]−W[0]−β∑aha¯¯¯¯¯¯¯m0 (15) −β22∑a(ha)2(1−qaa) = ln(′∑πeβ22∑abqabhπ(a)hπ(b)).

To obtain the -th cumulant, we need to consider for groups of replica with . Each group is subject to a different physical field , . This field is constant within a replica group. We remind that we use superscript indices to denote replicas, and subscript indices to label the replica groups. The resulting (and likewise for ) has the cumulant expansion

 Wp[h]=∑qβqq!nqp∑i1=1…p∑iq=1αi1…αiqS(q)(hi1,…hiq). (16)

The magnetization cumulants for can then be extracted as

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c=∂h1…∂hpS(p)(h1,…,hp) =limn→01npβp∏pi=1αi∂h1…∂hp~Wp[h] . (17)

This works, since the terms in (16) with vanish after the differentiation and the ones with vanish in the limit , leaving the desired term .

## Iii Calculation of moments

### iii.1 Generalized flow equation

To proceed, we decouple the ’s by a Hubbard-Stratonovich transformation,

 e~W[h]=⟨′∑πe∑ahπ(a)μa⟩μ, (18)

where are Gaussian random variables with variance , and denotes the average over them.

The sum over permutations in (18) is equivalent to a normalized sum (indicated by a prime) over assignments , describing the permutation :

 hπ(a)=hia. (19)

Since the permutation preserves the number of equivalent replica, we have the constraint . With this notation we obtain

 (20)

As we prove in App. C, this can be rewritten as

 e~Wp[h]=⟨∫∞−∞∏pi=1dyiδ(∑pi=1αiyi)∏na=1[∑pi=1exp(hiμa+yi)]∫∞−∞∏pi=1dyiδ(∑pi=1αiyi)[∑pi=1exp(yi)]n⟩μ, (21)

valid for any , and for any set of , with . This identity significantly generalizes the formula (D6) in Ref. BMP, .

In the case where is a hierarchical ultrametric matrix of Parisi type, parameterized by the Parisi function with , the average over of expression (21) can be performed extending the methods of Ref. Duplantier, . We recall that we use everywhere and rewrite

 e~Wp[h]=∫∞−∞∏pi=1dyiδ(∑iαiyi)g(x=n,{yi})∫∞−∞∏pi=1dyiδ(∑iαiyi)[∑iexp(yi)]n. (22)

We have defined

 g(x;{yi}) ≡ exϕ(x;{yi}) (23) ≡ ⟨x∏a=1(p∑i=1eyi+hiμ(x)a)⟩μ(x)

The auxiliary fields have Gaussian correlations . For convenience, we define .

Generalizing the method of Ref. Duplantier, to several groups, we find the flow equation for the function defined above:

 ∂ϕ∂x=−β22p∑i,j=1hihjdq(x)dx(∂2ϕ∂yi∂yj+x∂ϕ∂yi∂ϕ∂yj) . (24)

It must be solved with the boundary condition

 ϕ(x=1;{yi})=log(p∑i=1eyi)≡H({yi}) . (25)

Here and below we denote .

To simplify (22) we first evaluate the denominator. In the limit of one can show the general formula for any with constraint and :

 ∫dpyδ(∑iαiyi)enH(→y) =∫dpyδ(∑iαiyi)e−(−n)max(yi)[1+O(n)] =1∏iαi(−n)1−p[1+O(n)]. (26)

Expanding also the numerator and the exponential in (22) to lowest non-trivial order in , we find

 ∂h1…∂hp~Wp[h] (27) = n∫dpyδ(∑iαiyi)∂h1…∂hpϕ(0,→y)(−n)1−p/∏iαi[1+O(n)] = −(−n)pp∏i=1αi∫∞−∞dpyδ(∑iαiyi)∂h1…∂hpϕ(0,→y) ×[1+O(n)].

Inserting this into Eq. (17) we obtain the final formula for the -th cumulant of the reduced magnetization:

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c (28) =−(−T)p∫dpyδ(∑iαiyi)∂h1…∂hpϕ(0,→y).

This expression is independent of the choice of , as it must be. In case the order parameter function, , has a plateau for (as happens for the SK model in a magnetic field )

 ϕ(x=0,→y)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ϕ(xm,→y+zβ→h√q0)z, (29)

where denotes an average over , a unit-centered Gaussian variable,

 ¯¯¯¯¯¯¯¯¯¯f(z)z:=∫∞−∞dz√2πe−z2/2f(z) , (30)

see Eq. (96) in Ref. FRGRSB, .

### iii.2 TBL-shock expansion

We now solve the flow equation (24) perturbatively in the nonlinear term. This generates a low-temperature expansion which is well suited to study shocksFRGRSB (). We write

 ϕ(x,→y)=ϕ0(x,→y)+ϕ1(x,→y)+… . (31)

The successive terms satisfy

 ∂ϕ0∂x = −β22∑ijhihjdq(x)dx∂2ϕ0∂yi∂yj (32)

with initial condition .

 ∂ϕ1∂x=−β22∑ijhihjdq(x)dx(∂2ϕ1∂yi∂yj+x∂ϕ0∂yi∂ϕ0∂yj), (33)

with initial condition .

The leading-order equation (32) is a linear diffusion equation, and integrated as (for )

 ϕ0(x,→y) = ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯H(→y+zβ→h√q(1)−q(x))z. (34)

Taking into account (30), we find the contribution of to the magnetization cumulants

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c,(0) (35) = −(−T)p∫dpyδ(∑iαiyi) ×∂h1…∂hp¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯H(→y+zβ→h√q(1))z.

It is shown in appendix A that at this equals

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c,(0) (36) =

which is a constant independent of . In addition as . For one finds at any temperature the contribution

 ¯¯¯¯¯¯¯¯¯m2h1J,c,(0)=q(1). (37)

Even though at this is the full result for the sample-to-sample fluctuations of the magnetization, at finite there will be an additional piece from obtained below. Similarly, to obtain the full finite- expression of higher-order moments of , contributions from are needed. However, here we focus on .

We now turn to the calculation of the contributions which capture the information about jumps, which are of order in the limit . It is contained in the contribution of and only in that contribution, as was discussed in Ref. FRGRSB, . Higher-order functions contain contributions of order at encoding information of multi-shock correlations. To first order in the non-linear term we find, extending the calculation in Ref. FRGRSB, :

 ϕ1(x,→y) = ∫1xdx′β22∑ijdq(x′)dx′x′hihj¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂ϕ0∂yi(x′,→y+ηβ→hDx′x)∂ϕ0∂yj(x′,→y+ηβ→hDx′x)η (38) = ∫1xdx′β22∑ijdq(x′)dx′x′hihj¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂H∂yi(→y+β→h[ηDx′x+z1D1x′])∂H∂yj(→y+β→h[ηDx′x+z2D1x′])η,z1,z2.

As in Eq. (30), , and are independent unit-centered Gaussian random variables, and .

We now change integration variables from and define and , the “thermal boundary layer variable” FRGRSB () for the external field. Using Eq. (28) the contribution of the first-order term to the magnetization cumulant becomes, denoting , and , see Fig. 1,

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c,(1) (39) = (−1)p+1∂^h1…∂^hpT2∫qcqmdq^x(q)∫∞−∞p∏i=1dyiδ(∑iαiyi)¯∂A+∂A−H(→y+→^hA+)H(→y+→^hA−)A+,A− = (−1)p∂^h1…∂^hpT2∫qcqmdq^x(q)∫∞−∞p∏i=1dyiδ(∑iαiyi)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∂A+∂A−12[H(→y+→^hA+)−H(→y+→^hA−)]2A+,A− .

are centered Gaussian random variables with correlations defined from the above independent Gaussian variables as

 A+=η√q−qm+z1√qc−q (40) A−=η√q−qm+z2√qc−q . (41)

It is convenient to introduce

 F:=A++A− ,G:=A+−A− (42)

in terms of which one can integrate by part

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c,(1) = (−1)p∂^h1…∂^hpT2qc∫qmdq^x(q)−∞∫−∞dF−∞∫−∞dG(∂2F−∂2G)exp(−F24[qc+q−2qm]−G24[qc−q])2π√2[qc+q−2qm]2(qc−q) (43) ×p∏i=1∫∞−∞dyiδ(∑iαiyi)12[H(→y+12→^h(F+G))−H(→y+12→^h(F−G))]2 .

The differential operator acts only on the Gaussian measure. Note that its action is equivalent to . One can thus integrate by part over to get

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c,(1) = (−1)p+1T4∂^h1…∂^hpqc∫qmdqd^x(q)dq (44) ×p∏i=1∫∞−∞dyiδ(∑iαiyi)¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯[H(→y+12→^h(F+G))−H(→y+12→^h(F−G))]2F,G .

The measure over and is defined by

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯f(F,G)F,G:=−∞∫−∞dF−∞∫−∞dGexp(−F24[qc+q−2qm]−G24[qc−q])2π√2[qc+q−2qm]2(qc−q)f(F,G). (45)

Equivalently, one can write in terms of and

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1…mhpJ,c,(1)=(−1)p+1T4∂^h1…∂^hpqc∫qmdqd^x(q)dqp∏i=1∞∫−∞dyiδ(∑iαiyi)¯[H(→y+→^hA+)−H(→y+→^hA−)]2A+,A− (46)

with measure

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯f(A+,A−)A+,A−:=−∞∫−∞dA+−∞∫−∞dA−exp(−(A++A−)24[qc+q−2qm]−(A+−A−)24[qc−q])π√2[qc+q−2qm]2(qc−q)f(A+,A−). (47)

The boundary terms in the integration by part vanish, provided that whenever exhibits a plateau for it is included as a function.

Using the expression (25) for , the formula (44) allows us to compute the thermal boundary-layer form of the -th cumulant. We give here the result for :

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯mh1mh2J,c,(1) (48) = −14∫qcqmdqd^x(q)dq∫∞−∞dG[exp(−G24[qc−q])√4π[qc−q]] ×G3(h1−h2)coth((h1−h2)G2T),

recovering the form obtained in Ref FRGRSB (). For and one finds

 (49)

Added to Eq. (37), this gives the correct total sample-to-sample fluctuations of the magnetization. The fact that higher terms do not contribute to this variance can be verified by a direct expansion of Eq. (24) in .

For general we only study the limit . For convenience we introduce the notation ,