Anomalous finite size corrections in random field models

There is considerable interest in understanding the behaviour of ferromagnetic systems in random magnetic fields. Further, the most simple model one can conceive, the Random Field Ising Model (RFIM), has a plethora of applications to condensed matter problems [1–3], while a complete understanding of its static and dynamical properties is still lacking. Only in recent years have some long-standing issues such as the presence of a glassy phase [4, 5] and the universality of the zero-temperature fixed point [6] been partially settled. Many open problems remain; most noticeably the causes of the breakdown of dimensional reduction are still being actively investigated [7, 8].

Here we contribute to the literature on such random systems by discussing the convergence properties of the average free energy to its thermodynamic limit. Finite size corrections to the free energy have been investigated by the authors in two papers [9, 10] focusing on the paramagnetic phase of disorder Ising models on diluted graphs, in an effort to characterize the free energy correction induced by loops to the Bethe mean field approximation.

Perturbative expansions around mean field approximations are not yet numerically or analytically manageable, in spite of the large amount of work done on them [11–13]. It was suggested in [9, 10] that the study of finite size corrections for infinite range systems has many points in common with the computation of loop corrections in finite dimensional problems. In order to make progress the case of finite size corrections to the average free energy was carefully studied in the case of Erdös-Rényi and random regular graphs [9, 10]. The corrections were computed in the high temperature phase and they were divergent at the transition point. We show in this paper that the O (1/N) corrections are overshadowed in the whole low temperature phase, because $O(\sqrt{1/N})$ corrections arise. The computation of these corrections is the goal of this paper. We present both numerical simulations and analytical computations to support this claim.

These $O(\sqrt{1/N})$ corrections do not arise from loops, but from disorder induced fluctuations in the weight of the two macroscopically different states. A more complex version of this phenomenon should take place in the presence of many different equilibrium states. Some preliminary numerical analysis we conducted suggests that this phenomenon is also present in finite dimensional models. From the technical point of view the analytical computations in the replica framework are connected to the existence of multiples solutions of the mean field equations [14, 15].

In the following section we present the general analytical framework that can be employed to compute the $O(1/\sqrt{N})$ corrections to the average free energy density. Notice that the same corrections are of order $O(\sqrt{N})$ when referred to extensive quantities, and therefore the indistinct use in the text of $O(1/\sqrt{N})$ and $O(\sqrt{N})$ should not confuse the reader. In sections 3 and 4 we apply this formalism to the fully connected and to a diluted Ising model respectively, checking numerically the consistency of the results.

We consider for concreteness a system of N Ising spins, σ_i = ±1, and external random fields, though the following arguments apply to a general class of models possessing O(m) symmetry (once the average over disorder is taken). On a given graph G, the Hamiltonian of the RFIM is

$$\begin{equation} \mathcal{H}=-J \sum \limits_{(i,j)} \sigma_i \sigma_j - \sum \limits^N_{i=1} h_i\, \sigma_i, \end{equation} \tag{ 1 }$$

where the first sum is over adjacent spins, J ⩾ 0 is a ferromagnetic coupling and the the fields h_i are quenched i.i.d. random variables with zero mean and unit variance, i.e. $\mathbb{E}[h_i]=0$ and $\mathbb{E}[h_i\,h_j]=\delta_{ij}$.

It takes a simple argument to show that, at least at zero temperature and for J large enough, the subleading term in N to the average energy E(N) is of order $O(\sqrt{N})$. In fact in this case, for a given realization of the external fields {h_i}, the Gibbs measure is concentrated on the configuration with minimum energy, the candidates being the one with all the spins up and the one with all the spins down. For a given graph with M = O(N) edges, the energies of the two states, let us call them E₊ and are E₋ respectively, are given by

$$\begin{equation} E_{\pm} = - M J\mp \sum\limits^N_{i=1} h_i \qquad {\rm for}\quad J\gg 1. \end{equation} \tag{ 2 }$$

