Critical exponent $\nu$ of the Ising model in three dimensions with long-range correlated site disorder analyzed with Monte Carlo techniques

We study the critical behavior of the Ising model in three dimensions on a lattice with site disorder by using Monte Carlo simulations. The disorder is either uncorrelated or long-range correlated with correlation function that decays according to a power-law $r^{-a}$. We derive the critical exponent of the correlation length $\nu$ and the confluent correction exponent $\omega$ in dependence of $a$ by combining different concentrations of defects $0.05 \leq p_d \leq 0.4$ into one global fit ansatz and applying finite-size scaling techniques. We simulate and study a wide range of different correlation exponents $1.5 \leq a \leq 3.5$ as well as the uncorrelated case $a = \infty$ and are able to provide a global picture not yet known from previous works. Additionally, we perform a dedicated analysis of our long-range correlated disorder ensembles and provide estimates for the critical temperatures of the system in dependence of the correlation exponent $a$ and the concentrations of defects $p_d$. We compare our results to known results from other works and to the conjecture of Weinrib and Halperin: $\nu = 2/a$ and discuss the occurring deviations.


I. INTRODUCTION
The influence of quenched disorder on phase transition properties of a system is of great importance as many real-world materials show defects or impurities. The simplest way to introduce the disorder is by assuming it to be point-wise and uncorrelated. A prominent achievement in describing the critical behavior of such systems was done by Harris [1]. The result is known as Harris criterion. It states that if the system has a negative specific heat exponent in the pure case (without disorder, α pure < 0) the disorder does not influence the system's universality class. On the other hand, for α pure > 0 the disorder will change the system's universality class. This universality class will have new critical exponents which will not depend on the disorder concentration. Various studies [2][3][4][5][6] confirmed the change of the universality class of the three-dimensional Ising model for which α pure > 0 is true.
However, in nature the disorder usually comes with a certain structure. One possible way to introduce such disorder to a model is by adding a spatial correlation to the disorder. For a magnetic system this could be nonmagnetic lines or planes or clustered nonmagnetic impurities. Other interesting areas are magnetic foams and magnetic elements in porous media [7,8]. The correlated disorder in systems was intensively studied with the help of the renormalization group theory by Weinrib and Halperin [9] and the result is known as the extended Harris criterion. It states that a system with long-range correlated disorder where the spatial disorder correlation follows a power-law ∝ r −a will change its universality class if a < d and otherwise the standard Harris criterion will be recovered. Further, they claim that the critical * kazmin@mis.mpg.de exponent of the correlation length ν in the long-range correlated three-dimensional Ising model is given by They argue, but do not prove rigorously, that this result is exact. Several studies dealt with the Ising model with correlated disorder in two dimensions by applying Monte Carlo techniques [10,11] or renormalization group techniques [12]. In three dimensions Monte Carlo simulations were performed in Refs. [13][14][15][16][17][18][19] while renormalization group techniques were used in Refs. [9,20]. While it is generally accepted that the correlated disorder case belongs to a new universality class, the quantitative results and in particular the claim given in Eq. (1) are controversially discussed. One condition which is often overseen when assuming Eq. (1) is that d = 4 − ≈ 4 and a = 4 − δ ≈ 4 is a necessary condition in Ref. [9]. So it remains unclear which range of a values fulfills this requirement. As a further reinforcement of the prediction given in Eq. (1), Honkonen and Nalimov [21] claimed that Eq. (1) is exact to all orders in the -δ-expansion. This has been further analyzed in Refs. [22,23].
The results for the ν exponent obtained by different groups for the uncorrelated and the long-range correlated disordered three-dimensional Ising model are summarized in Table I. The ambiguity about the numerical values of the critical exponents and considerable differences in the literature motivated us to attack the problem once again.
We extensively analyzed a three-dimensional Ising lattice with power-law correlated site disorder by using Monte Carlo techniques. In contrast to previous works we performed simulations for various different correlation strengths a and a wide range of disorder concentrations p d . We focused on the critical exponent of the correlation length ν and the confluent correction exponent ω and obtained a global picture of their behaviors in the arXiv:2008.03169v1 [cond-mat.stat-mech] 7 Aug 2020 TABLE I. Various results of the critical exponent ν and the confluent correction exponent ω for the three-dimensional Ising model with uncorrelated and long-range power-law correlated disorder. We schematically denote the uncorrelated disorder case with a = ∞. The rest of the paper is structured as follows. In Section II we specify our model and the details of the performed simulations. In Section III we analyze the disorder realizations to confirm the desired power-law behavior. The main analysis of the Monte Carlo simulations of the Ising model and the obtained results are contained in Section IV. We present the derivation of the critical exponent ν as well as the correction exponent ω. We compare our results to the Weinrib and Halperin conjecture, ν = 2/a, and to the known results. Finally we obtain critical temperatures for different concentrations and correlation exponents. A conclusion in Section V completes this work.

