Diluted 3d-Random Field Ising Model at zero temperature with metastable dynamics

We study the influence of vacancy concentration on the behaviour of the three dimensional Random Field Ising model with metastable dynamics. We focus our analysis on the number of spanning avalanches which allows for a clean determination of the critical line where the hysteresis loops change from continuous to discontinuous. By a detailed finite size scaling analysis we determine the phase diagram and prove numerically the critical exponents along the whole critical line. Finally we discuss the origin of the curvature of the critical line at high vacancy concentration.


I. INTRODUCTION
Externally driven systems at low enough temperature often display rate independent hysteresis. This out-ofequilibrium phenomenon occurs because intrinsic disorder creates multiple energy barriers that the system cannot overcome due to the very weak thermal fluctuations.
The study of zero temperature models with metastable dynamics has been very succesful for understanding rate independent hysteresis. A prototype case is the Random Field Ising Model (RFIM) with single spin-flip relaxation dynamics [1,2]. Although the model is formulated in terms of magnetic variables (external field H and magnetization m) it can be applied to the study of many phenomena associated to low temperature first-order phase transitions in disordered systems, e.g. martensitic transformations [3], fluid adsorption in porous solids [4], ferroelectrics [5], etc.
Disorder is a intriguing concept: in the RFIM it is introduced via independent and quenched random fields on each lattice site, gaussian distributed with zero mean and standard deviation σ. In real materials, disorder is much more complicated and includes features at all length scales: vacancies, interstitials, composition fluctuations, dislocations, strain fields, grain boundaries, sample surfaces, edges and corners, etc. Thus it is interesting to add to the RFIM other sources of disorder in order to see how the non-equilibrium behaviour is modified.
The goal of this paper is to study the diluted RFIM at T = 0 with metastable dynamics and analyze the consequence of introducing a concentration c of quenched vacancies. The interplay between the two kinds of disorder (random fields and vacancies) will be at the origin of the properties of the σ − c phase diagram.
One of the striking results concerning the RFIM with metastable dynamics, as already pointed out in the seminal paper of Sethna et al. [1], is the occurence of a critical point when the amount of disorder σ is increased.
The m vs. H hysteresis loops change from discontinuous (like in a ferromagnet) when σ < σ c to continuous (like in a spin-glass) when σ > σ c . This result was demonstrated using mean-field analysis and numerical simulations in 3d systems. This problem was also studied within the Renormalization Group formalism [6,7]. Moreover, many properties of the critical point have been also studied analytically on Bethe lattices [8,9,10,11,12]. Experimental evidence for the occurence of such a critical point has been found in different magnetic systems [13,14].
Another interesting result of the RFIM with metastable dynamics is that it reproduces the experimental observation that the m(H) trajectories of such athermal systems are discontinuous at small scales. The evolution proceeds by avalanches from a mestastable state to another. In the RFIM the avalanche size distribution becomes a power-law at the critical point. Experimentally, scale free distributions of avalanche properties have been found in many systems [15,16,17,18,19,20,21,22]. A first attempt to study the influence of dilution in such avalanche size distributions was done some years ago [23]. The results of this work, however, should be considered as only qualitative, given the fact that the studied system was two-dimensional [28], the analysis focused only on the avalanche distributions and the results concerning the phase diagram were very approximate.
The order parameter that vanishes at the critical point is the size of the macroscopic discontinuity ∆m. The analysis of this quantity from simulations is very intrincate. In finite-size systems it is very difficult to make the distinction between a macroscopic jump and a microscopic avalanche. The measured order parameter only displays reasonable finite-size scaling (FSS) properties when the simulated systems are very large [24]. Recent studies [25,26], have shown how the critical point can be characterized in systems of moderate size. The key point is to detect the so-called "spanning" avalanches which are the magnetization jumps that involve a set of spins that spans the whole finite system (e.g. cubic lattice) from one face to the opposite one. By this method avalanches in finite systems can be classified as non-spanning, 1Dspanning, 2D-spanning, or 3D-spanning. The average numbers N 1 , N 2 of 1D and 2D-spanning avalanches dis-play a peak at a value of σ that shifts with system size L and tends to σ c when L → ∞. The numerical data can then be scaled according to the FSS hypothesis [25] where α = 1, 2. The exponent ν = 1.2 ± 0.1 characterizes the divergence of the correlation length (ξ ∼ (σ − σ c ) −ν ) whereas θ = 0.10 ± 0.02 characterizes the divergence of the number of critical avalanches. The scaling variable u(σ) is analytic and measures the distance to the critical point. It can be fitted by the second order expression with σ c = 2.21 and A = −0.2. The behaviour of the 3Dspanning avalanches is more complex because they are of two different kinds: (i) critical 3D-spanning avalanches that behave like the 1D and 2D ones and (ii) subcritical 3D-spanning avalanches that will corespond to the ∆m discontinuity in the thermodynamic limit. The analysis is more difficult and requires a double finite-size-scaling technique. This will not be used in the present paper. Instead we will focus only on the behaviour of the average numbers N 1 and N 2 in the presence of vacancies and propose a FSS hypothesis by using a bivariate scaling variable u(σ, c) that allows to study the full σ−c diagram.
In section II, we define the model and the dynamics. In section III, we present results of the numerical simulations. In section IV, we formulate the FSS hypothesis and determine the critical line. In section V, we propose approximations to the bivariate scaling variable u(σ, c). In section VI, we discuss the interplay between vacancies and avalanches and finally, in section VII, we summarize our main findings and conclude.

