On the critical exponent α of the 5D random-field Ising model

The random-field Ising model (RFIM) is one of the archetypal disordered systems [1–15], extensively studied due to its theoretical interest, as well as its close connection with experiments in condensed-matter physics [15–19]. Its beauty stems from the combination of random fields and the standard Ising model that creates rich and complicated physical phenomena, responsible for a great volume of research over the last 40 years and more. It is well established that the physically relevant dimensions of the RFIM lay between 2  <  D  <  6, where $D_{\rm l} = 2$ and $D_{\rm u} = 6$ are the lower and upper critical dimensions of the model, respectively. For $D \geqslant D_{\rm u}$ one expects the standard mean-field (MF) behavior [1, 8–10, 20, 21], whereas exactly at $D = D_{\rm u}$ the notoriously obscuring logarithmic corrections appear [22].

In the last few years, the development of a powerful panoply of simulation and statistical analysis methods [23] have set the basis for a fresh revision of the problem. 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 [24–26] and the restoration of supersymmetry and dimensional reduction at D  =  5 [27–29] (see also [30–33] for additional evidence in this respect).

In particular, the large-scale numerical simulations of the 5D RFIM reported in [27] have provided high-accuracy estimates for the spectrum of critical exponents and for several universal ratios (see table III in [27]), with one missing element: that of the direct computation of the critical exponent $\alpha$ of the specific heat. Let us point out that the specific heat of the RFIM is of experimental interest [18] and that the value of $\alpha$ has severe implications for the validity of the fundamental scaling relations, and in particular for the Rushbrooke relation, $\alpha+2\beta+\gamma = 2$, that has been the most controversial of all [34–38]. Therefore a strong command on this aspect of the model's critical behavior is necessary. In the current work we fill this gap by performing additional simulations and scaling analysis that allow us to directly compute $\alpha$ for the 5D RFIM and to therefore present a complete picture of the scaling behavior of the specific heat. Our final estimate, $\alpha = 0.12(2)$, agrees well with that of the 3D Ising universality class, $0.110\,087(12)$ [39], and therefore constitutes additional evidence in favor of our recently proposed restoration of dimensional reduction at D  =  5 [27–29].

The RFIM Hamiltonian is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{H} {{\mathcal H}} = - J \sum_{\langle xy\rangle} S_x S_y - \sum_{x} h_x S_x \; , \nonumber \end{align} \tag{ 1 }$$

with the spins $S_x = \pm 1$ occupying the nodes of a hyper-cubic lattice in space dimension D with nearest-neighbor ferromagnetic interactions J and h_x independent random magnetic fields with zero mean and dispersion $\sigma$. Here we consider the Hamiltonian (1) on a D  =  5 hyper-cubic lattice with periodic boundary conditions and energy units J  =  1. Our random fields h_x follow either a Gaussian $({{\mathcal P}}_G)$, or a Poissonian $({{\mathcal P}}_P)$ distribution of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{distribution} {{\mathcal P}}_G(h,\sigma) = {1\over \sqrt{2 \pi \sigma^2}} {\rm e}^{- {h^2 \over 2\sigma^2}};\ {{\mathcal P}}_P(h,\sigma) = {1\over 2 |\sigma| } {\rm e}^{- {|h| \over \sigma}}, \nonumber \end{align} \tag{ 2 }$$

where $-\infty< h < \infty$ and $\sigma$ the disorder-strength control parameter.

As it is well-established, in order to describe the critical behavior of the model one needs two correlation functions, namely the connected and disconnected propagators, $C^{\mathrm{(con)}}_{xy}$ and $C^{\mathrm{(dis)}}_{xy}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:anomalous} C^{\mathrm{(con)}}_{xy} \equiv \frac{\partial\overline{\langle S_x\rangle}}{\partial h_y} \; ; \, C^{\mathrm{(dis)}}_{xy} \equiv \overline{\langle S_x\rangle\langle S_y\rangle} , \nonumber \end{align} \tag{ 3 }$$

where the $\langle \ldots \rangle$ are thermal mean values as computed for a given realization, a sample, of the random fields {h_x}. Over-line refers to the average over the samples. Following the prescription of [23], for each of these two propagators we scrutinize the second-moment correlation lengths, denoted as $\xi^{\mathrm{(con)}}$ and $\xi^{\mathrm{(dis)}}$, respectively.

