Multiscaling in the 3D critical site-diluted Ising ferromagnet

We study numerically the appearance of multiscaling behavior in the 3D ferromagnetic Ising site-diluted model, in the form of a multifractal distribution of the decay exponents for the spatial correlation functions at the critical temperature. We have computed the exponents of the long-distance decay of higher moments of the correlation function, up to the 10th power, by studying three different quantities: global susceptibilities, local susceptibilities and correlation functions. We have found very clear evidence of multiscaling behavior.


Introduction
The extension of Statistical Mechanics to the study of systems with quenched disorder is of paramount importance.Many efforts have been devoted to this task in the last decades, and many important results have been established.Universality is here the main key.Critical exponents can take, for pure systems, few values, that only depend on crucial features of the system, like its dimensionality, the nature of its fields and the symmetries of the Hamiltonian.How does this idea extends to system where the Hamiltonian is characterized by quenched disorder has been discussed at length and is still, at least partially, an open problem.
The Diluted ferromagnetic Ising Model (DIM in the following) is one of the crucial prototypes on which the entrance of quenched disorder in Statistical Mechanics is analyzed.
Here even an infinitesimal amount of dilution (that can be a site or a link dilution, leading to models showing basically the same physical behavior) on a pure model, with a positive specific heat critical exponent, will change universality class and the critical exponents [1,2].
Here in order to further clarify crucial properties of the DIM we will be interested in the behavior of correlation functions and of susceptibilities.In a second order phase transition a correlation length ξ diverges at the critical point T c when the volume V of the system diverges.When the position of the critical point is known we can analyze correlation functions C there, where they are expected to scale as a power law of the distance, with an exponent, say, τ .An usual scaling law would suggest that the q-th moment of a correlation function C(r) would scale with an exponent τ (q) equal to qτ (1) (see later for details and precise definitions).
We will find here that in a three dimensional, D = 3 Diluted Ising Model, this is not true, and we get what is called a multi-fractal behavior [3,4].We will determine with high accuracy the values of these exponents, and show numerically that they do not obey the relation τ (q) = qτ (1).
This multiscaling behavior of the correlation functions has been thoroughly studied in the context of two-dimensional ferromagnetic spin models with quenched disorder (Potts models with more than two states and site or link dilution) using conformal field techniques and numerical simulations [5][6][7][8][9][10][11].In another context, we can also quote the multiscaling behavior which appears in the strong amplitude fluctuations of the wave functions in the Anderson localization [12].
To find such a results in a 3D "simple" disordered system (simple as opposed to spin glasses, where the same fact has been recently established numerically for all values of the temperature in the broken phase [13]) is in some sense unexpected, but indeed some important theoretical results by Davis and Cardy had clarified, already many years ago, the plausibility of this fact [14,15].Conformal invariance allows crucial developments about critical scaling, starting from two dimensional systems [5].By applying to this context renormalization group and an ϵ expansion Davis and Cardy were able to show that in the DIM one expects multifractality in 2 + ϵ dimensions [15].
Our numerical findings, that we support with theoretical arguments, are consistent with Davis and Cardy's computation [15] (see also A.3).As many results in the field explain, there are indeed many factors that have to be considered: logarithmic corrections also appear in these kinds of contexts, and we have been able to use numerics to clearly show that in this case we are observing a real, bona fide multifractality.We will discuss this in detail in the following.
The structure of the paper is the following.We shall first define our model (Sect.2.1), describe the numerical approach we follow (Sect.2.2), define the observable quantities that we compute (Sect 2.3) and then explain how we compute these quantities (Sect 2.4).Next, we shall show our numerical results in Sect.3. Finally, we shall discuss our results and draw some conclusions in Sect. 4. The paper is complemented with an Appendix, where we remind the reader of relevant theoretical results.

