text:
# The three-dimensional randomly dilute Ising model: Monte Carlo results

## I Introduction.

During the last few decades many theoretical and experimental studies have investigated the critical properties of statistical systems in the presence of quenched disorder. Typical examples are randomly dilute uniaxial antiferromagnets, for instance, FexZn1-xF2 and MnxZn1-xF2, obtained by mixing a uniaxial antiferromagnet with short-range interactions with a nonmagnetic material. Experiments find that, for sufficiently low impurity concentration 1–$`x`$, these systems undergo a second-order phase transition at $`T_c(x)<T_c(x=1)`$, with critical exponents independent of the impurity concentration. The experimental results have been summarized in Ref. , which reports $`\alpha =0.10(2)`$, $`\nu =0.69(1)`$, and $`\beta =0.350(9)`$. These estimates are definitely different from the values of the critical exponents of the pure Ising universality class, where, e.g., $`\alpha =0.1096(5)`$ (Ref. ), and thus indicate that the impurities change the nature of the transition that belongs to a new random universality class. In the presence of an external magnetic field, dilute uniaxial antiferromagnets show a different critical transition, belonging to the universality class of the random-field Ising model .

A simple model for dilute uniaxial systems is provided by the three-dimensional random Ising model (RIM) with Hamiltonian

$$_x=J\underset{<ij>}{}\rho _i\rho _js_is_j,$$
(1)
where the sum is extended over all nearest-neighbor sites, $`s_i`$ are Ising spin variables, and $`\rho _i`$ are uncorrelated quenched random variables, which are equal to one with probability $`x`$ (the spin concentration) and zero with probability $`1x`$ (the impurity concentration). For small $`1x`$, i.e. above the percolation threshold of the spins, this model shows a critical transition analogous to that observed in experiments and whose nature has been the object of many theoretical studies, see, e.g., Refs. .

Numerical Monte Carlo simulations  had long been inconclusive in setting the question of the critical behavior of the RIM. While the measured critical exponents were definitely different from the Ising ones, results apparently depended on the spin concentration, in disagreement with renormalization-group (RG) theory. Only recently has the question been clarified. Ref.  showed the presence of very strong concentration-dependent scaling corrections with exponent $`\omega =0.37(6)`$. Only if they are properly taken into account, the numerical estimates of the critical exponents become dilution independent as expected. Their final estimates are $`\nu =0.6837(53)`$ and $`\eta =0.0374(45)`$, from which one also derives $`\beta =0.3546(28)`$ and $`\alpha =0.051(16)`$ using scaling relations. These results are in good agreement with the experimental ones reported above, although the numerical estimate of $`\alpha `$ is slightly different.

Randomly dilute Ising systems can also be studied by using the field-theoretical (FT) approach . The starting point is the cubic-symmetric Hamiltonian

$$_{\mathrm{LGW}}=d^dx\left\{\frac{1}{2}\underset{i=1}{\overset{N}{}}\left[(_\mu \varphi _i)^2+r\varphi _i^2\right]+\frac{1}{4!}\underset{i,j=1}{\overset{N}{}}\left(u_0+v_0\delta _{ij}\right)\varphi _i^2\varphi _j^2\right\},$$
(2)
where $`\varphi _i`$ is an $`N`$-component field. By using the standard replica trick, it can be shown that in the formal limit $`N0`$ such a model corresponds to a system with quenched disorder effectively coupled to the energy density, as is the case of the RIM . As is well known, the limit $`N0`$ is a subtle one. In the standard perturbative approaches, the limit is taken naively—we simply set $`N=0`$ in the perturbative expansion—implicitly assuming that the replica symmetry is not broken. In recent years, however, this assumption has been questioned  on the ground that the RG approach may not take into account other local minimum configurations of the random Hamiltonian (2), which may cause the spontaneous breaking of the replica symmetry. However, a fixed-dimension perturbative two-loop calculation  in a perturbative approach proposed in Ref.  finds that the standard replica-symmetric fixed point is stable under any replica-symmetry breaking perturbation, thereby supporting the standard approach. In this paper, we do not further consider this issue and in the following we always assume that the standard approach is correct. Note that the good agreement among numerical and field-theoretical results supports this assumption, although one cannot exclude that replica-symmetry breaking effects can only be seen very close to the critical point.

In the FT approach one looks for stable fixed points in the region $`u_0<0`$ (or, equivalently $`u<0`$). If the pure fixed point at $`u=0`$ is stable, disorder is irrelevant, while the presence of a new stable fixed point with $`u<0`$, indicates that disorder is relevant and that dilute systems belong to a new universality class. Numerical and experimental results indicate that in dilute Ising systems the correct scenario is the second one and thus a new random fixed point should be present with $`u^{}<0`$.

The most precise FT estimates of critical quantities are presently obtained by using perturbative methods. However, in the case of random systems the perturbative approach faces new difficulties: the perturbative series are not only divergent, but are also non-Borel summable . This means that even the knowledge of the complete perturbative series does not allow the exact computation of the critical quantities. These difficulties are clearly present in the $`\sqrt{ϵ}`$-expansion and in the related minimal-subtraction scheme without $`ϵ`$-expansion . The expansion in $`\sqrt{ϵ}`$ is not well behaved and does not allow quantitative determinations of the critical exponents, while in the minimal-subtraction scheme results are very sensitive to the resummation method. If the Chisholm-Borel method is used , no random fixed point is found with the longest available series (five loops). Apparently, four-loop series provide the most accurate results and increasing the length of the expansion does not help improving the precision of the results. On the other hand, if a double Padé-Borel resummation is used as proposed in Ref. , a random fixed point is found also at five loops . The estimates of the critical exponents are in any case not very precise, and moreover, at variance with the fixed-dimension approach described below, the stability-matrix eigenvalues turn out to be complex.

The fixed-dimension perturbative expansion in powers of two independent zero-momentum quartic couplings $`u`$ and $`v`$ (directly related to $`u_0`$ and $`v_0`$ defined in Eq. (2)) is apparently better behaved. Up to six loops, a random fixed point $`u^{}`$, $`v^{}`$ is always found, although the estimates of $`u^{}`$ and $`v^{}`$ vary significantly with the order and the resummation method. In spite of that, the estimates of the critical exponents are quite precise, due to a relatively large insensitivity of the results to the position of the fixed point. The analysis of the six-loop series gives  $`\nu =0.678(10)`$, $`\eta =0.030(3)`$, $`\beta =0.349(5)`$, and $`\alpha =0.034(30)`$. The agreement with the experimental and numerical results is again quite satisfactory; only the estimate of $`\alpha `$ seems to be slightly larger than the experimental result.

In this paper we present a new numerical study of the RIM. The purpose is to extend and possibly improve the numerical results of Ref. . We estimate the critical exponents and, in particular, $`\alpha `$ in order to verify if the apparent small discrepancy between experiments and numerical results is really there. Moreover, we determine the four-point and the six-point zero-momentum couplings, and the universal ratio $`R_\xi ^+`$. As a byproduct we are able to check the accuracy of the FT approach by comparing Monte Carlo and FT estimates of the fixed-point values $`u^{}`$ and $`v^{}`$.

We have performed a high-precision Monte Carlo simulation of the model with Hamiltonian (1) at $`J=1`$ and density $`x=0.8`$. Such a value has been chosen on the basis of the results of Ref. , where it was shown that scaling corrections are particularly small for such a value of $`x`$. This is fully confirmed by our analysis: We do not observe scaling corrections with exponent $`\omega =0.37(6)`$, the correction-to-scaling exponent observed in Ref.  for generic values of the density. Note that the absence of corrections with exponent $`\omega `$ also implies the absence of corrections with exponents $`2\omega `$, $`3\omega `$, $`\mathrm{}`$ Therefore, we expect corrections to scaling with next-to-leading exponent $`\omega _2`$ ($`\omega _2=0.8(2)`$ according to field theory ). Unexpectedly, also these corrections are small. The RIM at $`x=0.8`$ is therefore an “improved” model , i.e. a model in which the leading correction to scaling is (approximately) absent in the expansion of any observable near the critical point.

First of all, we determine the critical exponents by using two different methods. A first estimate is obtained by employing the extrapolation method of Refs.  (similar methods have been discussed in Refs. ). It allows us to determine the critical exponents from the high-temperature behavior of the susceptibility and of the correlation length. We also use direct finite-size scaling (FSS) methods, obtaining consistent estimates. Our final results are

$`\nu =0.683(3),`$ (3)
$`\eta =0.035(2),`$ (4)
from which, using scaling and hyperscaling relations, we obtain

$`\gamma =\nu (2\eta )=1.342(6),`$ (5)
$`\beta ={\displaystyle \frac{\nu }{2}}(1+\eta )=0.3535(17),`$ (6)
$`\delta ={\displaystyle \frac{5\eta }{1+\eta }}=4.80(11),`$ (7)
$`\alpha =23\nu =0.049(9).`$ (8)
Our results are in good agreement with those of Ref.  and, in particular, confirm the discrepancy between the experimental and theoretical estimates of $`\alpha `$.

