# Self-organized criticality in the hysteresis of the Sherrington - Kirkpatrick model

\[

## Abstract

We study hysteretic phenomena in random ferromagnets. We argue that the angle dependent magnetostatic (dipolar) terms introduce frustration and long range interactions in these systems. This makes it plausible that the Sherrington - Kirkpatrick model may be able to capture some of the relevant physics of these systems. We use scaling arguments, replica calculations and large scale numerical simulations to characterize the hysteresis of the zero temperature SK model. By constructing the distribution functions of the avalanche sizes, magnetization jumps and local fields, we conclude that the system exhibits self-organized criticality everywhere on the hysteresis loop.

PACS numbers: 64.60.Lx, 75.60.Ej, 75.10.Nr

\]

Hysteresis in ferromagnetic systems is a century old physical problem. Efficient phenomenologies have already been developed , but an accepted microscopic theory is yet to be constructed. In soft magnets, where domain wall motion dominates the physics, considerable progress has been achieved recently . In hard magnets domain nucleation, domain wall motion and their interaction are all important. Hence they are better described on a more microsopic level as an assembly of strongly interacting spins or hysterons . Quantitative insight to such systems has been gained recently through studying the random field Ising model (RFIM) .

However a key aspect of the physics of real systems is missing from the RFIM: it does not include the long range dipolar (or magnetostatic) interactions. While these are negligible on atomic scales relative to the exchange term, they can dominate the collective behaviour of granular systems. This is so because the dipolar interaction is long ranged, so it involves every spin in the volume of the grains, whereas the exchange coupling scales only with the number of spins on the surface of the grain. These dipolar forces are important: they prevent the roughening of the domain walls  and determine the size of the domains . Crucially, the sign of these interactions changes with the angle. This introduces frustration into the system, which is not represented in the RFIM.

To capture the influence of frustration on hysteretic phenomena, we study the simplest system, containing long-range frustrated interactions, the Sherrington-Kirkpatrick (SK) model. Early numerical work demonstrated that this model exhibits hysteresis . However, in spite of its obvious importance, we could not find analytic studies of the hysteresis loop of the SK model. In this Letter we use scaling arguments, replica calculations and large scale numerical simulations to characterize the hysteresis of the zero temperature SK model. By constructing the distribution functions of the avalanche sizes, magnetization jumps and local fields, we conclude that the system exhibits self-organized criticality everywhere on the hysteresis loop.

The SK model consists of $`N`$ Ising spins ($`\sigma =\pm 1`$) on a fully connected lattice, described by the Hamiltonian

$$=\frac{1}{2}\underset{ij=1}{\overset{N}{}}J_{ij}\sigma _i\sigma _jh\underset{i=1}{\overset{N}{}}\sigma _i,$$
(1)   
where $`J_{ij}`$ is a random Gaussian number of zero mean and variance $`1/N`$. Throughout the paper we work at $`T=0`$.

First we summarize our numerical results. We start from a fully polarized state and change the external magnetic field $`h`$ adiabatically: for a given field we let all spins align according to their local field before varying $`h`$ again. During the avalanches we use sequential single spin flip updating to ensure the decrease of the total energy. The resulting hysteresis loop for the SK model is presented in Fig. 1. Finite size scaling analysis shows that the hysteretic trajectories are well-defined in the $`N\mathrm{}`$ limit, and the coercive field converges to a finite value.

We also analyzed the minor hysteresis loops of the SK model (inset Fig. 1). Within numerical accuracy they return to the major loop at the point of departure, exhibiting return point memory. This feature is present in many experimental systems, and it is also one of the criteria for the applicability of the Preisach phenomenology .

Next we establish some of the basic energy scales from elementary considerations. When spin $`\sigma _j`$ is flipped, the local field $`h_i`$ at another site changes by an amount proportional to $`2J_{ij}2/N^{1/2}`$. Thus the external field $`h`$ has to be changed by an amount $`dh1/N^{1/2}`$ to start a new avalanche. Now let $`S`$ be the change in the total magnetization during an avalanche, and $`dm=S/N`$ the jump of the magnetization $`m`$ during the avalanche. The average $`m(h)`$ curve is continuous and thus its derivative $`dm/dhS/N^{1/2}`$ is finite (Fig. 1), requiring: $`SN^{1/2}`$. This is possible only if the scale of the distribution of avalanches is set by $`N^{1/2}`$. This is characteristic of systems at criticality, whereas off-criticality the scale is set by some control parameter of the Hamiltonian. This leads to the central result of the paper: the SK model exhibits critical behaviour everywhere along its hysteresis loop. As this phenomenon is independent of the parameters of the Hamiltonian, it is a manifestation of self-organized criticality.

