Numerical signs for a transition in the 2d Random Field Ising Model at T=0

Intensive numerical studies of exact ground states of the 2-d ferromagnetic random field Ising model at T=0 with gaussian distribution of fields are presented. Standard finite size scaling analysis of the data suggests the existence of a transition at sigma_c= 0.64 +/- 0.08. Results are compared with existing theories and with the study of metastable avalanches in the same model.

The study of systems with quenched disorder has been a challenging problem since the last 15 years. The interplay between thermal fluctuations and disorder has a great influence on the existing phase transitions. Many systems are known to exhibit such phase diagrams highly determined by the degree of disorder (vacancies, impurities, dislocations, etc.) Among other, the most typical examples can be found in magnetism, superconductivity, structural phase transitions, etc. For such systems different models have been proposed. The Ising model with quenched disorder is one of the simplests and it has the advantage that the pure model is well known. The disorder can be of two types: (i) symmetry breaking terms like random-fields or random magnetic impurities, and (ii) non-symmetry breaking like randombonds, vacancies, etc. For all the cases different probability distributions of disorder have been studied. Here we will focus on the study of the Random Field Ising Model (RFIM) in two dimensions (2d) with a gaussian distribution of fields. Since many years ago there has been a discussion concerning the possibility that, for such a 2d model with symmetry breaking random fields, it exists or not order at low temperatures. The initial studies lead to a certain controversy: the Imry-Ma [1] argument suggestes that the lower critical dimension, below which ferromagnetic order is distroyed, is d l ≤ 2 with d = 2 being the limiting case. Renormalization group expansions [2] around d = 6 lead to the "dimensional reduction" argument suggesting that d l = 3, discarding the possibility for ordering in the 2d-RFIM. It has also been suggested [3] that there are new types of order for d > 1. This controversy is probably due to the difficulty in balancing the two ingredients of such models: disorder and thermal fluctuations.
More recently a new approach to disordered systems has been proposed, namely the study of disordered systems at T = 0, i.e. withouth thermal fluctuations. From a theoretical point of view this simplifies the problem withouth making it trivial. Moreover, several experimental systems exhibit phase transitions that can be catalogued into this "athermal" category: two examples are ferromagnetism at low temperature under an external magnetic field [4], and martensitic transformations [5]. Both systems present a first-order phase transition that can be crossed by sweeping a control parameter and are greatly affected by the presence of quenched disorder. We will concentrate on the study of the 2-d RFIM at T = 0 for different values of the standard deviation σ of the gaussian distribution of fields. Our goal is to look for signs for the existence of a phase transition at a certain σ c from a ferromagnetic ordered state for σ < σ c to a disordered state for σ > σ c . For the 3d-RFIM at T = 0, ground state studies [6,7] and renormalization group arguments [8] reveal the existence of such phase transition, but to our knowledge no results for the 2d case have been published. Figure 1 summarizes the finite size scaling study presented in this paper. Data corresponds to estimations of σ cL obtained using different methods, as a function of the linear system size L 1/ν (where 1/ν = 0.5 is the exponent characterizing the correlation length divergence). The standard extrapolation to L → ∞, as will be discussed, renders σ c = 0.64 ± 0.08 different from zero.
We consider the 2d-RFIM on a L × L square lattice with periodic boundary conditions and with the hamiltonian H = − n.n.
i,j S i S j − N i=1 S i h i , where i and j are indices sweeping the full lattice (i, j = 1, . . . , N = L × L), the sum refers to nearest-neighbours (n.n.) pairs, S i = ±1 are spin variables and h i are independent random fields distributed according to the gaussian probability density with h = 0 and h 2 = σ 2 . The advantage of using a continuous distribution is that, for almost any configuration of fields {h i } the ground-sate is not degenerated. The order parameter is the magnetization of the system defined as m = S i /N. Because the ground state is unique, thermal averages are meaningless. Since we are interested in the dependence of the system properties with the amount of disorder σ the only possible averages with physical meaning are the ensemble averages ( · · · (σ)) performed over different realizations of the random fields with a certain fixed degree of disorder σ. Experimentally this has to be understood as averaging measurements on different samples which have been prepared with the same amount of disorder.
The zeroth-order mean field (MF) theory was solved 15 years ago [9]. A solution with the order parameter m (σ) = 0 appears for σ < σ c = 8 √ 2π , and the phase transition is continuous. Of course this MF result cannot be expected to be correct, neither to reflect any dependence on dimensionality. Moreover, the MF studies can be extended to higher orders by exactly treating larger and larger clusters of spins. For thermal phase transitions this is known to extrapolate to the exact value of the critical temperature. The first-order approximation is the Bethe approximation which consists in solving exactly a cluster of a central spin and its four n.n. The method can be extended to larger clusters. We have found a continuous phase transition at: σ c = 2.76, 2.48, 2.12 and 1.98 for clusters of N = 5, 13, 25 and 41 spins respectively. These results are indicated with stars in Fig. 1 . A better approach consists in looking for exact ground states by using the max-flow min-cut theorem [6]. We have designed an algorithm which solves a set of different σ values with a minimization time that grows like L 4 . We have studied lattices with L = 4,8,16,32,64 and 128 and taken averages over 10 5 , 10 4 , 10 4 ,5 10 3 , 10 3 and 30 realizations of random fields respectivelly. The inset in Fig. 1 shows a typical example of a ground state for σ = 1.0 and L = 64. We have focused on the computation of different magnitudes. The order parameter has been estimated from |m| L (σ) and m 2 L (σ). The subscript L indicates that such quantities will, in general, depend on system size. We have also measured the susceptibility as χ L (σ) = N m 2 L − |m| 2 L (which, for large σ, tends to 1 independently of L) and the Binder's cumulant g L (σ) = 1 − m 4 L /(3 m 2 2 L ). Also the correlation length ξ L (σ) can be computed by fitting an exponential decay to the spin-spin correlation function. Figure 2 shows the behaviour of m 2 L as a function of σ for different system sizes. One can estimate σ cL as the inflection point of a fitted third order polinomium. Data can be scaled using the standard finite size scaling assumption:  Fig. 2. A similar scheme can be applied to the study of |m| L (σ) rendering β = −0.026±0.017 and 1/ν = 0.54±0.05. Susceptibility χ L (σ), shown in Fig. 3, exhibits a peak at σ = σ cL which shifts and increases when increasing L. Data can also be scaled using χ L ∼ L αχ L 1/ν (σ − σ cL ) . Power law fits to the heigth and curvature of the peak render α = 1.89 ±0.03 and 1/ν = 0.46 ±0.05. Scaled data is shown in the inset of Fig. 3. Figure 4 shows the behaviour of ξ L (σ). The peak gives an independent measure of σ cL . The continous line is an estimation of ξ(σ) for L → ∞ that will be discused later. Note that the behaviour of the curves is compatible with the finite size scaling hypothesis, i.e., ξ L follows the behaviour corresponding to the infinite system up to a certain ξ max = KL. (Data is compatible with K ≃ 0.08). A last estimation of σ cL can be obtained from g L (σ). For all the studied sizes g L (σ) takes a low value for large σ and reaches the value 2/3 at a certain σ cL . This estimation is independent of L as suggested in ref. [7]. We want to note that β ≃ 0, which means that the order parameter increases very fast to m ≃ 1 after the transition. A low-σ first-order expansion renders 1 − m ∼ 10 −11 for σ = 0.6. This fact may also explain why it is very difficult to measure g L (σ) with numerical accuracy enough to check for a crossing point which is the standard procedure to locate the transition. Note also that β ≃ 0 and α ≃ 2 would suggests that m exhibits a lack of self averaging [10]. The different estimations, as explained in the previous paragraph, of σ cL are plotted in Fig. 1 in front of L −1/ν . The open symbols correspond to the estimation from |m| L (σ) (•), m 2 L (✷), χ L (σ) (✸), ξ L (σ) (△) and g L (σ) (▽). In order to extrapolate the data to L → ∞, we have used the standard expansion for the divergence of ξ, up to second order: . Now, imposing that σ cL is determined by the condition ξ = KL one gets: σ cL = σ c + C 1 L −1/ν + C 2 L −2/ν . Such parabolic fits are also shown in Fig.  1 with continuous lines. The extrapolated σ c lay all within σ c = 0.64 ± 0.08. To get an idea of the error margins, we have also fitted the first order expansion (C 2 = 0) leaving ν free, rendering σ c = 0.65 ± 0.1 and ν = 1.8 ± 0.2; or fixing ν = 2 which renders σ c = 0.56 ± 0.06.
The existence of this phase transition is in apparent contradiction with previous results. It has been proved [11] that the RFIM has a unique Gibbs state in the thermodynamic limit, i.e. for a given configuration of the random fields the ground state is unique. This can be misunderstood [12] as a proof that the ordered phase cannot exist. When considering the ensemble of all possible realizations of the random fields corresponding to a certain value of σ it may well be that the distribution of magnetizations changes from a single peak one (for large σ) to a bimodal one for small values of σ. Thus the phase transition we are proposing should be understood as existing in this ensemble rather than for a single system for which the ground state is unique.It is true that there is an open question here concerning the size of this ensemble: in the thermodynamic limit, is there more than one realization of disorder compatible with a certain σ? We understand that given the discreteness of the Ising lattice and the continuity of the random fields one can still consider the existence of such an ensemble.
Another interesting point is the comparison of our results with studies suggesting an exponential divergence of the correlation length ξ ∼ exp{−B/(σ − σ c ) τ }. On the basis of the study of the interfaces separating regions with m > 0 and m < 0 Binder [13] derived a theory with τ = 2 and σ c = 0. We have tested the validity of this theory by studying the corresponding finite size scaling hypothesis [13,14]: 1/σ 2 cL ∝ ln(L −1 ). The inset in Fig. 4 shows the comparison of this behaviour with the standard one we propose in this work on top of the data corresponding to the estimations of σ cL from the inflection point in m 2 L (both fits have two free parameters). The standard theory works better. We have also tested that, assuming the standard hypothesis (ξ L ∼ σ−σ cL σ cL , the finite size scaling of ξ L is better than using Binder's hypothesis. Moreover, Binder's theory proposes that for large enough systems the configurations with total m ≃ 0 will be more and more frequent. We have not observed the existence of many of such "slab" configurations but found ground-states with closed domains like that in the inset of Fig. 1. Figure 5 shows the probability P (m) obtained from the computations of a very large number of exact ground-states for a system with L = 64. Clearly, for σ < σ c , the configurations with m ≃ 0 have much less probability than the configurations with m ≃ 1. The reason for the failure of Binder's theory could be that, in order to perform the thermodynamic limit he uses a very anisotropic system with open boundary conditions. Our data can be compared with the studies of the evolution of the RFIM at T = 0, obbeying a local relaxation dynamics. It has been found that when sweeping the external field the system evolves by avalanches between metastable states. At a certain degree of disorder σ c the distribution of avalanches becomes critical. In Fig. 1 we show the values of the σ cL corresponding to the 2d case from Ref. [15]. The behaviour is very similar to the equilibirum data. Different extrapolations to L → ∞ have been reported (σ c = 0.75 ± 0.03 [15], σ c = 0.54±0.04 [16]) but all are close to the equilibrium one. Concerning the exponents, for the metastable studies the exponent β has also been found to be very small, while previous reported values for ν are 1.6 ± 0.1 [15] and 5.3 ± 1.4 [16]. Therefore we suggest that the metastability phenomena found in the out-of-equilibrium studies might be associated with a real underlying equilibrium phase transition at σ < σ c for zero external field. It should be mentioned that in the context of these out of equilibrium phase transitions Sethna and collaborators have [16] proposed a theory with exponential divergence of ξ with τ = 1 and σ c = 0.42 ± 0.04. Our data is not consistent with such theory. If we perform a fit in the evolution of σ cL (L) leaving σ c and τ free we get τ = 0.6 and σ c = 0.25. We can obtain still a good fit (and good scalings) by taking τ = 1 and σ c = 0, although we cannot provide any physical explanation for such a behaviour. We finally want to point out that the phase transition we have found at T = 0 may also be related to the change in the type of growth found at σ = 0.33 in the studies of the depinning transition in the same model [17].
In conclusion, we have presented a finite size scaling analysis of numerical data for systems up to L = 128, that suggests that the RFIM with a gaussian distribution of fields, at T = 0 exhibits, a phase transition at σ c = 0.64 ± 0.08. The ensemble average of the magnetization changes from m = 0 for σ > σ c to a state with m = 0 for σ < σ c . The transition is characterized by the exponents 1/ν = 0.5 ± 0.05, β = −0.03 ± 0.02 and α = 1.89 ± 0.03. The possibility of a exponential divergence ξ ∼ exp(B/σ) cannot be excluded.
We acknowledge finantial support from CICyT project number (mat95-0504), CRAY and FCR. C.F. also acknowledges the Comissionat per a Universitats i Recerca (Generalitat de Catalunya) for finantial support.