Influence of the driving mechanism on the response of systems with athermal dynamics: the example of the random-field Ising model

We investigate the influence of the driving mechanism on the hysteretic response of systems with athermal dynamics. In the framework of local-mean field theory at finite temperature (but neglecting thermallly activated processes), we compare the rate-independent hysteresis loops obtained in the random field Ising model (RFIM) when controlling either the external magnetic field $H$ or the extensive magnetization $M$. Two distinct behaviors are observed, depending on disorder strength. At large disorder, the $H$-driven and $M$-driven protocols yield identical hysteresis loops in the thermodynamic limit. At low disorder, when the $H$-driven magnetization curve is discontinuous (due to the presence of a macroscopic avalanche), the $M$-driven loop is re-entrant while the induced field exhibits strong intermittent fluctuations and is only weakly self-averaging. The relevance of these results to the experimental observations in ferromagnetic materials, shape memory alloys, and other disordered systems is discussed.


I. INTRODUCTION
When slowly driven by an external force, systems with quenched-in disorder often respond via a sequence of random bursts or avalanches whose size distribution extends over several decades, ranging from microscopic to macroscopic. A classic example is the Barkhausen noise that occurs in soft magnetic materials as one varies the external magnetic field H. The random sequence of spikes observed in the induced voltage (that is, essentially, in the time-derivative of the induced magnetization M ) is then related to the fluctuations in the motion of the magnetic domain walls due to structural disorder [1]. Avalanche dynamics is also observed in other systems, as for instance in the acoustic emission associated to martensitic transformations [2] or during the drainage of superfluid 4 He from a mesoporous solid [3]. All these systems also exhibit rate-independent hysteresis when the external force is cycled up and down, thermal fluctuations being negligible on the experimental time scale due to the presence of very high energy barriers (the systems then explore a set of long-lived metastable states).
Avalanches are characteristic of a system driven by an external force. In some experimental situations, however, it may be more convenient to control the conjugated extensive variable (or its time-derivative) rather than the external field. In magnetic materials for instance, some suitable feedback control may be used so to impose the flux (or its time derivative) instead of the magnetic field. This allows to study materials that reverse their magnetization very rapidly [4]. In gas adsorption experiments, there are also cases where one controls the amount of gas introduced inside the porous material and measures the induced change in the pressure [5]. In shape memory alloys, stress-strain curves are routinely obtained by controlling the deformation of the sample and measur-ing the stress [6]. The following question then naturally arises: what is the relation between the force-response diagrams (magnetization curves, adsorption isotherms, strain-stress curves,etc...) observed in the two situations, say "field-driven" and "magnetization-driven" (hereafter, we shall adopt the language of magnetic systems). This question is especially intriguing when there is a macroscopic instability, e.g. a macroscopic jump in the magnetization curve obtained with the field-driven protocol.
In this case, one may force a system that would spontaneously jump to a stable (or metastable) situation to remain into some other metastable or maybe unstable state. This may lead to a re-entrant behavior in the hysteresis loop, as indeed observed experimentally [4,7,8].
To what extent this behavior reveals an intrinsic property of the free-energy landscape of the system or depends on the experimental protocol is not clear. To the best of our knowledge, this problem has not been seriously investigated in the literature, in spite of the intense theoretical activity focusing on out-of-equilibrium driven systems.
In order to address this issue, we consider here a microscopic model, the Random Field Ising Model (RFIM), which has been extensively studied at zero temperature as a prototype of disordered hysteretic system with avalanche behavior [9]. In particular, in dimension 3, the RFIM is known to exhibit an out-of-equilibrium phase transition between a strong-disorder regime where the magnetization hysteresis loop is smooth on the macroscopic scale and a low-disorder one where it displays a macroscopic jump. The nonequilibrium RFIM at T = 0 is probably the simplest model to study the connection between metastability and hysteresis, and we expect its behavior in the magnetization-driven case to be relevant to a wide class of systems. However, as is discussed elsewhere [10], the definition of a metastable M -driven dynamics for Ising spins at T = 0 is somewhat problem-atic, and it is more convenient to work with continuous variables. This can be done in various ways and here we choose to study the RFIM at finite temperatures, using a local mean-field approach in which each Ising spin is replaced by its thermal average (see e.g. Ref. [11]). One then follows the evolution of the system among the local minima of the free-energy landscape, discarding all thermally activated processes that would allow for energy-barrier crossings. All the important features of the behavior of the RFIM at zero temperature (hysteresis, avalanches, and out-of-equilibrium phase transition) are thus preserved. This approach is expected to be valid on an intermediate time scale during which the system can equilibrate locally but is not able to reach the global thermal equilibrium. In many disordered systems, the free-energy barriers are much larger than k B T and this intermediate time scale is the one probed by experiments.
In the next section, we describe the model and briefly present the local mean-field approach in the field-driven case. The algorithmic method for the magnetizationdriven situation is presented in section III and the numerical results in the strong and weak-disorder regimes are discussed in section IV. Some concluding remarks are given in section V.