A. Ising Model with Site Disorder
We will not discuss the standard Ising model here and refer to [24,25] as a good starting point for readers who need a deeper background. For the rest of the paper we will deal with the Ising model with site disorder which we will refer to as disordered Ising model 1 . The Hamiltonian of the Ising model with site disorder has a very similar form to the standard Ising model where the spins can take the values s x = ±1 and the defect variables can be η x = 1 when the spin is present 1 The Ising model with random couplings, i.e., bond disorder, is also called "disordered Ising model" in the literature. at site x and η x = 0 when the site x is empty (a defect). The sum runs over all next-neighbors denoted by xy . The coupling constant is set to J = 1 on the whole lattice and we work without an external magnetic field, i.e., h = 0. Schematically the Ising model with and without site disorder is presented in Fig. 1. We distinguish between two different disorder types. The first type is the uncorrelated disorder or random disorder. In this case the defects are chosen randomly according to the probability density where p s is the concentration of spins, p d = 1 − p s is the concentration of defects and δ is the Dirac-delta distribution. The second type is the correlated disorder. In this case the probability density for the defects is again given by Eq. (3). However, now additionally the spatial correlation between the defects decays according to a power-law where r(x, y) is the distance between sites x and y and a ≥ 0 is the correlation exponent. Note, that for both cases we work in the so-called grand-canonical approach where the desired concentration p d is a mean value over a large number of realizations while for each separate realization p d can vary. In Fig. 2 we show slices of a three-dimensional Ising model lattice with different concentrations of defects and different correlation exponents near the critical temperature. According to the Harris criterion and the extended Harris criterion the disordered three-dimensional Ising model falls into three different universality classes in dependence of the correlation exponent a and the concentration of defects p d . The pure case where no defects are present (p d = 0), the effectively uncorrelated case for a > d and the correlated case for a ≤ d. These cases are schematically shown in Fig. 3.

B. Monte Carlo Simulation Details
We performed Monte Carlo simulations of the disordered Ising model and used the Swendsen-Wang multicluster update algorithm [27]. The linear lattice sizes of our cubic lattices were in the range between L = 8 and L = 256 and we chose periodic boundary conditions in each direction. The correlation exponent values were a = 1.5, 2.0, 2.5, 3.0, 3.5 and ∞ which we will use symbolically for the uncorrelated case. For each a value we simulated eight defect concentrations p d = 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35 and 0.4. After a thermalization period of 500 sweeps we performed N = 10 000 measurement sweeps at each considered simulation point β sim = 1/(k B T sim ). Throughout the paper we will refer to the inverse temperature defined by β = 1/(k B T ) simply as "temperature". The temperatures were first chosen in a wide range and with larger spacing for small lattices. After the first analyses refined ranges (regions around the critical points for considered observables for finite lattice sizes L) were estimated and larger lattice sizes were simulated at less temperatures. For each parameter tuple (a, p d , L) we simulated N c = 1000 disorder realizations. After each sweep we measured and stored the total energy E and the total magnetization of the system M At the end we had two-dimensional arrays of values E c i and M c i where i = 1, . . . , N and c = 1, . . . , N c for each parameter tuple (a, p d , L, β sim ). This was needed in order to apply the reweighting technique in later analysis.