The sum on the right hand side is a random variable of variance N, converging to a Gaussian variable in the thermodynamic limit. Therefore it is easy to see that for the average energy $E(N)=\mathbb{E}[\min(E_+,E_-)]$ we have

$$\begin{eqnarray} &&E(N)= - M J -\sqrt{\frac{2}{\pi}}\sqrt{N}+o\big(\sqrt{N}\big) \quad {\rm for}\quad J\gg 1 . \end{eqnarray} \tag{ 3 }$$

It turns out that the $\sqrt{N}$ subleading behaviour we found in this limit case is present in the whole ferromagnetic region in the J − T (coupling-temperature) plane.

In fact in the ferromagnetic phase the statistical weight is concentrated on two disconnected regions of the configuration space, having positive and negative magnetization respectively, separated by a free energy barrier exponentially increasing in N. While the two regions are completely equivalent in the pure ferromagnetic system, once the disordered external field is turned on their free energies start to differ by a random quantity of order $O(\sqrt{N})$. The total free energy of the system is then simply given by the minimum of the two, except for some exponentially decaying terms. Let us first define, for a given realization of the disorder, the free energies F₊ and F₋ as the ones corresponding to configurations having positive and negative magnetization respectively (for simplicity we assume N to be odd). We assume, and verify a posteriori for the RFIM, that in the ferromagnetic phase the difference between the free energies of the two states, F₊−F₋, is of order $O(\sqrt{N})$. Then the average free energy of the system is given by

$$\begin{eqnarray} F(N) & = -\frac{1}{\beta}\,\mathbb{E}\left[ \log\left({\rm e}^{-\beta F_+}+{\rm e}^{-\beta F_-}\right)\right] \nonumber\\ & = \mathbb{E}\big[\min(F_+,F_-)\big]+ {\rm exp.~vanish.~terms.} \end{eqnarray} \tag{ 4 }$$

From these premises it follows simply that the average free energy has an expansion in N of the form

$$\begin{equation} F(N)= f_0 N +f_1\sqrt{N} +o(\sqrt{N}), \end{equation} \tag{ 5 }$$

where the coefficients f₀ and f₁ of the expansion are defined by

$$\begin{equation} f_0=\lim_{N\to+\infty}\frac{1}{N}\, \mathbb{E}\big[\,F_+\,\big]=\lim_{N\to+\infty}\frac{1}{N}\, \mathbb{E}\big[\,F_-\,\big], \end{equation} \tag{ 6 }$$

$$\begin{equation} f_1=\lim_{N\to+\infty}\frac{-1}{\sqrt{N}}\, \frac{\mathbb{E}\big[\,|F_+ - F_-|\,\big]}{2}. \end{equation} \tag{ 7 }$$

To compute f₁ we extend a method developed in [16] and [17–19]. At fixed N and for a given realization of the disorder, let us call Z₊ and Z₋ the partition functions constrained to the configurations with positive and negative magnetization. Considering n and m replicated systems respectively, one can easily see that for small n and m we have

$$\begin{eqnarray} &&\fl -\frac{1}{\beta}\log\,\mathbb{E}\big[\,Z_+^{n} Z_-^{m}\big]\sim\, n\, \mathbb E [F_+] + m\, \mathbb E [F_-]-\frac{n^2}{2} \beta {\rm Var}(F_+)\nonumber\\ &&-\frac{m^2}{2} \beta{\rm Var}(F_-) -n m\beta\{\rm Covar\}(F_+,F_-). \end{eqnarray} \tag{ 8 }$$

Obviously, for symmetry reasons, $\mathbb E [F_+]=\mathbb E [F_-]$ and Var [F₊] = Var [F₋]. The free energies F₊ and F₋ are jointly distributed random variables. If the limit