II. H-DRIVEN DYNAMICS FOR THE RFIM AT FINITE TEMPERATURE
We consider the ferromagnetic (Gaussian) RFIM defined by the energy function where S i (i = 1, ..., N ) are Ising spins placed on the sites of a 3d cubic lattice of size L × L × L and h i are quenched random fields drawn from a Gaussian distribution with zero mean and standard deviation σ. The first sum runs over all distinct pairs of nearest-neighbor sites.
In the "unconstrained" situation that was considered in all previous studies [9], the system is submitted to an external field H that is changed adiabatically. At zero temperature, it then tries to follow the field evolution by minimizing the magnetic enthalpy H = U − HM where M = i S i is the overall magnetization. The relaxation dynamics generally used in this case is the T = 0 version of the Glauber single-spin-flip dynamics where each spin is aligned to its local field f i + H, where and the sum is over the nearest neighbors of site i (this corresponds to a local, partial minimization of H [12]). A configuration {S i } is then (meta)stable according to the single-spin-flip dynamics when all the spins satisfy the stability condition When H is changed from an initially stable state, a spin may become unstable, which may initiate an avalanche of other spin flips in the neighborhood. The avalanche ends when a new stable state is reached. The resulting hysteresis loop (after averaging over disorder and extrapolating to the thermodynamic limit) is continuous for σ > σ c , with all avalanches being of microscopic size, and discontinuous for σ < σ c , the discontinuity corresponding to a macroscopic avalanche that involves a finite fraction of the spins. The critical disorder for the RFIM on the cubic lattice has been numerically located at σ c ≃ 2.2 [9,13]. Extension of this behavior to finite temperature is straightforward when using a local mean-field approach. In this approximation, the Helmholtz free-energy of the system for a given realization of the random fields is a function of the thermally averaged local magnetizations m i , where β = 1/(k B T ) [14]. In presence of an external field H, one considers the Gibbs free-energy where M = i m i . The minima of G define the metastable states and satisfy the stationary conditions δG/δm i = 0, i.e. where (in the following, the coupling constant J is taken as the energy unit and set to unity). The method then consists in following the evolution of the system among these local minima as they are displaced or destroyed when the temperature T or the field H is changed [11]. Since one forbids any barrier crossing, the system is forced to stay in the same miminum until the stability limit is reached (which is signaled by the vanishing of some eigenvalue of the Hessian matrix ∂ 2 F /∂m i ∂m j ). The system then jumps instantaneously to another nearby minimum. At constant temperature, this jump defines an H-driven avalanche. The dynamics has thus the same irreversible and hysteretic character as the T = 0 singlespin-flip dynamics [15]. The only effect of temperature is to change the shape of the free-energy landscape. In particular, there is only one minimum at high temperature and the evolution becomes reversible. In practice, Eqs. (6) and (7) are solved iteratively (i.e. m is the value of m i after n iterations), changing H by small increments and using as initial configuration the m i 's obtained with the preceding value of H. A crucial property of such an iterative scheme is that it only converges to minima of the free-energy (and not to maxima or saddle-points) [16,17]. Thanks to the convexity of the hyperbolic tangent function in Eqs. (6), this "dynamics" has also the same properties as the T = 0 dynamics: it is abelian [18] (i.e. the metastable state reached after an avalanche does not depend upon the order in which the local magnetizations are updated during the avalanche) and it satisfies return-point memory [15]. It is clear that the Gibbs free-energy obtained from Eq. (5) identifies with the exact enthalpy H = U − HM when T → 0 and m i → s i = ±1. The N coupled equations (6) then reduce to the stability conditions given by Eqs.
(3). The iterative dynamics at finite temperature thus reduces to the usual single-spin-flip metastable dynamics when T → 0. In particular, the hysteresis loop goes to the same T = 0 loop, as illustrated in Fig. 1 below and above the critical disorder (in this figure and in the following, the temperature unit is J/k B . The mean-field critical temperature of the pure system is thus T c = 6). On the other hand, since thermally activated processes are completely neglected, the T = 0 critical point where a macroscopic avalanche appears for the first time in the thermodynamic limit is now replaced by a full line σ c (T ) of out-of-equilibrium critical points [19], as shown schematically in Fig. 2 (note that a precise numerical determination of the critical line would require an extensive finite-size scaling study as in Ref. [19]). Moreover, since the theory fully accounts for the disorder-induced fluctuations that play a dominant role in random-field systems, the critical behavior along the line is expected to be non-classical and identical to the one observed at T = 0 [19]. We now consider the "constrained" situation where the magnetization of the system is changed externally and the field becomes an output variable. For simplicity, we assume that the magnetization m = M/N is monotonously increased from −1 to +1 (generalization to more complicated magnetization histories is straightforward). In this case, the system tries to minimize its Helmholtz free-energy F instead of G while satisfying the global constraint i m i = M .
Since there is no intrinsic dynamics for Ising spins, there is no obvious way to define an M -driven dynamics for the RFIM that would be the counterpart of the T = 0 single-spin-flip dynamics used in the H-driven case. In fact, it turns out that the simplest possible choice for a T = 0 dynamics leads to an hysteretic behavior of the system that does not allow for dissipation, which does not seem very realistic [10]. On the other hand, at finite temperature within the local mean-field approach, one is dealing with a standard problem of constrained optimization that can be solved in a natural way by the method of Lagrange mutipliers. One first introduces the function where λ is a Lagrange multiplier that has the meaning of a field coupled to the extensive variable i m i . Minimizing F ({m i }) with the constraint on the magnetization amounts to solve simultaneously the N + 1 coupled equations ∂L/∂m i = 0 and ∂L/∂λ = 0, and In a second step, one has to define an iterative scheme that tells the system how to go from one solution to the "nearest" one as the magnetization is slowly changed.
For that purpose, we rewrite Eqs. (9) as and sum over i so to express λ as a function of m and the local magnetizations, Eqs. (9), (7) and (12) define a set of nonlinear coupled equations that can be solved iteratively, as detailed below [20]. By construction, the configurations {m * i }, solutions of Eqs. (9), are extrema of the Gibbs free-energy for the special value of the field λ * (m) that satisfies the constraint equation, Eq.(10). In other words, the field λ adjusts itself in the course of the iterations so that the system eventually reaches an extremum of G with the correct magnetization. One can imagine that this is indeed the role of the feedback control used in magnetic materials [4]. It is worth noticing that the constraint i m i = M is not satisfied until convergence is reached, which of course is the very principle of the Lagrange method. It is therefore unclear if one can attribute a physical meaning to the intermediate stages of the iteration, although a similar problem may also occur in experiments [21].
In practice, to improve the convergence of the iteration procedure, we introduce a mixing parameter 0 < α ≤ 1 that controls how much of the previous iteration is retained. When changing the magnetization from m to m + δm, we thus search for a fixed point of the map applied to all m i in parallel. Convergence is assumed when max(|m The following general observations can be made: (i) In contrast with the H-driven situation, there is no obvious relation between the stability of the iterative scheme and the properties of the Hessian matrix [17]. However, we strongly believe that convergence to a fixed point of the map (13) only occurs when the final configuration {m * i } is a local minimum of G, i.e. when ∂ 2 F /∂m i ∂m j is positive definite. This was checked in small systems (up to N = 1000) by diagonalizing the Hessian for each solution {m * i (m), λ * (m)}, using a standard diagonalization routine [22]: in every case we found that all the eigenvalues were positive within numerical accuracy. In the following, we shall therefore assume that the configurations visited by the dynamics are always metastable states. This, however, does not mean that these states are the same as those reached with the field driving, as will be discussed below.
ii) Convergence does not occur at low temperatures (depending on σ and L). The iterative map either enters limit cycles that cannot be broken by further decreasing the value of α, or the constraint i m i (λ) = M starts to be violated, whatever the value of ǫ. This latter problem is obviously related to the fact that some m i become equal to ±1 within numerical accuracy. The problem could perhaps be cured by working with transformed variables (see e.g. Ref. [23]) but we have not tested such a modified procedure. Of course, another possibility is that there are no more solutions to Eqs. (9), (7), and (12) at low temperatures. We shall come back to this problem in section V.
iii) The Helmholtz free-energy F ({m i }) may increase during the iteration procedure. This is in contrast with what happens in the H-driven case where it is found empirically that G({m i }) always decreases during an avalanche (recall that H is kept constant during an avalanche). This is perhaps an indication that the transient behavior described by the map (13) does not correspond to any reasonable physical dynamics on the freeenergy landscape. But this may be a too pessimistic statement.
iv) The field λ varies in a complicated way during the iterations, as shown in Fig. 3. Accordingly, the evolution of the local magnetizations with m is not necessarily monotonous and the numerical results may depend on the choice of the increment δm and/or the mixing parameter α. In principle, this could be a serious shortcoming of the iterative scheme implying that the corresponding dynamics is not deterministic (or "adiabatic" in the sense of Ref. [15]), i.e., the final state depends on how the map (12) is applied. Note for instance that the final values of λ shown in the inset of Fig. 3 slightly differ for α = 0.95 and α = 0.15. However, there is good numerical evidence that an adiabatic behavior is recovered when δm is small enough. This is illustrated in Fig. 4 showing a portion of the curve H(m) computed in a single disorder realization of size L = 20 for two values of δm and α (in order to facilitate the comparison with the H-driven magnetization curve, we now use the same notation H for the field λ * (m) and plot the curve H(m) with m on the vertical axis). One can see that the trajectory is independent of α for δm = 10 −3 . Moreover, the same final states are obtained with δm = 10 −3 or 10 −4 . Unfortunately, one cannot compute a whole curve in a reasonable time using such small values of δm and most of the calculations presented in the following were computed with δm = 10 −2 . With this value, the dynamics is not fully adiabatic but the variations when changing α are reasonably small, as illustrated in Fig. 3.