III. CORRELATED DISORDER ANALYSIS
Before we move to the analysis of the Monte Carlo simulations of the Ising model we first take a look at the site disorder realization and analyze the generated ensembles. It is a necessary step to gain control over the correlation exponents a of the disorder ensembles on which we will perform the simulations later on.

A. Disorder Generation
In this work we mainly study the Ising model on a lattice with uncorrelated and long-range correlated site disorder. An important part is the generation of the site disorder for later Monte Carlo simulations. The uncorrelated disorder case is realized by setting the defect vari-ables η x for each site x of the lattice according to where 0 ≤ R x < 1 is a uniform random number drawn for each site x.
For the case of long-range correlated disorder let us first define the correlation function C η between two defects η at sites x and y at a distance r = |x − y| In this work we assume a power-law decay of the correlation function for large distances r 1 We used a modified Fourier method by Zierenberg et al. [26] for the generation of long-range correlated disorder. Initially the Fourier method was introduced by Makse et al. [28]. The code linked in [26] was used in this work. We will not discuss the details of the generation and only sketch the process: 1. Generate uncorrelated, normally distributed random variables.
2. Perform a Fourier transformation of these variables.
3. Correlate the transformed variables by multiplying with a spectral density generated from a chosen correlation function C 0 .
4. Fourier transform the correlated variables back to real space.
5. Truncate the final variables to {0, 1} such that the mean concentration of zeros equals the desired concentration of defects p d .
The resulting η variables are correlated and their correlation function is approximately given by C 0 . This approximation comes from the fact that the truncation in step 5 is not mathematically exact and introduces deviations from the desired function C 0 .
In order to overcome the infinity at r = 0 in Eq. (9) we used a slightly modified correlation function which asymptotically approaches Eq. (9) for large distances, We generated ensembles of disorder realizations by providing two parameters: the correlation decay exponent a and the concentration of defects p d . Because the step 5 is mathematically not exact and also because a modified correlation function was used, we verified both values a and p d for each ensemble numerically.

B. Mean Concentration of Defects
First we looked at the distribution of the concentrations of defects p d for each parameter tuple (a, p d , L). Examples of the distributions are shown in Fig. 4. We verified the normality of the distributions for each ensemble with the help of the Anderson-Darling test [29,30]. Apart from the strongest correlation with a = 1.5 at low p d ≤ 0.2 all distributions for L ≥ 24 were classified as normal with 95 % confidence. The results of the test for all parameter tuples (a, p d , L) are presented in Fig. 5. It can be seen that higher concentrations approach the normal distribution already for smaller L. The estimated concentrations p d as a mean over all lattice sizes for each ensemble are listed in Table II. They match the imposed concentrations p d perfectly in all cases.