The model
We have considered the D = 3 diluted Ising model defined on a cubic lattice with periodic boundary conditions, linear size L and volume V = L D .The Hamiltonian of the model, in zero magnetic field, has the form where s x are Ising variables and ϵ x are statistically independent quenched random variables, which take the value 1 with probability p and the value 0 with probability 1 − p.A set of {ϵ x } defines a sample.The sum in Eq. ( 1) runs over all pairs of lattice nearest-neighbors.As usual, we denote with ⟨(• • •)⟩ the thermal average for the {s x }, that is computed for a fixed set of disorder variables {ϵ x }.The average over the samples, (• • •), is only taken after thermal mean values have been computed for each sample.

Our numerical simulations
We have investigated the model defined in Eq. (1) through equilibrium numerical simulations.Specifically, we have brought to thermal equilibrium our samples by combining Wolff' singlecluster algorithm [16], with a local Metropolis method.The remarkable effectiveness of this combination of simulations algorithms in DIM simulations was demonstrated long ago [17,18].Specifically, our elementary Monte Carlo steps consisted on L single-cluster updates, followed by a sequential full-lattice Metropolis update.
In particular, we have fixed the (average) spin density to p = 0.8, because for this value the scaling corrections are very small, even negligible given our statistical errors: the model is almost governed by a perfect action [18,19].Furthermore, we have performed all our simulations at the infinite volume critical inverse temperature [19]: (4).For future use, we note that the anomalous dimension of the field is η = 0.036(1) [19].Further details about our simulations can be found in Table 1.

Definition of main observables
In the next three paragraphs we describe the different observables used to analyzed the multiscaling properties of our model (1).

Correlations
The relevant correlation functions are here As usual, we are assuming that rotational invariance is recovered in the scaling limit ξ ≫ 1, and we indicate only the dependence of C q (r) on the length r of the displacement vector r.
Finally, in order to compute the ζ exponent, minimizing the scaling corrections, we will analyze the ratio

Global susceptibilities
We define the χ q susceptibilities as which scale with L as provided that D > τ (q).The above relation allows us to to compute τ (q) and, from it, ζ(q).If D < τ (q) we have that Equations ( 8) and ( 9) provide a very direct test for multiscaling.Indeed, given that τ (1) = 1.036(1) in three dimensions, one notices that qτ (1) > D whenever q ≥ 3. Hence finding (as we shall do below) that χ q=3 scales with a positive power of L gives a direct confirmation of multiscaling behavior.

Local susceptibilities
Another strategy is based on computing a different "local susceptibility" χq defined as with The scaling of this observable, see A.1, is such that Notice that χ1 = χ 1 .
We have been able to compute numerically the first two derivatives of ζ(q) by computing (q ≥ 2) and (for the discrete proxy of the second derivative)

