Random-field Ising model: Insight from zero-temperature simulations

We enlighten some critical aspects of the three-dimensional ($d=3$) random-field Ising model from simulations performed at zero temperature. We consider two different, in terms of the field distribution, versions of model, namely a Gaussian random-field Ising model and an equal-weight trimodal random-field Ising model. By implementing a computational approach that maps the ground-state of the system to the maximum-flow optimization problem of a network, we employ the most up-to-date version of the push-relabel algorithm and simulate large ensembles of disorder realizations of both models for a broad range of random-field values and systems sizes $\mathcal{V}=L\times L\times L$, where $L$ denotes linear lattice size and $L_{\rm max}=156$. Using as finite-size measures the sample-to-sample fluctuations of various quantities of physical and technical origin, and the primitive operations of the push-relabel algorithm, we propose, for both types of distributions, estimates of the critical field $h_{\rm c}$ and the critical exponent $\nu$ of the correlation length, the latter clearly suggesting that both models share the same universality class. Additional simulations of the Gaussian random-field Ising model at the best-known value of the critical field provide the magnetic exponent ratio $\beta/\nu$ with high accuracy and clear out the controversial issue of the critical exponent $\alpha$ of the specific heat. Finally, we discuss the infinite-limit size extrapolation of energy- and order-parameter-based noise to signal ratios related to the self-averaging properties of the model, as well as the critical slowing down aspects of the algorithm.


Introduction
The random-field Ising model (RFIM) is one of the archetypal disordered systems [1][2][3], extensively studied due to its theoretical interest, as well as its close connection to experiments in hard [4,5] and soft condensed matter systems [6]. Its beauty is that the mixture of random fields and the standard Ising model creates rich physics and leaves many still unanswered problems. The Hamiltonian describing the model is where σ i = ±1 are Ising spins, J > 0 is the nearest-neighbor's ferromagnetic interaction, and h i are independent quenched random fields.
The existence of an ordered ferromagnetic phase for the RFIM, at low temperature and weak disorder, followed from the seminal discussion of Imry and Ma [1], when the space dimension is greater than two (d > 2) [7][8][9][10][11]. This has provided us with a general qualitative agreement on the sketch of the phase boundary, separating the ordered ferromagnetic phase from the high-temperature paramagnetic one. The phase-diagram line separates the two phases of the model and intersects the randomness axis at the critical value of the disorder strength h c , as shown in figure 1. Such qualitative sketching has been commonly used for the RFIM [12][13][14] and closed form quantitative expressions are also known from the early mean-field calculations [15][16][17]. However, it is generally true that the quantitative aspects of phase diagrams produced by mean-field treatments provide rather poor approximations. The criteria for determining the order of the low-temperature phase transition and its dependence on the form of the field distribution have been discussed throughout the years [15][16][17][18][19]. Although the view that the phase transition of the RFIM is nowadays considered to be of second order [20][21][22][23][24][25], the extremely small value of the exponent β continues to cast some doubts. Moreover, a rather strong debate with regard to the role of disorder, i.e., the dependence, or not, of the critical exponents on the particular choice of the distribution for the random fields and the value of the disorder strength, analogously to the mean-field theory predictions [15][16][17], was only recently put on a different basis [26]. Currently, even the well-known correspondence among the RFIM and its experimental analogue, the diluted antifferomagnet in a field (DAFF), has been severely questioned by extensive simulations performed on both models at positive-and zero-temperature [27]. In any case, the whole issue of the model's critical behavior is still under intense investigation [20][21][22][23][24][25][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].
Already from the work of Houghton et al. [44], the importance of the form of the distribution function in the determination of the critical properties of the RFIM has been emphasized. In fact, different results have been proposed for different field distributions, like the existence of a tricritical point at the strong disorder regime of the system, present only in the bimodal case [15][16][17]44]. Following the results of Houghton et al. [44], Mattis [45] reexamined the RFIM introducing a new type of a trimodal distribution where h defines the disorder (field) strength and p ∈ (0, 1). Clearly, for p = 1 one switches to the pure Ising model, whereas for p = 0 the well-known bimodal distribution is recovered. In general terms, the trimodal distribution (2) permits a physical interpretation as a diluted bimodal distribution, in which a fraction p of the spins are not exposed to the external field. Thus, it mimics the salient feature of the Gaussian distribution for which a significant fraction of the spins are in weak external fields. Mattis suggested that for a particular case, p = 1/3, equation (2) may be considered to a good approximation as the Gaussian distribution [45]. This in turn indicated that the two models should be in the same universality class. Further studies along these lines, using mean-field and renormalization-group approaches, provided contradicting evidence for the critical aspects of the p = 1/3 model and also proposed several approximations of its phase diagram for a range of values of p [46][47][48]. However, none of these predictions has been confirmed by numerical simulations up till now, thus remaining ambiguous, due to the approximate nature of the mean-field-type of the methods used. The scope of the present work is to shed some light towards this direction by examining several critical features of the phase diagram of the RFIM at d = 3, using both distributions described above in equa-

