# A statistical mechanical description of metastable states and hysteresis in the 3D soft-spin random-field model at

###### Abstract

We present a formalism for computing the complexity of metastable states and the zero-temperature magnetic hysteresis loop in the soft-spin random-field model in finite dimensions. The complexity is obtained as the Legendre transform of the free-energy associated to a certain action in replica space and the hysteresis loop above the critical disorder is defined as the curve in the field-magnetization plane where the complexity vanishes; the nonequilibrium magnetization is therefore obtained without having to follow the dynamical evolution. We use approximations borrowed from condensed-matter theory and based on assumptions on the structure of the direct correlation functions (or proper vertices), such as a local approximation for the self-energies, to calculate the hysteresis loop in three dimensions, the correlation functions along the loop, and the second moment of the avalanche-size distribution.

###### pacs:

75.10.Nr, 75.60.Ej, 64.60.av## I Introduction

The random-field Ising model (RFIM) at zero temperature is a prototype of many disordered systems which exhibit hysteretic and jerky behavior when an external parameter (e.g. magnetic field, pressure, strain) is slowly changedSDP2006 (). This behavior is related to the existence of a corrugated (free) energy landscape with many local minima (or metastable states) separated by barriers much larger than , therefore preventing relaxation towards equilibrium on experimental time scales. As a consequence, the response to an external driving field is a series of discontinuous jumps (called avalanches) from one metastable state to another, the number and size of these jumps varying with the amount of disorder. In the 3-dimensional RFIM with Gaussian random fields two regimes of avalanches are observedS1993 (); DS1996 (): at large disorder, there are many microscopic jumps resulting in a smooth magnetization curve macroscopically; at low disorder, there is a system-spanning avalanche resulting in a macroscopic jump in the magnetization curve. These two regimes are separated by a critical point at which avalanches of all sizes are observed. As shown recently, this type of nonequilibrium disorder-induced phase transition may underly the hysteretic behavior of He adsorbed in high porosity silica aerogelsDKRT2005 (); BLCGDPW2008 ().

Even though the RFIM and similar models have been extensively studied, there is currently no analytical tool to compute the saturation hysteresis loop in finite dimensions (even approximately) and no theory to describe the statistical properties of the metastable states, for instance their configurational entropy (also called ‘complexity’) as a function of magnetization. Furthermore, the behavior of the correlation functions along the hysteresis loop is unknown, although this is an issue of experimental interestDKRT2006 (). Of course, describing analytically such a nonequilibrium situation where the present state of the system depends on its past history is a difficult problem and at first sight there seems to be no other choice than to follow the dynamical evolution for given initial conditions (e.g. one of the two saturated states corresponding to infinitely positive or negative magnetic field). In this way, one can treat exactly ‘mean-field’ systems, i.e. fully connected latticesS1993 (); DS1996 () or lattices with a locally tree-like structure such as the Bethe latticeDSS1997 (); SSD2000 (), and obtain analytical expressions for the saturation hysteresis loop or the avalanche size distribution. To go further and build a field-theoretical description of these phenomena it then seems necessary to employ the Martin-Siggia-Rose formalismDS1996 (). In ferromagnetic systems, however, an alternative route is possible thanks to a remarkable property of the saturation hysteresis loop in the large-disorder regime: the loop is the convex hull of the metastable states in the two-dimensional field-magnetization planeDRT2005 (); PRT2008 (); RTP2009 (); RM2009 (). In other words, the complexity is zero outside the loop and positive insidenote1 (). Moreover, it exactly vanishes along the loop since there is a single ‘extremal’ metastable state at a given field (at least for a continuous distribution of the random field). Determining the hysteresis loop is then tantamount to counting the number of metastable states at fixed field as a function of their magnetization and finding the magnetization at which the complexity vanishes. The main difficulty is that one must count the typical (i.e. most likely) number of metastable states and compute the associated quenched complexity, which of course is a nontrivial task requiring the use of replicas. Nevertheless, at least in principle, one can thereby study hysteresis and avalanche statistics without referring to the dynamical evolution.

In the present paper we test this approach by studying the soft-spin version of the RFIM introduced in Ref.DS1996 (), where the spins take on continuous values between and and are confined by a bistable potential with a linear cusp. This model has the advantage over the standard RFIM to show hysteresis for all values of the disorder in mean-field theory and is anyhow expected to behave as the hard-spin model on long length scales. Our basic strategy is to express the complexity as the Legendre transform of the ‘free-energy’ associated to a certain action in replica space and then to use approximations borrowed from condensed-matter theory to calculate the corresponding correlation (or Green’s) functions and thermodynamic properties. Our focus on the correlation functions is also motivated by analytical and numerical calculations on the Bethe lattice which suggest that their spatial structure along the hysteresis loop is simpler than at equilibriumIR2010 ().

