We report the results of a Monte Carlo investigation of the (equilibrium) time-displaced correlation functions for the magnetization and energy of a simple cubic Ising model as a function of time, temperature, and lattice size. The simulations were carried out on a CDC CYBER 205 supercomputer employing a high-speed, vectorized multispin coding program and using a total of 5\ifmmode\times\else\texttimes\fi{}1012 Monte Carlo spin-flip trials. We used L\ifmmode\times\else\texttimes\fi{}L\ifmmode\times\else\texttimes\fi{}L lattices with periodic boundary conditions and L as large as 96. The short-time and long-time behaviors of the correlation functions are analyzed by fits to a sum of exponential decays, and the critical exponent z for the largest relaxation time is extracted using a finite-size-scaling analysis. Our estimate z=2.04\ifmmode\pm\else\textpm\fi{}0.03 resolves an intriguing contradiction in the literature; it satisfies the theoretical lower bound and is in agreement with the prediction obtained by \ensuremath{\epsilon} expansion. We also consider various small systematic errors that typically occur in the analysis of relaxation functions and show how they can lead to spurious results if sufficient care is not exercised.
This dataset and code repository accompanies the article of the same title (version 5, May 2026) and provides all simulation data, source code, and raw terminal output necessary to reproduce and verify the results. Study summary We investigate the critical dynamics of the two-dimensional Ising model at the exact critical point (Tc ≈ 2.269185) using the local Metropolis Monte Carlo algorithm. The central focus is the characterization and comparison of two widely used correlation time estimators: the exponential decay criterion τ₁/e and the integrated autocorrelation time τ_int, following the convention of Sokal (1997). Simulations were performed for seven system sizes L ∈ {16, 24, 32, 48, 64, 96, 128} with up to 50 independent replicas per system size, using both Python (NumPy/PCG64) and C++ (std::mt19937, compiled with -O2) implementations. The two implementations were developed independently, enabling cross-validation of results across languages. Main findings The coefficient of variation (CV = std/mean) of τ_int converges to a physical floor of approximately 0.44–0.52 regardless of the number of replicas (8 or 50) or programming language, confirming that this floor is determined by the physics of the estimator under local Metropolis dynamics at criticality — not by statistical sampling limitations. The distribution of τ₁/e is genuinely bimodal at the critical point: individual replicas alternate between low (~200–500 MC steps) and high (~3000–5000 MC steps) correlation times for the same system size, with mean-to-median ratios of up to 2.6× and bimodal ratios (mean of upper half / mean of lower half) of up to 23.5×. This bimodality does not disappear with more replicas, confirming it is a physical property of the estimator. No detectable power-law scaling τ ~ L^z is found for L ≤ 128 in any of the six independent runs (R² < 0.14 in all log-log fits using all 7 system sizes). The accepted asymptotic value z ≈ 2.17 (Li & Sokal, 1989) implies an expected τ_int of approximately 73,600 MC steps at L = 128 — 76× larger than the observed value of 972. This quantifies the dominance of finite-size corrections in the range studied, and justifies treating the effective dynamic exponent z_eff as a continuously varying quantity rather than a fixed constant. The C++ implementation achieves speedups of 98× to 782× relative to Python depending on system size, enabling the full 50-replica run across all seven system sizes to complete in approximately 33 minutes. Cross-validation caveat Cross-validation between Python and C++ is fully valid for τ₁/e. For τ_int, the Python scripts used a fixed equilibration of n_equil = 3000 steps for all L — documentably insufficient for L ≥ 24 where τ_int > 600 (n_equil/τ < 3–5τ). This generates systematic underestimation of up to 86% at L = 64. The valid cross-validation for τ_int is between the C++-8 and C++-50 runs, where CV agreement (~0.52 in both) confirms the physical floor is robust. Connection to the EIE/ZIE framework This study provides the empirical foundation for the EIE/ZIE adaptive meta-control framework (Patent BR 10 2025 024459 4), justifying three architectural decisions: (1) τ_int as the preferred estimator of the correlation time τ_corr; (2) the effective dynamic exponent z_eff treated as a continuously estimated dynamic variable rather than a fixed constant; and (3) the dynamic freezing time τ_freeze introduced as an observable independent of τ_corr, motivated by the bimodal behavior of τ₁/e which reveals that the system alternates between regimes of high and low correlation persistence — a phenomenon invisible to conventional autocorrelation-based diagnostics. Repository contents Cpp/ — C++ source code (3 scripts) and simulation data (6 CSV files + raw terminal output) Python/ — Python source code (2 scripts) and simulation data (4 CSV files + raw screen output) Paper/ — Article in English (PDF) and Portuguese (PDF) README.md — Full repository documentation including file descriptions, column definitions, reproduction instructions, and physical parameters LICENSE.txt — CC BY 4.0 for data; MIT for code; patent notice Joemerson da Silva Lima — Independent Researcher, Cuiabá, Mato Grosso, Brazil joemerson93lima@gmail.com Patent BR 10 2025 024459 4
Extending the Monte Carlo method to dynamic critical phenomena we investigated the time-dependent correlation functions in the two-dimensional one-spin-flip Ising model and the critical behavior of the associated relaxation times. These relaxation times are the following: ${\ensuremath{\tau}}_{\ensuremath{\delta}\ensuremath{\mu}}^{\ensuremath{\Delta}T}$, characterizing the approach of the order parameter to equilibrium after a change of temperature $\ensuremath{\Delta}T$ of the system; ${\ensuremath{\tau}}_{\ensuremath{\delta}\ensuremath{\mu}\ensuremath{\delta}\ensuremath{\mu}}$ and ${\ensuremath{\tau}}_{\ensuremath{\delta}\ensuremath{\mu}\ensuremath{\delta}\ensuremath{\mu}}^{A}$ characterizing the slowing down of the order-parameter correlation and autocorrelation functions, respectively; ${\ensuremath{\tau}}_{\ensuremath{\delta}\mathcal{H}\ensuremath{\delta}\mathcal{H}}$ and ${\ensuremath{\tau}}_{\ensuremath{\delta}\mathcal{H}\ensuremath{\delta}\mathcal{H}}^{A}$, characterizing the slowing down of the energy correlation and autocorrelation functions; and finally ${\ensuremath{\tau}}_{\ensuremath{\delta}\ensuremath{\mu}\ensuremath{\delta}\mathcal{H}}$, characterizing the cross-correlation function. We give estimates for the associated exponents ${\ensuremath{\Delta}}_{\ensuremath{\delta}\ensuremath{\mu}}^{\ensuremath{\Delta}T}\ensuremath{\approx}{\ensuremath{\Delta}}_{\ensuremath{\delta}\ensuremath{\mu}\ensuremath{\delta}\ensuremath{\mu}}\ensuremath{\approx}{\ensuremath{\Delta}}_{\ensuremath{\delta}\mathcal{H}\ensuremath{\delta}\mathcal{H}}\ensuremath{\approx}{\ensuremath{\Delta}}_{\ensuremath{\delta}\ensuremath{\mu}\ensuremath{\delta}\mathcal{H}}\ensuremath{\approx}1.90\ifmmode\pm\else\textpm\fi{}0.10$, and ${\ensuremath{\Delta}}_{\ensuremath{\delta}\ensuremath{\mu}\ensuremath{\delta}\ensuremath{\mu}}^{A}\ensuremath{\approx}1.60\ifmmode\pm\else\textpm\fi{}0.10$, ${\ensuremath{\Delta}}_{\ensuremath{\delta}\ensuremath{\mu}\ensuremath{\delta}\mathcal{H}}^{A}\ensuremath{\approx}0.95\ifmmode\pm\else\textpm\fi{}0.10$, ${\ensuremath{\Delta}}_{\ensuremath{\delta}\mathcal{H}\ensuremath{\delta}\mathcal{H}}^{A}\ensuremath{\approx}0$, which are consistent with the dynamic scaling hypothesis and with exact inequalities. A detailed comparison with recent high-temperature-expansion estimates is performed, and the reliability of the Monte Carlo results is carefully analyzed.
Critical dynamics in two-dimension Ising lattices up to 2048×2048 is simulated on field-programmable-gate-array- based computing devices. Linear relaxation times are measured from extremely long Monte Carlo simulations. The longest simulation has 7.1×10(16) spin updates, which would take over 37 years to simulate on a general purpose computer. The linear relaxation time of the Ising lattices is found to follow the dynamic scaling law for correlation lengths as long as 2048. The dynamic exponent z of the system is found to be 2.179(12), which is consistent with previous studies of Ising lattices with shorter correlation lengths. It is also found that Monte Carlo simulations of critical dynamics in Ising lattices larger than 512×512 are very sensitive to the statistical correlations between pseudorandom numbers, making it even more difficult to study such large systems.
Monte Carlo cluster algorithms are popular for their efficiency in studying the Ising model near its critical temperature. We might expect that this efficiency extends to the bond-diluted Ising model. We show, however, that this is not always the case by comparing how the correlation times τ_{w} and τ_{sw} of the Wolff and Swendsen-Wang cluster algorithms scale as a function of the system size L when applied to the two-dimensional bond-diluted Ising model. We demonstrate that the Wolff algorithm suffers from a much longer correlation time than in the pure Ising model, caused by isolated (groups of) spins which are infrequently visited by the algorithm. With a simple argument we prove that these cause the correlation time τ_{w} to be bounded from below by L^{z_{w}} with a dynamical exponent z_{w}=γ/ν≈1.75 for a bond concentration p<1. Furthermore, we numerically show that this lower bound is actually taken for several values of p in the range 0.5<p<1. Moreover, we show that the Swendsen-Wang algorithm does not suffer from the same problem. Consequently, it has a much shorter correlation time, shorter than in the pure Ising model even. Numerically at p=0.6, we find that its dynamical exponent is z_{sw}=0.09(4).