Multiscaling in the 3D critical site-diluted Ising ferromagnet

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,

$$\begin{align} {\cal H} = -\sum_{\langle\boldsymbol{x},\boldsymbol{y}\rangle} \epsilon_{\boldsymbol{x}} \epsilon_{\boldsymbol{y}} s_{\boldsymbol{x}} s_{\boldsymbol{y}}\,, \end{align} \tag{ 1 }$$

where $s_{\boldsymbol{x}}$ are Ising variables and $\epsilon_{\boldsymbol{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 $\{\epsilon_{\boldsymbol{x}}\}$ defines a sample. The sum in equation (1) runs over all pairs of lattice nearest-neighbors. As usual, we denote with $\langle (\cdots)\rangle$ the thermal average for the $\{s_{\boldsymbol{x}}\}$, which is computed for a fixed set of disorder variables $\{\epsilon_{\boldsymbol{x}}\}$. The average over the samples, $\overline{(\cdots)}$, is only taken after thermal mean values have been computed for each sample.

We have investigated the model defined in equation (1) through equilibrium numerical simulations. Specifically, we have brought our samples to thermal equilibrium by combining Wolff's single-cluster algorithm [16], with a local Metropolis method. The remarkable effectiveness of this combination of simulation 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]: $\beta_c^{(p = 0.8)} = 0.285\,7429(4)$. For future use, we note that the anomalous dimension of the field is $\eta = 0.036(1)$ [19]. Further details of our simulations can be found in table 1.

Table 1. Details of our numerical simulations. We report the number of samples (disorder realizations), $N_\mathrm{S}$, used in the different runs (χ_q susceptibilities, $\tilde{\chi}_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_\mathrm{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.

L	$N_\mathrm{S}$ (χq )	$N_\mathrm{S}$ ($\tilde{\chi}_q$)	$N_\mathrm{S}$ (Corr.)	$N_\mathrm{sweeps}$
8		110 592		16
12	98 304	108 902		32
16	98 304	109 666		32
24	98 304	109 428		64
32	81 920	796 367	5000	64
48	81 920	109 528	5030	128
64	81 920	84 774	9800	128
96	81 920	84 034		512
128	26 984	60 015		512

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

2.3.1. Correlations.

The relevant correlation functions are,

$$\begin{align} C_q\left(r\right) = \frac{1}{p V} \sum_{\boldsymbol{x}} \overline{\langle \epsilon_{\boldsymbol{x}} s_{\boldsymbol{x}} \epsilon_{\boldsymbol{x+r}} s_{\boldsymbol{x+r}}\rangle^q} \,. \end{align} \tag{ 2 }$$

As usual, we assume that rotational invariance is recovered in the scaling limit $\xi\gg 1$, and we indicate only the dependence of $C_q(r)$ on the length r of the displacement vector r .

In order to study multiscaling behavior, we introduce the $\zeta(q)$ and $\tau(q)$ exponents through the following relations:

$$\begin{align} C_q\left(r\right) \sim \frac{1}{r^{\tau\left(q\right)}} \sim \left[C_1\left(r\right)\right]^{\zeta\left(q\right)}\,. \end{align} \tag{ 3 }$$

It follows that $\tau(q)$ and $\zeta(q)$ are connected by the relation,

$$\begin{align} \tau\left(q\right) = \left(D-2+\eta\right) \zeta\left(q\right)\, ,\end{align} \tag{ 4 }$$

because,

$$\begin{align} C_1\left(r\right) \sim \frac{1}{r^{D-2+\eta}}\,. \end{align} \tag{ 5 }$$

By definition, $\zeta(1) = 1$ and $\tau(1) = D-2+\eta$. In the absence of multiscaling behavior one would have $\tau(q) = (D-2+\eta) q$ or, equivalently, $\zeta(q) = q$ (see A.1 for more details). Instead, we find that $\zeta(q)\lt q$ whenever q > 1.

Finally, in order to compute the ζ exponent, minimizing the scaling corrections, we analyze the ratio:

$$\begin{align} \frac{C_q}{C_{1}^q} \sim \left[C_{1}\right]^{\zeta\left(q\right)-q}\,.\end{align} \tag{ 6 }$$

2.3.2. Global susceptibilities.

We define the χ_q susceptibilities as,

$$\begin{align} \chi_q = \frac{1}{p V} \sum_{\boldsymbol{x} \boldsymbol{y}} \overline{\langle \epsilon_{\boldsymbol{x}} s_{\boldsymbol{x}} \epsilon_{\boldsymbol{y}} s_{\boldsymbol{y}}\rangle^q}\,, \end{align} \tag{ 7 }$$

which scale with L as,

$$\begin{align} \chi_q \sim \int^L{\textrm{d}}^D x ~C_q\left(r\right) \sim L^{D-\tau\left(q\right)}\,, \end{align} \tag{ 8 }$$

provided that $D\gt\tau(q)$. The above relation allows us to compute $\tau(q)$ and, from it, $\zeta(q)$. If $D\lt\tau(q),$ we have that,

$$\begin{align} \chi_q \sim L^0\,\,\,. \end{align} \tag{ 9 }$$

Equations (8) and (9) provide a very direct test for multiscaling. Indeed, given that $\tau(1) = 1.036(1)$ in 3D, one notes that $q\tau(1) \gt D$ whenever $q\unicode{x2A7E} 3$. Hence, finding (as we do below) that $\chi_{q = 3}$ scales with a positive power of L gives a direct confirmation of the multiscaling behavior.

2.3.3. Local susceptibilities.

Another strategy is based on computing a different 'local susceptibility' $\tilde{\chi}_q$ defined as,

$$\begin{align} \tilde{\chi}_q = \frac{1}{p V} \sum_{\boldsymbol{x}} \overline{\chi_{\boldsymbol{x}}^q}\,, \end{align} \tag{ 10 }$$

with,

$$\begin{align} \chi_{\boldsymbol{x}} = \sum_{\boldsymbol{y}} \langle \epsilon_{\boldsymbol{x}} s_{\boldsymbol{x}} \epsilon_{\boldsymbol{y}} s_{\boldsymbol{y}} \rangle\,. \end{align} \tag{ 11 }$$

The scaling of this observable, see A.1, is such that,

$$\begin{align} R^{\left(1\right)}_q \equiv \frac{\tilde{\chi}_q}{\tilde{\chi}_1^q}\sim L^{\left(D-2 +\eta\right)\left(q-\zeta\left(q\right)\right)/2}\,. \end{align} \tag{ 12 }$$

Note that $\tilde{\chi}_1 = \chi_1$.

We are able to compute numerically the first two derivatives of $\zeta(q)$ by computing ($q\unicode{x2A7E} 2$):

$$\begin{align} R^{\left(2\right)}_q\equiv \frac{\tilde{\chi}_q}{\tilde{\chi}_{q-1}} \sim L^{\left(D-2 +\eta\right)\left(\zeta\left(q-1\right)-\zeta\left(q\right)\right)/2}\end{align} \tag{ 13 }$$

and (for the discrete proxy of the second derivative),

$$\begin{align} R^{\left(3\right)}_q\equiv \frac{\tilde{\chi}_q \tilde{\chi}_{q-2} }{\tilde{\chi}_{q-1}^2} \sim L^{\left(D-2 +\eta\right)\left(2 \zeta\left(q-1\right)-\zeta\left(q\right)-\zeta\left(q-2\right)\right)/2}\ \,. \end{align} \tag{ 14 }$$

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 to the computation of the local susceptibility $\tilde{\chi_q}$), we have chosen to discuss these points in separate sections.

For the computation of χ_q and $\tilde{\chi}_q$ we shall be employing real replicas for each sample. Real replicas are statistically independent system copies that evolve under the same quenched disorder $\{\epsilon_{\boldsymbol{x}}\}$.

Furthermore, our computation of the two susceptibilities uses the enhanced cluster estimator [20]:

$$\begin{align} \langle \epsilon_{\boldsymbol{x}} s_{\boldsymbol{x}} \epsilon_{\boldsymbol{y}} s_{\boldsymbol{y}}\rangle = \langle \gamma_{\boldsymbol{y},\boldsymbol{x}}\rangle\,, \end{align} \tag{ 15 }$$

where $\gamma_{\boldsymbol{y},\boldsymbol{x}} = 1$ if, and only if,

• $\epsilon_{\boldsymbol{x}} = \epsilon_{\boldsymbol{y}} = 1$,
• When the lattice is decomposed on Fortuin–Kasteleyn clusters, both sites x and y happen to belong to the same cluster.

Otherwise, $\gamma_{\boldsymbol{y},\boldsymbol{x}} = 0$.

2.4.1. The computation of χ_q .

Our computation of the global susceptibilities χ_q uses 12 real replicas. Let $\gamma_{\boldsymbol{y},\boldsymbol{x}}^{(a)}$ be the enhanced estimator for correlation functions obtained from replica a, $(a = 1,2,\ldots,12)$.

Our estimator of χ_q is obtained in the following way. We start by randomly choosing a starting site x and trace the cluster to which x belongs in all 12 replicas. Next, we select a set of q distinct replicas $a_1, a_2, \ldots,a_q$, and compute,

$$\begin{align} M_{\boldsymbol{x}}^{\left(q\right)} = \sum_{\boldsymbol{y}}\, \left(\prod_{j = 1}^q \gamma_{\boldsymbol{y},\boldsymbol{x}}^{\left(a_j\right)}\right)\,. \end{align} \tag{ 16 }$$

Defining the number of spins as,

$$\begin{align} N_{\mathrm{spins}} = \sum_{\boldsymbol{x}} \epsilon_{\boldsymbol{x}}\,, \end{align} \tag{ 17 }$$

one easily finds that,

$$\begin{align} \frac{1}{N_\mathrm{spins}} \sum_{\boldsymbol{x}} \langle M_{\boldsymbol{x}}^{\left(q\right)}\rangle = \frac{1}{N_\mathrm{spins}} \sum_{\boldsymbol{x} \boldsymbol{y}} \overline{\langle \epsilon_{\boldsymbol{x}} s_{\boldsymbol{x}} \epsilon_{\boldsymbol{y}} s_{\boldsymbol{y}}\rangle^q}\,. \end{align} \tag{ 18 }$$

Of course, in order to improve the statistical accuracy, one can average over different choices of a subset of q distinct replica indices. Note that in practice we choose just one starting site x in equation (18). The starting site is chosen with uniform probability $1/N_\mathrm{spins}$.

Now, given that,

$$\begin{align} \overline{N_\mathrm{spins}} = p V \,, \end{align} \tag{ 19 }$$

and that the fluctuations of $N_\mathrm{spins}$ do not play any relevant role ⁴ , we may approximate,

$$\begin{align} \chi_q\approx \overline{\langle M_{\boldsymbol{x}}^{\left(q\right)}\rangle}\,. \end{align} \tag{ 20 }$$

A final note of warning is in order. In the update of the spin configuration, we cannot employ the same Fortuin–Kasteleyn clusters that we employ in the computation of the susceptibilities. Indeed, our 12 real replicas are to remain statistically independent at all times. Hence, we cannot start a single-cluster update from the same site x in all replicas. Here, we separate our simulations into two different phases:

• (i)  
  During the update phase, we choose 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 for all replicas. However, the spin configurations of the replicas are not modified after the $M_{\boldsymbol{x}}^{(q)}$ estimators are computed.

In both phases, the starting site for cluster tracing is chosen to be occupied (namely, $\epsilon_{\boldsymbol{x}} = 1$).

2.4.2. The computation of $\tilde{\chi}_q$.

Our computation of the local susceptibilities $\tilde{\chi}_q$ employs 16 real replicas. The separation between the update phase and the measuring phase, which we have discussed above also applies to this case.

During the measuring phase, we compute for each replica:

$$\begin{align} N_{\boldsymbol{x}}^{\left(a\right)} = \sum_{\boldsymbol{y}} \gamma_{\boldsymbol{y},\boldsymbol{x}}^{\left(a\right)}\,. \end{align} \tag{ 21 }$$

We proceed by selecting a subset of q distinct replica indices $a_1, a_2, \ldots,a_q$, and compute the estimator:

$$\begin{align} Y_{\boldsymbol{x},q} = \prod_{j = 1}^q N_{\boldsymbol{x}}^{\left(a_j\right)}\,, \end{align} \tag{ 22 }$$

which has the thermal expectation value:

$$\begin{align} \langle Y_{\boldsymbol{x},q}\rangle = \chi_{\boldsymbol{x}}^q\,, \end{align} \tag{ 23 }$$

where the local susceptibility $\chi_{\boldsymbol{x}}$ is defined in equation (11). Of course, one may average over different choices of a subset of q distinct replicas in order to increase the statistics.

As in the previous section, we choose just one starting point x during the measuring phase, with uniform probability $1/N_{\mathrm{spins}}$. Hence,

$$\begin{align} \tilde{\chi}_q\approx \overline{\langle Y_{\boldsymbol{x},q}\rangle}\, . \end{align} \tag{ 24 }$$

The above equation is not an equality due to the fluctuations on the number of spins $N_{\mathrm{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 on 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_\mathrm{m}$ ($N_\mathrm{m}$ being the number of measurements on a given sample), and problems arise when the magnitude of this bias is similar to the statistical error induced by the disorder (proportional to $1/\sqrt{N_\mathrm{S}}$ ($N_\mathrm{S}$ is the number of disorder realizations, samples).

To resolve this bias, we have followed the strategy introduced in [17]. Let us work on a given sample (i.e. fixed disorder) and at equilibrium. Our strategy consists, first, of computing the total average of a given observable O, denoted as $O^{(1)}$. Then, we use 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, we obtain the average of O on each quarter of the Monte Carlo history and compute their mean, obtaining $O^{(3)}$. Finally, the unbiased quadratic estimator is [17],

$$\begin{align} O^Q = \frac{8}{3} \, O^{\left(1\right)}-2 \, O^{\left(2\right)}+\frac{1}{3} \, O^{\left(3\right)}\,. \end{align} \tag{ 25 }$$

In table 1, we report the number of samples used in the different runs.

To quantify the $\zeta(q)$ exponent we have analyzed the dependence of C_q on C₁. In figure 1, we show the behavior of C_q as a function of C₁ 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₁ is very clear and accurate.

Zoom In Zoom Out Reset image size

In figure 2, we plot $C_q/C_1^q$ as a function of C₁ in order to assess the multiscaling behavior. In all the simulated cases ($q\unicode{x2A7D} 10$) we have found diverging behavior as $C_1 \to 0$, which marks the onset of multiscaling behavior in this model. Note again the good scaling of the data from different lattice sizes (L = 32, 48 and 64).

Zoom In Zoom Out Reset image size

We estimate the $\zeta(q)$ exponent by fitting the correlation function data to equation (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 $\zeta(q)$ using the jackknife method over the number of samples, following the method of [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 $\zeta(q)$ is not proportional to q. This clearly shows that the 3D diluted Ising model presents multiscaling behavior. Despite the fact that we are using a quasi-perfect action, we observe a small dependence of $\zeta(q)$ on the lattice size.

Table 2. Exponent $\zeta(q)$ computed using the L = 32, L = 48 and L = 64 correlation functions. This exponent has been computed using the jackknife method to address the high correlation of the data (different r in the correlations $C_q(r)$). We also report the C₁ interval used in the fits, $(C_1^{\mathrm{m}},C_1^{\mathrm{M}})$.

	L = 64		L = 48		L = 32
q	$(C_1^{\mathrm{m}},C_1^{\mathrm{M}})$	$\zeta(q)$	$(C_1^{\mathrm{m}},C_1^{\mathrm{M}})$	$\zeta(q)$	$(C_1^{\mathrm{m}},C_1^{\mathrm{M}})$	$\zeta(q)$
2	(0,0.3)	1.8109(3)	(0,0.3)	1.8094(7)	(0,0.3)	1.808(1)
3	(0,0.17)	2.457(1)	(0,0.3)	2.4525(19)	(0,0.3)	2.447(3)
4	(0,0.1)	3.026(4)	(0,0.16)	2.997(6)	(0,0.3)	2.993(8)
5	(0,0.1)	3.486(7)	(0,0.1)	3.474(13)	(0,0.3)	3.432(12)
6	(0,0.1)	3.849(9)	(0,0.1)	3.836(18)	(0,0.3)	3.803(18)
7	(0,0.1)	4.192(13)	(0,0.1)	4.17(2)	(0,0.3)	4.23(3)
8	(0,0.1)	4.483(17)	(0,0.1)	4.45(3)	(0,0.3)	4.54(5)
9	(0,0.1)	4.75(2)	(0,0.1)	4.72(4)	(0,0.3)	4.83(6)
10	(0,0.1)	5.00(3)	(0,0.1)	4.95(5)	(0,0.3)	5.10(7)

Figure 3 shows the behavior of the χ_q susceptibilities as a function of the lattice size.

Zoom In Zoom Out Reset image size

By fitting the susceptibilities presented in figure 3 we have estimated the $\tau(q)$ exponents using equation (8). The numerical estimates for the $\tau(q)$ exponents (and consequently those for the $\zeta(q)$ exponents) can be found in table 3.

Table 3.  $\tau(q)$ and $\zeta(q)$ exponents from the finite size scaling of the susceptibilities χ_q . For χ₄, we do not detect divergence with L. Thus, for $q\unicode{x2A7E} 4$, $\tau(q)\unicode{x2A7E} 3$ and $\theta(q) \unicode{x2A7E} 2.89(1)$. We have computed, as a test, the $1+\eta$ exponent in the q = 1 row, and it compares very well with the most accurate available estimate $1+\eta = 1.036(1)$ [19].

q	$\tau(q)$	$\zeta(q)$
1	1.039(2)	1 (by def.)
2	1.893(3)	1.833(5)
3	2.586(4)	2.503(6)
4	$\unicode{x2A7E} 3$	$\unicode{x2A7E} 2.896(3)$

In figure 4, we show the behavior of $R^{(1)}_q$ versus L for three different values of q, showing the power-law divergence implied by equation (12).

Zoom In Zoom Out Reset image size

In table 4, we present our final values for the ζ-exponents and the (discrete) first and second derivatives obtained analyzing the observables $R^{(1)}_q$, $R^{(2)}_q$ and $R^{(3)}_q$, computing their error bars with the bootstrap method, and using only the lattice sizes with $L\unicode{x2A7E} 48$ (as for the global susceptibility analysis) to avoid small scaling corrections ⁵ . We plot our results for $\zeta(q)$ in figure 5 where the strong departure from the linear regime is clear, providing strong indications of multifractal behavior.

Zoom In Zoom Out Reset image size

Table 4.  $\zeta(q)$ exponents and the (discrete) first and second derivatives from the finite size scaling of the ratios $R^{(1)}_q$, $R^{(2}_q$ and $R^{(3)}_q$ of the local susceptibilities $\tilde{\chi}_q$. Note that $\zeta(1) = 1$.

q	$\zeta(q)$	$\zeta(q)-\zeta(q-1)$	$\zeta(q)+\zeta(q-2)-2 \zeta(q-1)$
2	1.830(2)	0.830(2)
3	2.517(5)	0.688(3)	−0.143(2)
4	3.09(1)	0.573(6)	−0.115(3)
5	3.57(2)	0.481(9)	−0.092(4)
6	3.98(3)	0.40(1)	−0.077(6)
7	4.31(5)	0.33(2)	−0.07(1)
8	4.56(8)	0.25(4)	−0.08(2)
9	4.7(1)	0.16(6)	−0.09(2)
10	4.7(2)	0.05(9)	−0.10(3)

Furthermore, the results reported in table 4 show that the function $\zeta(q)$ is a non-decreasing function for the simulated values of q. In addition, it is a concave function of q, in agreement with the results from Ludwig [5] (see also A.2).

In the replicated theory the C_q correlation functions can be written as,

$$\begin{align} \langle \left(\phi^{a_1}\left(\boldsymbol{x}\right)\cdots \phi^{a_q}\left(\boldsymbol{x}\right)\right) ~ \left(\phi^{a_1}\left(\boldsymbol{y}\right)\cdots \phi^{a_q}\left(\boldsymbol{y}\right)\right) \rangle \,, \end{align} \tag{ A.1 }$$

with $1 \unicode{x2A7D} a_i\lt a_j \unicode{x2A7D} n$ for i < j, where $\phi^{a}(\boldsymbol{x})$ are the replicated fields with scaling dimension $[(D-2+\eta)/2]$ and n is the number of replicas (n → 0).

As discussed in [5], the composite operator $(\phi^{a_1}(\boldsymbol{x})\cdots \phi^{a_q}(\boldsymbol{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 $\Phi_q^{{\mu},\alpha}(\boldsymbol{x})$, and they each transform following the µ-irreducible representation (irrep) of S_n (α labels the multiplicity of the µ-irrep) [5], with scaling dimension $X_q^{{\mu},\alpha}$. Its asymptotic behavior is,

$$\begin{align} \langle \Phi_q^{{\mu},\alpha}\left(\boldsymbol{x}\right) \Phi_q^{{\mu},\alpha}\left(0\right) \rangle \sim \frac{1}{|\boldsymbol{x}|^{2 X_q^{\mu,\alpha}}}\,. \end{align} \tag{ A.2 }$$

Therefore, the behavior of C_q will be the sum (over the irreducible representations) of decays with powers $1/|\boldsymbol{x}|^{2 X_q^{\mu,\alpha}}$, and the smallest scaling dimension $X_q^{{\mu},\alpha}$ 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+\eta)/2] \zeta(q)$ is the scaling dimension of the full composite operator $(\phi^{a_1}(\boldsymbol{x})\cdots \phi^{a_q}(\boldsymbol{x}))$, which defines $\zeta(q)$ in such a way that $\zeta(1) = 1$.

Hence, the scaling behavior of C_q will be twice that of the composite operator $(\phi^1(\boldsymbol{x})\cdots \phi^q(\boldsymbol{x}))$, namely $[(D-2+\eta)/2] \zeta(q) \times 2$.

In the case of the local susceptibilities $\tilde{\chi}_q$ we need to compute,

$$\begin{align} \frac{1}{V} \sum_{\boldsymbol{x}} \sum_{\boldsymbol{y}_1}\cdots \sum_{\boldsymbol{y}_q} \overline{\langle S_{\boldsymbol{x}} S_{\boldsymbol{y}_1}\rangle \cdots \langle S_{\boldsymbol{x}} S_{\boldsymbol{y}_q}\rangle } \,, \end{align} \tag{ A.3 }$$

so that in terms of the replicated theory, the associated correlation function can be written as,

$$\begin{align} \langle \left(\phi^{a_1}\left(\boldsymbol{x}\right)\cdots \phi^{a_q}\left(\boldsymbol{x}\right)\right) ~\phi^{a_1}\left(\boldsymbol{y}_1\right) \cdots \phi^{a_q}\left(\boldsymbol{y}_q\right) \rangle\,, \end{align} \tag{ A.4 }$$

and the scaling dimension of the only composite operator $(\phi^{a_1}(\boldsymbol{x})\cdots \phi^{a_q}(\boldsymbol{x}))$, is again $[(D-2+\eta)/2] \zeta(q)$ and each of the q-factors $\phi^a(\boldsymbol{y}_a)$ have $(D-2+\eta)/2$ as their scaling dimension.

Note that a field φ with $\mathrm{dim}(\phi)$ as the scaling dimension transforms following the rule $\phi(b \boldsymbol{x}) = b^{-\mathrm{dim}(\phi)} \phi(\boldsymbol{x})$.

Equation (A.3) can be written in the continuum as,

$$\begin{eqnarray*} \tilde{\chi}_q &\sim \frac{1}{L^D}\int {\textrm{d}}^Dx~\int {\textrm{d}}^Dy_1 \cdots \int {\textrm{d}}^Dy_q\ \langle \left(\phi^{a_1}\left(\boldsymbol{x}\right)\cdots \phi^{a_q}\left(\boldsymbol{x}\right)\right) ~\phi^{a_1}\left(\boldsymbol{y}_1\right) \cdots \phi^{a_q}\left(\boldsymbol{y}_q\right) \rangle\,, \end{eqnarray*}$$

where the q + 1 integrals run in the volume $[a,L]^D$ (a is the lattice spacing). After the change of variables $\boldsymbol{\tilde{x}} = \boldsymbol{x}/L$ and $\boldsymbol{\tilde{y}}_i = \boldsymbol{y}_i/L$ ($i = 1,\dots, q$) we can write:

$$\begin{eqnarray*} \tilde{\chi}_q &\sim \frac{1}{L^D} L^{(q+1) D} \frac{1}{L^{\left((D-2+\eta)/2\right) \zeta(q)}} \frac{1}{L^{\left((D-2+\eta)/2\right) q}}\\ &\quad \times\int {\textrm{d}}^D \tilde{x}~\int {\textrm{d}}^D \tilde{y}_1 \cdots \int {\textrm{d}}^D \tilde{y}_q\ (\phi^{a_1}(\boldsymbol{\tilde{x}})\cdots \phi^{a_q}(\boldsymbol{\tilde{x}})) ~\phi^{a_1}(\boldsymbol{\tilde{y}}_1) \cdots \phi^{a_q}(\boldsymbol{\tilde{y}}_q) \rangle\\ &\sim L^{(q D - \left[(D-2+\eta)/2\right] (\zeta(q)+q)} \,, \end{eqnarray*}$$

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 obtain the behavior given in equations (12)–(14).

Ludwig showed that $X_q^{\mu,\alpha}$ is a concave function ⁸ of q [5]. The argument is simple, so we reproduce here for $\zeta(q)$. Note that if, for a given r, the correlation function for a given sample can take only values in the compact interval $[0,1]$, then its probability density function (over the samples) is determined by its (qth) positive moments, which we have denoted as C_q . Moreover, these moments are smooth functions of q. Therefore, the $\zeta(q)$ exponents can also be defined for real q. In addition, $C_q^{1/q}$ is a non-decreasing function of q, which implies that $\zeta(q)$ is a concave functions since $\zeta(q_1)/q_1 \unicode{x2A7D} \zeta(q_2)/q_2$ as $0\lt q_2\lt q_1$.

In [15], Davis and Cardy computed the scaling exponents $X_q^{\mu,\alpha}$ for the DIM in $2+\epsilon$ dimensions. In the notation used in this paper, their result can be written as,

$$\begin{align} q-\zeta\left(q\right) = \left( \frac{q \left(q-1\right)}{2}\right) y + O\left(y^2\right)\,, \end{align} \tag{ A.5 }$$

where y is the RG eigenvalue of the disorder strength, which is, at this order, $y = \alpha = O(\epsilon)$ (α is the specific heat exponent of the pure model and $\epsilon = D-2$). These results clearly show that the DIM in $2+\epsilon$ is expected to undergo multiscaling behavior. However, in order to have an accurate analytical estimate of the difference $q-\zeta(q)$ for the 3D model one would need to extend this computation to higher orders of y [2].