The paper is organized as follows. In section II, we define the model and the observables that we want to compute. In section III, we introduce the general formalism in the context of a replica method, focusing on the various correlation functions. In section IV, we first consider the random-phase approximation (RPA) which is equivalent to mean-field theory and becomes exact in infinite dimension. In section V, we go beyond the RPA by introducing a local self-energy approximation (LSEA) and we compare the predictions for the magnetization, the correlation functions, and the second-moment of the avalanche-size distribution along the hysteresis loop to simulation data. Concluding remarks and directions for future work are provided in Section VI. Additional details on the analytical calculations are provided in the appendix.

## Ii Model and observables

We consider a collection of soft spins placed on the sites of a -dimensional hypercubic lattice and interacting via the Hamitonian

(1) |

where is a ferromagnetic interaction that couples nearest-neighbor spins, is an external uniform field, and is a set of quenched random fields drawn independently from a Gaussian probability distribution with zero mean and variance . is a double-well potential for which we choose the same cuspy form as in Ref.DS1996 (); RM2009 ()

(2) |

so that all solutions of , i.e. all stationary points of the Hamiltonian, are local minima. This greatly facilitates the present study as there is no need to explicitly discard local maxima and saddle-points through the consideration of the Hessian of the energy function. These local minima are the so-called ‘metastable’ states in which each spin satisfies

(3) |

where is a neighbor of on the lattice.

At zero-temperature, one can define a local relaxation dynamics in which each spin is forced to satisfy Eq. 3 as the external field is changedDS1996 (). When adiabatically varying from to and backwards, the model exhibits hysteresis and the magnetization typically behaves as shown in Fig. 1. (The figure displays the results of a single simulation on a cubic lattice with and note5 ().) These values are arbitrarily chosen and will stay fixed in the rest of this work. The shape of the loop then changes with the coupling : one can see that and correspond to the large- and small-disorder regimes respectively, with a macroscopic jump in the latter case.

For a given realization of the disorder, i.e. a set of random fields , and a given value of the external uniform field , each metastable state is characterized at a macroscopic level by its magnetization, its energy, etc. In the present study, we only focus on the magnetization since this is sufficient to unambiguously determine the hysteresis loop. Our goal is then to compute properties averaged over all metastable states with a given magnetization per site at a given external magnetic field (with a flat measure). A central quantity is the quenched complexity , which is defined as

(4) |

where is the number of metastable states with magnetization at the field and the overbar denotes an average over the random-field distribution. We are also interested in the following two-point correlation (Green’s) functions:

(i) the spin-spin correlation functions,

(5) |

with the (disorder) connected and disconnected components defined as

(6) |

and

(7) |

the brackets denoting an averaged over all metastable states with magnetization at the field ,

(ii) the correlation function involving the spin variable and the random field,

(8) |

(iii) the spin-spin correlation function for two copies of the same disordered system coupled to different external fields and with different magnetizations,

(9) |

where the subscrit indicates that are fixed and similarly for the subscript .

Along the hysteresis loop, the complexity is zero (at least in the large disorder regime where the loop is continuous)DRT2005 (); PRT2008 (); RTP2009 (); RM2009 () and the connected spin-spin correlation is identically zero, as the fluctuations inside a metastable state vanish at zero temperature. On the other hand, the disconnected spin-spin correlation function and the spin-random-field one should be nontrivial.

Finally, we are also interested in the distribution of avalanche sizes along the loop (say, along the ascending branch). For a given disorder realization, the magnetization curve consists of smooth parts resulting from the harmonic response to the magnetic field and a series of jumps of size occurring at the fields (these jumps are not visible in Fig. 1 due to the scale of the figure) note2 (). The magnetization can thus be decomposed as

(10) |

where is the Heaviside step function. From this, one defines the jump (avalanche) density

(11) |

so that is the number of avalanches of size between and when the field is increased from to (note that is measured per site). In present work we will only compute the ‘unnormalized’ second moment (it is ‘unnormalized’ because as such is not a probability density and is the total number of avalanches between and ).

## Iii Formalism

In order to control the local magnetization, we introduce an additional source that is linearly coupled to the spins and we consider the following (disorder-dependent) ‘partition function’ in an external magnetic field which is momentarily taken as nonuniform, :

(12) |

where the symbol refers to the integration over all the spin variables, , and no Jacobian is needed (it would merely introduce a constant multiplicative factor). This partition function can be put into a more standard form by replacing the Dirac delta function by its Fourier representation,

(13) |