C. Mean Correlation Exponent
The correlation function C η (r) was calculated as a mean over all configurations for each parameter tuple (a, p d , L). It was measured for two different distance directions (along the x-axis and along the diagonal), and all possible distances in the corresponding direction. The correlation function was calculated by where C is the normalization constant such that C η (0) = 1 and N r is the number of possible realizations of the distance r on the lattice. From the chosenr i vectors and from periodic boundary conditions it follows that The sum in Eq. (13) runs over all site pairs x and y which have the vector distance rr i where i = 1, 2. The normalization constant turns out to be Once the correlation functions defined through Eq. (13) were measured for each disorder ensemble, we had to obtain the correlation exponent a. We performed a fit to the linearized ansatz on a logarithmic scale corresponding to the asymptotic behavior of Eq. (9) where a is the desired decay exponent. We had to find a minimal distance r min included into the fits in order to obtain the correlation exponent for r > r min 1 where the assumption of a power-law decay is valid. We used the condition that r min is the distance where the relative deviation between C 0 and C η became less than a threshold value of C = 0.05 for the first time, Note that the amplitudes for C 0 and C η were omitted as we assume them to be equal and to cancel in Eq. (17). Eq. (17) leads to the condition Furthermore we had to set a maximum distance r max to exclude the noisy tail of the correlation function and possible finite-size effects. Here we have chosen the distance r max where the absolute value of the measured correlation function C η was below a minimal threshold value of C min = 10 −5 for the first time, For small lattices with L ≤ 20 and weak correlations (large a) sometimes the found r min and r max where too close together or even r min > r max . Is such cases we reduced r min until a fit with 4 degrees of freedom was possible. The estimated a(p d , L) are shown in Fig. 6 and the final averages are summarized in Table III while in  TABLE III in Fig. 7 examples of the correlation function fits are presented. The final results a are means over all p d and L ≥ L min which were chosen for each a according to the quality of the fits. Please note that we will still refer to different ensembles by the imposed a for clarity.
As naturally follows from the described determination of r min and r max , smaller a(p d , L) have more degrees of freedom and therefore the estimated values a(p d , L) coincide better with the proposed a. For weak correlations with a ≥ 3.0 we exhibit poorer agreement and larger errors for lattice sizes L 128. Also a systematic underestimation of a can be seen in the results. It becomes more pronounced with larger a and smaller L. We have plotted the relative deviation of the estimates a to the ex- pected values a in Fig. 8. One can see a constant increase in the deviations for increasing a. For our largest a value the deviation reaches ≈ 5 %. Nevertheless, we can state that we achieve the desired a values within a precision of ≈ 5 %. A test involving more realizations considerably improved the results for the weak correlation cases but we wanted to stay with the number of disorder realizations for which the Monte Carlo simulations were performed later.

IV. FINITE-SIZE SCALING ANALYSIS
We will now discuss the extraction of the critical exponent of the correlation length ν and the confluent correction exponent ω. For the finite-size scaling analysis we chose the derivative with respect to the inverse temperature β = 1/(k B T ) of the logarithm of the magnetization ∂ β (ln [ |m| ]). It can be expressed in terms of expectation values as where · denotes the thermal average and [·] the disorder average and e = E/V , m = M/V are the normalized energy and magnetization, respectively. Note, that we use the common convention of taking the absolute value of m to avoid the trivial averaging to zero in the lowtemperature phase for finite lattice sizes. For the sake of clarity, we will omit the average brackets for the rest of this work and simply write ∂ β (ln |m|) = ∂ β (ln [ |m| ]).
The derivative of the logarithm of the magnetization ∂ β (ln |m|) is known to diverge at the critical temperature in the thermodynamic limit L → ∞. For finite system sizes it hence develops a minimum. The finite-size scaling behavior up to the first-order correction reads where ∂ β (ln |m|) min (L) is the (finite) minimum value of ∂ β (ln |m|)(β) for a given lattice size L. Fitting with this ansatz is difficult as it is a non-linear four-parameter fit. Therefore, we first determined the correction exponent ω separately and used it as a fixed parameter in the final estimation.
The whole finite-size scaling analysis can be split into three main steps. In the first step we derive the peaks of ∂ β (ln |m|). The second step is the extraction of the correction exponent ω which is needed for the fits in the last step. The last step is the fitting of ∂ β (ln |m|) min (L) with fixed ω and the extraction of ν.

A. Peaks of Observables
We start the analysis with the extraction of the peaks of the derivative of the logarithm of the magnetization ∂ β (ln |m|) min (L). Out of all simulated temperatures for each parameter tuple (a, p d , L) we chose three temperatures β i sim with i = 1, 2, 3 where the derivative of the logarithm of the magnetization calculated at these temperatures ∂ β (ln |m|)(β i sim ) was minimal. For these three β i sim we performed a single histogram reweighting of ∂ β (ln |m|) to find the minimum values ∂ β (ln |m|) i min and the corresponding temperatures β i min . The final ∂ β (ln |m|) min was chosen to be the minimum of all three ∂ β (ln |m|) i min values. A more detailed explanation of the reweighting and error estimation process through resampling is presented in Appendix A.
An important issue was to ensure that the histogram reweighting results lay within the reweighting range. This is an inevitable restriction coming from the limited statistics. We used the reweighting range approximation as defined in [31] ∆β rew = 1 We looked at the ratios of the differences between the simulation temperatures β sim and the found temperatures of the minimum values β min with respect to the reweighting range ∆β rew |β sim − β min | ∆β rew .
As can be seen in Fig. 9 all obtained β min were close enough to the corresponding β sim to assume that the use of the reweighting technique is valid.