Our numerical simulations for the 5D RFIM are described in [27]. We therefore outline here the very necessary details. We simulated lattice sizes from $L_{\rm min}=4$ to $L_{\rm max}=28$. For each pair of (L, $\sigma$) values we generated ground states for 10⁷ samples—for the additional simulations at the most accurate determinations of the critical points shown below in figures 4 and 5, 10⁶ samples were generated—exceeding previous relevant studies [22] by a factor of 10³ on average. The calculation of the ground states of the RFIM was based on the well-established mapping [34–38, 40–55] to the maximum-flow problem [56–58]. We used our own C version of the push-relabel algorithm of Tarjan and Goldberg [59], involving some technical modifications proposed by Middelton and collaborators for further efficiency [46, 47]. Suitable generalized fluctuation-dissipation formulas and reweighting extrapolations have facilitated our analysis, as exemplified in [23]. A comparative illustration in favor of the numerical accuracy of our scheme is shown in figure 1 for the universal ratio $\xi^{\rm (con)} / L$ of an L  =  10 Gaussian RFIM and four different simulation sets, as outlined in the panel.

Zoom In Zoom Out Reset image size

The specific heat of the RFIM can be estimated using ground-state calculations in two complementary frameworks, both based on the analysis of singularities of the bond-energy density E_J [60]. This bond-energy density is the first derivative $\partial E/\partial J$ of the ground-state energy with respect to the random-field strength $\sigma$ [34, 35]. The derivative of the sample averaged quantity $\overline{E}_{J}$ with respect to $\sigma$ then gives the second derivative with respect to $\sigma$ 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 $\overline{E}_{J}$, as $\overline{E}_{J}$ is just the integral of C with respect to $\sigma$. Thus, one may estimate $\alpha$ by studying the behavior of $\overline{E}_{J}$ at $\sigma = \sigma_{\rm c}$ [34], via the scaling form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:E_J_scaling} \overline{E}_{J}(L,\sigma_{\rm c}) = E_{J,\infty} + b L^{(\alpha-1)/\nu}(1+b'L^{-\omega}), \nonumber \end{align} \tag{ 4 }$$

where $E_{J,\infty}$, b, and $b'$ are non-universal constants, and $\omega$ is the universal corrections-to-scaling exponent.

Of course, the use of equation (4) for the application of standard finite-size scaling methods requires an a priori knowledge of the exact value of the critical random-field strength $\sigma_{\rm c}$ (see also the analysis below in figures 4 and 5). Although we currently have at hand such high-accuracy estimates of the critical fields for both types of the random-field distributions under study [27], we start our analysis with an alternative to this approach. In particular, we implement a three lattice-size variant of the original quotients method⁸, also known as phenomenological renormalization [61–63] that has been described in detail in [26] and already successfully applied to the D  =  3 [23] and D  =  4 [25] models. The main idea in this perspective, given that $\alpha - 1 < 0$, is the elimination of the non-divergent background term $E_{J,\infty}$ in equation (4) by considering three lattice sizes in the following sequence: $(L_1, L_2, L_3) = (L, 2L, 4L)$ (see figure 2 for an instructive illustration of the three-lattice variant of the quotients method based on the crossings of $\xi^{\rm (con)} / L$). Taking the quotient of the differences at the crossings of the pairs $(L,2L)$ and $(2L,4L)$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:FSS-3L} \hat Q_O=\frac{\left.(\overline{E}_{J,4L}-\overline{E}_{J,2L})\right|_{(\xi_{4L}/\xi_{2L})=2}}{\left(\overline{E}_{J,2L} - \overline{E}_{J,L})\right|_{(\xi_{2L}/\xi_{L})=2}}, \nonumber \end{align} \tag{ 5 }$$

one obtains the following scaling formula for the bond-energy density [26]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:QO_new} \hat Q_{\overline{E}_J}^{\rm (cross)}=2^{(\alpha - 1)/\nu}+\mathcal{O}(L^{-\omega}). \nonumber \end{align} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

Our results for the effective exponent ratio $(\alpha-1)/\nu$ as a function of $L^{-\omega}$—where $\omega$ is set to the 3D Ising value $0.829\,66$ [39]—are shown in figure 3. The dashed line marks the estimate $(\alpha-1)/\nu = -1.412\, 625\,34\ldots$ of the 3D Ising universality class, where we have used the values $\alpha = 0.110\,087(12)$ and $\nu = 0.629\,971(4)$ [39]. A few comments are in order: (i) Clearly, there exist large corrections to scaling for the sequence of smaller sizes $(2,4,8)$ and $(3,6,12)$ that obscure the application of any finite-size scaling approach. (ii) The remaining data points ($(4,8,16)$, $(5,10,20)$, and $(6,12,24)$) do not allow for a safe extrapolation of the ratio $(\alpha-1)/\nu$ to $L\rightarrow \infty$, although the general trend of the data appears to be on the right track and, in fact, joint polynomial fits with a shared constant term do approach the value $-1.45(6)$ but with a rather bad fitting quality. (iii) Larger system sizes would be needed to clarify this point, but are unfortunately out of reach with our current resources.

Zoom In Zoom Out Reset image size