We also carefully check the validity of the hyperscaling relation $`2\alpha =3\nu `$. We analyze the specific heat at the critical point obtaining $`\alpha /\nu =0.115(28)`$. Using the estimate (3) for $`1/\nu `$ we obtain

$$\frac{2}{\nu }\frac{\alpha }{\nu }=3.04(3),$$
(9)
which is fully consistent with 3. We also perform another check of hyperscaling, analyzing the specific heat and the energy at the critical point. We obtain

$$\frac{2}{\nu }\frac{\alpha }{\nu }=2.93(6),$$
(10)
again consistent with 3.

Beside the critical exponents we also measure the four-point couplings $`G_4^{}`$ and $`G_{22}^{}`$ defined in Eqs. (A19) and (A23), which can be directly related to the fixed-point values $`u^{}`$ and $`v^{}`$: $`G_4^{}=v^{}`$ and $`G_{22}^{}=u^{}/3`$. We obtain:

$`G_4^{}=43.3(2),`$ (11)
$`G_{22}^{}=6.2(1).`$ (12)
These estimates differ significantly from those reported in Ref. , which were obtained from the analysis of perturbative six-loop series: $`G_4^{}=38.0(1.5)`$ and $`G_{22}^{}=4.5(6)`$. Clearly, the non-Borel summability of the perturbative expansions gives rise to a large systematic error. It is also possible that the nonanalyticity of the RG functions  near the random fixed point plays an important role.

These discrepancies on the estimates of $`u^{}`$ and $`v^{}`$ call for a reanalysis of the perturbative expansions of the critical exponents. By using the Monte Carlo estimate of $`u^{}`$ and $`v^{}`$ we find $`\nu =0.686(4)`$, $`\eta =0.026(3)`$, and $`\gamma =1.355(8)`$. These estimates do not differ significantly from those obtained in Ref.  and are also in satisfactory agreement with the Monte Carlo results. Clearly exponents are quite insensitive to the exact location of the fixed point. We also try a different method for estimating critical quantities. It is based on an expansion around the Ising fixed point. Results are similar: $`\nu =0.690(8)`$, $`\eta =0.0345(20)`$, and $`\gamma =1.355(10)`$. Note that the estimate of $`\eta `$ is now in perfect agreement with the Monte Carlo result.

In this paper we also determine some other universal amplitude ratios that involve high-temperature quantities. First, we determine the six-point universal ratios $`r_6^{}`$, $`C_{42}^{}`$, and $`C_{222}^{}`$, defined in Eq. (A22). The coefficient $`r_6^{}`$ is particularly important since it parametrizes the small-magnetization expansion of the critical equation of state in the high-temperature phase . We find

$$r_6^{}=0.90(15).$$
(13)
Finally, we compute the universal ratio $`R_\xi ^+`$ defined by

$$R_\xi ^+(\alpha A^+)^{1/3}f^+,$$
(14)
where $`A^+`$ and $`f^+`$ are defined in terms of the singular behavior of the specific heat $`C`$ and of the correlation length $`\xi `$ , $`C_{\mathrm{sing}}A^+t^\alpha `$, $`\xi f^+t^\nu `$ for $`t\beta _c\beta 0^+`$. We obtain

$$R_\xi ^+=0.2885(15),$$
(15)
in good agreement with other theoretical results : $`R_\xi ^+=0.290(10)`$, obtained from the analysis of the corresponding six-loop perturbative series, and $`R_\xi ^+=0.282(3)`$, derived from a quite precise approximation of the equation of state.

The paper is organized as follows. In Sec. II we present the Monte Carlo results. In Sec. II B we determine the critical temperature by performing a careful analysis of the finite-size behavior of some RG invariant ratios near the critical point. In Sec. II C we determine the four-point and six-point couplings by using the extrapolation method of Refs. . In Sec. II D and II E we determine $`\nu `$ and $`\eta `$ by using again the extrapolation method and by also performing a more direct FSS analysis. Then, in Sec. II F we study the finite-size behavior of the energy and of the specific heat near the critical point. We obtain an independent estimate of $`\alpha `$, that allows us to check the validity of the hyperscaling relation $`2\alpha =3\nu `$. Finally, in Sec. II G we compute the universal ratio $`R_\xi ^+`$. For this purpose, we generalize the extrapolation method of Refs.  to the energy. In spite of the necessary subtractions, the method works quite well, providing a rather precise estimate. In Sec. III we reanalyze the six-loop perturbative series of Ref. , using the new Monte Carlo estimate of the fixed point. We employ the different resummation methods discussed in Ref.  and also a new method based on an expansion around the Ising fixed point. Finally, in the appendix we report the definitions of the quantities that are used throughout the paper.

## II Numerical results

### A The Monte Carlo simulation

We have performed a high-precision Monte Carlo simulation of the model with Hamiltonian (1) with $`J=1`$ at density $`x=0.8`$. Such a value has been chosen on the basis of the results of Ref. , who showed that for such a value of $`x`$ scaling corrections are particularly small. In the simulations we have considered cubic lattices of size $`L^3`$, $`L=16`$, 32, 64, 128, and 256, with periodic boundary conditions. Simulations have been performed for several values of $`\beta `$ between 0.275 and $`0.28578`$. Two thirds of the simulations refer to the interval $`0.275\beta 0.2856`$ (we will call the corresponding data the high-temperature results), while one third of the CPU time was used in simulations in a narrow interval around the critical point, $`0.28572\beta 0.28578`$. The average number of samples for each $`\beta `$ and $`L`$ has been approximately $`710^4`$ ($`L=16`$), $`3610^3`$ ($`L=32`$), $`2710^3`$ ($`L=64`$), $`1210^3`$ ($`L=128`$), and $`310^3`$ ($`L=256`$). The runs were performed on a cluster with Dual Athlon MP 1.2MHz processors. The total CPU time is approximately 17.4 CPU-years of a single processor. As random number generator we have used a combination of the Parisi-Rapuano generator  and of a congruential generator . Results for each sample have been obtained as follows. Starting from a random spin configuration, we perform 2000 iterations, each of them consisting alternatively of a Metropolis sweep and of a full Swendsen-Wang update. We use both a local and a nonlocal dynamics to guarantee equilibration of short-distance and long-distance modes. Then, we perform 2000 full Swendsen-Wang updates, measuring all quantities (see appendix for definitions) every 10 iterations. To estimate correlation functions we use improved estimators that significantly reduce the statistical errors. Note that we have been much more conservative than Ref. : There, only 200 Swendsen-Wang iterations were performed for equilibration. For quantities that involve the noise average of products of sample averages, there is a bias due to the finite length of the run for each sample. In order to take this bias into account we have performed a bias correction following Ref. .

### B Determination of the critical temperature

As a first step in our analysis we have determined the critical temperature $`\beta _c`$. For this purpose we consider the results of the simulations for $`0.28572\beta 0.28578`$, which is a small interval around $`\beta _c`$. We consider four invariant ratios, $`U_4`$, $`U_6`$, $`U_{22}`$, and $`\xi /L`$, see appendix for definitions. Standard FSS predicts that, in the FSS limit $`\beta \beta _c`$, $`L\mathrm{}`$ at $`(\beta \beta _c)L^{1/\nu }`$ fixed, each quantity $`R(\beta ,L)`$ behaves as

$$R(\beta ,L)=\widehat{R}[(\beta \beta _c)L^{1/\nu }],$$
(16)
where $`\widehat{R}(z)`$ is a universal function. Since $`\beta \beta _c`$ is particularly small for the data, we can expand $`\widehat{R}(z)`$ in powers of $`z`$, keeping only the first term (we checked that the addition of the term of order $`z^2`$ does not change the results). Thus, we fit each quantity $`R(\beta ,L)`$ by using

$$R(\beta ,L)=R^{}+a(\beta \beta _c)L^{1/\nu },$$
(17)
with $`R^{}`$, $`\beta _c`$, and $`\nu `$ free parameters. In each fit we include all data with $`LL_{\mathrm{min}}`$ and, in order to detect corrections to scaling, we use $`L_{\mathrm{min}}=16`$, 32, and 64. The results are reported in Table I.

There is a slight evidence of corrections to scaling, but it is interesting to note that they have opposite sign in $`U_4`$, $`U_6`$, and $`\xi /L`$. Conservatively, we would obtain

$$\beta _c=0.2857447(24),$$
(18)
where the central value is the average of the estimates obtained for $`U_4`$ and $`\xi /L`$ ($`L_{\mathrm{min}}=64`$), and the error is such to include one error bar for both $`U_4`$ and $`\xi /L`$.

These fits are not particularly sensitive to the exponent $`\nu `$, which is quite poorly determined. We obtain $`\nu =0.72(6)`$. One could imagine of improving the results by fixing $`\nu `$. However, the dependence of $`\beta _c`$ on $`\nu `$ is very small and no significant change is observed.

In order to include scaling corrections we also perform fits of the form