$$\begin{eqnarray} {\it\Delta}^2 &\equiv\lim_{N\to\infty}\frac{1}{2 N}\big[{\rm Var}(F_{+}) -{\rm Covar}(F_+,F_-) \big]\nonumber\\ &=\lim_{N\to\infty}\frac{1}{4 N}\mathbb{E}\big[(F_{+} - F_-)^2 \big] \end{eqnarray} \tag{ 9 }$$

is non-zero, the rescaled variables f₊ and f₋, defined by

$$\begin{equation} f_\pm =\frac{F_\pm - f_0 N}{\sqrt{N}}, \end{equation} \tag{ 10 }$$

become non-trivially jointly distributed Gaussian random variables for large N. In fact, while higher orders cumulants in the n and m expansion of equation (8) are of order O(N), they give vanishing contributions to f₊ and f₋.

Since asymptotically $\frac{f_+-f_-}{2}$ is itself a Gaussian random variable of zero mean and variance Δ², it follows that the coefficient f₁ of equations (5) and (7) is given by

$$\begin{equation} f_1 = -\int \frac{{\rm d} z}{\sqrt{2\pi {\it\Delta}^2}}\ |z|\ {\rm e}^{-\frac{z^2}{2{\it\Delta}^2}} =-{\it\Delta}\sqrt{\frac{2}{\pi}}. \end{equation} \tag{ 11 }$$

The computation of the n and m expansions in equation (8) can be performed using the standard replica techniques. Since it is impractical to work with the partial partition functions Z₊ and Z₋, we define the partition function Z_H of the system having an additional deterministic external field H acting on all the spins. Then we define the replicated free energy φ(n, m) as

$$\begin{equation} \phi(n,m) =\lim_{H\to 0^+}\lim_{N\to\infty} -\frac{1}{\beta N}\log\,\mathbb{E}\big[\,Z_H^{n}\, Z_{-H}^{m}\big]. \end{equation} \tag{ 12 }$$

This definition is completely consistent with the left hand side of equation (8) in the ferromagnetic region, and is even more physically sound in the paramagnetic region, allowing us to treat with a unified formalism the whole phase space.

We compute φ(n, m) for integer values of n and m using the saddle point technique for N → ∞. We obtain as usual a variational expression depending on a replicated order parameter and the particular structure of the field in equation (12) implies that we must look for a solution where the first n replicas are positively magnetized while the remaining m replicas are negatively magnetized. Then we make an analytical continuation of the solution to real value of n and m. By definition (12) the 0-th order term of the expansion is zero. Since we are interested only in the quadratic terms of the φ(n, m) expansion, and since φ(n, m) is variational in some order parameter, our calculations will involve only the order parameter computed at n = m = 0 and we do not have to consider its n and m dependence. As always happens with calculations involving replicas, the exchange of the n and N limits will be done without particular care.

It is interesting to note that the $O(\sqrt{N})$ corrections will be obtained evaluating expression (8) at the leading O(N) order, that is by taking the saddle point, while the more common O (1) corrections are typically controlled by Gaussian fluctuations around the saddle point.

In the following sections we will explicitly compute φ(n, m) given by equation (12) in two mean field models.

The first case we consider is that of the RFIM on the fully connected graph. Here the first sum in the Hamiltonian (1) runs over N(N − 1)/2 edges, and the coupling J has to be rescaled by a factor N⁻¹ in order to obtain an extensive thermodynamic behaviour. It is then easy to compute the replicated free energy of the system, defined in equation (12), through a standard Hubbard–Stratonovich transformation given by

$$\begin{eqnarray} &&\fl \phi_{FC}(n,m) = \min_{x\geq 0} \bigg\{\frac{1}{2} J (n+m)\, x^2\nonumber\\&&\ -\frac{1}{\beta}\log \mathbb{E}_h\left[ \Big(2\cosh\beta(J x+h)\big)^n \big(2\cosh\beta(-J x+h)\Big)^m\right]\bigg\}. \end{eqnarray} \tag{ 13 }$$