43003-2
tions (2) and (3). In particular, in the first part of our study we provide numerical evidence that clarify the matching between the trimodal (p = 1/3) and Gaussian models and we give estimates for the critical field h c and critical exponent ν that compare very well to the most accurate ones in the corresponding literature of the RFIM. In the second part of our study we concentrate on the most studied case of the Gaussian RFIM, for which we present a scaling analysis of critical data for the order parameter and the specific heat, i.e., data obtained at the best known estimate of the critical field h c . Our analysis points to a very small, but non-zero, value for the magnetic exponent ratio β/ν, and a critical exponent α → 0 − , in good agreement with experimental predictions [49,50]. We also discuss the infinite-limit size extrapolation of energy-and order-parameter-based noise to signal ratios related to the self-averaging properties of the model, as well as some technical aspects of the implemented numerical method. Our attempt benefits from: (i) the existence of robust computational methods of graph theory at zero temperature (T = 0) that allow us to simulate very large system sizes and disorder ensembles, necessary for an accurate investigation of the delicate properties discussed above, (ii) classical finite-size scaling (FSS) techniques, and (iii) a new scaling approach that involves the analysis of the sample-to-sample fluctuations of various well-defined quantities. In particular, sample-to-sample fluctuations and the relative issue of self-averaging have attracted much interest in the study of disordered systems [51]. Although it has been known for many years now that for (spin and regular) glasses there is no self-averaging in the ordered phase [52], for random ferromagnets such a behavior was first observed for the RFIM by Dayan et al. [53] and some years later for the random versions of the Ising and Ashkin-Teller models by Wiseman and Domany [54]. These latter authors suggested a FSS ansatz describing the absence of selfaveraging and the universal fluctuations of random systems near critical points that was refined on a more rigorous basis by Aharony and Harris [55]. Ever since, the subject of breakdown of self-averaging is an important aspect of several theoretical and numerical investigations of disordered spin systems [56][57][58][59][60][61][62][63][64][65][66][67][68]. In fact, Efrat and Schwartz [69] showed that the property of lack of self-averaging may be turned into a useful tool that can provide an independent measure to distinguish the ordered and disordered phases of the system. In view of this observation, we discuss here another useful application of the fluctuation properties of several quantities of the system to obtain information on the ground-state criticality of the RFIM. The rest of the paper is organized as follows: In the next section we describe the general framework behind the mapping of the RFIM to the corresponding network, outline the numerical approach, and provide all the necessary details of our implementation. The relevant FSS analysis that shows the equivalence of both distributions under study in terms of the critical exponent ν of the correlation length, using an approach based on the sample-to-sample fluctuations of the model, is presented in section 3. Then, in section 4 we focus our attention on the most studied case of the Gaussian model and we provide estimates for the magnetic exponent ratio β/ν and the critical exponent α of the specific heat, via the scaling of the order parameter and bond energy, respectively, at the best known estimate of the critical field value. We also investigate the self-averaging properties of the model at criticality, using properly defined noise to signal ratios and we provide an estimate for the exponent z that describes the critical slowing of the algorithm used. Finally, we synopsize our findings in section 5.

