Finite-size scaling of the random-field Ising model above the upper critical dimension

Finite-size scaling above the upper critical dimension is a long-standing puzzle in the field of Statistical Physics. Even for pure systems various scaling theories have been suggested, partially corroborated by numerical simulations. In the present manuscript we address this problem in the even more complicated case of disordered systems. In particular, we investigate the scaling behavior of the random-field Ising model at dimension $D = 7$, i.e., above its upper critical dimension $D_{\rm u} = 6$, by employing extensive ground-state numerical simulations. Our results confirm the hypothesis that at dimensions $D>D_{\rm u}$, linear length scale $L$ should be replaced in finite-size scaling expressions by the effective scale $L_{\rm eff} = L^{D / D_{\rm u}}$. Via a fitted version of the quotients method that takes this modification, but also subleading scaling corrections into account, we compute the critical point of the transition for Gaussian random fields and provide estimates for the full set of critical exponents. Thus, our analysis indicates that this modified version of finite-size scaling is successful also in the context of the random-field problem.


I. INTRODUCTION
The random-field Ising model (RFIM) represents one of the simplest models of cooperative behavior with quenched disorder [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16].Despite being seemingly simple in terms of definition, the combined presence of random fields and the standard Ising behavior accounts for a vast range of new physical phenomena, many of them remain unresolved even after 50 years of extensive research.Additionally, its direct relevance to two-and three-dimensionsal experimental analogues in condensedmatter physics, such as diluted antiferromagnets in a field, colloid-polymer mixtures, and others [14,[17][18][19] establishes the RFIM as one of the most prominent platform models for the designing and/or deciphering of experiments.Another asset is that a vast majority of non-equilibrium phenomena including critical hysteresis, avalanches, and the Barkhausen noise [20][21][22][23] can be studied through the RFIM.For a review but also a summary of most recent results we refer to Ref. [24].
It is well established that the physically relevant dimensions of the RFIM lay between 2 < D < 6, where D l = 2 and D u = 6 are the lower and upper critical dimensions of the model, respectively [1].Although the critical behavior of the RFIM at these dimensions has been scrutinized by a variety of methods, a consensus has not been reached for decades.Fortunately, over the last few years several ambiguities have been put at ease due to the development of a powerful panoply of simulation and statistical analysis methods, that have set the basis for a fresh revision of the problem [25].In fact, some of the main controversies have been resolved, the most notable being the illustration of critical universality in terms of different random-field distributions [26][27][28] see also Ref. [29] where it was shown that the diluted Ising model in a field belongs also to the same universality class with the RFIM as predicted by the perturbative renormalization group -and the restoration of supersymmetry and dimensional reduction at D = 5 [28,[30][31][32][33].We refer the reader to Refs.[34][35][36][37][38][39] for additional evidence supporting this latter respect.Furthermore the large-scale numerical simulations of Refs.[26-28, 30, 33] have provided high-accuracy estimates for the full spectrum of critical exponents, putting at rest previous fears of possible violations of fundamental scaling relations.
On the other hand for D ≥ D u the RFIM is expected to show dimension-independent mean-field behavior [1], with the critical exponents holding the mean-field values of the pure Ising ferromagnet (exactly at D = D u the well-known logarithmic corrections appear [40][41][42]).At this point we should emphasize that although the method of finite-size scaling has been successfully applied to the analysis of results by numerous numerical simulations for spin models at D < D u , the situation becomes more complicated when one considers the system above its D u , as discussed extensively for the 5D Ising model (note that D u = 4 for the pure Ising ferromagnet) [43][44][45][46][47][48][49][50][51][52].
For periodic boundary conditions a possible solution has been proposed.The key point in these studies [43][44][45][46][47][48][49][50][51][52] is that at dimensions D > D u the linear length scale L of the system should be replaced in finite-size scaling expressions by a new effective length scale of the form L eff = L D/Du , an ansatz originally proposed by Kenna and Lang in the framework of the ϕ 4 theory [40].In fact, the ratio D/D u is the so-called ϙ exponent introduced by Kenna and Berche in Ref. [50] and elaborated in several subsequent works -see Ref. [52] where an overview of the renormalization group as a successful framework to understand critical phenomena above the upper critical dimension is provided.The proposed scaling theory not only concerns the Ising model but it is believed to be more general.In particular, the finite-size scaling of percolation above its upper critical dimension has also been successfully analyzed in the same framework [53].However, we should point out that the problem is highly non-trivial as the selection of boundary conditions qualitatively changes the scaling [54], so that the case of free (or other type of) boundary conditions is not yet completely settled [49,52].
What is even more, for disordered systems, and in particular for the RFIM, not much has been achieved in this direction, with the exception of Ref. [42] where a qualitative picture of the transition has been provided at high dimensions.In the context of spin glasses, see Ref. [55].
To this end, we present in the current work an extensive numerical study of the RFIM at D = 7 using exact ground-state simulations and a suitable finite-size scaling method based on phenomenological renormalization that takes into account the effective length scale L eff .We locate the critical point of the transition for Gaussian fields and monitor the size evolution of effective critical exponents.Our final results are compatible up to a very good numerical accuracy with their mean-field expectations.Instrumental in our analysis is the use of a proper value for the corrections-to-scaling exponent ω.In this respect, we provide in Appendix A a detailed derivation of ω for the large-N limit of the O(N ) model, starting from Brézin's analysis [56].We find that the exponent ω corresponding to the O(N ) model plays a crucial role for a safe determination of the critical properties in the 7D RFIM.
The remainder of this manuscript is as follows: In Sec.II the model and methods employed are described shortly and in Sec.III our main results on the scaling aspects of the 7D RFIM are presented.We conclude in Sec.IV by providing a summary and an outlook for future work in this direction.

II. MODEL AND METHODS
The Hamiltonian of the RFIM is with the spins S x = ±1 on a D = 7 hypercubic lattice with periodic boundary conditions and energy units J = 1, and h x independent random magnetic fields with zero mean and variance σ 2 .Given our previous universality confirmations [26][27][28], we have restricted ourselves to Gaussian normal-distributed {h x }.We work directly at zero temperature [57][58][59][60][61] because the relevant fixed point of the model lies there [8][9][10].The system has a ferromagnetic phase at small σ, that, upon increasing the disorder, becomes paramagnetic at the critical point σ c .Obviously, the only relevant spin configurations are ground states, which are non-degenerate for continuous random-field distributions.An instance of random fields {h x } is named a sample and thermal mean values are denoted as ⟨• • • ⟩.The subsequent average over samples is indicated by an overline, (e.g., for the magnetization density m = x S x /L D , we consider both ⟨m⟩ and ⟨m⟩).
The scaling theory of the RFIM entails an analysis of two correlation functions, namely the connected and disconnected propagators C For each of these two propagators we scrutinize the second-moment correlation lengths [62], denoted as ξ (con) and ξ (dis) , respectively.Hereafter, we shall indicate with the superscript "(con)", e.g., ξ (con) , quantities computed from the connected propagator.Similarly, the superscript "(dis)", e.g., ξ (dis) , will refer to the propagator C (dis) .We also compute the corresponding connected susceptibility χ (con) to obtain the anomalous dimension η, as well the dimensionless Binder ratio As it is well-known, the random field is a relevant perturbation at the pure fixed point, and the random-field fixed point is at T = 0 [8][9][10].The main assumption leading to this result is that the scale of variation of the effective free energy in a correlation volume scales as ξ θ , where θ is so-called violation of hyperscaling exponent [10].This is a consequence of the observation that the important competition yielding the phase transition is between the exchange interactions and the random field which implies that the controlling critical fixed point is at zero temperature with the temperature irrelevant.Hence, the critical behavior is the same everywhere along the phase boundary and we can predict it simply by staying at T = 0 and crossing the phase boundary at the critical field point.This is a convenient approach because we can determine the ground states of the system exactly using efficient optimization algorithms [16,26,42,57,[59][60][61][63][64][65][66][67][68][69][70][71][72][73] through an existing mapping of the ground state to the maximum-flow optimization problem [58,74,75].A clear advantage of this approach is the ability to simulate large system sizes and disorder ensembles in rather moderate computational times.The application of maximum-flow algorithms to the RFIM is nowadays well established [72].One of the most efficient network flow algorithms used to solve the RFIM is the push-relabel algorithm of Tarjan and Goldberg [76].In the present study we prepared our own C version of the algorithm that involves a modification proposed by Middleton et al. [59][60][61] that removes the source and sink nodes, reducing memory usage and also clarifying the physical connection [60,61].Further details on the numerical implementation can be found in Ref. [25].
One big advantage of our numerical toolkit is that it allows us from simulations at a given σ to compute σderivatives and extrapolate to neighboring σ values by means of a reweighting method -see Ref. [25] for full mathematical derivations of fluctuation-dissipation and reweighting formulas.In the present work we consider lattice sizes within the range L min = 2 to L max = 10.For each pair of {L, σ} values we compute exact ground states for 10 6 samples (initial exploratory runs were performed using 10 5 samples), outperforming previous studies -For comparison, 5000 samples with L max = 8 were used in Ref. [42].
We follow the quotients method for finite-size scaling [62,77,78], taking into account the modification L → L eff = L 7/6 , as we work above the upper critical dimension with periodic boundary conditions.As mentioned above, Kenna and Berche identify this ratio 7/6 as the new critical exponent ϙ, giving extensive discussions in Refs.[50,52].In practice, we focus on three dimensionless quantities g(σ, L eff ) that, barring correction to scaling, are independent of the system size at the critical point, namely ξ (con) /L eff , ξ (dis) /L eff , and U 4 .Given a dimensionless quantity g, we consider a pair of lattices sizes (L eff , 2L eff ) and determine the crossing σ c,L eff , where g(σ c,L eff , L eff ) = g(σ c,L eff , 2L eff ), see Fig. 1.This allows us to compute three such σ c,L eff , a first for ξ (con) /L eff , another for ξ (dis) /L eff , and a third for U 4 .
Dimensionful quantities O scale with ξ in the thermodynamic limit as ξ x O /ν , where x O is the scaling dimension of O and ν the critical exponent of the correlation length.At finite system sizes we consider the quotient Q cross O,L eff can be evaluated at the crossings of ξ (con) /L eff , ξ (dis) /L eff , and U 4 .Renormalization group tells us that x O , ν, and the leading corrections-to-scaling exponent ω are universal.Instances of dimensionful quantities used in this work are the derivatives of correlation lengths ξ (con) and ξ (dis) [x Dσξ (con) = x Dσξ (dis) = 1 + ν] and the connected susceptibility , etc.Note that as we applied the quotients method at the crossings of ξ (con) /L eff , ξ (dis) /L eff , and U 4 , the data sets of our simulations were tripled for each pair of system sizes used and thus our practice was to use joint fits imposing a common extrapolation to the thermodynamic limit.Finally, the exponent ω is fixed to the value ω = 1/2 throughout the analysis below, see Appendix A.
Finally, some comments on the fitting procedure: We restrict ourselves to data with L ≥ L min and to determine an acceptable L min we employ the standard χ 2 /DOF-test for goodness of fit, where χ 2 is computed using the complete covariance matrix and DOF denotes the number of degrees of freedom.Specifically, we consider a fit as being fair only if 10% < Q < 90%, where Q denotes the probability of finding a χ 2 value which is even larger than the one actually found from our data [79].
Figure 3. Effective critical exponent ν (eff) vs. L −ω eff .Results are shown obtained from the derivatives of ξ (con) , ξ (dis) , and U4 for all data sets at hand.

III. RESULTS
We start the presentation of our results in Fig. 2 where a joint fit of the form provides the estimate σ c = 9.48391 (50) for the critical field, in excellent agreement (but higher numerical accuracy) with the earlier result 9.48(3) of Ref. [42].The coefficients b k with k = 1, 2, 3 are just scaling amplitudes and the quality is quite good (Q ∼ 45%).Figures 3 and 4 document the infinite-limit size extrapolations of the main critical exponents ν and η using also joint fits of the form (3) in linear and quadratic L −ω eff order and with cutoff sizes L min = 2 and 3, respectively.In both cases a fair fit quality is obtained, namely Q ∼ 25% and 18%, respectively.Evidently, the obtained estimates ν = 0.516 (18) and η = 0.014 (23) are compatible to the mean-field (MF) values ν (MF) = 1/2 and η (MF) = 0.
Obtaining the critical exponent α of the specific heat is much more trickier in most cases, and the random-field problem is no exception [16,25,30,33].The specific heat of the RFIM can be computed via ground-state calculations and the bond-energy density E J [80].This is the first derivative ∂E/∂J of the ground-state energy with respect to the random-field strength σ [16,59].The σderivative of the sample averaged quantity E J then gives the second derivative with respect to σ of the total energy and thus the sample-averaged specific heat C. The singularities in C can also be studied by computing the singular part of E J , as E J is just the integral of C over σ.Thus, one may estimate α from E J at σ = σ c [80] via the scaling form where E J,∞ , b, and b ′ are non-universal constants.Since α (MF) = 0 and ν (MF) = 1/2 above the upper critical dimension as already noted above, it is expected that Obviously, the use of Eq. ( 5) for the application of standard finite-size scaling methods requires an a priori knowledge of the exact value of the critical random-field strength σ c [81].Fortunately, we currently have at hand such a high-accuracy estimate of the critical field, see Fig. 2. Thus, we have performed additional simulations exactly at the critical point σ c = 9.48391 for all range of the accessible system sizes using the standard averaging of 10 6 samples.Data for the bond-energy density are shown in the main panel of Fig. 5 as a function of 1/L eff .The solid line is a fair fit (Q ∼ 23%) of the form (5) excluding the smaller system sizes (L min = 5) while fixing the exponents α, ν, and ω to their expected values.
As an additional consistency check we present in the inset of Fig. 5 the scaling behavior of a "specificheatlike" quantity C obtained from the bond-energy density derivative with respect to the random-field strength σ at the critical point σ c = 9.48391 and using again 10 6 samples.For C the following scaling ansatz is expected since α/ν = 0 at the mean-field level.As it is evident from the plot, the data become rather noisy with increasing system size.This is typical of all derivatives obtained from a fluctuation-dissipation formula.There is no bias, but errors are large because the quantities involved in the fluctuation-dissipation formula are not self-averaging themselves -see also the discussion Ref.[33].Therefore we exclude from our fitting attempt the largest system size L = 10 where statistical errors are larger than 30%.
The solid line shows a simple linear fit of the form (6) excluding the smaller sizes (L min = 4) with an acceptable fitting quality (Q ∼ 89%).Finite-size scaling behavior of the bond-energy density at the critical random-field strength σc (main panel) and the "specific-heatlike" quantity C (inset).

IV. SUMMARY
We have presented a finite-size scaling analysis of the 7D random-field Ising model with a Gaussian field distribution and periodic boundary conditions.Indeed, above the upper critical dimension the choice of boundary conditions remains crucial [54].Ground-state simulations in combination with recent advancements in finite-size scaling and reweighting methods for disordered systems [25] allowed us to provide a high-accuracy confirmation of the mean-field behavior of the model.A major point has been the numerical verification for the use of an effective length-scale L eff = L D/Du (where D/D u = ϙ in the formulation of Ref. [52]) in all finite-size scaling relations as has been proposed for the pure Ising ferromagnet [40,[48][49][50][51][52] and also the clarification with respect to the corrections-to-scaling exponent ω in Ising systems above the upper critical dimension.Currently, we are working exactly at D u , where characteristic logarithmic scaling violations have been reported [42] but still await for a detailed confirmation.
Let us start by recalling the basic definitions from the original work by Brézin [56].We consider a ferromagnetic system with an O(N )-symmetric, nearest-neighbor Hamiltonian on a hypercybic lattice of linear size L with periodic boundary conditions.From this point and on we shall be using the dimensionless inverse temperature β = J/T .The model greatly simplifies in the limit N → ∞.In the paramagnetic phase, where Note that the dispersion relation λ(q) depends crucially on our choice of the nearest-neighbor lattice interaction.
In fact, the only feature shared by all local-interaction Hamiltonians is λ(q → 0) = q 2 + O(q 4 i ).As it is well-known, the problem becomes much simpler in the thermodynamic limit (where anyway the choice of boundary conditions becomes inconsequential) where B.Z. stands for the first Brillouin zone and −π < q i < π for i = 1, 2, . . ., D. Note that the integral in Eq. (A4) is convergent for D > 2 even if we plug m 2 ∞ = 0.The problem we shall be dealing here is the precise connection between Eqs. (A3) and (A4) as L grows, for D > 4. The alert reader will note that this connection cannot be smooth because of the singular behavior at m 2 L = 0 and q = 0 (the strong singularity is characteristic of the periodic boundary conditions) The analysis by Brézin [56] puts the above observation in a sound mathematical footing.

The (finite) Poisson summation formula
Let H(q) be a smooth, periodic function H(q) = H(q + 2π).One starts by recalling the (finite) Poisson summation formula If the function H depends on a D-dimensional argument, H(q), and if it is periodic (with period 2π) along every one of the D axes in the q space, then one can use Eq.(A6) in a nested way Let us now use the notation n = (n 1 , n 2 , ..., n D ) and the short hand n to refer to the multi-dimensional series in the r.h.s. of Eq. (A7) [ ′ n will be the series in which the term n = (0, . . ., 0) has been excluded].Hence, the gap equation (A3) can be rewritten as Let us now introduce the notation and analyze the remainder term , m L = y/L .
(A10) On the view of Eq. (A5), one may expect for small m L that Our analysis is based on the above asymptotic estimate (that we shall now derive).However, because we are interested in corrections to scaling, we shall need to extend this analysis by obtaining as well the next-to-leading term in Eq. (A11).
Brézin did the following simplification that is only valid at small q, and which is, fortunately, the regime of interest  .(A30) (Note that for D > 6, corrections of the order of y D−4 /L D−4 become corrections of order y 2 /L 2 ).Now, recalling Eq. (A28), we see that y 2 ∼ L (D−4)/2 .On the other hand, (y/L) D−4 ∼ 1/L [(D−4)D]/4 (that becomes 1/L (D/2) for D > 6).Therefore, in the regime 4 < D < 6 we identify a dominant exponent ω 1 and a subleading one ω 2 , as follows And, of course, one should expect all kind of sub-leading corrections terms, such as L −2ω1 , L −(ω1+ω2) , etc. Relating the result of Eq. (A31) to the random-field problem (where D u = 6 rather than 4) leads to our main result Hence, for the present case of D = 7 we obtain ω = 1/2.

Figure 1 .
Figure1.Connected correlation length in units of the effective system size as a function of the random-field strength σ.Lines join data obtained from reweighting extrapolation.The inset presents typical illustrations with respect to the sample-to-sample fluctuations and the errors induced by the reweighting extrapolation for a size L = 5.N samples denotes the number of disorder realizations and σ (sim) the value of random field at which the simulation was performed.The comparative data for the different values of N samples refer to the simulation value σ (sim) = 9.47716 and serve in favor of our numerical scheme.
Figure 5.Finite-size scaling behavior of the bond-energy density at the critical random-field strength σc (main panel) and the "specific-heatlike" quantity C (inset).

2 L
and the mass term m 2 L is the inversesquared correlation length m 2 L = 1/ξ 2 L .One relates m and β through the gap equation which simply codes the constraint G(r = 0) = 1