B. Confluent Correction Exponent ω
The quotient method which we used for the determination of the confluent correction exponent ω was successfully used in other works, e.g., [2,32,33]. Starting from an observable O which has a peak at the critical temperature we build quotients of O at different lattice sizes L and sL where the observables are taken at the critical temperatures for the given lattice sizes L and sL, respectively, and s is an arbitrary positive (integer) factor. For a dimensional observable, e.g., ∂ β (ln |m|), the finite-size scaling of Q O in leading order reads [33] Q where x O is the critical exponent of O.
We calculated the quotients defined through Eq. (24) for O = ∂ β (ln |m|) with x O = 1 and for s = 4. This allowed us to have 8 independent Q values without using the same lattice size twice. We used the peak values ∂ β (ln |m|) min (L) and performed a global fit to Q ∂ β (ln|m|) (L, p d ) according to Eq. (25) but using all p d simultaneously where we explicitly denote the dependence of the amplitudes A p d on the concentrations of defects with the index p d and relate the constant C to Eq. (25) with In Fig. 10 we present the ω results and the qualities of the fits χ 2 red for p max  Table IV.
From Fig. 11 we can clearly see a distinction between the uncorrelated and correlated cases. The correction exponent for the uncorrelated case ω = 0.373(53) matches the prediction ω = 0.37(6) made by Ballesteros et al. [2]. The correction exponent ω = 1.047(90) for the case a = 2.0 is in good agreement with the value ω = 1.01 (13) obtained by Ballesteros and Parisi [13]. A value around ω ≈ 0.95 (10) is also found for all other a parameters. As the errors (ω) are quite large for all correlated cases, a = ∞, chances are that the correction exponent ω does not depend on a and has a value of roughly ω ≈ 1. Visually it can be verified in Fig. 11.
where we explicitly denoted the dependence of A p d and B p d on p d . We performed least squares fits to Eq. (28) with A p d , B p d and 1/ν as parameters and used fixed correction exponents ω(a) listed in Table IV. Examples of the resulting fits are shown in Fig. 12.
We performed the fits for various minimal lattice sizes 20 ≤ L min ≤ 64. We also varied the smallest concentration p min   Table V. Additionally, the ν values are shown in Fig. 14.
The obtained value for the uncorrelated case ν = 0.6875(47) is in very good agreement with the results from other groups listed in Table I. Please note that in most works the ν exponent was concentration dependent  in contrast to the present work. Therefore the comparison must be done with care. Altogether we can conclude that our extraction method and in particular the global fit ansatz work well for the uncorrelated case which can be seen as a good verification.
For the correlated disorder cases we first compare our results to the prediction of the extended Harris criterion.
All obtained values lie about 10 % above the prediction of ν = 2/a. Nevertheless we see the right tendency of the ν values in being proportional to 1/a and in approaching the uncorrelated case somewhere around a ≈ 3.0. The crossover region around a ≈ 3.0 shows the largest deviations from the extended Harris criterion estimate as well as from the uncorrelated case. This behavior is expected for finite systems. The estimate for a = 1.5 has a huge error and therefore is not very representative. Probably more realizations are needed to get a better result for such strongly correlated case.
Considering the ν values for the correlated case with a = 2.0 we see a discrepancy between our results and results from other groups listed in Table I (see also the  summary plot in Fig. 17). There are several possible reasons for such deviations. Comparing to the work of Ballesteros and Parisi [13] and Ivaneyko et al. [15], we used much more finer lattice size stepping; 18 lattice sizes in the range of 8 ≤ L ≤ 256 versus 5 sizes in the range 8 ≤ L ≤ 128. The number of realizations in our case was smaller by a factor of 10 but we measured 10   Table III and scaled to 1/a. The uncorrelated disorder case critical exponent was set to ν ∞ = 0.683 as an average value from other works listed in Table I. Expected values ν = 2/a according to the prediction of the extended Harris criterion are shown for comparison for all a ≤ d where the extended Harris criterion is assumed to be valid.
times longer time series on each realization. Further, we used the derivative of the logarithm of the magnetization ∂ β (ln |m|) as our primary observable whereas in the other works the derivatives of Binder cumulants ∂ β U 2 and ∂ β U 4 were used. Additionally, the concrete methods of generating the correlated disorder and extracting the critical exponent ν were very different. Finally, but probably most importantly, the method used in this work included all p d values in the critical exponent ν estimation. Comparing our ν exponents to the results of Prudnikov et al. [20] and [14] we do not see any agreement. The reason for this remains unclear to us.