Guided by these qualitative results of the phenomenological-renormalization approach, we have performed, at a second stage, additional simulations at the critical points $\sigma_{\rm c}(G) = 6.023\,95$ and $\sigma_{\rm c}(P) = 5.590\,38$ of the Gaussian and Poissonian models, respectively [27]. In figures 4 and 5 we report on the finite-size scaling behavior of the bond-energy density at these critical points for the whole spectrum of system sizes studied, alongside with the resulting estimates for the ratio $(\alpha - 1)/\nu$. In both panels the solid lines are fits of the form (4), where the different colors correspond to different fixed values of $\omega$. Black curves correspond to the value $0.829\,66$ of the 3D Ising universality class [39], whereas red curves to the value 0.66 estimated in [27]. The fitting quality, measured in terms of $\chi^2$/dof, where dof measures the number of degrees of freedom, and the minimum system size, $L_{\rm min}$, used in the fits are as follows: $\chi^2/{\rm dof} = 1.8/3$, $L_{\rm min} = 8$ for the Gaussian model (figure 4) and $\chi^2/{\rm dof} = 6.1/4$, $L_{\rm min} = 6$ for the Poissonian model (figure 5). Note that there was practically no variation in the fitting quality moving from $\omega = 0.829\,66$ down to 0.66⁹. Using now the estimate $\nu = 0.629\,971(4)$ for the critical exponent of the correlation length, simple algebra and error propagation produces values for $\alpha$ within the range 0.10–0.13. Taking an average over the values of $\alpha$ obtained from the black curves with $\omega = 0.829\,66$, we give our final estimate for the critical exponent $\alpha$ to be

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:alpha} \alpha = 0.12(2). \nonumber \end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This is compatible to the value $0.12(5)$ obtained in [27] via the modified hyperscaling relation $\newcommand{\e}{{\rm e}} \alpha = 2 - \nu(D - 2 + \bar\eta - \eta)$, where $\newcommand{\e}{{\rm e}} \eta$ and $\newcommand{\e}{{\rm e}} \bar\eta$ are the corresponding anomalous dimensions of the connected and disconnected correlation functions (see equation (3)) and also agrees nicely with the 3D Ising universality benchmark $\alpha = 0.110\,087(12)$ [39].

As an additional consistency check of our results shown in figures 4 and 5, we depict in figure 6 the scaling behavior of the specific heat C, obtained from the derivative of the bond-energy density with respect to the random-field strength $\sigma$, at the critical point. Note that the horizontal axis has been rescaled to $L^{\alpha/\nu}$ (remember that as in the standard case $C \sim L^{\alpha/\nu}$), and $\alpha/\nu$ has been set to the value $0.174\,749\ldots$ via $\alpha = 0.110\,087$ and $\nu = 0.629\,971$ of the 3D Ising universality class [39]. As expected the data become rather noisy with increasing system size, forcing us to exclude from our fittings the larger system sizes L  =  20 and L  =  24, where statistical errors are larger than 30%. Although we illustrate for the benefit of the reader data for the complete spectrum of system sizes studied, the solid lines are simple linear fits within the range L  =  4–16 with a very good fitting quality indeed: $\chi^2/{\rm dof} = 4.16/6$ and 2.03/6 for the Gaussian and Poissonian models, respectively.

Zoom In Zoom Out Reset image size

To summarize, using extensive numerical simulations at zero temperature we provided a high-precision estimate of the specific-heat's critical exponent of the 5D RFIM. Our final result $\alpha=0.12(2)$ is fully consistent with the estimation coming from the modified hyperscaling relation given in [27], and also supports the recent results of [29] for the restoration of supersymmetry and dimensional reduction in the RFIM at D  =  5. We close this contribution with figure 7 and an overview of the critical exponent $\alpha$ of the RFIM at all physically relevant dimensions. Two sets of data points are shown, as outlined in the caption, corroborated by a graphical validation of the Rushbrooke relation in the corresponding inset. Whilst the collative results of figure 7 are reassuring and settle down previous controversies in the random-field problem originating from defective estimations of the critical exponent $\alpha$, for reasons of clarity we sould like to point out that the large error at D  =  3 stems from the joint fits of $[(\alpha - 1)/\nu]^{\rm (ef\,\!f)}$ performed over several random-field distributions (including the double Gaussian distribution) and the large scaling corrections via $\omega(D = 3) = 0.52$ [23, 24]—for further details and graphical explanations on this aspect we refer the interested reader to figures 6 and 7 of [23].

Zoom In Zoom Out Reset image size

We acknowledge partial financial support from Ministerio de Economía, Industria y Competitividad (MINECO, Spain) through Grants No. FIS2015-65078-C2 and PGC2018-094684-B-C21, and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant No. 694925).