IV. RESULTS
The main results of our study are summarized in Fig. 5 where we compare the m(H) and H(m) hysteresis loops obtained with the H-driven and M -driven dynamics, respectively, for three different values of σ at T = 2(well below the mean-field critical temperature of the pure system). These calculations were performed in samples of linear size L = 60 or L = 80 [24]. One can readily notice two remarkable features: i) the M -driven hysteresis loops display strong fluctuations with m, especially at low σ, and ii) the behavior is quite different in the strong and weak-disorder regimes characterized respectively by the absence or presence of a macroscopic (Hdriven) avalanche. (Although we do not know the exact critical line σ c (T ), the changes in the shape of the m(H) curve in Fig. 5 clearly indicate that 2 < σ c (T = 2) < 3). Before discussing these results in more detail, let us first comment on the mechanism that underlies the behavior of the induced field as one slowly increases the magnetization. At the beginning, when m is close to −1, H has a large negative value and follows m smoothly because there exists only one stable state, the equilibrium one. Hence the two curves H(m) and m(H) are identical. Then, for a certain value of m, there is a discontinuous jump in the field towards a lower value. This corresponds to the first avalanche in the m(H) curve. One can thus associate this jump to the disappearance of the initial minimum: the system has to find another metastable state with the new magnetization. One then again follows this mimimum smoothly, until it reaches its own stability limit. Then, there is a new jump to the left and the evolution proceeds like this until m = +1. A similar behavior is observed along the descending branch. As can be seen in Fig. 4, when increasing m, the excursions of the field are always to the left of the H-driven magnetization curve and there may be several discontinuities in H(m) corresponding to a single avalanche in m(H). Indeed, because of the ferromagnetic nature of the interactions, one can prove (at T = 0 [25], but also at finite temperature in the framework of local mean-field theory) that all the metastable states are inside the Hdriven saturation loop. If this interpretation is correct, the magnitude of the horizontal jumps in H(m) is thus related to the organization and the number of metastable states inside the H-driven hysteresis loop. This explains the two different behaviors observed above and below the critical disorder σ c (T ).