$$R(\beta ,L)=R^{}+a(\beta \beta _c)L^{1/\nu }+bL^\omega ,$$
(19)
where we include the leading scaling correction. These fits are not sensitive to the value of $`\nu `$ and thus we fix it, taking $`\nu =0.69`$. We keep $`\omega `$ as a free parameter, since we do not know which is the most important correction to scaling for our data. Indeed, the leading correction has exponent $`\omega =0.37(6)`$, but there is evidence that for $`p=0.8`$ leading corrections have a very small amplitude . In order to be able to keep $`\omega `$ as a free parameter, we analyzed at the same time two different observables. We restrict our attention to $`U_4`$, $`U_6`$, and $`\xi /L`$, since $`U_{22}`$ is too noisy. Using all data with $`L16`$ (if the data with $`L=16`$ are discarded the fit is unstable) and taking properly into account the covariance between the two observables, we obtain:

* analysis of $`\xi /L`$ and $`U_4`$: $`\omega =0.70(11)`$, $`\beta _c=0.2857435(8)`$, $`(\xi /L)^{}=0.5943(8)`$, $`U_4^{}=1.6502(24)`$, $`b(\xi /L)=0.017(7)`$, $`b(U_4)=0.13(3)`$; $`\chi ^2=47.5`$, DOF$`=50`$.
* analysis of $`\xi /L`$ and $`U_6`$: $`\omega =0.71(11)`$, $`\beta _c=0.2857435(8)`$, $`(\xi /L)^{}=0.5942(8)`$, $`U_6^{}=3.318(11)`$, $`b(\xi /L)=0.017(7)`$, $`b(U_4)=0.67(16)`$; $`\chi ^2=46.7`$, DOF$`=50`$.

Here DOF is the number of degrees of freedom. The results are quite stable, indicating the presence of corrections with exponent $`\omega =0.7(1)`$, in agreement with the idea that scaling corrections with exponent $`\omega 0.4`$ are very small. The effective exponent $`\omega =0.7(1)`$ is quite close to the next-to-leading exponent predicted by perturbative field theory , i.e. $`\omega _2=0.8(2)`$. Thus, we mainly observe next-to-leading corrections, which in any case are quite small. In particular, they are of little relevance for $`\xi /L`$. The coefficient $`b`$, cf. Eq. (19), is very small and the estimates of $`\beta _c`$ obtained from the combined fits are fully compatible with those obtained for $`\xi /L`$ without scaling corrections.

These analyses that keep into account scaling corrections hint at values of $`\beta _c`$ lower than the estimate (18). We are thus led to consider

$$\beta _c=0.285744(2)$$
(20)
as our final estimate.

From the above-reported analyses, we also obtain estimates of the invariant ratios $`R^{}`$ at the critical point. We obtain

$`\left({\displaystyle \frac{\xi }{L}}\right)^{}`$ $`=`$ $`0.5943(9),`$ (21)
$`U_4^{}`$ $`=`$ $`1.650(9),`$ (22)
$`U_6^{}`$ $`=`$ $`3.32(5),`$ (23)
$`U_{22}^{}`$ $`=`$ $`0.1480(10).`$ (24)
We quote the results obtained in the fits (19), while the error is such to include also the result of the fit without scaling corrections and $`L_{\mathrm{min}}=64`$. For $`U_{22}`$ we only consider the fits without scaling corrections. Note that $`U_{22}^{}0`$, indicating the absence of self-averaging at the critical point, in agreement with the theoretical arguments of Ref. .

We can compare our results with those of Ref. . They found $`\beta _c=0.2857421(52)`$, $`U_4^{}=1.653(20)`$, $`U_{22}^{}=0.145(7)`$, $`(\xi /L)^{}=0.598(7)`$, which are in full agreement with our final results.

### C Determination of the four-point and six-point couplings

In this Section we wish to determine the four-point couplings $`G_4^{}`$ and $`G_{22}^{}`$, and the six-point universal ratios $`r_6^{}`$, $`C_{42}^{}`$, and $`C_{222}^{}`$, see appendix for definitions. Note that these quantities are defined in the high-temperature phase and one should take first the infinite-volume limit and then the limit $`\beta \beta _c`$, cf. Eq. (A23). In order to perform this task we have applied the extrapolation method of Refs.  to our high-temperature data, i.e. to the results with $`\beta 0.2856`$ (corresponding to $`\xi _{\mathrm{}}(\beta )89`$). This method is extremely powerful in order to compute the infinite-volume behavior of critical quantities and it has been applied to several models, including spin glasses .

The idea is the following. Given a long-distance quantity $`S(\beta ,L)`$, in the FSS limit we can write

$$\frac{S(\beta ,sL)}{S(\beta ,L)}=F_S\left[\xi (\beta ,L)/L\right]+O(L^\omega ,\xi ^\omega ),$$
(25)
where $`s`$ is an arbitrary (rational) number (here we always consider $`s=2`$). Here $`F_S(z)`$ is a universal function defined for $`0zz^{}(\xi /L)^{}=0.5943(9)`$ such that $`F_S(0)=1`$ and $`F_S(z^{})=s^\sigma `$, where $`\sigma `$ is the exponent characterizing the critical behavior of $`S(\beta ,L)`$ at the critical point, i.e., $`S(\beta _c,L)L^\sigma `$. Eq. (25) is the basis of the extrapolation technique since, in the absence of scaling corrections, it allows us to compute $`S(\beta ,sL)`$ on a lattice of size $`sL`$ in terms of quantities defined on a lattice of size $`L`$ and of the function $`F_S(z)`$. In practice, one works as follows. First, one performs several runs, determining $`S(\beta ,sL)`$, $`S(\beta ,L)`$, $`\xi (\beta ,sL)`$, and $`\xi (\beta ,L)`$. By means of a suitable interpolation, this provides the function $`F(z)`$ for $`S`$ and $`\xi `$. Then, $`S_{\mathrm{}}(\beta )`$ and $`\xi _{\mathrm{}}(\beta )`$ are obtained from $`S(\beta ,L)`$ and $`\xi (\beta ,L)`$ by iterating Eq. (25) and the corresponding equation for $`\xi (\beta ,L)`$. Of course, one must be very careful about scaling corrections, discarding systematically lattices with small values of $`L`$ till the results become independent of $`L`$ within error bars.

Let us first discuss the four-point couplings for which we will obtain quite precise estimates. In Figs. 1, 2, and 3 we report the data for $`S(\beta ,2L)/S(\beta ,L)`$ for $`G_4`$, $`G_{22}`$, and $`\xi `$, together with a fit of the data with $`L64`$. As discussed in Refs. , we parametrize $`F_S(z)`$ with a polynomial  in $`e^{1/z}`$ of order $`n`$, increasing $`n`$ until the $`\chi ^2`$ of the fit changes by less than one by going from $`n`$ to $`n+1`$. In these analyses we have taken $`n=7`$. The parametrization of $`F_S(z)`$ as a polynomial in $`e^{1/z}`$ is theoretically motivated: indeed, for zero-momentum quantities, $`F_S(z)`$ approaches 1 with corrections of order $`e^{a/z}`$, $`a1`$, as $`z0`$. This choice is not strictly correct for $`\xi `$ since in this case  $`F_\xi (z)=1+O(z^2)`$. However, these power corrections are expected to be very small for our definition of finite-volume correlation length , and therefore the systematic error due to our choice of parametrization should be small.

Looking at the figures, it is quite difficult to distinguish any correction to scaling, i.e. systematic deviations from the fitted curve. However, at a closer look one may see that some points with $`L=16`$ are out of the curve (in all cases by less than three error bars, so that these differences are barely significant). Conservatively, we have decided to discard all $`L=16`$ data. In order to check further for corrections to scaling we have computed infinite-volume estimates $`S_{\mathrm{}}(L)`$ using only data with $`LL_{\mathrm{min}}`$, $`L_{\mathrm{min}}=32`$, 64. Additionally, we have also systematically discarded points that are far from the critical point by including only data with $`\beta \beta _{\mathrm{min}}`$ for several values of $`\beta _{\mathrm{min}}`$.

Using the extrapolation procedure we have outlined above, for each $`L_{\mathrm{min}}`$ and $`\beta _{\mathrm{min}}`$ we obtain infinite-volume estimates $`\xi _{\mathrm{}}(\beta )`$, $`G_{4,\mathrm{}}(\beta )`$, and $`G_{22,\mathrm{}}(\beta )`$. To obtain the estimate at the critical point, the extrapolated values for the coupling constants have been fitted by using

$$S_{\mathrm{}}(\beta )=S^{}+a(\beta _c\beta ),$$
(26)
with $`\beta _c=0.285744(2)`$. The results are reported in Tables II and III for $`0.275\beta _{\mathrm{min}}0.284`$ (corresponding to $`4.45\xi _{\mathrm{}}15.86`$).

To check for corrections to scaling, we have also performed a different analysis. First we fit $`S(\beta ,2L)/S(\beta ,L)`$ taking into account a scaling correction with exponent $`\omega `$, i.e. assuming