II. MODEL AND SIMULATIONS
The diluted 3D-RFIM on a cubic lattice with N sites (N = L × L × L) is defined by the following Hamiltonian (magnetic enthalpy): where S i = ±1 are Ising spin variables, c i = 0, 1 indicates the presence of a vacancy (c i = 0) or not (c i = 1) on each site, h i are quenched random fields gaussian distributed with zero mean and standard deviation σ and H is the driving field. The first sum extends over all distinct nearest-neighbour (n.n.) pairs. Vacancies are quenched and randomly distributed over the lattice. Their concentration is measured by The metastable dynamics is implemented as follows: the system is externally driven by the field H which is adiabatically swept from −∞ where the system is fully negatively magnetized (S i = −1) to +∞. (S i = +1). The spins flip according to a local relaxation dynamical rule, where the sum extends over all the n.n. of S i . When a spin flips, it may trigger an avalanche. The unstable spins are flipped synchronously until a new stable situation is reached.
The hysteresis loop is obtained by computing the magnetization as a function of the applied field H. Magnetization avalanches are recorded along the whole increasing field branch and their spanning properties are analyzed by using "mask" vectors (as explained in Ref. 25) that allow to classify them as non-spanning, 1D-spanning, 2Dspanning and 3D-spanning. In this work we shall mainly study the number of spanning avalanches of each kind which are recorded in the full upwards branch. These numbers, N 1 , N 2 and N 3 which depend on L, σ and c corresponds to averages over more than 10 4 realizations with different random fields and random vacancy positions. The disorder averages are denoted by the symbol · . We study sytems of sizes ranging from L = 8 to L = 64 in a number of points on the σ − c diagram, as indicated schematically in Fig. 1.