Details about the computation of the observables
Since the computation of the correlation function C q (r) technically differs from the computation of the global susceptibility χ q (that, on its side, is also very different from the computation of the local susceptibility χq ), we have chosen to discuss these points in separate paragraphs.
For the computation of χ q and χq we shall be employing real replicas for every sample.Real replicas are statistically independent system copies that evolve under the same quenched disorder {ϵ x }.
Furthermore, our computation of the two susceptibilities uses the enhanced cluster estimator [20] where γ y,x = 1 if, and only if: • When the lattices is decomposed on Fortuin-Kasteleyn clusters, both sites x and y happen to belong to the same cluster.
Our estimator of χ q is obtained in the following way.We start by choosing randomly a starting site x and trace the cluster to which x belongs in all the 12 replicas.Next, we select a set of q distinct replicas a 1 , a 2 , . . ., a q , and compute Defining the number of spins as one easily finds that Of course, in order to improve the statistical accuracy, one can average over different choices of the subset of q distinct replica indices.Note that in practice we choose just one starting site x in Eq. ( 18).This starting site is picked with uniform probability 1/N spins .Now, given that and that the fluctuations of N spins do not play any relevant role, ‡ we may approximate A final note of warning is in order.We cannot employ in the update of the spin configuration the same Fortuin-Kasteleyn clusters that we employ in the computation of the susceptibilities.Indeed, our twelve real-replicas are to remain statistically independent at all times.Hence we cannot start the single-cluster update from the same site x in all the replicas.This is why we separate our simulations in two different phases: (i) During the update phase, we chose independently the starting site x for the single-cluster update of each replica.The spin-configuration is updated after the cluster is traced.
(ii) During the measuring phase, instead, the starting point for the cluster tracing is the same in all replicas.However, the spin configurations of the replicas are not modified after the M (q) x estimators are computed.In both phases, the starting site for cluster tracing is chosen to be occupied (namely, ϵ x = 1).‡ More quantitatively, trading 1/pV with 1/N spins only induces additional corrections to scaling, of size L −a with a = D − 1 ν ≥ D/2.This estimate can be obtained by combining the following three observations: (i) 2.4.2.The computation of χq Our computation of the local susceptibilities χq employs 16 real replicas.The separation between update phase and measuring phase that we have discussed above also applies to this case.
During the measuring phase, we compute for each replica We proceed by selecting a subset of q distinct replica indices a 1 , a 2 , . . ., a q , and compute the estimator which has the thermal expectation value where the local susceptibility χ x was defined in Eq. (11).Of course, one may average over the different choices of the subset of q distinct replicas in order to increase the statistics.
As in the previous paragraph, we pick just one starting point x during the measuring phase, with uniform probability 1/N spins .Hence The above equation is not an equality because of the fluctuations on the number of spins N spins .
2.4.3.The computation of C q (r) As opposed to the strategy we have followed in the computation of the susceptibilities based in the simulation of a high number of replicas (in the same realization of the disorder), in the computation of the correlation function C q we have chosen to simulate only one replica.This greatly simplifies the simulation but induces strong bias in the statistical estimator of C q .In particular, it is possible to show [17] that this bias is proportional to 1/N m (N m being the number of measurements on a given sample), and the problems arise when the magnitude of this bias is similar to the statistical error induced by the disorder (proportional to 1/ √ N S (N S is the number of disorder realizations, samples).To cure this bias, we have followed the strategy introduced in Ref. [17].Let us work on a given sample (i.e.fixed disorder) and at equilibrium.Our strategy consists, firstly, in computing the total average of a given observable O, denoted as O (1) .Then we take two halves of the Monte Carlo simulation, and compute the average value of O in these two halves and compute their mean, denoted as O (2) .Finally, obtain the average of O on each quarter of the Monte Carlo history and compute the mean of them, getting O (3) .The final, unbiased, quadratic estimator is [17] Finally, in Table 1 we report the number of samples used in the different runs.

Numerical Results
In this section we present our numerical results regarding multiscaling by analyzing the behavior of the correlation functions and the global and local susceptibilities.
Table 1: Details of our numerical simulations.We report the number of samples (disorder realizations), N S , used in the different runs (χ q susceptibilities, χq local susceptibilities and correlation functions).We also report the number of sweeps (as defined in the text) used for thermalizing the different systems, N sweeps .In the χ q susceptibility runs we have simulated 12 replicas per sample while for computing the local susceptibilities we have simulated 16 replicas per sample.In the correlation (Corr.)runs we have simulated only one replica per sample.

Correlation functions
To quantify the ζ(q) exponent we have analyzed the dependence of C q on C 1 .We show in Fig. 1 the behavior of C q as a function of C 1 for three different lattice sizes and three values of q.
The scaling of the data (for L = 32, 64 and 128) is very good and the power law dependence of C q on C 1 is very clear and accurate.In Fig. 2 we plot C q /C q 1 as a function of C 1 in order to asses the multiscaling behavior.In all the simulated cases (q ≤ 10) we have found a diverging behavior as C 1 → 0 which marks the onset of multiscaling behavior in this model.Notice again the good scaling of the data from different lattice sizes (L = 32, 48 and 64).
We can estimate the ζ(q) exponent by fitting the correlation function data to Eq. ( 6).We fit C q (r) computed on a given lattice size against C 1 (r) computed on the same lattice size.We have computed the error bars of the ζ(q) using the jackknife method over the number of samples, following the recipe of Ref. [22]: we perform independent fits on all the jackknife blocks using the diagonal covariance matrix, and then, the error bars in the parameters of the fit can be computed from the fluctuations among all the jackknife blocks.
We report our results in Table 2.The exponent ζ(q) is not proportional to q.This clearly shows that the three-dimensional diluted Ising model presents multiscaling behavior.Despite the fact that we are using a quasi-perfect action, we can observe a small dependence of ζ(q) on the lattice size.

Global Susceptibilities
Figure 3 shows the behavior of the χ q susceptibilities as a function of the lattice size.
By fitting the susceptibilities presented in Fig. 3 we have estimated the τ (q) exponents using Eq. ( 8).The numerical estimates for the τ (q) exponents (and consequently those for the ζ(q) exponents) can be found in Table 3.  Table 2: Exponent ζ(q) computed using the L = 32, L = 48 and L = 64 correlation functions.This exponent has been computed using the jackknife method to tackle the highly correlation of the data (different r in the correlations C q (r)).We also report the C 1 -interval used in the fits, (C m 1 , C M 1 ).
Figure 2: C q (r)/C q 1 (r) versus C 1 (r) for q = 2, 5 and 10, and for L = 32 (blue triangles), L = 48 (red circles) and L = 64 (black squares).Notice the potential growth in all the cases which is a clear signal of the multiscaling behavior.All the curves take the value 1 as C 1 = 1 (the rightmost point of the plot).The error bars of the ratio have been computed using the jackknife method.The size of the error bars is smaller than the ones of the symbols.Table 3: τ (q) and ζ(q) exponents from the finite size scaling of the susceptibilities χ q .For χ 4 we do not detect divergence with L, so for q ≥ 4, τ (q) ≥ 3 and θ(q) ≥ 2.89 (1).We have computed, as a test, the 1 + η exponent in the q = 1 row, and it compares very well with the most accurate available estimate 1 + η = 1.036(1) [19].

Local susceptibilities
In Fig. 4 we show the behavior of R (1) q versus L for three different values of q, showing the power law divergence implied by Eq. (12).
derivatives obtained analyzing the observables R q and R q , computing their error bars with the bootstrap method, and using only the lattice sizes with L ≥ 48 (as for the global susceptibility analysis) to avoid the small scaling corrections §.We plot our results for ζ(q) in Fig. 5 where the strong departure of the linear regime is clear, providing strong indications of multi-fractal behavior.
Furthermore, the results reported in Table 4 show that the function ζ(q) is a non-decreasing function for the values of q simulated.In addition, it is a concave function of q, in agreement with the results from Ludwig [5] (see also A.2).

Discussion and Conclusions
Along the years, many disordered ferromagnetic systems in two dimensions have been shown to exhibit multiscaling behavior (see e.g.Ref. [5]), with the surprising exception of the DIM.This is however a marginal behavior.The DIM is the Q = 2 instance of the diluted Potts model with Q states (which does display multiscaling for Q > 2).Furthermore, it suffices to consider § An analysis using all the sizes and assuming the leading scaling correction term gives essentially the same exponents.
Table 4: ζ(q) exponents and the (discrete) first and second derivatives from the finite size scaling of the ratios R q versus L for q = 2, 4, and 7 in a double logarithm scale.The straight lines are fits to Eq. ( 12) considering scaling corrections.The error bars of all three ratios R (i) q have been computed using the bootstrap method.long range interactions to find multiscaling in the two-dimensional DIM (see [23] and references therein).The behavior of the DIM (with short range interactions) in D = 2 is marginal in a different way, as well.Although there is no multiscaling in D = 2, in a pioneering work, Davis and Cardy [15] found multifractality in the DIM in space dimensions D = 2 + ϵ.Here, we have extended Davis and Cardy's [15] result to D = 3 through extensive numerical simulations for ζ(q) q Figure 5: The exponent ζ(q) versus q (see results of Table 4).By definition ζ(1) = 1.We have also plotted the line ζ = q to show the appearance of the multi-fractal behavior in the model (ζ(q) ̸ = q).
the site-diluted Ising model, that is studied under equilibrium conditions.We have computed the exponents τ (q) for integer values q = 1, 2, . . ., 10.The q-th power of the correlation function, Eq. ( 2), decays at long distances as 1/r τ (q) (3).We have found a very clear numerical evidence for τ (q) < qτ (1) for q > 1, which is the hallmark of multifractality.
Let us emphasize that multiscaling in no way contradicts the general Renormalization Group expectation of a universal large-L limit for relative cumulants of the order-parameter [24].Specifically, for the very same model studied here, it was shown long ago that the relative variance of the squared order parameter, reaches its expected universal limit [18] (this behavior is sometimes named no selfaveraging [24]).One could think that having a finite large-L limit for g 2 does somehow contradict multiscaling.The correct conclusion is indeed different, as we now explain starting from the results discussed in Appendix A.1.Consider the second moment, q = 2, which is one relevant in the computation of g 2 .The key point regards the spatial sites one integrates over: • If one integrates ϵ x ϵ y ⟨S x S y ⟩ 2 over y the global susceptibility χ q=2 is obtained and the full multifractal signal is recovered.In summary, the extreme self-averaging violations evinced by the multiscaling analysis can only be identified when studying correlations at the local level.Carrying out too many spatial integrations erases the multifractal signal.In fact, to unearth multiscaling we have pursued here the same local approach that has been used in Ref. [13] in the context of the out-equilibrium dynamics of spin-glasses ¶.
As discussed in Ref. [5], the composite operator (ϕ a 1 (x) • • • ϕ aq (x)) transforms following a reducible representation of the symmetric group (S n ) and, therefore, it is not a scaling operator at the random critical point.However, it can be expressed as a linear combination of scaling fields Φ µ,α q (x), and each of them transforms following the µirreducible representation (irrep) of S n (α labels the multiplicity of the µ-irrep) [5], with scaling dimension X µ,α q .Its asymptotic behavior is ⟨Φ µ,α q (x)Φ µ,α q (0)⟩ ∼ 1 |x| 2X µ,α q . (A.2) Therefore, the behavior of C q will be a sum (over the irreducible representations) of decays with powers 1/|x| 2X µ,α q , and the smallest scaling dimension X µ,α q will determine its large scale behavior.
It is possible to show that for the spin operator (which is our case), all the different representations become degenerate as opposed to the energy operator [5,15].Therefore, [(D − 2 + η)/2]ζ(q) is the scaling dimension of the full composite operator (ϕ a 1 (x) • • • ϕ aq (x)), which defines ζ(q) in such a way that ζ(1) = 1.
Equation (A.3) can be written in the continuum as where the q + 1 integrals run in the volume [a, L] D (a is the lattice spacing).After the change of variables x = x/L and ỹi = y i /L (i = 1, . . ., q) we can write ∼ L (qD−[(D−2+η)/2](ζ(q)+q) , since the q + 1 integrals in the tilde variables provided a constant in the large L-limit: their integration limits are constrained to the volume [a/L, 1] D which is [0, 1] D for large L.
Taking into account this result, it is easy to get the behaviors given in Eqs. ( 12), ( 13) and ( 14).

Figure 1 :
Figure1: C q (r) versus C 1 (r) for q = 2, 5 and 10, and for L = 32 (blue triangles), L = 48 (red circles) and L = 64 (black squares).The size of the error bars is smaller than the ones of the symbols.

•
If one integrates instead ϵ x ϵ y 1 ϵ y 2 ⟨S x S y 1 ⟩⟨S x S y 2 ⟩ over y 1 and y 2 the resulting quantity is the local susceptibility χq=2 which has only half the multifractal signal.∥• Yet, the needed quantity in the computation of g 2 is ⟨M 2 ⟩ 2 .Hence, we are to integrate ϵ x 1 ϵ x 2 ϵ y 1 ϵ y 2 ⟨S x 1 S y 1 ⟩⟨S x 1 S y 2 ⟩ over x 1 , x 2 , y 1 and y 2 .But the integrals over x 2 , y 1 and y 2 already suffice to completely suppress multifractal scaling.In fact, neglecting scaling corrections, ⟨M 2 ⟩ 2 scales with L in the same way that ⟨M 2 ⟩ 2 does.