Notice that the order parameter x, the magnetization of the up-state, enters in the last term of equation (13) both with a plus and a minus sign. As a technical note, in equation (13) the minimization condition does not ever turn to a maximization when the number of replicas is small, as happens in glassy models [20], since the dimension of the order parameter in the full replica space is n+m and is always greater than zero in our calculations.

To obtain the $O(\sqrt{N})$ contribution to the free energy we have to compute the small and n and m expansion of φ_FC(n, m). As already discussed at the end of the previous section, only the n = m = 0 saddle point solution, given by the non-negative solution of

$$\begin{equation} x=\mathbb{E}_h\tanh \beta (J x + h), \end{equation} \tag{ 14 }$$

appears in the second order expansion of the replicated action. Therefore the coefficient Δ, defined in equation (9) and characterizing the $O(\sqrt{N})$ correction though (11), is given by

$$\begin{eqnarray} &&\fl {\it\Delta}_{FC}^2 = \frac{1}{2\beta^2}\big\{\mathbb{E}_h[(\log\cosh \beta(J x +h))^2] -\mathbb{E}_h[\log\cosh \beta(J x +h)\log\cosh \beta(J x - h)]\big\}. \end{eqnarray} \tag{ 15 }$$

The last expression can be easily computed for any value of J and β. To compare the analytical prediction with exact results, it is easier to consider the zero temperature limit of (15). This is given by

$$\begin{equation} \lim_{T\to 0}{\it\Delta}_{FC}^2= \frac{1}{4}\,\mathbb{E}_h\big(|J x +h|-|J x -h|\big)^2 \end{equation} \tag{ 16 }$$

with x non-negative solution of

$$\begin{equation} x=\mathbb{E}_h\,{\rm sign}(J x + h). \end{equation} \tag{ 17 }$$

We focus on the zero temperature limit for two reasons: the zero temperature fixed point is the one controlling the flow of the renormalization group also starting from finite temperature [21]; and the ground state of the RFIM can be computed efficiently using exact numerical algorithms.

We implemented an exact and very efficient algorithm to compute the ground state of a fully connected RFIM that takes advantage of the topological equivalence of all the spins. In fact it is easy to realize that among all the configurations having a certain total magnetization M = ∑_iσ_i, the one with lowest energy is the one where only the first $\frac{N-M}{2}$ spins with the lowest external field are down.

Therefore, once the spins are sorted according to their external fields (an operation of time complexity Θ(NlogN)), we have to look only to these N + 1 configurations characterized by M = −N, −N + 2, ⋅, N to find the ground state (an Θ(N) operation). We have thus produced an algorithm of time complexity Θ(NlogN) and with Θ(N) memory requirements. This is a great improvement over the min-cut algorithm that we use on diluted graphs (as we shall explain in next section), that has time and memory complexity Θ(N³) for fully connected graphs. With this algorithm we were able to perform highly precise averages of systems up to ∼10⁶ spins. We can then subtract the leading order (in N) term f₀N, that can be computed exactly, to obtain the numerical estimate of the coefficient f₁ up to subleading finite size effects. In figure 1 we show the perfect agreement between the results of our exact algorithm and the analytical prediction $f_1 =-{\it\Delta}\sqrt{\frac{2}{\pi}}$, with Δ given by equation (16), for the RFIM at T = 0.

Zoom In Zoom Out Reset image size

The analytical computation of the $O\big(\sqrt{N}\big)$ correction to the average free energy in diluted models is slightly more involved. Here we focus on the random regular graph (RRG) ensemble. Each element of the ensemble is chosen uniformly at random from all the graphs where each node has fixed connectivity z. The RRG and other diluted graph ensembles have the property of being locally tree-like, in the sense that each finite neighborhood of a randomly chosen node is with high probability a tree in the large N limit, and the density of finite loops goes to zero [22]. This property allows for the analytical solution of such models, at least at the leading order in N.

Following the lead of [19], we apply the replica technique to equation (12), to obtain the replicated free energy