where the action is defined by

(14) |

and are auxiliary (imaginary) variables; for conciseness, the factor associated to the integration of along the imaginary axis is adsorbed into .

This defines a ‘free energy’ which is a random object whose cumulants give access to full information about the system, including the complexity and the correlation functions. The complexity is obtained from the first cumulant,

(15) |

via a Legendre transform, where

(16) |

and which for uniform sources takes the form

(17) |

(Here and below we use square brackets when the arguments are locally varying and parenthesis when they are uniform.) As discussed in Refs.DRT2005 (); PRT2008 (); RTP2009 (); RM2009 (), the hysteresis loop in the large-disorder regime identifies with the curve in the limit whereas the typical properties of the metastable states are obtained for which corresponds to the maximum of the complexity.

The information about the distribution of avalanche sizes is contained in the higher-order cumulants ,…, where

(18) |

etc…

The correlation (Green’s) functions are obtained by derivation of the cumulants with respect to the sources. For instance, the physical two-point spin-spin correlation functions introduced above are given by

(19) |

(20) |

where both right-hand sides are evaluated for uniform sources and is considered as a function of and through the Legendre transform in Eq. (17). The spin-random-field correlation function requires a little more thought. By using the property of Gaussian distributions,

(21) |

one finds

(22) |

where

(23) |

When all the sources are uniform, this yields

(24) |

where is the Fourier transform of . (In particular, this equation is valid in the limit , i.e. along the hysteresis loop: surprisingly, it seems that this extension of the ‘susceptibily sum-rule’ to the nonequilibrium magnetization curve has not been noticed before.)

To compute the average over disorder, the common procedure is to replicate the system times and take the limit at the end. (This is in contrast with the Martin-Siggia-Rose formalism in which the partition function is directly averaged over disorder, making the use of replicas unnecessaryDS1996 (); here indeed, the partition function in Eq. (13) is nontrivial so that one must average and not simply .) If one is interested in computing the cumulants of the random free-energy for generic arguments, one must introduce sources that act separately on each replicaLW2003 (); TT2008 (); MT2010 (). After performing the average over the random-field distribution, we obtain a ‘replica partition function’,

(25) |

with the replicated action given by

(26) |

The ‘thermodynamic potential’ can then be expanded in increasing number of free replica sumsLW2003 (); TT2008 (); MT2010 (),

(27) |

where the ’s are continuous and symmetric functions of their arguments. This coincides with the cumulant expansionTT2008 (). The thermodynamic potential generates the ‘magnetizations’ and

(28) |

and the correlation (or Green’s) functions, e.g. at the pair level,

(29) |

where denotes an average over the replicated action. As above, the magnetizations and the correlation functions can be expanded in increasing number of free replica sums:

(30) |

and similarly for , whereas after decomposing the two-point functions as (where does not contain any explicit Kronecker delta), one has TT2008 (); MT2010 ()

(31) |

(32) |

and similarly for and .

In the present approach, however, the central object is not but the ‘effective action’ which is the Legendre transform of with respect to the two sets of sources and ,

(33) |

so that

(34) |

is the generating functional of the ‘direct’ correlation functions [or one-particle irreducible (1PI) functions, or else proper vertices, in field-theoretic language]. At the pair level,

(35) |

As above, , , , etc… can be expanded in increasing number of free replica sums.

In the following, all the two-point functions will be put together as components of a matrix with both replica indices and spatial coordinates

(36) |

where , , have for elements , , ; a similar notation is used for the direct correlation functions, all collected in . The matrix is then just the inverse of , i.e.

(37) |

The matrix inversion with respect to spatial coordinates is easily realized by Fourier transformation when all the sources are taken as uniform. The inversion with respect to the replica indices can be performed by using the expansion in free replica sums for and (see Eqs. (31,III)) and proceeding to a term-by-term identification. The zeroth-order terms are then given by

(38) | |||

(39) |

where are matrices containing the components , , as in Eq. (36), and similarly for . Note that the zeroth-order functions obtained when fixing the replica sources and fixing the replica ‘magnetizations’ are generically related through with and being the zeroth-order expressions as in Eq. (28). Since the zeroth-order contributions already contain the physics of the problem, we will not consider higher-order terms and will drop the superscript in the following.

When all replica sources or magnetizations are taken as equal, and provided that this limiting process is regular enough (this will be discussed in more detail below), replica symmetry is recovered and one obtains a set of ‘Ornstein-Zernike’ (OZ) equations (in the language of liquid-state theoryHM1986 ()) which takes the following explicit form:

(40) |

and

(41) |

where all functions depend on and .