D. Critical Temperature
Once we have derived the peaks of ∂ β (ln |m|) in Section IV A, we also had the corresponding temperatures β min . This allowed us to study the critical temperatures for for all correlation exponents a and concentrations of defects p d . Note, that unlike for the critical exponent ν we need to attend each p d separately and cannot perform a global fit as the critical temperature depends on it. To obtain the critical temperatures β c for all a and p d values we used the finite-size scaling relation in the leading order where β min are the temperatures corresponding to the minimal values of the derivative of the logarithm of magnetization ∂ β (ln |m|) min at different L and β c is the desired critical temperature at L → ∞. We performed the fits to the ansatz given in Eq. (29) by using the extracted exponents ν for the corresponding a values listed in Table V. The quality of the fits was moderate and varied for different p d and a significantly. Therefore we set L min = 32 for all parameter tuples. Finally, to incorporate the uncertainties in the ν estimates we performed the fits in a bootstrapped way by randomly choosing a ν i = Normal(ν, (ν)) according to a normal distribution and performing 10 000 fits. All final quantities were averages over these bootstrapped fits. The resulting temperatures and the qualities of the fits are presented in Fig. 15 and Table VI. Examples of the fits for different a and p d can be found in Fig. 16.
The qualitative behavior of the temperature curves is in strong agreement with the expectations. When   Table V. the concentration of defects vanishes, p d → 0, the inverse temperature goes to the pure Ising model case with β c = 0.221 654 628(2) [34]. On the other hand, when the concentration approaches the percolation threshold concentration, p d →p d (a), Ref. [26], the inverse temperature becomes infinity, β c → ∞. In contrast to the mini-mal values ∂ β (ln |m|) min which were obtained with a high accuracy, it was not possible to get such precise temperatures β min . The main difficulty was the large width of the peaks of ∂ β (ln |m|)(β) for stronger correlations. Additionally, in some cases the reweighting range was not large enough to cover the temperature of the peak sufficiently. Nevertheless, the estimates provide a consistent picture and can serve as a good starting point for later analyses.