A. Strong-disorder regime
We first consider the strong-disorder regime σ > σ c (T ). As can be seen in Fig. 5(a), the fluctuations in H(m) are small, which may be traced back to the fact that there are many metastable states in the vicinity of the H-driven magnetization curve. More precisely, one expects the number of metastable states inside the hysteresis loop to scale exponentially with N when σ > σ c (T ). This was proven analytically in Ref. [25] for the 1-d RFIM at T = 0 and is probably always true in the strong-disorder regime where the saturation loop is continuous. In particular, Fig. 4 in Ref. [25] shows that the complexity, i.e. the logarithm of the number of metastable states divided by N , increases very rapidly in the vicinity of the boundary (whereas it is zero outside the loop, indicating that there are no metastable states in this region when N → ∞). Since, moreover, the H-driven avalanches remain of microscopic size when σ > σ c (T ), we expect the small wiggles in the H(m) curve to become infinitesimally small in the thermodynamic limit. This is indeed what is observed in Fig. 6(a) where the curves obtained for different system sizes are compared. Therefore, we predict that the two curves H(m) and m(H) should become identical in the thermodynamic limit, as suggested by Fig. 5(a). Note that these heuristic and numerical arguments imply that H(m) is a self-averaging quantity for σ > σ c (T ). H(m) has a different value for each realization {h i } of the random fields and the issue of self-averaging concerns the behavior of the width of the probability distribution P L (H(m)) as L increases. If H(m) is self-averaging, "most" realizations (in a sense that can be made precise) lead to the same value of H(m) in the thermodynamic limit. If not, the measurement of the field in a single sample, no matter how large, does not give a meaningful result and must be repeated on many samples. The situation is clearly different from the usual H-driven case where the magnetization is the output variable: H(m) is indeed an intensive quantity and imposing the magnetization M of a large sample does not impose the magnetization in large subsamples. Therefore, the standard Brout argument [26]  the thermodynamic limit, H(m) is a self-averaging quantity, like m(H). This is confirmed by Fig. 7 that shows that the variance V H (m) does vary like L −3 (note that this calculation has only been performed in the initial part of the curve, in the range −1 < m < −0.75, so to reduce the computational effort).