When the replica sources are not equal, the zeroth-order terms have the capability to capture a nonanalytic behavior in the dependence on the replica sources or replica magnetizations. This important feature is already displayed by the noninteracting system corresponding to (hereafter called the ‘reference’ system) whose properties are listed in the appendix (Eq. (111) and below with and ). For instance, the 2-replica function behaves like

(42) |

when . We stress that this is not a spurious behavior due to the presence of a cusp in . Indeed, even if were a smooth double-well potential, it would be necessary to distinguish between the two minima in order to count the metastable states and, whatever the method, this would introduce a nonanalyticity in and lead to a cusp in the 2-replica function . This behavior is thus an intrinsic feature of the problem under consideration and is also intimately connected to the magnetization discontinuities along the hysteresis loop, i.e. to the avalanches. A similar connection is discussed in Ref. LW2009 () in the context of random elastic systems where the statistics of static avalanches (or ‘shocks’) is studied via the functional renormalization group. (In this case, however, the cusp only shows up in the course of the renormalization flow whereas it is always present here, even in the large disorder regime.)

To show that the cusp is related to the (unnormalized) second moment of the avalanche size distribution, we consider the quantity

(43) |

where it is implicit that the limit corresponding to the hysteresis loop (in the large-disorder regime) has been taken. Then, using Eqs. (10) and (42), taking the second derivative with respect to and in the limit , and identifying the singular contributions proportional to on both sides of the equation lead to

(44) |

The nonanalyticity in also implies that the 2-replica correlation function contains a singular contribution proportional to . This induces singular contributions in the three disconnected direct correlation functions via the OZ equations and leads to formally diverging terms proportional to ‘’ when the replica fields are equal. The function , however, has no obvious physical meaning and this problem is in principle harmless, although it may be a source of difficulties in an approximate treatment, as will be discussed below in section V.

Finally, it is instructive to examine the behavior of the correlation functions in the limit , i.e. along the two branches of the hysteresis loop in the large-disorder regime. Assuming that the general behavior is the same as in the reference system (which seems quite reasonable, at least in the regime under consideration), we find

(45) |

and

(46) |

where the notation indicates that both the regular and the singular contributions of and are of order . Note that we have considered the dependence on the source . A similar dependence is found on when the latter goes to plus or minus infinity (recall that we work at the zeroth-order level of the expansion in free replica sums).

It is not surprising that vanishes as . This is due to the already mentioned fact that there is only one metastable state along the loop at a given field so that the statistical fluctuations only come from the quenched disorder and are thus contained in the disconnected functions. (As a consequence, also vanishes.) The important feature is that and vanish exponentially fast with . Indeed, this implies that the OZ equations for the physically relevant functions and take in this limit the much simpler form:

(47) |

Note that does not contain any diverging part when (or ), as it must be.

## Iv Random phase approximation (RPA)

Equipped with the formalism of the previous section, we wish to introduce approximations based on assumptions on the structure of the direct correlation functions (proper vertices). Formally, these functions may be written as

(48) |

where is the ‘bare’ inverse propagator matrix (not taking into account the local potential ), whose components satisfy

(49) |

and is a ‘self-energy’ matrix (the minus sign is chosen so to match the usual definition of this quantity in field theory and many-body physics). In the above equations, is if and are nearest neighbors on the lattice, zero otherwise.

From now on in this section, except if explicitly stated, we only consider uniform sources that are equal for all replicas. Our first approximation, akin to the usual random phase approximation (RPA), consists in assuming that the self-energies are identical to those of the noninteracting () reference system. In Fourier space, this then leads to

(50) |

where is the connectivity of the hypercubic lattice and its characteristic function. The direct correlation functions of the reference system are obtained from the Green’s functions computed in the appendix [Eqs. (116)-(118) and (A)-(135) with and ]. The RPA effective action is immediately obtained by integrating the ‘susceptibility sum-rule’ , which yields

(51) |

From this expression, one obtains

(52) |

It can be easily checked that is identical to the quenched complexity (not to be confused with a self-energy) obtained in the mean-field model of Ref. RM2009 () (with replaced by ). (On the other hand, the auxiliary field does not coincide with the parameter introduced in this reference. In particular, satisfies the Legendre relation , in contrast with . In this respect, the RPA is a nicer way of obtaining the mean-field limit.)

One expects the RPA to be valid when the coupling is sufficiently weak, or, equivalently, when and are sufficiently large (all ‘thermodynamic’ quantities can be expressed in terms of the reduced variables , , and ). This is indeed what is observed in Fig. 2 where the predictions of the RPA for the magnetization along the ascending branch of the hysteresis loop are compared to numerical simulations (the descending branch is obtained using the symmetry , ). The theoretical value is given by