To elucidate this point, in Fig. 2 we show the distribution functions of $`S`$, and the number of spin flips in an avalanche (its “size”), $`n`$; $`𝒫(S)`$ and $`𝒟(n)`$, respectively, measured in the interval $`m[0.3,0.3]`$ for various system sizes. Both distributions exhibit power law behavior and can be well described by the finite size scaling forms:

$`𝒟(n)=(B/\mathrm{ln}N)n^\varrho d(n/N^\sigma ),`$ (2)   
$`𝒫(S)=(A/\mathrm{ln}N)S^\tau p(S/N^\beta ),`$ (3)   
with $`\tau ,\varrho =1\pm 0.1`$, $`\sigma =0.9\pm 0.1`$, and $`\beta =0.6\pm 0.1`$. The logarithmic prefactors were necessary to achieve satisfactory scaling collapse. Since such terms are needed only to keep distributions with an exponent $`1`$ normalized, this strongly suggests that $`\tau =\varrho =1`$ exactly. Unfortunately, because the cutoffs of the distributions $`𝒫(S)`$ and $`𝒟(n)`$ scale with different powers of $`N`$ the attractive picture of a diffusive motion of the local fields due to the randomness in $`J_{ij}`$  would lead to an infinite diffusion constant $`Dn/N^{1/2}`$ and is thus inapplicable.

Adopting the $`\tau =1`$ equality and combining it with $`SN^{1/2}`$ immediately yields the relation $`\beta =1/2`$, with logarithmic corrections, in good agreement with the above measured value. Also, because the $`J_{ij}`$’s take negative values as well, spins of both signs are destabilized in an avalanche. Therefore the number of participating spins is only bounded from below by $`S/2`$, yielding the exponent - bound $`\sigma \beta =1/2`$. An upper bound for $`\sigma `$ can be obtained from estimating the dissipated energy, $`E_d`$, during a finite but small sweep of the external-field $`h_1h_2=h_1+\mathrm{\Delta }h`$: $`E_d=Nm\mathrm{\Delta }hN`$. Also, since the average energy dissipation per spin is at least $`2dh1/N^{1/2}`$, $`E_d`$ can be estimated as $`E_d>2n_{total}/N^{1/2}`$, where $`n_{total}`$ is the number of flips during all avalanches from $`h_1`$ to $`h_2`$. But the number of avalanches during this sweep is proportional to $`\mathrm{\Delta }h/dhN^{1/2}`$, i.e. $`n_{total}N^{1/2}n`$. Combining all the above gives $`NE_d>nN^\sigma `$ implying the upper bound $`\sigma 1`$, which is nearly saturated according to our numerics.

The above distributions imply that the average value of $`\chi dm/dh`$ is dominated by a few very large avalanches, whereas its typical value scales to zero as $`1/N^{1/2}`$, which we confirmed independently numerically. Therefore the hysteresis loop for a specific disorder realization has a slope zero with unit probability, interrupted by a few macroscopically large avalanches. This feature is characteristic of the Barkhausen noise and establishes the frustrated spin glasses as possible candidates to describe certain classes of hysteretic magnets.

We also studied the correlations of consecutive avalanches. We measured the Hausdorf dimensions of the numerically determined hysteresis loop and that of a sequence of independent avalanches, generated with the above distributions. Having found the two Hausdorf dimensions equal suggests that avalanches are uncorrelated.

These results, in particular the size $`N`$ as the sole cutoff of the different distribution functions, which all exhibit power law behaviour, confirm the above - stated self-organized criticality of the entire hysteresis loop of the SK model. To shed more light on the underlying physics we explore the local-fields, $`h_i=_jJ_{ij}\sigma _j+h`$ by studying the local stabilities, $`\lambda _i=\sigma _ih_i`$, which are all positive for stable spin configurations. Their distribution, $`P(\lambda )`$ is shown in Fig. 3.

Remarkably – unlike the local field distribution $`P(\lambda )`$ is essentially the same at any point of the hysteresis loop. This suggests that the avalanche dynamics of the SK model organizes the system into special states with similar properties everywhere along the hystersis loop. A careful finite size analysis shows that $`P(\lambda =0)1/\sqrt{N}`$ and $`P(\lambda )C\lambda ^\alpha `$ with $`C\alpha 1`$ for small $`\lambda `$’s. As we now show, this latter result establishes once again that these special states are critical. To prove this let us flip $`n_{flip}`$ arbitrary spins starting from a given stable spin configuration $`\{\sigma _i\}`$ with $`\lambda _i>0`$ and calculate the average number of new unstable spins, $`n_{unst}`$, distinguished by negative stabilities $`\lambda _i^{}=\lambda _i+\mathrm{\Delta }\lambda _i<0`$:

$$\lambda _i^{}=\lambda _i2\underset{j\mathrm{flipped}}{}\sigma _iJ_{ij}\sigma _j.$$
(4)   
The system is critical if $`n_{unst}=n_{flip}`$, as for $`n_{unst}<n_{flip}`$ the avalanches die out exponentially fast while in the opposite case they explode . Assuming that the $`n_{flip}`$ random terms at the rhs. of Eq. (4) are independent, the probability $`P_d`$ of destabilizing a given spin is:

$$P_d=_0^{\mathrm{}}𝑑\lambda P(\lambda )_{\mathrm{}}^\lambda d(\mathrm{\Delta }\lambda )Q(\mathrm{\Delta }\lambda ),$$
(5)   
where $`Q(\mathrm{\Delta }\lambda )=\mathrm{exp}\{N\mathrm{\Delta }\lambda ^2/8n_{flip}\}\sqrt{N/8\pi n_{flip}}`$ is the probability distribution of the $`\mathrm{\Delta }\lambda `$ term in Eq. (4), and $`P(\lambda )`$ is approximated by its asymptotic form, $`P(\lambda )=C\lambda ^\alpha `$. The average number of destabilized spins is then $`n_{unst}=NP_d=\stackrel{~}{C}(\alpha )N(n_{flip}/N)^{(\alpha +1)/2}`$, with $`\stackrel{~}{C}(\alpha )`$ an $`\alpha `$-dependent constant, $`\stackrel{~}{C}(1)=C`$. For $`\alpha >1`$ (or $`\alpha =1`$ and $`C<1`$) $`n_{unst}<n_{flip}`$ and the system cannot give rise to large avalanches. On the other hand, for $`\alpha <1`$ (or $`\alpha =1`$ and $`C>1`$) $`n_{unst}>n_{flip}`$, and the state is unstable. Thus the criticality condition is characterized by $`\alpha =1`$ and $`C=1`$. These are exactly the values found in our numerical simulations, once again underlining the criticality of the system.

The physical mechanism of self-organized criticality can be qualitatively understood as follows. As the avalanche rolls, at any given time step $`t`$ the stabilities of the spins are shifted only by those spins, which changed sign at step $`t1`$. These spins have flipped because the second term of Eq. (4) for their stabilities was negative, pulling their $`\lambda _i`$’s downward. However once $`\lambda _i`$ changed sign, the very same term now enhances this stability. More importantly, this term being symmetric, it also pushes upward the stabilities of the other spins of the avalanche, which pulled spin $`i`$ down and flipped it in the first place. This effect is suppressing the density of states with low local fields, reminescent of the formation of the Coulomb gap in the disordered electron problem . The stabilities of the spins which did not participate in the avalanche will be shifted by a random amount by the just-flipped spins. However in the presence of a slope in their distribution $`P(\lambda )`$, this will have a net effect, moving the stabilities of more spins downward than upward. In short, correlations between the spins of an avalanche move the stabilities of the already flipped spins upward; at the same time the random couplings between all spins drive a net downward drift. The competition of these two forces keeps the system critical.

To understand the shape of the measured major hysteresis loop more in detail we first observe that the states where an avalanche stops must always be single spin-flip stable (SSS). Let us therefore define the average number of SSS states,

$$V(m,h)=\mathrm{Tr}\left\{\underset{i=1}{\overset{N}{}}\mathrm{\Theta }(\lambda _i)\delta (mN\underset{i}{}\sigma _i)\right\},$$
(6)   
where the angular bracket indicates annealed disorder averaging and the trace stands for the summation over all spin configurations. The product of the theta functions in Eq. (6) projects out those states where all the spins have positive stabilities, while the delta function selects states with a given magnetization.

Using the integral representations $`\delta (y)=_{\mathrm{}}^{\mathrm{}}\frac{dx}{2\pi }e^{ixy}`$ and $`\mathrm{\Theta }(\lambda _i)=_{\mathrm{}}^{\mathrm{}}\frac{idz_i}{2\pi z_i}e^{iz_i\lambda _i}`$ the function $`V(m,h)`$ can be rewritten in an exponential form, $`V(m,h)=\frac{dx}{2\pi }_i\left[\frac{idz_i}{2\pi z_i}\mathrm{exp}\{i(\lambda _iz_ix\sigma _i+mx)\}\right]`$, and the disorder average and the spin summation can be easily carried out. After the disorder averaging the effective action contains a term proportional to $`(_iz_i)^2`$. Decoupling this term with a new Hubbard-Stratonovich field $`R`$, one finally arrives at the following expression:

$`V=\sqrt{N}{\displaystyle \frac{dx}{2\pi }\frac{dR}{\sqrt{2\pi }}\mathrm{exp}\left\{N[\mathrm{ln}Q\frac{R^2}{2}imx]\right\}},`$ (7)   
$`Q(h,x,p)={\displaystyle \frac{idz}{2\pi z}\mathrm{\hspace{0.33em}2}\mathrm{cos}(xzh)e^{\frac{1}{2}z^2iRz}}.`$ (8)   
The above integral can readily be evaluated in the $`N\mathrm{}`$ limit by the saddle point method. The saddle point equation, $`Q/x=imQ`$ can be solved analytically, and the variable $`x`$ can be completely eliminated resulting in the following expression for $`V(m,h)`$:

$`V(h,m)\mathrm{exp}\left(N\mathrm{\Omega }_{sp}(m,h)\right),`$ (9)   
$`\mathrm{\Omega }_{sp}={\displaystyle \frac{1m}{2}}\mathrm{ln}{\displaystyle \frac{1\varphi ^{}}{1m}}+{\displaystyle \frac{1+m}{2}}\mathrm{ln}{\displaystyle \frac{1+\varphi ^+}{1+m}}{\displaystyle \frac{R^2}{2}},`$ (10)   
where $`R`$ is determined from $`\mathrm{\Omega }_{sp}(m,h,R)/R=0`$, and $`\varphi ^\pm =\varphi (\frac{h\pm R}{\sqrt{2}})`$ with $`\varphi (x)=\frac{2}{\sqrt{\pi }}_0^xe^{t^2}𝑑t`$.

In Fig.4. we plotted the contour of $`\mathrm{\Omega }_{sp}=0`$. Outside this line the density of SSS states scales to $`0`$ exponentially, thus they are definitely unable to arrest the avalanches. Inside this line the number of SSS states is exponentially large, and is thus comparable to the total number of states, themselves expoential in $`N`$. Therefore avalanches get trapped with a higher probability in one of the SSS states. Hence the $`\mathrm{\Omega }_{sp}=0`$ contour constitutes a strict outer bound for the true hysteresis loop.

Comparing Figs.1 and 4 shows that the $`\mathrm{\Omega }_{sp}=0`$ contour considerably overestimates the size of the hysteresis loop. We pursued two refinements of this calculation. We developed a replica symmetric description as well as a de Almeida - Thouless type replicon instability analysis: these will be reported separately .

Finally, we briefly discuss the effect of a finite ferromagnetic coupling, $`J_0>0`$. $`J_0`$ simply shifts the value of the magnetic field $`hh+J_0m`$ in Eq. (8), and results in a shear  of the entire contour $`\mathrm{\Omega }_{sp}=0`$. For a branch of the hysteresis loop $`m`$ must be a monotonic function of $`h`$. Since the $`\mathrm{\Omega }_{sp}=0`$ contour is an outer bound of the hysteresis loop, therefore, when this loop would force a non-monotonic $`m(h)`$ relation (Fig.4), the major hysteresis loop must develop a finite jump. Since the slope of the major hysteresis loop for $`J_0=0`$ is finite, one expects this transition to occur at a finite critical coupling, $`J_0=J_c`$. Our numerical data agree with this picture .

We end with a comparison to the random field Ising model. In that model our initial simple scaling considerations yield $`S𝒪(1)`$, i.e. a non-critical avalanche distribution. As shown in Ref., there is only a single critical point at the coercive field at some specific value of the disorder. Therefore the distribution functions exhibit a scaling behaviour, with the cutoff set by the distance from this critical point, rather than by the system size, as happens for the SK model. In short, the RFIM exhibits “plain old criticality” , while the SK model exhibits self-organized criticality. Also, the avalanche distribution exponents $`\tau `$ in the two models are different: on the mean field level $`\tau =1.5`$ for the RFIM and $`1.0`$ for the SK model. In finite dimensions the exponents typically increase: numerical simulations of the 3D RFIM found $`\tau =1.6`$ . In contrast, numerical studies of realistic 3D models , as well as experimental  works report $`\tau `$ in the $`1.1\mathrm{}1.4`$ regime. This raises the possibility that the finite dimensional extensions of the frustrated models might provide a $`\tau `$ closer to the experimental values.

In sum, we studied the hysteretic behaviour of the SK model. We determined numerically the distribution functions of the avalanches, the magnetization jumps, and the local fields. The model exhibits self-organized criticality everywhere along the hysteresis loop. We recalculated the loop with analytic methods as the location of one-spin-flip stable states, and found satisfactory agreement with the numerical results.

We acknowledge useful discussions with R. Scalettar, G. Bertotti, E. Della Torre, G. Kádár, M. Mézard and J. Sethna. This research has been supported by the U.S - Hungarian Joint Fund 587, Hungarian Grants OTKA T026327 and OTKA D29236 and by NSF DMR 95-28535. GTZ acknowledges the kind hospitality of the George Washington University.