$$\begin{eqnarray} &&\phi_{\rm RRG}(n,m) = \min_{\rho} \left\{\frac{z}{2}\, \phi_{\rm edge}[\rho] - (z-1)\, \phi_{\rm site}[\rho]\right\}, \end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray} &&\fl \phi_{\rm edge}[\rho]=\log\left[ \mathbb{E}_{h}\sum_{\sigma_1, \sigma_2} \rho^{z-1}( \sigma_1) \times {\rm e}^{\beta (h_1\sum_{a} \sigma_1^a+J\sum_{a} \sigma^a_1 \sigma^a_2 +h_2\sum_{a} \sigma^a_2)}\,\rho^{z-1}( \sigma_2)\right]. \end{eqnarray} \tag{ 19 }$$

and

$$\begin{equation} \phi_{\rm site}[\rho] = \log\left[\mathbb{E}_h \sum_ \sigma {\rm e}^{\beta h\,\sum_{a=1}^{n+m} \sigma^a}\,\rho^z( \sigma)\right]\!. \end{equation} \tag{ 20 }$$

A similar expression, for the replicated free energy of spin glasses on RRG, was derived in [23] (see also equation (7) of [24] for a more general formulation) and relies on the hypothesis that the graph contains few short loops.

Equation (18) is a straightforward generalization of these results, the difference being that the replicas are divided into two blocks of size n and m respectively. The minimization condition in equation (18) is imposed over all the functions ρ(σ) = ρ(σ¹, ..., σ^n+m) of an n+m replicated Ising spin, taking particular attention to constrain the first n and last m replicas to be in the up and down state respectively. Here, as in the fully connected case, the dimension of the order parameter is always positive, since it is 2^n+m. Therefore we do not incur that common peculiarity of replica calculations, the exchange between minimization and maximization conditions [20].

Under replica symmetry assumptions the order parameter ρ(σ) can be parametrized in the form

$$\begin{equation} \rho( \sigma)=\int \frac{{\rm d} P(u^+,u^-)}{(2\cosh \beta u^+)^n (2\cosh \beta u^-)^m}\, {\rm e}^{\beta(u^+\sum_{a=1}^n \sigma^a +u^-\sum_{a=n+1}^m \sigma^a)}. \end{equation} \tag{ 21 }$$

In the small n and m limit the minimum condition gives

$$\begin{equation} \fl P(u^+,u^-)=\mathbb{E}_h\int \prod^{z-1}_{k=1} {\rm d} P(u_{k}^+,u_{k}^-)\ \delta\bigg(u^+ -g(h+\sum_k u_{k}^+)\bigg)\,\delta\bigg(u^- -g(h+\sum_k u_{k}^-)\bigg), \end{equation} \tag{ 22 }$$

where the function g(x) is the usual cavity iteration rule $g(x)\,{=}\,\frac{1}{\beta}{\rm atanh}(\tanh(\beta J)$ tan h(β x)). The vanishing auxiliary external field H of equation (12) selects the solution of equation (22), supposed to be unique, such that $\overline{u^+}\,{\geq}\,0\,{\geq}\,\overline{u^-}$ (denoting with $\overline{\,\cdot\,}$ the expectation over P(u⁺, u⁻)). The algorithmic equivalent of the distributional equation (22) is nothing else than the standard Belief Propagation algorithm applied two times with opposite initializations on the same sample.

Expanding the replicated free energy φ_RRG(n, m) to the second order in n and m we can than derive the $O(\sqrt{N})$ coefficient of the free energy according to equations (8) and (11). We note that only the fixed point distribution P(u⁺, u⁻) computed at n = m = 0, that is the solution of (22), contributes to the quadratic order of φ(n, m). The coefficient Δ appearing in $f_1=-{\it\Delta}\sqrt{2/\pi}$ then reads

