{'fraction_non_alphanumeric': 0.05371257309082361, 'fraction_numerical': 0.030382276026209853, 'mean_word_length': 3.9713701431492843, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Corrections to finite-size scaling in the 3D Ising model based on non-perturbative approaches and Monte Carlo simulations

Corrections to scaling in the 3D Ising model are studied based on non-perturbative analytical arguments and Monte Carlo (MC) simulation data for different lattice sizes L. Analytical arguments show the existence of corrections with the exponent (gamma-1)/nu (approximately 0.38), the leading correction-to-scaling exponent being omega =<(gamma-1)/nu. A numerical estimation of omega from the susceptibility data within 40 =<L =<2048 yields omega=0.25(33). It is consistent with the statement omega =<(gamma-1)/nu, as well as with the value omega = 1/8 of the GFD theory. We reconsider the MC estimation of omega from smaller lattice sizes to show that it does not lead to conclusive results, since the obtained values of omega depend on the particular method chosen. In particular, estimates ranging from omega =1.274(72) to omega=0.18(37) are obtained by four different finite-size scaling methods, using MC data for thermodynamic average quantities, as well as for partition function zeros. We discuss the influence of omega on the estimation of exponents eta and nu.


Introduction
The critical exponents of the three-dimensional (3D) Ising universality class have been a subject of extensive analytical as well as Monte Carlo (MC) studies during many years. The results of the standard perturbative renormalization group (RG) methods are well known [1][2][3][4][5]. An alternative analytical approach has been proposed in [13] and further analyzed in [6], where this approach is called the GFD (Grouping of Feynman Diagrams) theory. A review of MC work till 2001 is provided in [7]. More recent papers are [8][9][10][11][12].
In this paper we will focus on the exponent ω, which describes the leading corrections to scaling. A particular interest in this subject is caused by recent challenging nonperturbative results reported in [6], showing that ω ≤ (γ−1)/ν holds in the ϕ 4 model based on a rigorous proof of certain theorem. The scalar 3D ϕ 4 model belongs to the 3D Ising universality class with (γ − 1)/ν ≈ 0.38. Therefore, ω is expected to be essentially smaller than the values of about 0.8 predicted by standard perturbative methods and currently available MC estimations. The results in [6] are fully consistent with the predictions of the alternative theoretical approach of [13], from which ω = 1/8 is expected. We have performed a Monte Carlo analysis of the standard 3D Ising model, using our data for very large lattice sizes L up to L = 2048, to clarify whether ω, extracted from such data, can be consistent with the results of [6] and [13]. Since our analysis supports this possibility, we have further addressed a related question how a decrease in ω influences the MC estimation of critical exponents η and ν. We have also tested different finite-size scaling methods of estimation ω from smaller lattice sizes to check whether such methods always give ω consistent with 0.832 (6), as one can be expected from the references in [8].
Models with the so-called improved Hamiltonians are often considered instead of the standard Ising model for a better estimation of the critical exponents [8,9]. The basic idea of this approach is to find such Hamiltonian parameters, for which the leading correction to scaling vanishes. However, this correction term has to be large enough and well detectable for the estimation of ω. So, this idea is not very useful in our case.

Analytical arguments
In [6], the ϕ 4 model in the thermodynamic limit has been considered, for which the leading singular part of specific heat C sing V can be expressed as assuming the power-law singularity ξ ∼ t −ν of the correlation length ξ at small reduced temperature t → 0. Here G(k) is the Fourier-transformed two-point correlation function, and G * (k) is its value at the critical point. This expression is valid for any positive Λ ′ < Λ, where Λ is the upper cut-off parameter of the model, since the leading singularity is provided by small wave vectors with the magnitude k =| k |→ 0 and not by the region Λ ′ ≤ k ≤ Λ. In other words, C sing V is independent of the constant Λ ′ . The leading singularity of specific heat in the form of C sing V ∝ (ln ξ) λ ξ α/ν and the two-point correlation function in the asymptotic form of G(k) = ℓ≥0 ξ (γ−θ ℓ )/ν g ℓ (kξ), G * (k) = ℓ≥0 b ℓ k (−γ+θ ℓ )/ν with θ 0 = 0 and θ ℓ > 0 for ℓ ≥ 1 have been considered in [6]. These expressions are consistent with the conventional scaling hypothesis, g ℓ (kξ) being the scaling functions. The exponent λ is responsible for possible logarithmic correction in specific heat, whereas the usual power-law singularity is recovered at λ = 0.
According to the theorem proven in [6], the two-point correlation function of the ϕ 4 model contains a correction with the exponent θ ℓ = γ+1−α−dν, if C sing V can be calculated from (1), applying the considered here scaling forms, if the result is Λ ′ -independent, and if the condition γ + 1 − α − dν > 0 is satisfied for the critical exponents. Applying the known hyperscaling hypothesis α + dν = 2, it yields θ ℓ = γ − 1 for γ > 1. Apparently, the listed here conditions of the theorem are satisfied for the scalar 3D ϕ 4 model. Since the critical singularities are provided by long-wave fluctuations, the condition of Λ ′ -independence is generally meaningful. The assumption ξ ∼ t −ν (with no logarithmic correction) and the considered here scaling forms (with λ = 0), as well as the relation γ + 1 − α − dν > 0 (or γ > 1 according to the hyperscaling hypothesis) are correct for the scalar 3D ϕ 4 model, according to the current knowledge about the critical phenomena.
The correction with the exponent θ ℓ corresponds to the one with ω ℓ = θ ℓ /ν in the finite-size scaling. In such a way, the discussed here analytical arguments predict the existence of a finite-size correction with the exponent (γ − 1)/ν in the scalar 3D ϕ 4 model. As discussed in [6], nontrivial corrections tend to be cancelled in the 2D Ising model, in such a way that only trivial ones with integer θ ℓ are usually observed. However, there is no reason to assume such a scenario in the 3D case. Therefore, the existence of corrections with the exponent (γ − 1)/ν is expected in the 3D Ising model, since it belongs to the same universality class as the 3D ϕ 4 model. Because this correction is not necessarily the leading one, the prediction is ω ≤ ω max , where ω max = (γ − 1)/ν is the upper bond for the leading correction-to-scaling exponent ω. Using the widely accepted estimates γ ≈ 1.24 and ν ≈ 0.63 [3] for the 3D Ising model, we obtain ω max ≈ 0.38. The prediction of the GFD theory [13] is γ = 5/4, ν = 2/3 and, therefore, ω max = 0.375. Thus, we can state that in any case ω max is about 0.38. The value of ω is expected to be 1/8 according to the GFD theory considered in [6,13].
3 MC estimation of ω from finite-size scaling 3.1 The case of very large lattice sizes L ≤ 2048 We have simulated the 3D Ising model on simple cubic lattice with periodic boundary conditions. The Hamiltonian H of the model is given by where T is the temperature measured in energy units, β is the coupling constant and ij denotes the pairs of neighboring spins σ i = ±1. The MC simulations have been performed with the Wolff single cluster algorithm [14], using its parallel implementation described in [15]. An iterative method, introduced in [15], has been used here to find pseudocritical couplings β c (L) corresponding to certain value U = 1.6 of the ratio U = m 4 / m 2 2 , where m is the magnetization per spin. We have evaluated by this method the susceptibility χ = L 3 m 2 an the derivative ∂Q/∂β at β = β c (L), where Q = 1/U . The results for 16 ≤ L ≤ 1536 are already reported in Tab. 1 of our earlier paper [11]. We have extended the simulations to lattice sizes L = 1728 and L = 2048, using approximately the same number of MC sweeps as for L = 1536 in [11]. Thus, Tab. 1 of [11] can be now completed with the new results presented in Tab. 1 here.
The exponent ω describes corrections to the asymptotic finite-size scaling. In particular, for the susceptibility at β = β c (L) we have We define the effective exponent η eff (L) as the mean slope of the − ln χ vs ln L plot, evaluated by fitting the data within [L/2, 2L]. It behaves asymptotically as η eff (L) = Table 1: The values of β c , as well as χ/L 2 , and 10 −3 ∂Q/∂β at β = β c depending on L. η + O (L −ω ). It has been mentioned in [11] that ω might be as small as 1/8, since the plot of the effective exponent η eff vs L −1/8 looks rather linear for large lattice sizes (see Fig. 6 in [11]). This observation is confirmed also by the extended here data, as it can be seen from Fig. 1 (left). An estimate of ω can be obtained by fitting the η eff (L) data. Here we use a more direct method, which gives similar, but slightly more accurate results. We consider the The correction amplitude B is larger for a larger b value, whereas a smaller b value allows us to obtain more data points for Φ b (L). The actual choice b = 2 is found to be optimal for our data. Like the η eff (L) vs L −1/8 plot, also the Φ 2 (L) vs L −1/8 plot can be well approximated by a straight line for large enough lattice sizes, as shown in Fig. 1. Thus, ω could be as small as 1/8. We have fit the quantity Φ 2 (L) to (4) within L ∈ [L min , 1024] (estimated from the χ/L 2 data within L ∈ [L min /2, 2048]) to evaluate ω. The results are collected in Tab. 2. The estimated ω values are essentially decreased for L min ≥ 80 as compared to smaller L min values. Moreover, the quality of fits is remarkably improved in this case, i. e., the values of χ 2 of the fit per degree of freedom (χ 2 /d.o.f.) become smaller. Note that L min = 80 corresponds to the fit interval for Φ 2 (L) in Fig. 1, where the data are well consistent with ω = 1/8. From a formal point of view, ω = 0.25(33) at L min = 80 can be considered as the best estimate from our data, since it perfectly agrees with the results for L min > 80 and has the minimal statistical error within L min ≥ 80. The estimate ω = 0.06(38) at L min = 96 most clearly shows the deviation below the usually accepted values at about 0.8, e. g., ω = 0.832(6) reported in [8]. Our estimation is fully consistent with the analytical arguments in Sec. 2, since all our ω values for L min ≥ 80 are smaller than ω max ≈ 0.38 and also well agree with 1/8. Unfortunately, the statistical accuracy of this estimation is too low to rule out a possibility that the dropping of ω to smaller values at L min ≥ 80 is caused by statistical errors in the data. However, the decrease in ω for large enough lattice sizes is strongly supported by the theorem discussed in Sec. 2. Note also that the recent MC analysis of the 2D ϕ 4 model [16] is consistent with this theorem. These facts make our MC estimation plausible.
Note that there exist many quantities, which scale asymptotically as A + BL −ω with different values of coefficients A and B -see, e. g., [8,17], as well as the examples in the next section. In principle, all of them can be used to estimate ω. However, it is possible that the leading correction term BL −ω for a subset of such quantities is too small as compared to statistical errors and, therefore, it is not well detectable at large lattice sizes. It means that a correction with small ω of about 1/8, probably, will not be detected by MC analysis in many cases, but this still does not imply that such a correction does not exist. Thus, it is sufficient to demonstrate clearly that such a correction exists in one of the cases. Our MC analysis shows that Φ 2 (L) is an appropriate quantity where corrections, i. e., variations in Φ 2 (L), are well detectable even for very large values of L. Moreover, it suggests that a correction with such small ω as 1/8, very likely, exists here.

Different estimates from the data for smaller lattice sizes
The quality of fits with L min < 80 is remarkably improved if 5 data points for the largest lattice sizes are discarded, i. e., if Φ 2 (L) is fit within L ∈ [L min , 432]. Choosing also not too large values of L min , we obtain formally quite good (provided by good fits with sufficiently small χ 2 /d.o.f. values) and stable estimates from remarkably smaller lattice sizes than those in Sec. 3.1. These are presented in Tab. 3. The estimate ω = 1.171(96) at L min = 32 is accepted as the best one from this reduced data set, since it perfectly agrees with the results for L min > 32 and has the smallest statistical error. We have tested another finite-size scaling method. Based on our simulations discussed in [12], we have evaluated U = U (L) at the pseudocritical couplingβ c (L), corresponding to the maximum of specific heat C V . It scales as at large L. This method is similar in spirit to the one used by Hasenbusch for the 3D ϕ 4 model in [17]. The only difference is that another pseudocritical coupling (corresponding to certain value of Z a /Z p , where Z p and Z a are partition functions for the lattice with periodic and antiperiodic boundary conditions) has been used in [17]. We have found that our U (L) data provide a good fit to (5) within 8 ≤ L ≤ 384, where L = 384 is similar to the maximal size L = 360 simulated in [8]. These data are listed in Tab. 4, and the fit results are presented in Tab. 5. The estimate ω = 1.247(73) at L min = 8 seems to be the best one, as it has the smallest statistical error, a good fit quality, and it perfectly agrees with the results for L min > 8. This value disagrees (the discrepancy is 5.7 standard deviations) with the best estimate ω = 0.832(6) of [8], obtained by a different finite-size scaling method. It well agrees with the other value ω = 1.171(96) reported here. Searching for a different method, we have evaluated the Fisher zeros of partition function from MC simulations by the Wolff single cluster algorithm, following the method described in [10]. The results for 4 ≤ L ≤ 72 have been reported in [10]. We have performed high statistics simulations (with MC measurements after each max{2, L/4} Wolff clusters, omitting 10 6 measurements from the beginning of each simulation run, and totally 5 × 10 8 measurements used in the analysis for each L) for 4 ≤ L ≤ 128. Two different pseudo-random number generators, discussed and tested in [11], have been used to verify that the results agree within error bars of about one or, sometimes, two standard deviations. Considering β = η + iξ as a complex number, the results for the first Fisher zero Re u (1) + i Im u (1) in terms of u = exp(−4β) are reported in Tab. 6. Our values are obtained, evaluating R = cos(ξE) η + i sin(ξE) η (where E is energy) by the histogram reweighting method and minimizing | R | (see [10]). Reliable results are ensured by the fact that, for each L, the simulation is performed at the coupling β sim which is close to Reβ (1) -see Tab. 6. We have reached it by using the results of [10] and finite-size extrapolations. We have also estimated the second zeros for L = 4, 32, 64 from different simulation runs -see Tab. 7.
Our results in Tabs. 6 and 7 are reasonably consistent with those of [10], but are more accurate and include larger lattice sizes. Like in [10], the results for the second zeros are much less accurate than those for the first zeros. Therefore only the latter ones are used here in the analysis, considering the ratios       The ratios Ψ 1 (L) = Im u (1) (L)/(Re u (1) (L) − u c ) and Ψ 2 (L) = Im u (1) (L)/Im u (1) (L/2) depending on the lattice size L.
Here u c = exp(−4β c ) is the critical u value, corresponding to the critical coupling β c . The estimation of correction-to-scaling exponent ω from fits of Ψ 1 (L) to A + BL −ω has been considered in [10], assuming the known approximate value 0.2216546 of β c . The use of Ψ 2 (L) instead of Ψ 1 (L) is another method, which has an advantage that it does not require the knowledge of the critical coupling β c . However, a disadvantage is that the data for two sizes, L and L/2, are necessary for one value of Ψ 2 (L). The values of Ψ 1 (L) and Ψ 2 (L) are listed in Tab. 8. The standard errors of Ψ 1 (L) are calculated by the jackknife method [18], thus taking into account the statistical correlations between Re u (1) and Im u (1) . As in [10], the errors due to the uncertainty in β c are ignored, assuming that β c = 0.2216546 holds with a high enough accuracy. According to [11], this β c value, likely, is correct within error bars of about ±3 × 10 −8 . It justifies the actual estimation.
The values of exponent ω, extracted from the fits of Ψ 1 (L) to A + BL −ω within L ∈ [L min , L max ], are collected in Tab. 9. The results ω = 0.903(60) for L ∈ [12,64] and ω = 0.84(10) for L ∈ [16, 64] agree within error bars with the results for similar fit intervals in [10], i. e., ω = 0.77(9) for L ∈ [12,72] and ω = 0.63 (16) for L ∈ [16,72]. However, the fits with L max = 128 are preferable for a reasonable estimation of the asymptotic exponent It turns out that the estimated value of ω becomes larger when L max is increased from 64 to 128. One of possible explanations, which is consistent with the data in Tab. 8, is such that the Ψ 1 (L) plot has a minimum near L = 128 or at somewhat larger L values. In this case the actual method is really valid only for remarkably larger lattice sizes. The results of the other method, using the ratio Ψ 2 (L) instead of Ψ 1 (L), are collected in Tab. 10. Here L max = 128 is fixed and only L min is varied. The standard errors of ω are calculated, taking into account that fluctuations in Im u (1) are the statistically independent quantities. As we can see, the estimated exponent ω decreases with increasing of L min in the considered range. Since χ 2 /d.o.f. is about unity for moderately good fits, the estimates ω = 0.61 (19) at L min = 16 and ω = 0.18(37) at L min = 24 are acceptable.
Summarizing the results of this section, we conclude that three of the considered here methods give larger values of ω (1.171(96), 1.274(72) and 0.997(42)) than ω = 0.832 (6) reported in [8], whereas the fourth method tends to give smaller values (0.61 (19) and 0.18(37)). Thus, it is evident that the estimation of ω from finite-size scaling, using the data for not too large lattice sizes (comparable with L ≤ 360 in [8]

Influence of ω on the estimation of exponents η and ν
Allowing a possibility that the correction-to-scaling exponent ω of the 3D Ising model is, indeed, essentially smaller than the commonly accepted values of about 0.8, we have tested the influence of ω on the estimation of critical exponents η and ν (or 1/ν). We have fit our susceptibility data at β = β c (L) to the ansatz with m = 1 and m = 2 to estimate η at three fixed values of the exponent ω, i. e., ω = 0.8, ω = 0.38 and ω = 1/8. The first one is very close to the known RG value ω = 0.799 ± 0.011 [3] and also is quite similar to a more recent RG estimate ω = 0.782(5) [19] and the MC estimate ω = 0.832(6) of [8]. The second value corresponds to the upper bound ω max ≈ 0.38 stated in Sec. 2, and the third value 1/8 is extracted from the GFD theory [6,13]. There exist different corrections to scaling, but the two correction terms in (6) are the most relevant ones at L → ∞, as it can be seen from the analysis in [8,12]. The results of the fit within L ∈ [L min , 2048] depending on ω, m and L min are shown in Tab. 11. Similarly, we have fit our ∂Q/∂β data at β = β c (L) to the ansatz and have presented the results in Tab. 12.
Considering the fits with only the leading correction to scaling included (m = 1), one can conclude from Tab. 11 that the estimated critical exponent η increases with decreasing of ω, whereas the exponent 1/ν in Tab. 12 is rather stable. The sub-leading correction to scaling (m = 2) makes the estimated exponents η and 1/ν remarkably less stable for small ω values, such as ω = 1/8. The latter value is expected from the GFD theory [6,13], so that the estimation at ω = 1/8 is self-consistent within this approach. In this case, the estimation of η appears to be compatible with the theoretical value η = 1/8 = 0.125 of [13], taking into account that the evaluated η increases with L min . Moreover, the selfconsistent estimation of 1/ν is even very well consistent with ν = 2/3 predicted in [13]. In particular, 1/ν = 1.525(50) can be considered as the best estimate at ω = 1/8 and m = 2 (it has the smallest χ 2 /d.o.f. value and much smaller statistical error than the estimate   (5) at L min = 64), which well agrees with 1.5. Thus, contrary to the statements in [10], the value ν = 2/3 of the GFD theory is not ruled out, since it is possible that ω has a much smaller value than 0.832(6) assumed in [10].
A question can arise about the influence of ω value on the estimation of critical exponents in the case of improved Hamiltonians [8,9]. It is expected that the leading corrections to scaling vanish in this case, and therefore the influence of ω is small. However, in the case if the asymptotic corrections to scaling are described by the exponent ω ≤ ω max ≈ 0.38, as it is strongly suggested by the theorem discussed in Sec. 2, the vanishing of leading corrections cannot be supported by the existing MC analyses of such models. Indeed, in these analyses the asymptotic corrections to scaling are not correctly identified (probably, because of too small lattice sizes) if ω ≤ ω max ≈ 0.38, since one finds that ω ≈ 0.8.

Comparison of recent results
It is interesting to compare our MC estimates and those of [8] with the most recent RG (3D expansion) values of [19] cited in [8]. Note that the estimates of ω in [8] and [19], i. e., ω = 0.832 (6) and ω = 0.782 (5), are clearly inconsistent within the claimed error bars. This discrepancy, however, can be understood from the point of view of our MC analysis, suggesting that the real uncertainty in the MC estimation of ω can be rather large.
The comparison of critical exponents η and ν is provided in Tab. 13. This comparison includes only some recent or relatively new results, since older ones are extensively discussed in [7][8][9]. Our values correspond to the fits within L ∈ [64, 2048] at m = 2 and ω = 0.8. The choice of ω = 0.8 is reasonable here, since this value is close enough to the above mentioned estimates ω = 0.832 (6) and ω = 0.782 (5), and practically the same results are obtained if ω = 0.8 is replaced by any of these two values. According to the claimed statistical error bars, the estimates of [8] seem to be extremely accurate. Note, however, that these estimates are extracted from much smaller lattice sizes (L ≤ 360) as compared to ours (L ≤ 2048).
The values of ν in Tab. 13 are consistent with each other. The MC estimates of η are consistent, as well. However, the recent RG value of [19] appears to be somewhat smaller and not consistent within the error bars with the actual MC estimations, even if the assumed values of ω are about 0.8, as predicted by the perturbative RG theory. In particular, the discrepancy with the MC value of [8] is about 45 standard deviations of the MC estimation or about 15 error bars of the RG estimation.
Recently, the conformald field theory (CFT) has been applied to the 3D Ising model [20] to obtain very accurate values of the critical exponents, using the numerical conformal bootstrap method. The conformal-symmetry relations for the correlation functions, like (2.1) in [20], are known to hold asymptotically in two dimensions, whereas their validity in 3D case can be questioned. Here "asymptotically" means that the limit L/x → ∞, ξ/x → ∞ and a/x → 0 is considered, where x is the actual distance, a is the lattice spacing and ξ is the correlation length. In other words, the conformal symmetry is expected to hold exactly for the asymptotic correlation functions on an infinite lattice (L = ∞) at the critical point (β = β c ). These asymptotic correlation functions are obtained by subtracting from the exact correlations functions (at L = ∞ and β = β c ) the corrections to scaling, containing powers of a/x. The existence of the conformal symmetry in the 3D Ising model has been supported by a non-trivial MC test in [21].
Apart from the assumption of the validity of (2.1) in [20] for the 3D Ising model, the following hypotheses have been proposed: (i) There exists a sharp kink on the border of the two-dimensional region of the allowed values of the operator dimensions ∆ σ = (1 + η)/2 and ∆ ǫ = 3 − 1/ν; (ii) Critical exponents of the 3D Ising model correspond just to this kink.
These hypotheses have been supported by the MC estimates of the exponents η, ν and ω in [8]. However, the obtained in [20] exponent ω = 0.8303 (18) is not supported by our MC value ω = 0.25(33), obtained from the susceptibility data for very large lattice sizes L ≤ 2048. Moreover, it does not satisfy the inequality ω ≤ (γ − 1)/ν, following from the theorem discussed in Sec. 2. This apparent contradiction can be understood from the point of view that corrections to scaling are not fully controlled in the CFT. Indeed, the prediction for ω in this CFT is based on the assumption that ω = ∆ ǫ ′ − 3 holds, where ∆ ǫ ′ is the dimension of an irrelevant operator in the conformal analysis of the asymptotic four-point correlation function. It means that corrections to scaling of the exact fourpoint correlation function are discarded (to obtain the asymptotic correlation function, as discussed before) and not included into the analysis.
More recently, a modified conformal bootstrap analysis has been performed in [22], where the mentioned two hypotheses have been replaced with the hypothesis that the operator product expansion (OPE) contains only two relevant scalar operators. The results for the exponents ∆ σ and ∆ ǫ (or η and ν) are consistent with those of [20]. This consistency is not surprising, since both methods agree with the idea that the operator spectrum of the 3D Ising model is relatively simple, so that the true values of ∆ σ and ∆ ǫ are located inside of a certain narrow region (as in [22]) or on its border (as in [20]), where many operators are decoupled from the spectrum. Apparently, the analysis in [22] does not lead to a contradiction with the two relations ω ≤ (γ − 1)/ν (the theorem) and ω = ∆ ǫ ′ − 3, since only ∆ ǫ ′ > 3 is assumed for the dimension ∆ ǫ ′ . Thus, both relations can be satisfied simultaneously, if the 3D Ising point in Fig. 1 of [20] is located inside of the allowed region, rather than on its border. This possibility is supported by the behavior of the effective exponent η eff in our Fig. 1. It suggests that the asymptotic exponent η could be larger than it is usually expected from MC simulations for relatively small lattice sizes, as in [8]. Thus, it is important to make further refined estimations, based on MC data for very large lattice sizes, in order to verify the hypotheses proposed in [20,22].
If the hypothesis (i) about the existence of a sharp kink is true, then this kink, probably, has a special meaning for the 3D Ising model. Its existence, however, is not evident. According to the conjectures of [20], such a kink is formed at N → ∞, where N is the number of derivatives included into the analysis. As discussed in [20], it implies that the minimum in the plots of Fig. 7 in [20] should be merged with the apparent "kink" at N → ∞. This "kink" is not really sharp at a finite N . Nevertheless, its location can be identified with the value of ∆(σ), at which the second derivative of the plot has a local maximum. The minimum of the plot is slightly varied with N , whereas the "kink" is barely moving [20]. Apparently, the convergence to a certain asymptotic curve is remarkably faster than 1/N , as it can be expected from Fig. 7 and other similar figures in [20]. In particular, we have found that the location of the minimum in Fig. 7 of [20] is varied almost linearly with N −3 . We have shown it in Fig. 2 by solid circles, the position of the "kink" being indicated by empty circles. The error bars of ±0.000001 correspond to the symbol size. The results for N = 153, 190, 231 are presented, skipping the estimate for the location of the "kink" at N = 153, which cannot be well determined from the corresponding plot in Fig. 7 of [20]. The linear extrapolations (dashed lines) suggest that the minimum, very likely, is moved only slightly closer to the "kink" when N is varied from N = 231 to N = ∞. The linear extrapolation might be too inaccurate. Only in this case a refined numerical analysis for larger N values can possibly confirm the hypothesis about the formation of a sharp kink at N → ∞.
The results of both [20] and [22] strongly support the commonly accepted 3D Ising values of the critical exponents η and ν. In particular, the estimates η = 0.03631(3) and ν = 0.62999 (5) have been reported in [20]. However, these estimates are obtained, based on certain hypotheses. If these hypotheses are not used, then the conformal bootstrap analysis appears to be consistent even with the discussed here GFD values η = 1/8 and ν = 2/3. Indeed, the corresponding operator dimensions ∆ σ = (1 + η)/2 and ∆ ǫ = 3 − 1/ν lie inside of the allowed region in Fig. 1 of [20].
The hypotheses (i) and (ii) can be questioned in view of the observations summarized in Fig. 2. The hypothesis of [22] about the existence of just two relevant scalar operators might be supported by some physically-intuitive arguments. In particular, one needs to adjust two scalar parameters P (pressure) and T (temperature) to reach the critical point of a liquid-vapor system. A real support for this hypothesis is provided by the already known estimations of the critical exponents. Taking into account the non-perturbative nature of the critical phenomena, the most reliable estimates are based on non-perturbative methods, such as the Monte Carlo simulation. An essential point in this discussion is that the MC estimates can be remarkably changed, if unusually large lattices are considered, as it is shown in our current study.

Summary and conclusions
Analytical as well as Monte Carlo arguments are provided in this paper, showing that corrections to scaling in the 3D Ising model are described by a remarkably smaller exponent ω than the usually accepted values of about 0.8. The analytical arguments in Sec. 2, which are based on a rigorous proof of certain theorem, suggest that ω ≤ (γ−1)/ν holds, implying that ω cannot be larger than ω max = (γ−1)/ν ≈ 0.38 in the 3D Ising model. The analytical prediction of the GFD theory [6,13] is ω = 1/8 in this case. Our MC estimation of ω from the susceptibility (χ) data of very large lattices (Sec. 3.1) is well consistent with these analytical results. Numerical values, extracted from the χ data within 40 ≤ L ≤ 2048 and 48 ≤ L ≤ 2048 (or Φ 2 (L) = 2 −4 χ(2L)/χ(L/2) data within 80 ≤ L ≤ 1024 and 96 ≤ L ≤ 1024) are ω = 0.25(33) and ω = 0.06(38), respectively. Unfortunately, the statistical errors in ω are rather large.
As discussed in [16], our analytical predictions generally refer to a subset of n-vector models, where spin is an n-component vector with n = 1 in two dimensions and n ≥ 1 in three dimensions. Our recent MC analysis agrees with these predictions for the scalar (n = 1) 2D ϕ 4 model [16], where statistical errors are small enough. The 3D case with n = 2 has been tested in [6], based on accurate experimental data for specific heat in zero gravity conditions very close to the λ-transition point in liquid helium. The test in Sec. 4 of [6] reveals some inconsistency of the data with corrections to scaling proposed by the perturbative RG treatments, indicating that these corrections decay slower, i. e., θ = νω is smaller than usually expected. This finding is consistent with the theorem discussed in Sec. 2. The mentioned here facts emphasize the importance of our MC analysis.
Our proposed values of ω may seem to be incredible in view of a series of known results, yielding ω at about 0.8 for the 3D Ising model. However, it is meaningful to reconsider these results from several aspects.
• First of all, they disagree with non-perturbative arguments in the form of the rigorously proven theorem, discussed in Sec. 2.
• This theorem states that ω ≤ (γ − 1)/ν, whereas the perturbative RG estimates are essentially larger. In view of the recent analysis in [23] (see also the discussions in [6]), this discrepancy can be understood as a failure of the standard perturbative RG methods. The actually discussed (Sec. 5) discrepancy between the recent RG and MC estimates of the critical exponent η also points to problems in the perturbative approach. Moreover, any perturbative method, also the high-and low-temperature series expansions, can fail to give correct results in critical phenomena, since it is not the natural domain of validity of the perturbation theory.
• The previous MC estimations of ω are based on simulations of lattices not larger than L ≤ 360 in [8]. We have clearly demonstrated in Sec. 3.2 that the values obtained from finite-size scaling with such relatively small (as compared to L ≤ 2048 in our study) lattice sizes depend on the particular method chosen. For example, different estimates ranging from ω = 1.274(72) to ω = 0.18(37) are obtained here, which substantially deviate from the usually reported values between 0.82 and 0.87 (see [7,8]).
Taking into account the possibility that the correction-to-scaling exponent ω can be remarkably smaller than the usually accepted values at about 0.8, we have tested in Sec. 4 the influence of ω on the estimation of critical exponents η and ν. We have concluded that the effect is remarkable if ω is changed from 0.8 to a much smaller value, such as ω = 1/8 of the GFD theory. In this case, the error bars strongly increase, and the estimation becomes compatible, or even well consistent, with the predictions of the GFD theory. In particular, the estimate 1/ν = 1.525(50) agrees with the GFD theoretical value 1.5.