B. Weak-disorder regime
We now consider the case σ < σ c (T ) illustrated in Figs. 5(b) and 5(c). There are two striking differences with Fig.  5(a): i) the curve H(m) shows a pronounced re-entrant behavior, and ii) the field fluctuates much strongly with m. In line with the preceding discussion, we attribute the re-entrant behavior to the peculiar distribution of the metastable states in the field-magnetization plane in the weak-disorder regime. This issue was discussed in Ref. [25] where it was suggested (but not proven) that the macroscopic jump in the m(H) curve is related to the absence of metastable states in a certain region of the field-magnetization plane (see Fig. 10 in that reference). If this scenario is correct, the field, in the M -driven situation, has to jump to the other side of this region so to find metastable states with the correct magnetization. We thus suspect that the H(m) curve is close to (or may be, in the thermodynamic limit, identical to) the boundary of this region. We want to stress that the re-entrant behavior discussed here has nothing to do with the Sshape predicted by the global mean-field theory where all the m i 's are taken equal to the average magnetization m . In Figs. 8(a) and 8(b), the lower part of the curves H(m) are compared to the mean-field prediction where z = 6 is the lattice coordination number and P (h) is the Gaussian probability distribution. As can be seen in Fig. 8(a) for σ = 2, when increasing m, the field H(m) deviates from the mean-field curve much before the spinodal H sp (σ) is reached and it starts to strongly fluctuate with m just at the onset of the macroscopic avalanche in m(H) (note that the intermediate part of the meanfield curve is unstable whereas all the states visited by the M -driven dynamics are metastable). Clearly, an average description cannot properly describe the metastable behavior associated to the existence of local inhomogeneities in the system (it is only when m is close to ±1 and the system is quasi-homogeneous that the local and global descriptions become similar). The same is true for σ = 1 in Figs. 5(c) or 8(b), but the re-entrant behavior is much more pronounced and the field almost drops to zero. Indeed, one expects the metastable states to concentrate in the vicinity of H = 0 as the disorder decreases (see for instance Fig. 1 in Ref. [25]). When σ = 0, the only metastable states that remain in non-zero field are the two extremal ones that correspond to the metastable branches of the mean-field curve for the pure system, solution of m = tanh[β(cm + H)]. Therefore, in the limit σ → 0, as m is increased from m = −1, the field H(m) follows the lower metastable branch up to H sp (σ = 0), then drops to zero and stays to this value until the upper metastable branch is recovered and followed up to m = +1. The intermediate metastable states in zero field then correspond to bubbles of the minority phase that grow progressively as m is increased (controlling m is like controlling the spinodal decomposition of the system). When σ is finite but small, one expects a similar scenario: initially, the system is able to follow the imposed magnetization while remaining quasi-homogeneous; then, for a certain value of m (and a corresponding value of H that is smaller than H sp (σ) as can be seen in Figs. 8(a) and 8(b)), this situation becomes unstable and a nonequilibrium multidomain state is formed while the field decreases suddenly; this multi-domain state then evolves progressively until the system becomes homogeneous again. This scenario may be compared with the one presented in Ref. [4] (see p. 321 in that reference) to explain the re-entrant behavior observed in picture-frame single crystals when the rate of change of the magnetic flux is imposed by using a feedback control [28]. In this case, the motion of a single domain wall is the dominant process. However, one can notice the similarity of the curve plotted in Fig. 5(c) with the experimental curve [28]. A similar behavior is also observed in strain-stress curves measured in shape memory alloys under displacement control loading condition (see e.g. Refs. [7,8]). In this case, the physical origin of the sudden drop in the stress after it has reached a critical value is still debated (in the recent simulations of Ref. [29], it is associated to the formation of a multi-domain state, like in the present work). Our results suggest that this is probably a generic behavior of the hysteresis loops obtained in weakly disordered systems when controlling the extensive variable.
We now move to the discussion of the discontinuous fluctuations of the field H with the magnetization. As shown in Figs. 6(b) and 6(c), the magnitude of these fluctuations decrease with system size but they are still very significant for L = 100, the largest size investigated in the present work (we note however that the strongest discontinuities observed in Fig. 6(c) are displaced to larger and larger values of m as L increases, which suggests that they disappear as L → ∞). One may wonder whether one should observe such a sporadic, discontinuous behavior of the induced field in a macroscopic sample. This is again related to the issue of self-averaging. As shown in Figs. 9(a) and 9(b) (for selected values of m in the lower part of the H(m) curve), we find that the variance V H (m) in the weak-disorder regime still decreases with L, but the finite-size scaling exponent becomes significantly smaller than 3 (ρ ≃ 1.75) as soon as m is in the range of the macroscopic avalanche. If this behavior does not change for larger L, this seems to indicate that H(m) is only weakly self-averaging [30], a rather remarkable behavior for a non-critical system. On the other hand, one may note that the corresponding histograms shown in Fig. 10 are still very well fitted by a Gaussian distribution.
What happens just at the onset of the avalanche for  A similar behavior is observed for the Helmholtz freeenergy F as a function of the magnetization: it is strongly self-averaging at large disorder (with a variance decaying like 1/L 3 ) and only weakly self-averaging when there is a macroscopic H-driven avalanche. In this latter case, however, the exponent ρ seems to be larger (ρ ≃ 2.45) than the one associated to the fluctuations of the field.

