{'fraction_non_alphanumeric': 0.03975076443777765, 'fraction_numerical': 0.026596665320486934, 'mean_word_length': 4.369888475836431, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
92 Years of the Ising Model: A High Resolution Monte Carlo Study

Using extensive Monte Carlo simulations that employ the Wolff cluster flipping and data analysis with histogram reweighting and quadruple precision arithmetic, we have investigated the critical behavior of the simple cubic Ising model with lattice sizes ranging from 163 to 10243. By analyzing data with cross correlations between various thermodynamic quantities obtained from the same data pool, we obtained the critical inverse temperature Kc = 0.221 654 626(5) and the critical exponent of the correlation length ν = 0.629 912(86) with precision that improves upon previous Monte Carlo estimates.


Introduction
The Ising Model [1] has long played a special role in the theory of phase transitions and has served as a "fruit fly" model for testing new numerical and theoretical approaches. While Ising published the exact solution in 1-dim in 1925 [1] and Onsager the 2-dim version in 1944 [2], 3dim Ising models have defied solution. Monte Carlo simulations [3], non-equilibrium relaxation Monte Carlo [4], Monte Carlo renormalization group [5,6], field theoretic methods [7,8] and high-temperature series expansions [9] have provided precise values for the location of the phase transition [10] and critical exponents, although the results did not all agree within the error bars. Simulations were insufficient to test Rosengren's "exact conjecture" for the critical temperature of the simple cubic Ising model [11]; moreover, Fisher showed that other "exact conjectures" gave quite similar values [12]. (For a nice review of results prior to 2002 see Ref. [10].) Recent developments have renewed interest in the critical behavior of the 3d Ising model. The conformal bootstrap method provided precise estimates for the critical exponent ν [13,14,15], and Monte Carlo simulations [16,17,18] and tensor renormalization group theory with highorder singularity value decomposition [19] also yielded very precise results. Although Zhang claimed to have solved the 3d Ising model exactly [20], convincing arguments showed that his solution was incorrect [21,22,23].

Three-dimensional Ising model
We have considered the simple cubic, ferromagnetic Ising model with nearest-neighbor interactions on L × L × L lattices with periodic boundary conditions. Each of the lattice sites i has a spin, σ i = ±1 which interact with Hamiltonian where i, j denotes distinct pairs of nearest-neighbors and J is the interaction constant. The dimensionless energy E is The dimensionless coupling constant K = J/k B T will be useful for discussing critical behavior.

Monte Carlo algorithm
We sampled states using the Wolff cluster flipping algorithm [32] in which single clusters are grown and flipped sequentially. Clusters are created by drawing bonds to all nearest neighbors of the growing cluster with probability Random numbers were produced using the Mersenne Twister generator [33]. No difference was discernible between data obtained using 32-bit and 53-bit versions. Simulations were performed at K 0 = 0.221 654 (the estimate for the critical inverse temperature K c by MCRG analysis [6], also used in an earlier, high resolution Monte Carlo study [3]). Lattice sizes were in the range 16 ≤ L ≤ 1024, and measurements were taken for each run after 2 × 10 5 Wolff steps were discarded for equilibrium. Even for the largest lattice, L = 1024, the energy had reached equilibrium by 130, 000 cluster steps, and the simulation was then run at least another ten times the equilibrium relaxation time before data accumulation began. We performed 6,000 to 12,000 runs of 5 × 10 6 measurements for each L using a total of around 2 × 10 7 CPU core hours and producing more than 5TB data.

Methods of analysis
To locate the peaks in response functions precisely we used histogram reweighting [24,25] to transform a probability distribution measured at one coupling to one at an adjacent value of coupling. After N configurations are generated at inverse temperature K 0 with a probability proportional to the Boltzmann weight, exp(−K 0 E), the probability of simultaneously observing the system with total (dimensionless) energy E and total magnetization M is estimated, where Z(K 0 ) is the partition function and H(E, M ) is the measured histogram of states. The (estimated) probability distribution for arbitrary K is then where ΔK = K 0 − K and the average value of any function of E and M , f (E, M ),  [3] showed that the critical exponent ν of the correlation length can be estimated more precisely from Monte Carlo simulation data if multiple quantities, including traditional quantities which still have the same critical properties, are included. The logarithmic derivative of any power of the magnetization for i = 1, 2, ..., can yield an estimate for ν. We also included (reduced) magnetization cumulants U 2i [34] defined by whose derivatives with respect to K can also be used to estimate ν. The inverse critical temperature K c (L) can be estimated from the locations of the peaks in the above quantities as well as the specific heat the derivative of |m| with respect to K, the finite-lattice susceptibility, and the zero of the fourth-order energy cumulant Note that Eq. (11) is for the finite-lattice susceptibility, not the "true" susceptibility calculated from the variance of m, χ = KL d ( m 2 − m 2 ) which has no peak for finite L. (For long runs, m = 0 for H = 0 so that any peak in χ is merely due to finite run length.) All quantities were calculated using the GCC Quad-Precision Math Library (quadruple, 128 bit) precision. The critical behavior of a system near a 2nd order transition can be extracted from the size dependence of the singular part of the free energy density [28,29,30,31]. As a function of L (linear dimension) and T (temperature) as variables, the free energy of a system is described by the finite size scaling (FSS) ansatz, where ε = (T − T c )/T c (T c is the infinite-lattice critical temperature), h is the magnetic field, and the critical exponents α, β, γ and ν are the infinite lattice values. Various thermodynamic properties can be determined from Eq. (13) and have corresponding scaling forms, e.g., where M 0 (x), χ 0 (x) and C 0 (x) are scaling functions. Because we are interested in zero-field properties (h = 0), x is the only relevant thermodynamic variable. Then, for finite L, e.g. the specific heat peak occurs where the scaling function C 0 is maximum, i.e., when The temperature corresponding to the peak is the finite-lattice "transition" temperature T c (L), on the condition x = x * which varies with L asymptotically as The FSS ansatz is valid only for sufficiently large L and T sufficiently close to T c . Corrections to scaling and finite-size scaling terms appear for smaller systems and temperatures away from T c due to: irrelevant scaling fields which can be expressed in terms of an exponent θ, e.g. a 1 ε θ + a 2 ε 2θ + · · · , and others due to the non-linear scaling fields, e.g. b 1 ε 1 + b 2 ε 2 + · · · . The correction terms can then be expressed by: We first estimate the critical exponent ν and then insert it into Eq. (19), so that there is one less unknown parameter to do the non-linear fit to Eq. (19). To do this we can use the critical scaling form without the prior knowledge of the transition coupling K c Once ν is determined from the fit of Eq. (20) and Eq. (21), we can estimate the critical inverse temperature K c with a fixed value of ν. K c can also be estimated from Binder's 4th order cumulant crossing technique [34]. As L → ∞, the fourth-order magnetization cumulant U 4 → 0 for K < K c and U 4 → 2/3 for K > K c . U 4 can be plotted as a function of K for different lattice sizes, and the location of the intersections between curves for the two lattice sizes is given by where L is the size of the smaller lattice, b = L /L is the ratio of two lattice sizes, and ω 1 , ω 2 are correction exponents in the finite-size scaling formulation. Ideally, each configuration only depends on the previous configuration, but, in practice, it is likely to be correlated to earlier configurations. The correlation decays with separation in Monte Carlo time, but fluctuations appear smaller than they should be. To deal with this issue, we can consider blocks of the original data, and use jackknife resampling [35].
An important advance by Weigel and Janke [26,27] was the seminal observation that cross correlation between different quantities could lead to systematic bias in the estimates of critical quantities extracted from the data. For a set of n measurements of a random variable x = (x 1 , x 2 , · · · , x n ), and an estimatorθ = f (x), the Jackknife estimate of the value and error ofθ is found by leaving out one measurement at a time. We define the jackknife average, x J i by, where i = 1, 2, · · · , n, so x J i is the average of all the x values except x i . Similarly, we definê The jackknife estimate ofθ = f (x) is the average ofθ J i , i.e.
and the jackknife error σ(θ), is given by, In general there can be multiple adjacent elements in each block, and significant cross correlations may exist between estimatesθ (k) andθ (l) from the same original time-series data, To reduce these effects, we calculated the jackknife covariance matrix G ∈ R m×m [35]. For a number of estimatesθ (k) , the r th row, c th column entry of matrix G is given by, The m different estimatesθ (k) (k = 1, 2, ..., m) for the same parameterθ, should have the same expectation value. So the estimated value forθ can be determined by a linear combination, where k α k = 1. Based on the cross correlation analysis from Ref. [26,27], with the optimal choice for the weights, the variance can be expressed by, At our level of resolution the finite precision of the pseudorandom number generator and the numbers that can be stored in memory become relevant. When the Wolff algorithm probability of adding a spin to the cluster is converted to a 32-bit unsigned number for comparison with a pseudorandom number generated in the simulation, it is truncated due to finite precision. When reconverted back into a value of K the result differs from 0.221 654 in the 10th decimal place. For the largest L, this is only a factor of 20 smaller than the statistical error. Simulations with a 53bit pseudorandom number generator showed that this is not significant for our analysis, but for future studies of larger systems and/or higher precision, a 32-bit random number generator might not suffice. We also used the corrected, effective K 0 instead of 0.221654 and for L = 1024 we used the multiple-histogram method [25] to combine results for the 32 and 53-bit pseudorandom number generators.
Comparisons of results from Wolff cluster flipping runs and Metropolis single spin-flip runs for L = 32 showed no differences. The Wolff cluster simulation for L=32 was repeated using the MRG32K3A random number generator from Pierre L'Ecuyer, "Combined Multiple Recursive Random Number Generators", Operations Research, 47, 1 (1999), 159-164. The locations and values of the maxima in all quantities agreed, to within the error bars, as those from the Metropolis simulations and the Wolff simulations with the Mersenne Twister. Hence the problems found by Ferrenberg et al. [36] using other random number generators were not noticeable here.  Table 1. Estimates for ν with: one fixed correction exponent ω 1 = 0.83, two fixed exponents ω 1 = 0.83, ω 2 = 4, and three fixed exponents ω 1 = 0.83, ω 2 = 4, ω 3 = 1.6 as a function of L min .

Finite-size scaling analysis
The conformal bootstrap [13,14] gives a theoretical prediction for the confluent correction exponents, ω 1 = 0.8303 (18), ω 2 ≈ 4. (29) and the correction term corresponding to the non-linear scaling fields is [37], Thus, the fitting model for the peak value of quantity X is, We analyzed the data fitting model Eq. (31) using between one and three correction terms with fixed powers ω 1 = 0.83, ω 2 = 4, ω 3 = 1.6. By calculating the covariance matrix and doing the cross-correlation analysis, we give estimates for ν in Table 1 where the minimum lattice size included in the analysis, L min , is eliminated one by one. With three correction exponents, two for confluent corrections (ω 1 = 0.83, ω 2 = 4) and one for the non-linear scaling fields (ω 3 = 1.6), the estimate for the critical exponent ν shows statistical fluctuations; but the value for ν was estimated by taking the average of ν obtained from fits for different L min and the error estimated with the jackknife method on estimates of ν from the analysis using different ranges of relatively large L min . Using values of L min from 80 to 144, we find Since higher order terms like a 1,2 L −2ω 1 + · · · could be present and larger than the higher order confluent corrections, we repeated the analysis using these higher order corrections instead. The results for the asymptotic values of the critical temperature and exponent were essentially identical. Thus, either the higher order terms have small prefactors or we do not have sufficient resolution to differentiate their contribution from those of higher order confluent corrections.
The fluctuation of K c when L min ≤ 80 is larger than when L min ≥ 80 and finite-size effects diminish for larger L. From the average of K c for L min = 80 to 144 we estimate K c = 0.221 654 626 2. Likewise, a jackknife analysis has been done on the estimates for K c which are obtained from the three correction terms analysis. From values of L min from 80 to 144 we estimate K c = 0.221 654 626 2(23) whereas for L min = 16 to 144, K c = 0.221 654 625 5(42). Therefore, we conclude from the finite size scaling analysis, with conservative error bars, that K c = 0.221 654 626(5).
As seen in Fig. 2, effects of higher order correction terms are clearly visible.

Crossing technique of the 4th order magnetization cumulant
As the lattice size L → ∞, the fourth-order magnetization cumulant U 4 → 0 for K < K c and U 4 → 2/3 for K > K c . For large enough lattice sizes, curves for U 4 cross as a function of inverse temperature at a "fixed point" U * , and the location of the crossing "fixed point" is K c . Because of finite-size correction terms, not all curves cross at a common intersection. The locations of the cumulant crossings have been fitted to Eq. (22) with two correction terms. The critical coupling appears to be stable if L min ≥ 96 and the average of K c values for L min ≥ 96 yields

Discussion
The combination of an efficient, cluster-flipping Monte Carlo algorithm, high statistics simulations, histogram reweighting, and a cross-correlation jackknife analysis enables the high resolution of the results presented above. In Fig. 3, we show the results of our Monte Carlo study as well as other high-resolution simulations estimates for ν and K c . The boxes represent the quoted error bars in both ν and K c assuming independent errors. Our estimate ν = 0.629 912(86) is consistent (i.e. within the error bars) with the recent conformal bootstrap result of Kos et al. [15], as well as that from an older work by El-Showk et al. [14]. Also, our result agrees well with the Monte Carlo result of Hasenbusch [16] but is lower      . High-resolution estimates for K c and ν obtained using different methods than that of Weigel and Janke [27]; however, within the respective error bars there is agreement although we have substantially higher precision than either of these previous studies. Our system sizes and statistics are substantially greater than those used by Weigel and Janke and Hasenbusch examined the behavior of the ratio of partition functions Z a /Z p , and the second moment correlation length over the linear lattice size ξ 2 /L so the methodologies differ. Our estimate for K c differs from that obtained by Kaupuzs et al [17] using a parallel Wolff algorithm by an amount that is slightly outside the error bars. Somewhat perplexingly, they were able to fit their data to two rather different values of ν, so no comparison of critical exponents is possible. The recent tensor renormalization group result for K c [19] disagrees with our result by many times the respective error bars.

Conclusion
We have studied a 3d Ising model with the Wolff cluster flipping algorithm, histogram reweighting, and finite size scaling including cross-correlations using quadruple precision arithmetic for the analysis. Using a wide range of system sizes, with the largest containing more than 10 9 spins, and including corrections to scaling, we have obtained results for K c , ν, and γ that are comparable in precision to those from the latest theoretical predictions and can provide independent verification of the predictions from those methods. Our values provide further numerical evidence that none of the purported "exact" values are correct. To within error bars we obtain the same value for the critical exponent ν as that predicted by the conformal bootstrap; however, our estimate for the critical temperature K c does not agree with that from the tensor renormalization group.
As computer power increases, new sources of error must be taken into account. To improve precision more stringent tests of the random number generator and much greater statistics for much larger lattice sizes will be needed.