$$\begin{eqnarray} &&\fl {\it\Delta}_{\rm RRG}^2 = \frac{1}{2\beta^2}\left\{\frac{z}{2}\left(\mathbb{E}[A^2_+]-\mathbb{E}[A_+ A_-]\right) -(z-1)\left(\mathbb{E}[B^2_+]-\mathbb{E}[B_+ B_-]\right)\right\}. \end{eqnarray} \tag{ 23 }$$

The terms A_± and B_± stem from the edge and site replicated free energies (φ_site and φ_edge) respectively. The expectations $\mathbb{E}[\,\cdot\,]$ are over both the distribution of the external random field and of the cavity fields. The site terms B₊ and B₋ are defined by

$$\begin{equation} B_{\pm}=\log \frac{2 \cosh \beta (h +\sum_{k=1}^{z} u_{k}^\pm)}{\prod_{k=1}^{z} 2\,\cosh \beta u_{k}^\pm}. \end{equation} \tag{ 24 }$$

Here the fields $u_{k}^+$ and $u_{k}^-$ are distributed according to $P(u_{k}^+,u_{k}^-)$, a solution of equation (22), and h is distributed as the external random fields. The edge terms A₊ and A₋, appearing in equation (23), are defined by

$$\begin{equation} A_{\pm}=\log\sum_{\sigma_1, \sigma_2} \frac{\exp{\beta( h^{\pm}_1 \sigma_1 +J \sigma_1 \sigma_2 + h^{\pm}_2 \sigma_2)}}{\prod_{k=1}^{z-1}4\,\cosh \beta u_{1k}^\pm\,\cosh \beta u_{2k}^\pm}. \end{equation} \tag{ 25 }$$

The random field $h_1^+$ is distributed as $h+\sum_{k=1}^{z-1} u_{1k}^+$, where h is an external random field, and analogous definitions follow for the other cavity fields.

We notice that only in the terms $\mathbb{E}[A_+ A_-]$ and $\mathbb{E}[B_+ B_-]$ of equation (23) is the full joint distribution P(u⁺, u⁻) needed, not just its marginals. The computation of Δ², and therefore of the analytical finite size correction f₁, to a high level of precision through equation (23), is a computationally easy task. We solved numerically the fixed point condition equation (22) through a population dynamic algorithm. In this case each element of the population is a couple of cavity messages, u⁺ and u⁻, each of them encountering the same external random fields h during the iterations of the algorithm. As an initial condition, in each couple the message u⁺ is set to a high positive value, while the message u⁻ is set to a low negative value.

In the paramagnetic phase the stable solution of equation (22) takes the trivial form P(u⁺, u⁻) = P(u⁺)δ(u⁺ −u⁻). The $O(\sqrt{N})$ finite size correction is thus zero. In the ferromagnetic phase instead, the messages u⁺ and u⁻ become non-trivially correlated.

We computed with the population dynamics algorithm the solution of the fixed point equation (22) at temperature T = 0, for many values of the coupling J and for connectivity z = 4. The expectations we find in the expression of Δ² given equation (23) are then computed sampling from the population.

As in the case of the fully connected model, it is easier to verify the analytical predictions working at zero temperature, since the free energy and the energy coincides and the ground state of the system can be obtained through exact polynomial algorithms.

Therefore the analytical result is compared with an exact numerical algorithm that exploits the equivalence between the problem of finding the ground state of the RFIM on an arbitrary graph and the minimum cut optimization problem [25]. We used the implementation of Goldberg–Tarjan's preflow push-relabel algorithm provided by the open source Lemon Graph Library [26], whose worst case complexity is $\Theta (N^{\frac{5}{2}})$ for instances of the RRG ensemble. The numerical estimate of the $O(\sqrt{N})$ coefficient is obtained from a linear combination of the free energy of systems of different sizes, given by

$$\begin{equation} \tilde f_1(N) = \frac{\mathbb{E}[F(2N)-2F(N)]}{c \sqrt{N}} \end{equation} \tag{ 26 }$$