V. CONCLUDING REMARKS
To summarize, we have presented a detailed study of the nonequilibrium response of the random field Ising model when one cycles adiabatically the magnetization (M -driven protocol) instead of the magnetic field (Hdriven protocol). The study has been performed in the framework of the local mean-field approach at finite temperature, neglecting all thermally activated processes. We have shown that two regimes can be observed, depending on the disorder strength. In the strong disorder regime, where all H-driven avalanches are of microscopic size, the two protocols yield the same hysteresis loop in the thermodynamic limit. On the other hand, in the weak-disorder regime, where the H-driven magnetization curve is discontinuous, the M -driven hysteresis loop is re-entrant and the output field displays strong fluctuations with the magnetization. There are also considerable sample-to-sample fluctuations that decrease only weakly with system size.
From the results that have been presented above, a question naturally arises: what is the behavior of the system as the temperature is decreased to zero ? As noted in the introduction, temperature was introduced in the model so to work with continuous variables instead of Ising spins and to use the Lagrange constrained optimization method. The free-energy landscape changes with T (which may be very important for the actual physical behavior, as shown by the study of fluid adsorption in porous solids [19]) but the dynamics is essentially a zero-temperature dynamics since all metastable states have an infinite life-time and free-energy barriers cannot be crossed. All the preceding calculations were thus performed rather arbitrarily at the reduced temperature T = 2 (below the mean-field critical temperature of the pure system), as we were mainly interested in the qualitative behavior of the system as a function of disorder strength. The problem is that the limit T → 0 seems to be non-trivial. On the one hand, one would expect a gradual evolution with T , as is the case with the Hdriven protocol (see Fig. 1). As an illustration, we show in Fig. 11 the evolution of the ascending branch of H(m) as T is decreased from 2 to 1. Not surprisingly, the fluctuations with m increase as T decreases but one could imagine that it is possible to decrease the temperature down to zero. This optimistic view is reinforced by the hope that the envelope of the metastable states in the H − M plane (which we think plays an important role in understanding the behavior of the M -driven system) has a well-defined limit at T = 0, as discussed in Ref. [25]. Unfortunately, the iteration procedure does not converge any more at lower T for reasons that are not com-pletely clear. As noted earlier, this may be a problem of numerical accuracy in the solution of the iterative map. This may also be related to the fact that the magnetization M becomes a discrete variable when T → 0 and can only change by an increment of ±2. It may be then very difficult for the algorithm to find a path from one metastable state to another due to the existence of these gaps in the magnetization. Clearly, it would be much more convenient to define a proper dynamics directly at T = 0. However, as discussed elsewhere [10], this does not seem to be an easy task either.
Finally, we would like to remark that some of the features found in the present calculations, in particular the re-entrant behavior and the intermittent fluctuations in the induced field as a function of magnetization, are observed in a variety of systems characterized by an athermal dynamics. For instance, many amorphous materials, including polymers and metallic alloys, exhibit a mechanical response to deformation (e.g. the stress induced by a quasistatic shear strain) that closely resemble the fieldmagnetization curves computed in this work. Similar behavior is also observed in granular materials and foams.
The relation of this behavior to changes in the underlying energy landscape has been discussed in the literature (see e.g. Ref [31]). However, it seems to be the first time that a direct comparison between the different driving mechanisms is performed.