III. NUMERICAL RESULTS
The general evolution of the average hysteresis loops as a function of σ and c is shown in Fig. 2. One can observe the transition from discontinuous loops to smooth loops when σ or c are increased. It can also be seen that the saturation magnetization decreases with increasing c.    The histograms include all avalanches irrespective of their spanning properties. The qualitative picture is that power law distributions are obtained along a critical line with an exponent that seems to be the same for all values of c. Apparently, no differences can be observed when comparing the transition induced by changing σ from the transition induced by changing c. Below the critical line the distributions show a peak for large values of s which correspond to the 1D, 2D and 3D spanning avalanches. Above, the distributions have an exponentially damped character. Figure 5 shows the average number of 1D, 2D and 3D spanning avalanches as a function of σ for increasing values of the vacancy concentration c ranging from 0 to 0.5. Data corresponds to a system with size L = 16.
The same information is displayed in Fig. 6 for a system with size L = 48.
The behaviour for small and intermediate vacancy concentration is qualitively similar to that found for the nondiluted model [25]. The average numbers N 1 and N 2 display peaks, whereas N 3 shows a peak on the edge of a step function. Note that for L = 16 the peak height in N 1 (σ, c, L) and N 2 (σ, c, L) seem to decrease with increasing c. This behaviour, however, is much less apparent for larger systems (L = 48). Thus, it is possibly due to a finite size effect. At higher concentrations (c > 0.4) N 1 and N 2 begin to develop a flat plateau at low σ. The reason for this plateau can be well understood by looking at the 3d-plot in Fig. 7, which represents the average number N 1 (σ, c, L) for L = 32. The plateau in the constant c cuts of Figs. 5 and 6 is due to the fact that the crest of the N 1 and N 2 functions bends towards the σ = 0 axis.