with $c= \frac{1}{\sqrt{2}}-1$. The value of $\tilde f_1(N)$, obtained by averaging the minimum cut results over many samples of the system, should converge for large N to the analytical value computed through the population dynamic algorithm applied to equation (23). The data plotted in figure 2 show a very good agreement between experiments and predictions, although the convergence is slow due to the presence of subleading $O(\frac{1}{\sqrt{N}})$ finite size effects.

Zoom In Zoom Out Reset image size

Moreover the analytic value of f₁, when plotted as a function of the absolute magnetization of the system, $m \equiv \lim_{N\to\infty} N^{-1}\mathbb E[ |\sum_i < \sigma_i> |]$, see figure 3, shows a weak dependence on the connectivity of the system and an almost linear behaviour. A simple estimate of f₁ on both the fully connected graph and RRGs is then given by $f_1 \approx -\sqrt{2/ \pi}\ m$.

Zoom In Zoom Out Reset image size

We are arguing for the existence of an anomalous $O(1/\sqrt{N})$ subleading correction to the thermodynamic average free energy f₀ of systems with a zero-mean external random field. This correction is limited to the ferromagnetic phase, since it is caused by the difference of free energy between the two pure states. If one is interested only in asymptotic quantities such as f₀, the slowing down of the convergence can be avoided in numerical simulations such as Markov Chain Monte Carlo methods by choosing an initial condition uncorrelated to the realization of the disorder (as is indeed usually done), such that the dynamics gets trapped with equal probability in the state with lowest free energy or in the other one. On the other hand estimations of some observables on systems of finite size through exact algorithms, such as the minimum cut algorithm, are bound to follow a behaviour of the type $a(N)\sim a_0 + \frac{a_1}{\sqrt{N}}$ ; therefore the convergence is much slower than the usual $O(\frac{1}{N})$ behaviour one finds in the paramagnetic phase or in pure systems. Constraining the sampling to an observable state symmetrically fluctuating around its mean asymptotic value is in this sense more convenient than sampling from the pure state with greater statistical weight.

Using a formalism inspired by some recent works [16–19] we have presented a general framework for the computation of the coefficient f₁ of the $O(\sqrt{N})$ term in the average free energy. The computation has been carried out, using a variant of the replica trick, in two solvable mean field systems, the RFIM on the fully connected graph and on the RRG ensemble. The analytical results obtained, equations (11), (15) and (23), are found to be in strong agreement with the numerical simulations. The coefficient of the $O(\sqrt{N})$ correction goes to zero at the critical point. By studying the way it vanishes, it could be possible to relate anomalous finite size corrections to critical avalanches, which have been found in [27] to play a fundamental role in the critical behaviour of the RFIM.

A different but equivalent approach to the problem can be taken, based on a peculiar replica symmetry breaking scheme, using the techniques of [14, 15]. Instead of computing the large deviation function φ(n, m), where n and m replicas are constrained to be in the up and down state respectively, and expanding to small values of n and m, the same conclusions can be obtained summing over all the saddle point contributions obtained from the partitioning of the replicas in the two sets. We preferred the approach of [16–19] because it is conceptually more clear and analytically less involving (although they share many similarities).

While the formalism we have developed in section 2 is completely general, the exact computation of the coefficient Δ² (thus of f₁) can be achieved only in mean field models. In finite dimensions one has to resort to a perturbative diagrammatic expansion of a replicated field theory. We did not take this path, but our numerical simulations with a min-cut algorithm, using the same procedure described in section 4 for the RRG ensemble show, qualitatively and also quantitatively, the scenario depicted in figure 2. Summing it up, as a general feature of models with zero-mean random external fields, the average free energy density has a first finite size correction of order $O(\frac{1}{N})$ in the paramagnetic phase, and of order $O(\frac{1}{\sqrt{N}})$ in the ferromagnetic phase.

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement No. 247328 and from the Italian Research Minister through the FIRB project No. RBFR086NN1.