V. CONCLUSIONS
We applied Monte Carlo simulation techniques to the three-dimensional Ising model on a lattice with longranged correlated site disorder. The correlation of the disorder was assumed to be proportional to a power-law ∝ r −a with a correlation exponent a. We provided a decent analysis of the disorder correlation in our disorder ensembles verifying the correlation exponent a numerically.
We found the critical exponents of the correlation length ν and the confluent correction exponents ω as well as the critical temperatures β c of the system for various correlation exponents 1.5 ≤ a ≤ 3.5 as well as for the uncorrelated case a = ∞. Contrarily to other works we performed a global fit where we included different disorder concentrations into one simultaneous fit. Such a study was not possible before because all known works only considered one particular correlation exponent choice a = 2.0 and only one or two different concentrations p d whereas in this work we used a wide range of a and p d values.
In Fig. 17 we give a visual comparison of the critical exponents ν obtained in this work, results known from other works and predictions by the extended Harris criterion. We obtain a value ν = 0.6875(47) for the uncorrelated case which matches the results from other groups listed in Table I and plotted in Fig. 17. Also the correction exponent ω = 0.373(53) coincides with Ref. [13].
The estimated ν values for the correlated disorder cases show the 1/a behavior predicted by the extended Harris criterion qualitatively but are approximately 10 % higher than the prediction 2/a. On the other hand, we strongly disagree with the renormalization group predictions made by Prudnikov et al. [20] and their Monte Carlo simulation result for a = 2.0 in Ref. [14]. The correction exponent ω = 1.047(90) for the case a = 2.0 is in good agreement with Ballesteros and Parisi [13]. For all correlated cases we measure a value which is compatible with ω = 0.95(10) ≈ 1.
Our estimation of the critical temperatures β c provides a global picture of the system for different a and p d parameters and can serve as a good starting point for further analyses and simulations.
In upcoming studies we will consider other critical exponents like β and γ and hopefully tackle down the problem even more.  [14]. The inset shows a close up of the uncorrelated case a = ∞. The uncorrelated disorder case critical exponent was set to ν ∞ = 0.683 as an average value from other works listed in Table I. The results of this work lie about 10 % above the prediction of the extended Harris criterion ν = 2/a. On the other hand, they also do not coincide with other works. The main reason for this discrepancy is probably the global fit ansatz of the present work which combines all p d into one single fit. disorder realization c The average over the disorder realizations is denoted by [·] and the final estimate O reads For variables of the type of Eq. (A1) a histogram reweighting technique can be used to reweight the observable from the simulated temperature β sim to a different temperature β. We used the form given in [31] Rew where the reweighting is performed separately for each disorder realization c and the final estimate at the temperature β is the disorder average Not every observable of interest, in particular the derivative of the logarithm of the magnetization ∂ β (ln |m|) has the form of Eq. (A1). Let P denote a composite observable of the following form Let us summarize what we have achieved so far. Starting with the arrays of raw observables E and M we are able to use the histogram reweighting technique to obtain practically any observable calculable from E and M as a function of β.
Let us now assume that the finite-size scaling analysis of P (β) predicts a minimumP at a certain temperaturě β. Without loss of generality we assume a minimum of P (β), otherwise we transform P → −P . In the thermodynamic limit L → ∞ we expectβ → β c . We can apply an optimization routine by plugging in Rew(P )(β) as the target function and obtainP andβ, P = min β (Rew(P )(β)) . (A8) However, we will not be able to estimate the errors (P ) and (β) as only one final value is calculated through Eq. (A8) from all simulated data. In order to overcome this problem, we can use a resampling technique. We have chosen the jackknife resampling technique which is described, e.g., in [35] in full length. We will only sketch the main steps applied in this work. As our measurements were two-dimensional arrays consisting of time series i = 1, . . . , N and disorder realizations c = 1, . . . , N c , we apply the resampling in both directions separately and combine the estimates at the end. For each jackknife resampling step j in the time series direction we leave out a block J j ⊂ {1, . . . , N } of measurements for each disorder realization c so that the thermal average defined through Eq. (A2) becomes where |J j | is the number of left-out samples. Analogously, for each resampling step k in the disorder direction we leave out a block J k ⊂ {1, . . . , N c } of disorder realizations so that the disorder average defined through Eq. (A3) becomes where |J k | is the number of left-out realizations. Starting from the modified thermal averages (O c ) j and disorder averages (O) k respectively, all steps in the following analysis remain the same. Let A be a final estimate coming from a certain analysis, e.g., minimum search as in Eq. (A8). By repeating the analysis for N j different jackknife blocks in the time direction and N k blocks in the disorder direction we are given two arrays of estimates (A) j and (A) k respectively. We calculate two jackknife means A j and A k A a = 1 N a N a a=1 (A) a with a = j, k , and two corresponding jackknife errors (A) j and (A) with a = j, k .

(A12)
As the last step we combine the two means and errors in a standard (uncorrelated) manner The mean A and the corresponding error (A) are the final results for a given analysis after applying jackknife resampling.