IV. FINITE SIZE SCALING HYPOTHESIS
The hypothesis that we want to check numerically is that, in the presence of vacancies, the critical point found at c = 0 transforms into a critical line for a wide range of concentrations. Thus the critical exponents found previously should be equally valid for the description of the behaviour of the average numbers N 1 and N 2 with c > 0. According to this hypothesis we shall propose the following corresponding FSS behaviour: where α = 1, 2 and u(σ, c) is a bivariate scaling variable measuring the distance to the critical line. The exponents θ and ν as well as the functionsÑ α were already   found in previous works [25]. Therefore, the hypothesis is quite strong and indicates that all the N 1 and N 2 data corresponding to different sizes L, different vacancy concentrations c and different amounts of disorder σ must collapse into an already known function. The only freedom is in the determination of the bivariate scaling variable u that should be analytic. Before constructing it in the next section, we can make a first test of Eq. 6, by cheking the scaling on the critical line. Note that, setting u = 0, Eq. 6 becomes: where σ and c should be on the critical line.   cuts we have constructed it. The result is shown in Fig.  9. Note that the process can be repeated independently with N 1 and N 2 . The two independent lines overlap almost perfectly.
It is remarkable that the finite size scaling hypothesis allows the collapse of the data up to large values of c, far from the point c = 0 where the scaling function and the exponents where determined. It is also remarkable that scaling works even after the bend that is observed for c > 0.3. (Note that the crossing point shown in Fig.  8(a)  For small values of σ, the critical line displays a vertical behaviour. The critical value of the vacancy concentration c c above which the hysteresis loops do not display a discontinuity can be fitted to c c = 0.426 ± 0.003

V. BIVARIATE SCALING VARIABLE
In general the bivariate scaling variable is a function that can be expanded as: u(σ, c) = a 0 + a 1 σ + a 2 c + a 3 σc + a 4 σ 2 + a 5 c 2 + . . . (8) Since c and σ are not necessarily very small along the critical line, it is difficult to know a priori how many terms in the expansion will be needed in order to find a good scaling collapse. The direct determination of a large number of coefficients from the numerical data is difficult. Therefore we shall take a different strategy taking into account, as much as possible, the previously known data.
As a first step we will concentrate in the region c ≤ 0.3 where the coexistence line shows a linear behaviour and we will try to use an expansion up to quadratic terms only. By forcing that the condition u = 0 is satisfied on the fitted coexistence line, we deduce that u satisfies: We should also consider the fact that the scaling variable is known to be well described by a second order expansion (up to σ 2 ) for c = 0 as indicated in Eq. 2. After some algebra one can determine the two parameters b 0 and b 1 : Therefore we are left with a single free parameter b 2 that should allow to collapse into a single curve all the data in the scaling region, for different values of σ, c and L. We have considered all the available data in region I of Fig.  1. Note also that the same b 2 parameter must be used to scale both N 1 and N 2 data. The best two collapses are shown in Fig. 10 and Fig. 11 for b 2 = −0.13.
Note that the data for c = 0 are also included in this plot. Therefore we are obtaining two scaling functionsÑ 1 and N 2 compatible with those of Ref. 25. In that reference the scaling functions were approximated by Gaussians, although it was also shown that there were systematic deviations. In this work we have tried to fit the data with more complex functions (with three free parameters). We have found a very good χ 2 by using the following modified Lorentzians, which are represented by a continuous line As a second step we try to build u(σ, c) for the data very close to the σ = 0 axis. In this region II (see Fig. 9) the transition line is again quite linear and in fact is almost vertical. This means that to measure the distance to the critical line it should be sufficient to use the variable (c − c c ). We have considered the following second-order expansion: Note that k ′ is not a free parameter. It can be fixed by imposing that the definitions of the scaling variables (9) and (14)  The continuous lines in both figures correspond to the same lines as in Figs. 10 and 11. We can thus conclude that we have built two good approximations (given by Eq. 9 and 14) in regions I and II, to the unique scaling variable which will display a more complex behaviour in the intermediate region were the critical line bends, probably with higher order terms.

VI. DISCUSSION
In this section we try to understand why the critical line exhibits such a curvature. There must be a physical reason that goes beyond the mere effect of the dilution of the system and unstabilizes even more the phase with the ferromagnetic-like discontinuity. We propose that the effect is related to the percolation of vacancies above c p = 0.3116, a value which is, indeed, very close to the limit where the critical line loses its linearity. To justify this hypothesis numerically we have studied, for each particular realization of disorder the distribution of the clusters of vacancies and the position of the avalanches. In particular we have determined the spatial position of the largest vacancy cluster (which above c p will correspond to the percolating cluster in the thermodynamic limit). It is clear that the neigbouring sites of this percolating cluster of vacancies are an easy path for the propagation of an avalanche since these sites have a smaller number of neighbours. To distinghuish such sites we have defined a local flag that takes values b i = 1 when a site belongs to the border of the largest cluster of vacancies or b i = 0 otherwise. We have also recorded the largest avalanche during the H-scan (which will correspond to the spanning avalanche below the critical line in the thermodynamic limit) and marked its position with a flag ǫ i = 1. With these two variables we have defined the correlation between the border of the largest cluster of vacancies and the largest avalanche as: Note that since ǫ i and b i take values 1, 0 only, the power 2 in the first bracket inside the square roots can be suppressed. This correlation is equal to 1 when the spanning avalanche sits exactly on the border of the spanning cluster of vacancies. The behaviour of ρ ǫ,b as a function of c is shown in Fig. 14 for three different values of σ that correspond to the dashed lines indicated in Fig. 9 and for increasing system sizes as indicated by the legend. The important observation is that the curves for σ = 0.5 and σ = 0.1 exhibit two crossing points. One is located at c p and the other on the critical line (it thus shifts with σ). For a concentration of vacancies below c p or above the critical line, the behaviour of the curves with increasing L indicates that the correlation vanishes in the thermodynamic limit, whereas in the region between the two crossing points the correlation increases with increasing system size. A value ρ = 1 is probably not reached since the spanning avalanche is larger than the border of the percolating cluster of vacancies. By this analysis, we have thus identified the origin of the curvature of the critical line: when vacancies percolate the spanning avalanche propagates along the border of the percolating cluster of vacancies. The propagation in such a constrained environment decreases the amount of disorder needed to break the infinite macroscopic avalanche into small microscopic jumps. However, as shown in the previous section this mechanism does not changes the values of the critical exponents.