$$\frac{S(\beta ,2L)}{S(\beta ,L)}=F_S\left(\xi (\beta ,L)/L\right)+\frac{1}{L^\omega }G_S\left(\xi (\beta ,L)/L\right),$$
(27)
and use both functions $`F_S(z)`$ and $`G_S(z)`$ to perform the infinite-volume extrapolation. Then, the infinite-volume results for the coupling constants are fitted by using

$$S_{\mathrm{}}(\beta )=S^{}+a(\beta _c\beta )^{\omega \nu }.$$
(28)
In order to perform the analysis we should fix the exponents $`\omega `$ and $`\nu `$. We use $`\nu =0.69`$, and repeat the analysis using $`\omega =0.8`$ (the next-to-leading exponent predicted by field theory ) and $`\omega =0.4`$, which is the leading exponent determined in Ref. . The results corresponding to $`\omega =0.8`$ are reported in Table IV. It is essential to include the results with $`L=16`$ in the analysis, otherwise the data do not show FSS corrections and the fit is unstable. Therefore, we cannot check the goodness of the Ansatz by discarding data with small $`L`$, i.e. present results for different values of $`L_{\mathrm{min}}`$ as done before.

Let us first discuss the results for $`G_4^{}`$. The fits without corrections to scaling show a significant (at the level of the reported errors) decrease as $`\beta _{\mathrm{min}}`$ is increased and also a slight dependence on $`L_{\mathrm{min}}`$. Corrections to scaling are positive and the estimate decreases with increasing $`\beta _{\mathrm{min}}`$, so that one only obtains an upper bound $`G_4^{}43.4`$. On the other hand, the fit with $`\omega =0.80`$ gives results independent of $`\beta _{\mathrm{min}}`$ within error bars; moreover the $`\chi ^2`$ of the fit (28) is systematically lower than that of the fit (26). Clearly the data are very well fitted by assuming a correction-to-scaling exponent $`\omega =0.80`$. For $`\omega =0.40`$ the results strongly depend on $`\beta _{\mathrm{min}}`$, varying from $`G_4^{}=42.49(14)`$ for $`\beta _{\mathrm{min}}=0.275`$ to $`G_4^{}=43.47(40)`$ for $`\beta _{\mathrm{min}}=0.2835`$. Also, the $`\chi ^2`$ is larger than the $`\chi ^2`$ obtained using $`\omega =0.8`$. There is therefore little evidence for such a small correction-to-scaling exponent, confirming again that for density $`x`$=0.8 the leading scaling corrections are very small. As final estimate we take

$$G_4^{}=43.3(2),$$
(29)
which is consistent with all results.

Let us now consider the results for $`G_{22}^{}`$. The results of Table III show no dependence on $`\beta _{\mathrm{min}}`$ and a tiny dependence on $`L_{\mathrm{min}}`$ which could be of purely statistical origin. With the present error bars there is no evidence for nonanalytic scaling corrections and indeed the results obtained using $`\omega =0.80`$ (see Table IV) are perfectly consistent with those of the fits with purely analytic corrections. Our final estimate is

$$G_{22}^{}=6.2(1).$$
(30)
The error is rather conservative and is such to include all estimates.

We should note that our results are compatible with the FT predictions of Ref. , where it is shown that in infinite volume $`G_{22,\mathrm{}}(\beta )`$ and $`G_{4,\mathrm{}}(\beta )`$ have scaling corrections with next-to-leading exponent $`\omega =0.8(2)`$ of similar relative size. In particular, the results of Ref.  give $`a_{G_{22}}=0.23(10)a_{G_4}`$ where $`a`$ is the coefficient defined in Eq. (28). From the fits we obtain instead $`a_{G_4}=20(5)`$ and $`a_{G_{22}}=1.5(1.0)`$. The errors should be taken with caution: They simply give the variation of the parameter $`a`$ with $`\beta _{\mathrm{min}}`$ for $`0.275\beta _{\mathrm{min}}0.281`$ (for larger values the statistical error becomes larger than $`|a|`$) and do not include any possible systematic effect. Assuming these values with their errors, we estimate $`a_{G_{22}}=0.08(5)a_{G_4}`$, which is in reasonable agreement with the FT result.

As a check we have repeated the analysis for $`H_4G_4+3G_{22}`$. Since the procedure is nonlinear, this represents an important consistency check. We obtain

$$H_4^{}=24.7(2),$$
(31)
which is in full agreement with the estimates reported above.

Finally, we consider the six-point couplings $`r_6`$, $`C_{42}`$, and $`C_{222}`$. We apply again the extrapolation procedure we have used for $`G_4`$ and $`G_{22}`$. However, in this case there are larger systematic errors. The extrapolation curve $`F(z)`$, cf. Eq. (25), is poorly determined for $`z0.3`$, since the six-point couplings have large statistical errors when $`\xi (\beta ,L)/L`$ is small. A large error on the curve gives a large systematic error on the extrapolations and induces correlations among the results for different $`\beta `$ (such correlations are instead small for $`G_4`$ and $`G_{22}`$). Because of them, it is difficult to set correct error bars and to determine the final results. For these reasons we have taken a conservative attitude. We have generated different sets of extrapolated data by changing the degree $`n`$ of the interpolation polynomial and $`\beta _{\mathrm{min}}`$ (in all cases we set $`L_{\mathrm{min}}=16`$ in order to have a sufficiently large number of points in the small-$`z`$ region). Then, we determine the range that includes most (approximately 2/3) of the extrapolated data with their error bars. The central value of such an interval gives the final result, while its half width gives the error. We obtain

$`r_6^{}`$ $`=`$ $`0.90(15),`$ (32)
$`C_{42}^{}`$ $`=`$ $`0.12(5),`$ (33)
$`C_{222}^{}`$ $`=`$ $`0.45(15).`$ (34)
In Fig. 4 we report the final results together with the results for a single extrapolation, so that the reader can see the quality of our numerical results and how much conservative our final error bars are. Note that in this analysis we have not made any attempt to evaluate corrections to scaling. In any case, we expect them to be small compared with the large errors we quote.

The numerical estimates can be compared with the FT results. We shall discuss the four-point couplings in Sec. III. For what concerns the six-point couplings, the analysis of the available four-loop series in the fixed-dimension scheme  is not very precise. We only mention the estimate $`r_6^{}=1.1_{0.5}^{+0.1}`$ reported in Ref. , which agrees with the more precise result (34).

### D The exponent $`\nu `$

In order to determine the exponent $`\nu `$ we have followed two different strategies: (a) we have determined extrapolated values $`\xi _{\mathrm{}}(\beta )`$ using our extrapolation method and then we have performed a fit $`\xi _{\mathrm{}}(\beta )(\beta _c\beta )^\nu `$; (b) we have performed a purely FSS analysis without extrapolations. The advantage of method (b) is that we can keep both $`\beta _c`$ and $`\omega `$ as free parameters, and thus check the previously determined value for $`\beta _c`$ and our claim that scaling corrections with small values of $`\omega `$ are tiny.

In the first case we proceed as before. If we neglect scaling corrections, we can generate infinite-volume estimates $`\xi _{\mathrm{}}(\beta )`$ by using the extrapolation precedure based on Eq. (25) and then we can determine $`\nu `$ from fits of the form

$$\mathrm{ln}\xi _{\mathrm{}}(\beta )=\nu \mathrm{ln}(\beta _c\beta )+a+b(\beta _c\beta ).$$
(35)
As a second possibility we can include corrections with exponent $`\omega `$ in the extrapolation, see Eq. (27), and fit the extrapolated data with

$$\mathrm{ln}\xi _{\mathrm{}}(\beta )=\nu \mathrm{ln}(\beta _c\beta )+a+b(\beta _c\beta )^{\omega \nu }.$$
(36)
We report in Fig. 5 the results for the analytic fit with $`L_{\mathrm{min}}=32`$ and $`L_{\mathrm{min}}=64`$ and for the nonanalytic fit with $`\omega =0.80`$ and $`L_{\mathrm{min}}=16`$ for several values of $`\beta _{\mathrm{min}}`$. We observe that results with and without nonanalytic corrections differ significantly, much more than the statistical errors. However, corrections have opposite trends. Analytic fits give results that decrease as $`\beta _{\mathrm{min}}`$ is increased while nonanalytic fits have the opposite behavior. Compatible results are obtained for $`\beta _{\mathrm{min}}0.283`$. For $`\beta =0.283`$ we obtain $`\nu =0.6847(15)`$ from the first analysis with $`L_{\mathrm{min}}=64`$ and $`\nu =0.6818(24)`$ from the second analysis. A reasonable estimate would be $`0.683(3)`$ which includes all results. It should be noted that the results depend strongly on $`\beta _c`$ that has been fixed to $`\beta _c=0.285744(2)`$. By varying $`\beta _c`$ within one error bar, $`\nu `$ varies approximately by $`\pm 0.003`$. Thus, collecting everything together, this method gives the final result

$$\nu =0.683(3+3).$$
(37)
The estimate (37) has been obtained by using only the high-temperature data, i.e. those with $`\beta 0.2856`$. A more precise estimate can be obtained by performing a direct FSS analysis that allows us to use both the high-temperature and the critical-point results for the correlation length. In this way, we do not need to fix $`\beta _c`$ and we can also keep $`\omega `$ as a free parameter. We start from the general FSS expression