Simulation protocol
As already discussed extensively in the literature (see reference [70] and references therein), the RFIM captures essential features of models in statistical physics that are controlled by disorder and have frustration. Such systems show complex energy landscapes due to the presence of large barriers that separate several meta-stable states. When such models are studied using simulations mimicking the local dynamics of physical processes, it takes an extremely long time to encounter the exact ground state. However, there are cases where efficient methods for finding the ground state can be utilized and, fortunately, the RFIM is one such case. These methods escape from the typical direct physical representation of the system, in a way that extra degrees of freedom are introduced and an expanded problem is finally solved. By expanding the configuration space and choosing proper dynamics, the algorithm practically avoids the need of overcoming large barriers that exist in the original physical configuration space. An attractor state in the expended space is found in time polynomial in the size of the system and when the algo-

43003-3
rithm terminates, the relevant auxiliary fields can be projected onto a physical configuration, which is the guaranteed ground state.
The random field is a relevant perturbation at the pure fixed-point, and the random-field fixed-point is at T = 0 [7][8][9][10]. Hence, the critical behavior is the same everywhere along the phase boundary of figure 1, and we can predict it simply by staying at T = 0 and crossing the phase boundary at h = h c . This is a convenient approach, because we can determine the ground states of the system exactly using efficient optimization algorithms [20,21,25,65,66,[71][72][73][74][75][76] through an existing mapping of the ground state to the maximum-flow optimization problem [77]. A clear advantage of this approach is the ability to simulate large system sizes and disorder ensembles in rather moderate computational times. We should underline here that, even the most efficient T > 0 Monte Carlo schemes exhibit extremely slow dynamics in the lowtemperature phase of these systems and are upper bounded by linear sizes of the order of L max 32 [70]. Further advantages of the T = 0 approach are the absence of statistical errors and equilibration problems, which, on the contrary, are the two major drawbacks encountered in the T > 0 simulation of systems with rough free-energy landscapes [5].
A short direct sketching of how this mapping may in principle occur through some simple considerations is as follows: Let G = (V, E ) be a directed, weighted graph consisting of a set V of nodes and a set E of edges, each of the latter connecting two nodes. In a directed graph, for each edge a direction is specified. The property of being weighted means that to each edge from node i to node j a capacity c i j is assigned. Let the number of nodes be n + 2. We enumerate the nodes V = {0, 1, 2, . . . , n, n + 1} and define the first node 0 as source s and the last node n + 1 as the sink t . The remaining nodes will be associated to the lattice sites of the RFIM. We call a directed, weighted graph G with source s, sink t , and capacities c, of the set of nodes V into two disjoint sets S and S (S ∩ S = ∅ and S ∪ S = V ) with s ∈ S and t ∈ S. In other words, one can imagine a cut as a partition that divides the network into two parts, one part belonging to the source and the other to the sink. Generally, there are many different possible cuts in a network. We can assign to each of them a capacity C (S, S), defined as the sum of the capacities of the edges that the cut Note that edges are directed, that is why only edges that start at the source side of the cut and end at the sink side contribute to the capacity of the cut. Now, the central idea that allows us to map the RFIM into a network defined above, consists of describing a cut by a vector X with the property: X i = 1 if i ∈ S and X i = 0 otherwise, i.e., if i ∈ S. Then, by definition, X 0 = 1 and X n+1 = 0. Using this representation, the formula for the cut capacity may be written in the following form An expansion of equation (5) leads to and already a structural similarity to the fundamental Hamiltonian definition of the RFIM [equation (1)] is clearly seen. Further information on this structural similarities, including a detailed algebra, may be found for the interested reader in the relevant literature (see for instance reference [70] and references therein). The application of maximum-flow algorithms to the RFIM is nowadays well established [75]. The most efficient network flow algorithm used to solve the RFIM is the push-relabel (PR) algorithm of Tarjan and Goldberg [78]. For the interested reader, general proofs and theorems on the PR algorithm can be found in standard textbooks [77]. The version of the algorithm implemented in our study involves a modification proposed by Middleton et al. [21,79,80] that removes the source and sink nodes, reducing memory usage and also clarifying the physical connection [79,80].
The algorithm starts by assigning an excess x i to each lattice site i , with x i = h i . Residual capacity variables r i j between neighboring sites are initially set to J . A height variable u i is then assigned to each node via a global update step. In this global update, the value of u i at each site in the set T = j |x j < 0 of negative excess sites is set to zero. Sites with x i 0 have u i set to the length of the shortest path, via edges with positive capacity, from i to T . The ground state is found by successively rearranging the excesses x i , via push operations, and updating the heights, via relabel operations. The order in which sites are considered is given by a queue. In this paper, we consider a first-in-first-out (FIFO) queue. The FIFO structure executes a PR step for the site i at the front of a list. If any neighboring site is made active by the PR step, it is added to the end of the list. If i is still active after the PR step, it is also added to the end of the list. This structure maintains the cycles through the set of active sites.
When no more pushes or relabels are possible, a final global update determines the ground state, so that sites which are path connected by bonds with r i j > 0 to T have σ i = −1, while those which are disconnected from T have σ i = 1. A push operation moves excess from a site i to a lower height neighbor j , if possible, that is, whenever x i > 0, r i j > 0, and u j = u i − 1. In a push, the working variables are modified according to Push operations tend to move the positive excess towards sites in T . When x i > 0 and no push is possible, the site is relabelled, with u i increased to 1+min { j |r i j >0} u j . In addition, if a set of highest sites U become isolated, with u i > u j +1, for all i ∈ U and all j ∉ U , the height u i for all i ∈ U is increased to its maximum value, V , as these sites will always be isolated from the negative excess nodes.
Periodic global updates are often crucial to the practical speed of the algorithm [79,80]. Following the suggestions of references [21,79,80], we have also applied global updates here every V relabels, a practice found to be computationally optimum [25,76,79,80]. Using this scheme we performed large-scale simulations of the RFIM with both type of distributions discussed above in section 1. Let us note here that prior to the commencement of these large-scale simulations, a set of preliminary runs with smaller system sizes revealed the critical h-regime that we should work on. In particular, for the trimodal (p = 1/3) In both cases a disorder-strength step of δh = 0.02 was used. Regarding the disorder averaging procedure, which is of paramount importance in the study of the RFIM, for each pair (L, h) an extensive averaging over N s = 50×10 3 independent random-field realizations has been undertaken, much larger than in previous relevant studies of the model [21,65,66,72,73]. Additionally, for the Gaussian RFIM we performed some further and even more extensive simulations, at the best known estimate of the critical field h c , using in this case an ensemble of N s = 200 × 10 3 random realizations.

Universality aspects
As the outcome of the PR algorithm is the spin configuration of the ground state, we can calculate for a given sample of a lattice with linear size L the magnetization via m = V −1 i σ i . Taking the average over different disorder configurations we may define the order parameter of the system M = [|m|], where [· · · ] denotes disorder averaging. Another physical parameter of interest is the bond energy per spin that corresponds to the first term of the Hamiltonian (1), i.e. e J = −V −1 〈i,j 〉 σ i σ j , and its disorder average, defined hereafter as E J = [e J ]. Our analysis in the sequel will be mainly based on these three thermodynamic quantities, as well as a relevant algorithmic quantity, namely the number of primitive operations of the PR algorithm, that is the number of relabels per spin R.
simultaneous meaning that the values of h c and ν for all data sets in the fitting procedure are shared during the fit. The quality of the fit is fair enough, with a value of χ 2 /dof of the order of 0.6, where dof refers to the degrees of freedom, and produces the estimates h c = 2.745(7) and ν = 1.37(2) for the critical disorder strength and the correlation length's exponent, in agreement with recent estimates in the literature [76]. We now turn our discussion on the Gaussian RFIM. For this case we show in figure 3 (a) the number of relabels per spin R as a function of the disorder strength for various lattice sizes in the range L = 24−156. Again, we observe that for every lattice size L, R has a maximum at a certain value of h, denoted as before with h * L , that may be considered now as a relevant pseudo-critical disorder strength. Following a similar procedure, we extracted the values of the peak-locations (h * L ) as well as the corresponding error bars, whose shift-behavior is now plotted in panel (b) of figure 3. The straight line is power-law fitting attempt of the same form (7) and the outcome for h c and ν is 2.274(4) and 1.37(1), respectively. The quality of the fit is also in this case good, with a value of χ 2 /dof of the order of 0.4.
A few comments on the scaling analysis are now in order: Having simulated more than five lattice size-points in each case, we also tried to perform the above analysis including higher-order scaling corrections of the form (1 + b ′ L −ω ), where ω is the well-known correction-to-scaling exponent, obtained very recently to be ω = 0.52 (11) for this model [26], using the scaling behavior of universal quantities.
However, no improvement has been observed in the quality of the fit of our data. On the contrary, the

43003-6
Random-field Ising model: Insight from zero-temperature simulations Our suggestion of choosing these newly defined pseudo-critical disorder strengths h * L as a proper measure for performing FSS, closely follows the analogous considerations of Hartmann and Young [72] and Dukovski and Machta [73], also for the Gaussian RFIM. The first authors [72] considered pseudocritical disorder strengths at the values of h at which a specific-heat-like quantity obtained by numerically differentiating the bond energy with respect to h attains its maximum. On the other hand, the authors of reference [73] identified the pseudo-critical points as those in the H − h plane (with H a uniform external field), where three degenerate ground states of the system show the largest discontinuities in the magnetization.
Respectively, Middleton and Fisher [21] using similar reasoning on the Gaussian RFIM, characterized the distribution of the order parameter by the average over samples of the square of the magnetization per spin and the root-mean-square sample-to-sample variations of the square of the magnetization. They identified a similar behavior to that of figures 2 (a) and (b), i.e., with increasing L, the peak magnitude of this quantity moved its location to smaller values of h, defining another relevant pseudo-critical disorder strength. However, in reference [21] the authors were only interested in the scaling behavior of the height of these peaks. The practice followed in the current paper, employing the FSS behavior of the peaks of the sample-to-sample fluctuations of several quantities of physical (M and E J ) and technical (R) origin, was inspired by the intriguing analysis of Efrat and Schwartz [69]. These authors, studying also the d = 3 RFIM, showed that the behavior of the sample-to-sample fluctuations in a disordered system may be turned into a useful tool that can provide an independent measure to distinguish between the ordered and disordered phases of the system. The analysis of figures 2 (a) and (b) above verifies their prediction, and the accuracy in the estimation of relevant phase diagram features, like the critical field h c and the critical exponent ν, turns out to be a clear test in favor of the overall scheme.
Let us make at this point a small comment concerning the errors inherent in these types of approximations. The errors induced in the scheme based on the sample-to-sample fluctuations of figures 2 (a), 2 (b), or the primitive operations of the PR algorithm shown in figure 3 (a), have their origin in the application of some polynomial, or peak-like, function in order to extract the relevant position of the maximum in the h-axis. On the contrary, in similar definitions of pseudo-critical points, such as through the use of some properly defined specific-heat-like quantity at T = 0 [72], one should numerically differentiate the data of the bond energy E J , and then consider a smoothing function to locate the position of the maximum. This scheme is subjected to two successive fitting approximations, thus increasing the errors in the estimation of the pseudo-critical points.
To summarize, in this section, we have investigated the matching between the trimodal, p = 1/3, RFIM and the Gaussian RFIM. Clearly enough, our estimates for the critical exponent ν of both models

43003-7
indicate an equivalence among both distributions within error bars, justifying the original prediction of Mattis [45]. Furthermore, we have suggested the values for the critical field h c which compare very well to the most accurate estimations of the literature. For instance, the best known estimate for the Gaussian RFIM is h c = 2.27205 [26], very close to the value 2.274(4) of the present work. This is also true for the reported values of the correlation-length's exponent, as for the Gaussian RFIM, previous high-accuracy estimates suggest a value of ν = 1.37 [21,26,72]. An interesting aspect of our analysis that led to the above results was the illustration that quantities related to the sample-to-sample fluctuations of several quantities of the system or simply the, originally technical, number of primitive operations of the PR algorithm, constitute a useful alternative to investigate criticality.

Gaussian RFIM
In this last part of our work, we concentrate on the Gaussian distribution, which is the most studied case in the literature of the RFIM, and present further results on important aspects of its critical behavior. As already mentioned above, we have performed additional runs at the best-known value of the critical field, that is the value h c = 2.27205 [26]. Thus, the data and analysis of this section are based on extensive simulations performed at this value of the critical field.
In principal, we are interested in the extraction of an accurate estimate for the magnetic exponent ratio β/ν, whose small value casts some doubts on the order of the transition of the RFIM. The route we follow here is via the scaling of the order parameter M at the critical field. This is shown in figure 4, and the solid line is a power-law fitting of the form M (h=h c ) ∼ L −β/ν . The resulting estimate of the magnetic exponent ratio, given also in the figure, is β/ν = 0.0131 (3), a rather small, but non-zero value, also in agreement with some of the most accurate estimations in the literature [21]. The next part of our FSS analysis concerns the controversial issue of the specific heat of the RFIM. The specific heat of the RFIM can be experimentally measured [49,50] and is, for sure, of great theoretical importance. Yet, it is well known that it is one of the most intricate thermodynamic quantities to deal with in numerical simulations, even when it comes to pure systems. For the RFIM, Monte Carlo methods at T > 0 have been used to estimate the value of its critical exponent α, but were restricted to rather small systems sizes and have also revealed many serious problems, i.e., severe violations of selfaveraging [61,64]. A better picture emerged throughout the years from T = 0 computations, proposing estimates of α ≈ 0. However, even by using the same numerical techniques, but different scaling approaches, some inconsistencies were recorded in the literature. The most prominent was that of reference [72], where a strongly negative value of the critical exponent α was estimated. On the other hand, experiments on random field and diluted antiferromagnetic systems suggest a clear logarithmic divergence of the specific heat [49,50].
In general, one expects that the finite-temperature definition of the specific heat C can be extended to T = 0, with the second derivative of 〈E 〉 with respect to temperature being replaced by the second derivative of the ground-state energy density E gs with respect to the random field h [21,72]. The first derivative ∂E gs /∂J is the bond energy E J , already defined above. The general FSS form assumed is that the singular part of the specific heat C s behaves as C s ∼ L α/νC (h − h c )L 1/ν . Thus, one may estimate α by studying the behavior of E J at h = h c [21]. The computation from the behavior of E J is based on integrating the above scaling equation up to h c , which gives a dependence with c i constants. Alternatively, following the prescription of [72], one may calculate the second derivative by finite differences of E J (h) for values of h near h c and determine α by fitting to the maximum of the peaks in C s , which occur at h * L − h c ≈ L −1/ν . However, as already noted in [21], this latter approach may be more strongly affected by finite-size corrections, since the peaks in C s found by numerical differentiation are somewhat above h c , and furthermore it is computationally more demanding, since one must have the values of E J in a wide and very dense range of h-values.
In the present case, where the critical value h c is known with good accuracy, the first approach seems to be more suitable to follow. The numerical data of the critical bond energy and the relevant scaling analysis are presented in figure 5. The solid line is a power-law fitting of the form (8) and the estimate for the exponent ratio (α− 1)/ν is −0.799 (28), as also given in the figure. Using now our estimate ν = 1.37(1), we calculate the critical exponent α of the specific heat, resulting in an estimate α = −0.095 (37), which is fairly compatible to the experimental scenario of a logarithmic divergence (α = 0) [49,50].
Following the discussion above in section 1, our numerical studies of disordered systems are carried out near their critical points using finite samples; each sample is a particular random realization of the quenched disorder. A measurement of a thermodynamic property, say Z , yields a different value for every sample. In an ensemble of disordered samples of linear size L, the values of Z are distributed according to a probability distribution. The behavior of this distribution is directly related to the issue of self-averaging. In particular, by studying the behavior of the width of this distribution, one may qualitatively address the issue of self-averaging, as has already been stressed by previous authors [54 , 55, 58]. In general, we characterize the distribution by its average [Z ] and also by the relative variance that we employ here to investigate the self-averaging properties of the RFIM.
In particular, we study the behavior of the ratio R Z , in the framework of the two main quantities typically used, the order parameter M and the bond energy E J of the model. In figure 6 we plot the ratio R Z , estimated at the critical field h c , for both quantities, as a function of the inverse linear size. The solid lines are simple linear extrapolations to the infinite-limit size L → ∞. As it is straightforward from the extrapolations, the order-parameter that carries the effect of the disorder -we remind here that the 43003-9 random field couples to the local spins -is a strongly non-self-averaging quantity. On the other hand, as expected, the bond energy restores self-averaging in the thermodynamic limit.
Closing, we present some computational aspects of the implemented PR algorithm and its performance on the study of the Gaussian RFIM. Although its generic implementation has a polynomial time bound, its actual performance depends on the order in which operations are performed and which heuristics are used to maintain auxiliary fields for the algorithm. Even within this polynomial time bound, there is a power-law critical slowing down of the PR algorithm at the zero-temperature transition [21,71]. This critical slowing down is certainly reminiscent of the slowing down seen in local algorithms for statistical mechanics at finite temperature, such as Metropolis, and even for cluster algorithms [81]. In fact, Ogielski was the first to note that the PR algorithm takes more time to find the ground state near the transition in three dimensions from the ferromagnetic to paramagnetic phase [71]. This has already been qualitatively seen in figure 3 (a), where, indeed, the number of primitive operations R of the PR algorithm is maximized in the suitably defined pseudo-critical fields h * L . Assuming the standard scaling of the form R ≈ L z w (h − h c ) −1/ν L , where the dynamic exponent z describes the divergence in the running time at h = h c , and w(x) ∼ x −z at large x and w(x) ∼ |x| −z ln |x| as x → −∞, to be consistent with convergence to constant R for h > h c and R ∼ ln L for small h. Our fitting attempt of this scaling form is plotted in figure 7 and the obtained estimate for the dynamic critical exponent z is 0.487(7).

Summary and outlook
To summarize, we have investigated the ground-state criticality of the d = 3 RFIM with two types of the random-field distribution, namely a uniform trimodal and a Gaussian distribution. In particular, we have estimated for both cases the critical disorder strength h c and the critical exponent ν of the correlation length. These values, compare well enough to the most accurate estimates of the literature, with the values of ν placing the trimodal (p = 1/3) RFIM into the universality class of the Gaussian model, thus verifying a scenario suggested many years ago by Mattis [45]. Technically, our effort became feasible through the implementation of a modified version of the PR algorithm that enabled us to simulate very large system sizes, up to 156 3 spins, and disorder ensembles of the order of up to 200 × 10 3 , for several values of the random-field strength.
On physical grounds, we have implemented a FSS approach based on the sample-to-sample fluctuations of various quantities of physical and technical origin and the primitive operations of the PR algorithm. The outcome of this analysis indicated that the fluctuations of the system may be used as an alternative successful approach to criticality, paving the way to even more sophisticated studies of disordered systems under this perspective. Furthermore, we have provided high-accuracy estimates for the controversial issues of the magnetic-exponent ratio of the order parameter β/ν and the critical exponent α of the specific heat. In particular, the magnetic exponent ratio β/ν was found to be very small, but clearly non zero, ruling out the possibility of a first-order phase transition, whereas the exponent α was found to be compatible with zero, in agreement with the experiments [49,50]. Particular interest has been paid to the self-averaging properties of the model, by studying the infinite-limit size extrapolation of energy-and order-parameter-based noise to signal ratios, as well as the critical slowing down aspects of the PR algorithm.
A future challenge emerging out of the current work, is the production of the phase diagram of the trimodal RFIM in the h c −p plane and the verification, or challenge, of the originally proposed for p = 1/3 universal behavior of the trimodal and Gaussian models in higher dimensions, below the upper critical dimension of the RFIM d u . Preliminary simulations for various values of the probability p in the spectrum 0.1 − 0.5 of the trimodal RFIM indicate a smooth scaling behavior of the sample-to-sample fluctuations of the order parameter, as illustrated in figure 8 for a system size of L = 48 and ensembles of the order of N s = 5 × 10 3 realizations. We expect this task and analysis to bring forward new results on the RFIM that will be useful to the community of disordered systems.