$`{\displaystyle \frac{\xi }{L}}=F[u_1(t)L^{1/\nu },u_2(t)L^\omega ],`$ (38)
where only one correction-to-scaling operator has been taken into account. Here $`t(\beta _c\beta )`$, and $`u_1(t)`$ and $`u_2(t)`$ are nonlinear scaling fields satisfying $`u_1(0)=0`$ and $`u_2(0)0`$. Equation (38) can be rewritten in the scaling limit as

$`{\displaystyle \frac{\xi }{L}}=f_1(tL^{1/\nu })+L^{1/\nu }f_2(tL^{1/\nu })+L^\omega f_3(tL^{1/\nu }),`$ (39)
with $`f_1(0)0`$ and $`f_3(0)0`$. The three scaling functions represent the three types of contributions we expect: $`f_1(x)`$ is the leading FSS curve, $`f_2(x)`$ corresponds to the analytic corrections, and $`f_3(x)`$ is the nonanalytic FSS correction. The function $`f_2(x)`$ is related to $`f_1(x)`$ by $`f_2(x)x^2f_1^{}(x)`$, a relation that follows from the fact that this correction is due to the expansion of the scaling field $`u_1(t)=t+at^2+O(t^3)`$. The existence of the infinite-volume limit fixes the behavior of these functions for $`x\mathrm{}`$: $`f_1(x)x^\nu `$, $`f_2(x)x^{1\nu }`$ and $`f_3(x)x^{(\omega 1)\nu }`$.

Now, we wish to determine $`\beta _c`$, $`\nu `$, and $`\omega `$ by fitting the numerical data for $`\xi /L`$ with Eq. (39). For this purpose, we must parametrize the scaling functions reported above with expressions containing a finite number of parameters. Thus, we consider an even integer $`n`$ and parametrize:

$`{\displaystyle \frac{\xi }{L}}`$ $`=`$ $`\left[P_n(tL^{1/\nu })+L^\omega (1+ptL^{1/\nu })^{\omega \nu }Q_n(tL^{1/\nu })\right]^{\nu /n},`$ (40)
$`P_n(x)`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{n}{}}}a_ix^i,Q_n(x)={\displaystyle \underset{i=0}{\overset{n/2}{}}}b_ix^{2i},`$ (41)
where $`t(\beta _c\beta )`$. This parametrization has the correct behavior both for small and large $`x`$ for any $`n`$. Note that our choice of writing $`Q_n(x)`$ as a polynomial in $`x^2`$ is due only to practical considerations: since corrections are small, $`Q_n(x)`$ cannot be parametrized with as many parameters as $`P_n(x)`$; on the other hand, to guarantee the correct large-$`x`$ behavior $`Q_n(x)`$ and $`P_n(x)`$ must have the same degree. We also made some analyses writing $`Q_n(x)`$ as a polynomial in $`x^3`$—in this case $`n`$ must be a multiple of 3—without finding significant differences. The value of $`n`$ has been chosen on the basis of the $`\chi ^2`$ of the fit: $`n`$ has been increased until the $`\chi ^2`$ does not change significantly as $`n`$ is increased by 2. The results of the analyses we present correspond to $`n=12`$.

In the analyses we have kept $`\nu `$, $`\omega `$, $`p`$, $`\beta _c`$, $`\{a_i\}`$, and $`\{b_i\}`$ as free parameters. The fit is stable only if we include the data with $`L=16`$, otherwise there are no scaling corrections at the level of our error bars and therefore it is impossible to determine $`\omega `$. To check for corrections, we have systematically discarded points far from $`\beta _c`$, including in the analysis only data with $`\beta \beta _{\mathrm{min}}`$ for several values of $`\beta _{\mathrm{min}}`$. In order to have a mostly linear fit, we have fitted $`(\frac{\xi }{L})^{n/\nu }`$. The results are reported in Table V.

The results for $`\beta _c`$ are quite stable, in full agreement with the analyses reported in Sec. II B for $`\xi /L`$ and with the final estimate (20). The estimates of $`\omega `$ vary significantly and have a quite large error. There is however very little evidence of scaling corrections with $`\omega 1`$, as already discussed in Sec. II B. Finally, let us consider $`\nu `$. As it can be seen from Fig. 5 the results are in rough agreement with those found before and apparently become independent of $`\beta _{\mathrm{min}}`$ for $`\beta 0.282`$. Our final estimate is

$$\nu =0.683(3),$$
(42)
which is compatible with the results of all analyses. The estimate (42) is in very good agreement with the result $`\nu =0.6837(53)`$ of Ref. .

### E The exponent $`\eta `$

In order to compute the exponent $`\eta `$ and correspondingly $`\gamma =(2\eta )\nu `$, one can analyze the critical behavior of the susceptibility $`\chi `$. However, we have found much more convenient to analyze $`\zeta \chi /\xi ^2`$. Indeed, because of statistical correlations between $`\chi `$ and $`\xi `$, the relative error on $`\zeta `$ is significantly smaller than that on $`\chi `$. Moreover, the FT analysis of Ref.  indicates that $`\zeta `$ has much smaller scaling corrections than $`\chi `$.

In order to determine $`\eta `$ we perform two different analyses. The first one uses the high-temperature data and follows closely the analysis of the correlation length presented in Sec. II D. Given $`\zeta (\beta ,L)`$ we determine the curve $`F_\zeta (z)`$, cf. Eq. (25). The ratios $`\zeta (\beta ,2L)/\zeta (\beta ,L)`$ are reported in Fig. 6, together with a fit of the points with $`L64`$. There are clear scaling corrections, especially when $`\xi (\beta ,L)/L0.3`$: the data with $`L=16`$, and also some points with $`L=32`$, are systematically above the curve. As before, we discard all points with $`L=16`$ and perform the extrapolation using data with $`LL_{\mathrm{min}}`$, with $`L_{\mathrm{min}}=32,64`$, in order to detect scaling corrections. The extrapolated values are fitted by using

$$\mathrm{ln}\zeta _{\mathrm{}}(\beta )=\eta \nu \mathrm{ln}(\beta _c\beta )+a+b(\beta _c\beta )$$
(43)
and fixing $`\beta _c=0.285744(2)`$. The results are reported in Table VI. They show a systematic dependence on $`L_{\mathrm{min}}`$ and $`\beta _{\mathrm{min}}`$ and seem to indicate $`\eta \nu 0.0240`$, but it is difficult to set an error bar. Conservatively, we write

$$\eta \nu =0.0240(15),$$
(44)
that includes the estimates with $`\beta _{\mathrm{min}}=0.284`$ with their error bars. Using the estimate of $`\nu `$ reported in Sec. II D, $`\nu =0.683(3)`$, we obtain

$$\eta =0.035(2).$$
(45)
We have also repeated the analysis including scaling corrections with $`\omega =0.8`$. In the extrapolation procedure we use Eq. (27), and then we fit the extrapolated data by using

$$\mathrm{ln}\zeta _{\mathrm{}}(\beta )=\eta \nu \mathrm{ln}(\beta _c\beta )+a+b(\beta _c\beta )^{\omega \nu }.$$
(46)
The results are reported in Table VII. They show a systematic upward trend, indicating that scaling corrections with $`\omega =0.8`$ do not correctly describe the data. This is also evident from the $`\chi ^2`$ of the fit (46) which is quite large for small $`\beta _{\mathrm{min}}`$. A lower value of $`\omega `$, say $`\omega =0.4`$, gives even worse results, confirming that corrections with $`\omega 1`$ are very small, in agreement with the FT analysis of Ref. . In any case the results of Table VII indicate that $`\eta \nu 0.0220`$, in agreement with Eq. (44).

The second analysis we consider uses the data at the critical point. For $`\beta \beta _c`$ we have

$$\zeta (\beta ,L)=L^\eta f_1(tL^{1/\nu })+L^{\eta \omega }f_2(tL^{1/\nu }),$$
(47)
where $`t\beta _c\beta `$. If we expand for $`t0`$ we can write

$$\mathrm{ln}\zeta (\beta ,L)=\eta \mathrm{ln}L+a+btL^{1/\nu }+cL^\omega .$$
(48)
We first perform an analysis neglecting scaling corrections ($`c=0`$) for several values of $`L_{\mathrm{min}}`$. If we fix $`\beta _c=0.285744(2)`$ and $`\nu =0.683(3)`$ we obtain:

* $`L_{\mathrm{min}}=16`$: $`\eta =0.0331(2+6)`$, $`\chi ^2`$/DOF = 76/26;
* $`L_{\mathrm{min}}=32`$: $`\eta =0.0347(3+8)`$, $`\chi ^2`$/DOF = 27/21;
* $`L_{\mathrm{min}}=64`$: $`\eta =0.0354(6+14)`$, $`\chi ^2`$/DOF = 22/13.

The first error is statistical, while the second one is due to the error on $`\beta _c`$. The error due to the error on $`\nu `$ can be neglected. There are apparently some scaling corrections, but it should be noted that the difference between the results with $`L_{\mathrm{min}}=32`$ and 64 is not significant given the statistical error bars. We also perform an analysis with scaling corrections with exponent $`\omega =0.80`$. Using all data, i.e. taking $`L_{\mathrm{min}}=16`$, we obtain $`\eta =0.0382(8+16)`$, $`c=0.111(16+23)`$, $`\chi ^2`$/DOF = 29/25. This result is higher than the estimates obtained before, but still compatible with the quoted errors.

The analysis at the critical point gives therefore results that are in full agreement with those obtained before, confirming the correctness of the estimate (45), and with the result of Ref. , $`\eta =0.0374(45)`$.

### F The exponent $`\alpha `$ and hyperscaling

In this Section we wish to determine the exponent $`\alpha `$ and check the hyperscaling relation $`2\alpha =d\nu `$. Unfortunately, we have measured the specific heat $`C(\beta ,L)`$ only near the critical point and thus we can determine $`\alpha `$ only from the behavior for $`\beta \beta _c`$.

RG predicts that, in the FSS limit, the energy scales as

$$E(\beta ,L)E_{\mathrm{bulk}}(\beta )+L^{(\alpha 1)/\nu }g_1(tL^{1/\nu })+L^{(\alpha 1)/\nu \omega }g_2(tL^{1/\nu }),$$
(49)
where $`t\beta _c\beta `$ and we have included one scaling correction. Note that $`E_{\mathrm{bulk}}(\beta )`$ is expected to have an exponentially small dependence on $`L`$ which can be neglected for all practical purposes. Near the critical point we can expand this expression in powers of $`tL^{1/\nu }`$ obtaining for the energy and the specific heat the expressions

$`E(\beta ,L)`$ $``$ $`a_E+b_Et+c_EL^{(\alpha 1)/\nu }+d_EtL^{\alpha /\nu }+e_EL^{(\alpha 1)/\nu \omega },`$ (50)
$`C(\beta ,L)`$ $``$ $`a_C+b_Ct+c_CL^{\alpha /\nu }+d_CtL^{(\alpha +1)/\nu }+e_CL^{\alpha /\nu \omega },`$ (51)
with $`a_C=b_E`$ and $`c_C=d_E`$.

We wish now to determine $`\alpha /\nu `$. We consider the specific heat, since in this case the singular behavior is stronger. In Eq. (51) we have not used hyperscaling, so that the expression depends on two independent exponents $`\alpha `$ and $`\nu `$, or more precisely $`\alpha /\nu `$ and $`1/\nu `$. In order to simplify the analysis, we now analyze the specific heat using hyperscaling to rewrite the correction term. Thus, Eq. (51) becomes

$`C(\beta ,L)`$ $`=`$ $`a_C+b_Ct+c_CL^{\alpha /\nu }+d_CtL^{3/2(\alpha /\nu +1)}+e_CL^{\alpha /\nu \omega },`$ (52)
that depends only on $`\alpha /\nu `$. We begin by neglecting the correction-to-scaling term, i.e. we set $`e_C=0`$. Fixing $`\beta _c=0.285744(2)`$, a five-parameter fit of the data with $`LL_{\mathrm{min}}`$ gives:

* $`L_{\mathrm{min}}=16`$, $`\alpha /\nu =0.119(7+6)`$, $`\chi ^2`$/DOF = 12.0/24;
* $`L_{\mathrm{min}}=32`$, $`\alpha /\nu =0.115(17+11)`$, $`\chi ^2`$/DOF = 11.6/19.

The first error is the statistical one, while the second one gives the variation of the estimate as $`\beta _c`$ is varied by one error bar. Two things should be noticed: the analytic term proportional to $`t`$ is very small compared with the statistical errors. Indeed, $`b_C=23(83+8)`$ and $`b_C=30(167+24)`$ in the two fits. Moreover, corrections to scaling are apparently small, since the results do not depend on $`L_{\mathrm{min}}`$. As a check, we have also performed an analysis with a correction term ($`e_C0`$). We find $`\alpha /\nu =0.109(50+6)`$, $`e_C=0.3(1.4)`$, and $`\chi ^2`$/DOF = 11.9/23, for $`\omega =0.8`$ and $`L_{\mathrm{min}}=16`$. The coefficient $`e_C`$ is compatible with zero, while the $`\chi ^2`$ of the fit changes only by 0.1 in spite of the fact that we have added an additional parameter. There is no evidence of scaling corrections. As final estimate of $`\alpha /\nu `$ we take the result with $`L_{\mathrm{min}}=32`$ and $`e_C=0`$,

$$\frac{\alpha }{\nu }=0.115(28).$$
(53)
If we use hyperscaling and the estimate of $`\nu `$ of Sec. II D, $`\nu =0.683(3)`$, we obtain $`\alpha /\nu =0.072(13)`$. This estimate is in reasonable agreement with (53) confirming the validity of hyperscaling. Using Eq. (53) and $`\nu =0.683(3)`$ we obtain $`\alpha =0.079(19)`$, which is in reasonable agreement with the estimate obtained using hyperscaling $`\alpha =0.049(9)`$.

We can also check hyperscaling directly by comparing the results for the energy and the specific heat. If we define $`\theta _1\alpha /\nu `$, $`\theta _2(2\alpha )/\nu `$, we can rewrite Eqs. (50) and (51) as

$`E(\beta ,L)`$ $`=`$ $`a_E+b_Et+c_EL^{(\theta _1\theta _2)/2}+d_EtL^{\theta _1},`$ (54)
$`C(\beta ,L)`$ $`=`$ $`b_E+d_EL^{\theta _1}+d_CtL^{(3\theta _1+\theta _2)/2},`$ (55)
where we have neglected scaling corrections and we have set $`b_C=0`$ in agreement with the previous analysis. If hyperscaling is satisfied we should find $`\theta _2=3`$. A combined analysis of $`E(\beta ,L)`$ and $`C(\beta ,L)`$ fixing $`\beta _c`$ gives

* $`L_{\mathrm{min}}=16`$: $`\alpha /\nu =0.119(6+6)`$, $`\theta _2=2.952(17+8)`$, $`\chi ^2`$/DOF = 33.4/51;
* $`L_{\mathrm{min}}=32`$: $`\alpha /\nu =0.112(16+9)`$, $`\theta _2=2.930(43+12)`$, $`\chi ^2`$/DOF = 29.3/41.

Here, as before, the first error is the statistical one, while the second gives the error due to $`\beta _c`$. The estimates of $`\alpha /\nu `$ are compatible with Eq. (53), while $`\theta _2`$ is in reasonable agreement with 3, confirming again the validity of hyperscaling.

### G The universal ratio $`R_\xi ^+`$

We now wish to compute the universal ratio $`R_\xi ^+`$ that is related to the universality of the singular part of the free energy per correlation volume. More precisely, let us consider the infinite-volume free energy density $`_{\mathrm{}}(\beta )`$

$$_{\mathrm{}}=\frac{1}{V}\overline{\mathrm{ln}Z},$$
(56)
and the infinite-volume specific heat $`C_{\mathrm{}}(\beta )=^2_{\mathrm{}}(\beta )/\beta ^2`$. In the critical limit $`t\beta _c\beta 0`$ we can write

$`_{\mathrm{}}(\beta )=_{\mathrm{bulk}}(\beta )+F^+t^{2\alpha }+\text{corrections},`$ (57)
$`C_{\mathrm{}}(\beta )=C_{\mathrm{bulk}}(\beta )+A^+t^\alpha +\text{corrections},`$ (58)
where $`A^+=(2\alpha )(1\alpha )F^+`$. If $`\xi _{\mathrm{}}(\beta )f^+t^\nu `$ in the same limit $`t0^+`$, RG predicts the universality of $`F^+(f^+)^3`$, or equivalently of

$`R_\xi ^+\left(\alpha A^+\right)^{1/3}f^+.`$ (59)
We now compute $`R_\xi ^+`$ using our high-temperature results for the energy $`E(\beta ,L)`$ and for the correlation length $`\xi (\beta ,L)`$. We define a quantity $`R(\beta ,L)`$,

$$R(\beta ,L)\left[\frac{\alpha }{1\alpha }(E(\beta ,L)a_Eb_Et)t\right]^{1/3}\xi (\beta ,L),$$
(60)
where $`a_E`$ and $`b_E`$ are defined in terms of the expansion of $`_{\mathrm{bulk}}(\beta )`$ for $`t0`$:

$$\frac{d_{\mathrm{bulk}}(\beta )}{d\beta }E_{\mathrm{bulk}}(\beta )=a_E+b_Et+O(t^2).$$
(61)
It is easy to check that

$$\underset{t0}{lim}\underset{L\mathrm{}}{lim}R(\beta ,L)=R_\xi ^+.$$
(62)
In order to compute $`R(\beta ,L)`$, we must specify the values of the two constants $`a_E`$ and $`b_E`$ in Eq. (60). For this purpose we exploit the fact that $`E_{\mathrm{bulk}}(\beta )`$ is the same function in Eq. (49) and in Eq. (61), so that $`a_E`$ and $`b_E`$ coincide with the constants defined in Eq. (50). Thus, $`a_E`$ and $`b_E`$ can be determined independently by using the critical-point data for the energy and the specific heat. We thus proceed as follows. We consider Eq. (55), fix $`\beta _c=0.285744(2)`$, $`\theta _2=3`$, $`\alpha /\nu =2/\nu 3`$, $`\nu =0.683(3)`$, and compute $`a_E`$ and $`b_E`$ by analyzing $`E(\beta ,L)`$ and $`C(\beta ,L)`$ near the critical point. Then, we determine $`R(\beta ,L)`$. The error on $`R(\beta ,L)`$ takes into account the error on $`E(\beta ,L)`$, $`\xi (\beta ,L)`$, $`a_E`$, and $`b_E`$, and also the variation of the estimates as $`\nu `$ and $`\beta _c`$ vary within one error bar. In order to be conservative, we use a worst-error estimate summing all errors together. Once $`R(\beta ,L)`$ has been computed we use the extrapolation method presented in Sec. II C.

In Fig. 7 we report the ratios $`R(\beta ,2L)/R(\beta ,L)`$ together with a fit of the data with $`L64`$ (a good fit is obtained by using a polynomial with $`n=10`$). Apparently, there are no scaling corrections, but at a closer look one finds systematic deviations for $`L=16`$. As before, these points will be discarded in the analytic fits.

The results of the fits with analytic corrections are reported in Table VIII for $`L_{\mathrm{min}}=32`$ and $`L_{\mathrm{min}}=64`$. For small $`\beta _{\mathrm{min}}`$ they show an upward trend and then apparently stabilize around $`0.2874(3)`$. In order to check the role of the corrections to scaling, we have repeated the analysis by adding scaling corrections with exponent $`\omega =0.80`$. The results are presented in Table IX. They are now substantially independent of $`\beta _{\mathrm{min}}`$ confirming that the data are very well fitted by assuming such an exponent. The final estimate is somewhat higher than that obtained in the analytic fits, indicating that in this case nonanalytic scaling correction may play an important role. We do not know which of the two fits is the most reliable one and thus we have taken as final estimate

$$R_\xi ^+=0.2885(15),$$
(63)
which is compatible with the results of both analyses.

The estimate (63) is in good agreement with the results of other methods. A six-loop computation in the fixed-dimension FT approach gives  $`R_\xi ^+=0.290(10)`$, while, by using approximate parametric representations of the equation of state, one obtains  $`R_\xi ^+=0.282(3)`$.

## III Comparison with field-theory results

The critical behavior of the RIM has been extensively studied using the FT approach. Quantitative predictions can be obtained by using different techniques: perturbative methods in the four-point renormalized couplings in fixed dimension $`d=3`$ or in $`\sqrt{ϵ}`$, $`ϵ4d`$, or nonperturbative methods based on approximate RG equations, see Refs.  The most accurate results have been obtained in the first approach: six-loop expansions for the $`\beta `$-functions and the critical exponents have been derived and analyzed in Refs. . The corresponding estimates of the critical exponents, e.g., $`\nu =0.678(10)`$ and $`\eta =0.030(3)`$, are in satisfactory agreement with the Monte Carlo results presented before.

The main problem of the perturbative approach is the non-Borel summability of the series . This fact makes the analysis more subtle and less precise than in the case of the pure Ising model. The difficulties of the perturbative approach appear in the determinations of the fixed-point values $`u^{}`$ and $`v^{}`$ of the renormalized couplings (they are normalized so that at tree level $`u=u_0/m`$, $`v=v_0/m`$, $`m`$ being the renormalized mass), which are directly related to the quantities we have measured in the Monte Carlo simulation (see Ref.  for a derivation of these relations): $`G_4^{}=v^{}`$, $`G_{22}^{}=u^{}/3`$, and $`H_4^{}G_4^{}+3G_{22}^{}=u^{}+v^{}`$. The analysis of the six-loop series provided results somewhat dependent on the resummation method . Indeed, we found : $`G_4^{}=37.7(2)`$, $`G_{22}^{}=4.3(6)`$, $`H_4^{}=24.8(1.8)`$ (double Padé-Borel method); $`G_4^{}=36.8(3)`$, $`G_{22}^{}=4.0(1)`$, $`H_4^{}=24.8(6)`$ (conformal Padé-Borel method); $`G_4^{}=38.6(7)`$, $`G_{22}^{}=4.8(2)`$, $`H_4^{}=24.2(9)`$ (direct conformal method). The estimates of $`H_4^{}`$ are in good agreement with the Monte Carlo result (31), $`H_4^{}=24.7(2)`$. On the other hand, the estimates of $`G_4^{}`$ and $`G_{22}^{}`$—combining the results we would have guessed $`G_4^{}=38.0(1.5)`$ and $`G_{22}^{}=4.5(6)`$ with errors that are, at first sight, quite conservative—differ significantly from the Monte Carlo estimates (29) and (30).

These discrepancies call for a reanalysis of the perturbative series of the exponents, verifying if the use of the Monte Carlo results for $`G_4^{}`$ and $`G_{22}^{}`$ leads to significantly different estimates. We have thus repeated the analysis, using the different resummation methods outlined in Refs. . We find

$$\nu =0.686(4),\eta =0.026(3),\gamma =1.355(8),$$
(64)
where the errors include the results of the different resummation methods. It is reassuring that these estimates are close to those found in Ref. , $`\nu =0.678(10)`$, $`\eta =0.030(3)`$, and $`\gamma =1.330(17)`$, and also reasonably close to the Monte Carlo estimates. The small variation of the estimates of the critical exponents is due to the particular structure of the perturbative series: if they are rewritten in terms of $`yu+v`$ and $`u`$, the resummations depend mostly on $`y^{}=H_4^{}`$, which is correctly determined by FT methods, and only slightly on $`u^{}`$ that is instead poorly known. We should also observe that the new estimate of $`\nu `$ is closer to the Monte Carlo result, while the estimate of $`\eta `$ is slightly worse. Therefore, the FT estimates do not become more accurate if more precise results for $`u^{}`$ and $`v^{}`$ are used. This is an indication that, at least for the critical exponents, the location of the fixed point is not the main source of error on the results.

We also tried an alternative procedure based on an expansion of the RG functions around the unstable Ising fixed point $`u=0`$, $`v=g_I^{}`$, where  $`g_I^{}=23.56(2)`$. The analysis of the Ising-to-RIM RG flow reported in Ref.  and the discussion reported above show that it is convenient to introduce new variables $`yu+v`$ and $`zu`$. In terms of $`y`$ and $`z`$, the RIM fixed point is located in $`y^{}=H_4^{}=24.7(2)`$ and $`z^{}=3G_{22}^{}=18.6(3)`$, while the Ising fixed point is at $`y_I=g_I^{}`$, $`z=0`$. Then, we write $`y=g_I+\delta y`$, obtaining for any RG function $`f(y,z)`$,

$`f(y,z)={\displaystyle \underset{i,j}{}}c_{ij}(g_I)\delta y^iz^j,`$ (65)
$`c_{ij}(g_I)={\displaystyle \underset{k}{}}f_{ijk}g_I^k.`$ (66)
The value of $`f(y,z)`$ at the RIM fixed point is then obtained as follows. First, we compute the coefficients $`c_{ij}(g_I^{})`$ at the Ising fixed point, by using the conformal-mapping method and exploiting the known large-order behavior of the expansion of $`c_{ij}(g_I)`$ that is determined by the Ising fixed point. Then, we evaluate the double series at $`\delta y^{}=H_4^{}g_I^{}=1.14(20)`$ and $`z^{}=3G_{22}^{}=18.6(3)`$. Here, we are neglecting the fact that the RG functions are nonanalytic at the Ising fixed point . Note that $`f(g_I^{},0)`$ is the value of the same quantity for the Ising universality class, so that the expansion (65) provides the differences between RIM and Ising critical exponents, i.e. $`\mathrm{\Delta }f=f_{\mathrm{RIM}}f_{\mathrm{Ising}}`$, which are expected to be rather small. Of course, this expansion is, at most, asymptotic. But one may hope that the RIM and Ising fixed points are sufficiently close, so that the first few terms of the expansion around the Ising fixed point allow us to obtain more accurate estimates of the exponents of the RIM. In Table X we report the estimates of the first coefficients $`c_{ij}(g_I^{})`$ for $`\eta \eta _I`$, $`\nu \nu _I`$, $`\gamma \gamma _I`$, $`1/\nu 1/\nu _I`$, and $`1/\gamma 1/\gamma _I`$. The results for the critical exponents are reported in Table XI as a function of the order $`o\mathrm{Max}[i+j]`$ of the expansion. We thus obtain the following estimates: $`\eta \eta _I=0.0017(13)`$, $`\nu \nu _I=0.060(5)`$, $`1/\nu 1/\nu _I=0.135(15)`$, $`\gamma \gamma _I=0.12(1)`$, $`1/\gamma 1/\gamma _I=0.07(1)`$. The estimate we quote corresponds to $`o=3`$, while the error is such to include the results with $`o=2`$ and $`o=4`$. Then, using the estimates  $`\eta _I=0.0364(5)`$, $`\nu _I=0.6301(4)`$, $`\gamma _I=1.2372(5)`$, we finally obtain

$$\nu =0.690(5),\eta =0.035(2),\gamma =1.357(10),$$
(67)
which are again in substantial agreement with the MC results.

The expansion around the Ising fixed point can also be performed along the Ising-to-RIM RG trajectory , which is obtained as the limit $`u_00^{}`$ of the RG trajectories in the $`u`$, $`v`$ plane. At least in principle, this expansion is expected to be better behaved than the previous one, since RG functions should be analytic near the Ising fixed point only along this trajectory, up to the random fixed point where nonanalyticites are again present .

An effective parametrization of the curve is given by the first few terms of its expansion around $`z=0`$, which is given by

$$yy_I=T(z)=c_2z^2+c_3z^3+\mathrm{}$$
(68)
where $`c_2=0.0033(1)`$ and $`c_3=1(2)\times 10^5`$. The fact that $`yy_I`$ is of order $`z^2`$ is the main reason why we introduced the variable $`y`$ and is due to the identity

$`{\displaystyle \frac{\beta _v}{u}}|_{u=0}+{\displaystyle \frac{\beta _u}{u}}|_{u=0}{\displaystyle \frac{\beta _v}{v}}|_{u=0}=0.`$ (69)
Substituting the expansion (68) into the double expansion (65), we obtain an expansion in powers of $`z`$

$$f(y,z)=f(g_I+T(z),z)=\underset{i}{}e_i(g_I)z^i,$$
(70)
which should be then evaluated at $`g_I=g_I^{}`$ and $`z=z^{}`$. The values of the coefficients $`e_i(g_I^{})`$ at the Ising fixed point have been computed by using a conformal mapping and a Borel transform. The results for $`i=1,2,3`$ are reported in Table XII. The estimates of the difference between the critical exponents of the Ising and RIM universality classes are $`\eta \eta _I=0.0020(18)`$, $`\nu \nu _I=0.060(5)`$, $`1/\nu 1/\nu _I=0.136(20)`$, $`\gamma \gamma _I=0.119(10)`$, and $`1/\gamma 1/\gamma _I=0.070(10)`$. These results are obtained by truncating the expansion (70) to third order, while the error is the sum of the uncertainty due to the resummation, due to the truncation of the series (the difference between the second-order and the third-order result), and due to the uncertainty on $`z^{}`$ (we used the Monte Carlo result). Then, by using the estimates  $`\eta _I=0.0364(5)`$, $`\nu _I=0.6301(4)`$, $`\gamma _I=1.2372(5)`$, we finally obtain

$$\nu =0.690(8),\eta =0.0345(20),\gamma =1.355(10),$$
(71)
that do not differ significantly from the estimates (67).

Note that the estimate of $`\eta `$ obtained by using the expansion around the Ising fixed point is now in perfect agreement with the Monte Carlo result, at variance with the direct estimate (64). The estimate of $`\nu `$ is also in substantial agreement with the numerical estimate $`\nu =0.683(3)`$. Therefore, the expansion around the Ising fixed point appears to be a useful alternative method to compute the critical properties of the RIM.

## Acknowledgments

V.M.-M. is a Ramon y Cajal research fellow and is partly supported by MCyT (Spain), project Nos. FPA2001-1813 and FPA2000-0956. The numerical computations were performed at the PC Cluster at the University of Pisa. We thank Maurizio Davini for his indispensable technical assistance.

## A Notations

We consider the Hamiltonian (1) with $`J=1`$ on a finite lattice $`L^3`$ with periodic boundary conditions. Given a quantity $`𝒪`$ depending on the spins $`\{s\}`$ and on the random variable $`\{\rho \}`$ we define the sample average at fixed distribution $`\{\rho \}`$

$$𝒪(\beta ,\{\rho \})\frac{1}{Z(\{\rho \})}\underset{\{s_i\}}{}𝒪e^{\beta [s,\rho ]},$$
(A1)
where $`Z(\{\rho \})`$ is the sample partition function. Of course, we are interested in averaging over the random dilution and thus we consider disorder-averaged quantities

$$\overline{𝒪}(\beta )=[d\rho ]𝒪(\beta ,\{\rho \}),$$
(A2)
where

$$[d\rho ]=\underset{i}{}[x\delta (\rho _i1)+(1x)\delta (\rho _i)].$$
(A3)
We define the two-point correlation function and the susceptibility $`\chi (\beta ,L)`$

$`G(x;\beta ,L)`$ $``$ $`\overline{\rho _0s_0\rho _xs_x}`$ (A4)
$`\chi (\beta ,L)`$ $``$ $`{\displaystyle \underset{x}{}}G(x;\beta ,L)`$ (A5)
We also consider the second-moment correlation length $`\xi `$. In infinite volume it is defined as

$$\xi _{\mathrm{}}^2(\beta )\frac{1}{6\chi _{\mathrm{}}(\beta )}\underset{x}{}|x|^2G_{\mathrm{}}(x;\beta ).$$
(A6)
The finite-volume generalization is by no means unique. We use

$$\xi ^2(\beta ,L)\frac{\widehat{G}(0;\beta ,L)\widehat{G}(q_{\mathrm{min}};\beta ,L)}{\widehat{q}_{\mathrm{min}}^2\widehat{G}(q_{\mathrm{min}};\beta ,L)},$$
(A7)
where $`q_{\mathrm{min}}(2\pi /L,0,0)`$, $`\widehat{q}2\mathrm{sin}q/2`$, and $`\widehat{G}(q;\beta ,L)`$ is the Fourier transform of $`G(x;\beta ,L)`$. This finite-volume definition has the correct infinite-volume limit and shows a fast convergence as $`L\mathrm{}`$ .

We also define the energy $`E(\beta ,L)`$ and the specific heat $`C(\beta ,L)`$:

$`E(\beta ,L)`$ $``$ $`3G(e;\beta ,L),`$ (A8)
$`C(\beta ,L)`$ $``$ $`{\displaystyle \frac{E(\beta ,L)}{\beta }},`$ (A9)
where $`e=(1,0,0)`$. We also compute higher-order couplings. Setting

$$\mu _k(\underset{i}{}\rho _is_i)^k,m_{k_1k_2\mathrm{}k_n}\overline{\mu _{k_1}\mu _{k_2}\mathrm{}\mu _{k_n}},$$
(A10)
we define the connected $`n`$-point susceptibilities $`\chi _n`$ averaged over random dilution by

$`V\chi _2(\beta ,L)V\chi (\beta ,L)=m_2,`$ (A11)
$`V\chi _4(\beta ,L)m_43m_{22},`$ (A12)
$`V\chi _6(\beta ,L)m_615m_{42}+30m_{222},`$ (A13)
where $`VL^3`$ is the volume. Moreover, we also define

$`V\chi _{22}(\beta ,L)m_{22}m_2^2,`$ (A14)
$`V\chi _{42}(\beta ,L)m_{42}m_4m_23m_{222}+3m_{22}m_2,`$ (A15)
$`V\chi _{222}(\beta ,L)m_{222}3m_{22}m_2+2m_2^3.`$ (A16)
Then, we define the four-point couplings

$`G_4(\beta ,L){\displaystyle \frac{\chi _4}{\xi ^3\chi _2^2}},`$ (A17)
$`G_{22}(\beta ,L){\displaystyle \frac{\chi _{22}}{\xi ^3\chi _2^2}},`$ (A18)
$`H_4(\beta ,L)G_4+3G_{22},`$ (A19)
and the six-point universal ratios

$`r_6(\beta ,L)10{\displaystyle \frac{\chi _6\chi _2}{\chi _4^2}},`$ (A20)
$`C_{42}(\beta ,L)4{\displaystyle \frac{\chi _{42}\chi _2}{\chi _4\chi _{22}}},`$ (A21)
$`C_{222}(\beta ,L)6{\displaystyle \frac{\chi _{222}\chi _2}{\chi _{22}^2}}.`$ (A22)
We will be interested in the critical value of these quantities. If $`S(\beta ,L)`$ is any of them, we compute its fixed-point value (note that the order of the limits cannot be interchanged)

$$S^{}=\underset{\beta \beta _c}{lim}\underset{L\mathrm{}}{lim}S(\beta ,L).$$
(A23)
Finally, we define the Binder parameters

$`U_n(\beta ,L){\displaystyle \frac{m_n}{m_2^{n/2}}},`$ (A24)
$`U_{22}(\beta ,L){\displaystyle \frac{m_{22}}{m_2^2}},`$ (A25)
and the corresponding critical-point values

$$U^{}=\underset{L\mathrm{}}{lim}\underset{\beta \beta _c}{lim}U(\beta ,L).$$
(A26)
Note that the order of limits is reversed with respect to Eq. (A23).