XY Model: Temperature Sweep and BKT Crossover
The two-dimensional XY model places a continuous planar spin \(\vec{s}_i = (\cos\theta_i, \sin\theta_i)\) on every site of a square lattice, coupled through \(H = -J \sum_{\langle i,j \rangle} \cos(\theta_i - \theta_j)\). Unlike the Ising model, the continuous symmetry forbids spontaneous symmetry breaking at any finite temperature (Mermin-Wagner theorem), so there is no conventional ordered phase with a nonzero magnetization [1]. Instead, the system undergoes a Berezinskii-Kosterlitz-Thouless (BKT) transition at \(T_{\mathrm{BKT}} \approx 0.893\) between a low-temperature phase with bound vortex-antivortex pairs and algebraically decaying correlations and a high-temperature disordered phase with exponential correlation decay and free vortices [2] [3].
This notebook runs a thermodynamic temperature sweep through the BKT crossover region and records the same set of observables used in the Ising sweep: magnetization, energy, susceptibility, specific heat, entropy, and integrated autocorrelation time. The BKT transition is not a standard second-order transition, so the thermodynamic signatures are broader and less singular than those of the Ising model. The susceptibility and specific heat show rounded humps rather than sharp divergences, and the entropy curve changes slope gently across a wider temperature interval [3].
[ ]:
from __future__ import annotations
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from models.xy_model import XYSimulation
from utils.exceptions import ZeroVarianceAutocorrelationError
from utils.statistics import (
DEFAULT_CONFIDENCE_LEVEL,
calculate_autocorr,
summarize_derived_observable,
summarize_entropy_observable,
summarize_primary_observable,
)
from utils.equilibration import (
convergence_equilibrate,
estimate_relaxation_time_two_start,
_detect_quasi_steady_stuck,
)
def _find_repo_root() -> Path:
cwd = Path.cwd().resolve()
for candidate in (cwd, *cwd.parents):
if (candidate / 'pyproject.toml').exists() and (candidate / 'results').exists():
return candidate
return cwd
REPO_ROOT = _find_repo_root()
T_BKT = 0.893
CACHE_PATH = REPO_ROOT / 'results' / 'xy' / 'temperature_sweep_data.npz'
PALETTE = {
'mag': '#4878CF',
'eng': '#6ACC65',
'chi': '#EE854A',
'cv': '#D65F5F',
'entropy': '#956CB4',
'tau': '#8C613C',
}
plt.rcParams['figure.dpi'] = 120
print(f'T_BKT (theoretical) = {T_BKT}')
print(f'Repository root: {REPO_ROOT}')
print(f'Primary sweep path: {CACHE_PATH}')
Ordered and Disordered Configurations
The XY model’s phases look quite different from the Ising case. Below \(T_{\mathrm{BKT}}\) the spins form large coherent patches with slowly varying orientation, while above the transition the angle field becomes uncorrelated over short distances. The snapshots below show the local spin angle \(\theta_i\) after equilibration with checkerboard Metropolis updates on a \(48 \times 48\) lattice.
[ ]:
def _equilibrated_state(*, temp: float, size: int = 48, seed: int = 0) -> np.ndarray:
sim_random = XYSimulation(
size=size,
temp=temp,
update='checkerboard',
init_state='random',
seed=seed,
)
sim_ordered = XYSimulation(
size=size,
temp=temp,
update='checkerboard',
init_state='ordered',
seed=seed,
)
convergence_equilibrate(
sim_random,
sim_ordered,
chunk_size=120,
max_steps=6000,
)
sim_random.run(n_steps=200)
return np.array(sim_random.spins, copy=True)
snapshot_temperatures = [0.5 * T_BKT, 1.5 * T_BKT]
snapshot_titles = [r'Below $T_{\mathrm{BKT}}$', r'Above $T_{\mathrm{BKT}}$']
snapshot_states = [
_equilibrated_state(temp=float(temp), seed=200 + idx)
for idx, temp in enumerate(snapshot_temperatures)
]
fig, axes = plt.subplots(1, 2, figsize=(9, 4.5))
for ax, spins, title, temp in zip(axes, snapshot_states, snapshot_titles, snapshot_temperatures):
angles = np.arctan2(spins[:, :, 1], spins[:, :, 0])
ax.imshow(angles, cmap='hsv', vmin=-np.pi, vmax=np.pi, interpolation='nearest')
ax.set_title(f'{title}\n$T = {temp:.3f}$')
ax.set_xticks([])
ax.set_yticks([])
fig.suptitle('Equilibrated checkerboard-Metropolis XY snapshots (spin angle)', y=0.98)
fig.tight_layout()
plt.show()
Load or Compute Sweep Data
The notebook first looks for the project cache at results/xy/temperature_sweep_data.npz. If that file contains a sufficiently dense temperature grid, the plots are built directly from the precomputed data. The cache is produced by scripts/xy/temperature_sweep.py, which uses two-start convergence with quasi-steady stuck detection at low temperature and records standardized uncertainty fields for every observable.
If no valid cache is available, a lightweight fallback sweep runs inline on a small lattice with a reduced number of measurement steps. The fallback results are noisier and cover fewer temperature points, but they are sufficient to show the qualitative shape of each observable curve across the BKT crossover.
[ ]:
REQUIRED_KEYS = {
'temperatures',
'avg_m_value', 'avg_m_ci_low', 'avg_m_ci_high', 'avg_m_n_eff',
'avg_e_value', 'avg_e_ci_low', 'avg_e_ci_high',
'susc_value', 'susc_ci_low', 'susc_ci_high',
'spec_h_value', 'spec_h_err', 'spec_h_ci_low', 'spec_h_ci_high',
'tau_int',
'undefined_autocorr_flag', 'tau_interval_unstable_flag',
}
fallback_meas_steps = 1200
def _run_local_sweep_point(
*,
temp: float,
size: int,
meas_steps: int,
eq_probe_steps: int,
eq_max_steps: int,
seed: int,
) -> dict[str, float]:
sim_random = XYSimulation(
size=size,
temp=temp,
update='checkerboard',
init_state='random',
seed=seed,
)
sim_ordered = XYSimulation(
size=size,
temp=temp,
update='checkerboard',
init_state='ordered',
seed=seed,
)
convergence_equilibrate(
sim_random,
sim_ordered,
chunk_size=eq_probe_steps,
max_steps=eq_max_steps,
)
mags, engs = sim_random.run(n_steps=meas_steps)
mags_arr = np.asarray(mags, dtype=float)
engs_arr = np.asarray(engs, dtype=float)
mag = summarize_primary_observable(
time_series=mags_arr,
confidence=DEFAULT_CONFIDENCE_LEVEL,
)
eng = summarize_primary_observable(
time_series=engs_arr,
confidence=DEFAULT_CONFIDENCE_LEVEL,
)
chi = summarize_derived_observable(
magnetization_series=mags_arr,
temperature=temp,
L=size,
observable='chi',
confidence=DEFAULT_CONFIDENCE_LEVEL,
)
cv = summarize_derived_observable(
energy_series=engs_arr,
temperature=temp,
L=size,
observable='cv',
confidence=DEFAULT_CONFIDENCE_LEVEL,
)
return {
'avg_m_value': float(mag['value']),
'avg_m_ci_low': float(mag['ci_low']),
'avg_m_ci_high': float(mag['ci_high']),
'avg_m_n_eff': float(mag['n_eff']),
'avg_e_value': float(eng['value']),
'avg_e_ci_low': float(eng['ci_low']),
'avg_e_ci_high': float(eng['ci_high']),
'susc_value': float(chi['value']),
'susc_ci_low': float(chi['ci_low']),
'susc_ci_high': float(chi['ci_high']),
'spec_h_value': float(cv['value']),
'spec_h_err': float(cv['err']),
'spec_h_ci_low': float(cv['ci_low']),
'spec_h_ci_high': float(cv['ci_high']),
'tau_int': float(mag['tau_int']),
'tau_int_ci_low': float(mag['ci_low']) if not np.isfinite(mag['tau_int']) else float(max(0.0, mag['tau_int'] - mag['err'])),
'tau_int_ci_high': float(mag['ci_high']) if not np.isfinite(mag['tau_int']) else float(mag['tau_int'] + mag['err']),
}
cache_loaded = False
cache_reason = ''
requested_n_seeds = 1
retained_n_seeds = 1
if CACHE_PATH.exists():
cache = np.load(CACHE_PATH)
if REQUIRED_KEYS.issubset(set(cache.files)) and int(cache['temperatures'].size) >= 10:
temperatures = np.asarray(cache['temperatures'], dtype=float)
avg_m = np.asarray(cache['avg_m_value'], dtype=float)
avg_m_ci_low = np.asarray(cache['avg_m_ci_low'], dtype=float)
avg_m_ci_high = np.asarray(cache['avg_m_ci_high'], dtype=float)
avg_m_n_eff = np.asarray(cache['avg_m_n_eff'], dtype=float)
avg_e = np.asarray(cache['avg_e_value'], dtype=float)
avg_e_ci_low = np.asarray(cache['avg_e_ci_low'], dtype=float)
avg_e_ci_high = np.asarray(cache['avg_e_ci_high'], dtype=float)
susc = np.asarray(cache['susc_value'], dtype=float)
susc_ci_low = np.asarray(cache['susc_ci_low'], dtype=float)
susc_ci_high = np.asarray(cache['susc_ci_high'], dtype=float)
spec_h = np.asarray(cache['spec_h_value'], dtype=float)
spec_h_err = np.asarray(cache['spec_h_err'], dtype=float)
spec_h_ci_low = np.asarray(cache['spec_h_ci_low'], dtype=float)
spec_h_ci_high = np.asarray(cache['spec_h_ci_high'], dtype=float)
tau_int = np.asarray(cache['tau_int'], dtype=float)
tau_int_ci_low = (
np.asarray(cache['tau_int_ci_low'], dtype=float)
if 'tau_int_ci_low' in cache.files
else tau_int
)
tau_int_ci_high = (
np.asarray(cache['tau_int_ci_high'], dtype=float)
if 'tau_int_ci_high' in cache.files
else tau_int
)
undefined_autocorr_flag = np.asarray(cache['undefined_autocorr_flag'], dtype=np.uint8)
tau_interval_unstable_flag = np.asarray(cache['tau_interval_unstable_flag'], dtype=np.uint8)
low_effective_sample_flag = (
np.asarray(cache['low_effective_sample_flag'], dtype=np.uint8)
if 'low_effective_sample_flag' in cache.files
else np.zeros_like(temperatures, dtype=np.uint8)
)
entropy = (
np.asarray(cache['entropy_value'], dtype=float)
if 'entropy_value' in cache.files
else np.asarray(cache['entropy'], dtype=float)
)
entropy_ci_low = (
np.asarray(cache['entropy_ci_low'], dtype=float)
if 'entropy_ci_low' in cache.files
else np.asarray(cache['entropy'], dtype=float)
)
entropy_ci_high = (
np.asarray(cache['entropy_ci_high'], dtype=float)
if 'entropy_ci_high' in cache.files
else np.asarray(cache['entropy'], dtype=float)
)
confidence_level = float(cache['confidence_level']) if 'confidence_level' in cache.files else DEFAULT_CONFIDENCE_LEVEL
requested_n_seeds = int(cache['requested_n_seeds']) if 'requested_n_seeds' in cache.files else int(cache['n_seeds']) if 'n_seeds' in cache.files else 1
retained_n_seeds = int(cache['retained_n_seeds']) if 'retained_n_seeds' in cache.files else int(cache['n_seeds']) if 'n_seeds' in cache.files else 1
n_seeds = retained_n_seeds
data_source = (
f'Loaded cached sweep from {CACHE_PATH} '
f'using {retained_n_seeds} retained seeds '
f'out of {requested_n_seeds} requested.'
)
cache_loaded = True
else:
n_cached = int(cache['temperatures'].size) if 'temperatures' in cache.files else 0
cache_reason = f'cache exists but only provides {n_cached} temperature points'
else:
cache_reason = 'cache file not found'
if not cache_loaded:
temperatures = np.linspace(0.2, 1.8, 15)
avg_m = np.empty_like(temperatures)
avg_m_ci_low = np.empty_like(temperatures)
avg_m_ci_high = np.empty_like(temperatures)
avg_m_n_eff = np.empty_like(temperatures)
avg_e = np.empty_like(temperatures)
avg_e_ci_low = np.empty_like(temperatures)
avg_e_ci_high = np.empty_like(temperatures)
susc = np.empty_like(temperatures)
susc_ci_low = np.empty_like(temperatures)
susc_ci_high = np.empty_like(temperatures)
spec_h = np.empty_like(temperatures)
spec_h_err = np.empty_like(temperatures)
spec_h_ci_low = np.empty_like(temperatures)
spec_h_ci_high = np.empty_like(temperatures)
tau_int = np.empty_like(temperatures)
tau_int_ci_low = np.empty_like(temperatures)
tau_int_ci_high = np.empty_like(temperatures)
spec_h_samples = np.empty((temperatures.size, 1), dtype=float)
for idx, temp in enumerate(temperatures):
point = _run_local_sweep_point(
temp=float(temp),
size=32,
meas_steps=fallback_meas_steps,
eq_probe_steps=150,
eq_max_steps=12_000,
seed=1_000 + idx,
)
avg_m[idx] = point['avg_m_value']
avg_m_ci_low[idx] = point['avg_m_ci_low']
avg_m_ci_high[idx] = point['avg_m_ci_high']
avg_m_n_eff[idx] = point['avg_m_n_eff']
avg_e[idx] = point['avg_e_value']
avg_e_ci_low[idx] = point['avg_e_ci_low']
avg_e_ci_high[idx] = point['avg_e_ci_high']
susc[idx] = point['susc_value']
susc_ci_low[idx] = point['susc_ci_low']
susc_ci_high[idx] = point['susc_ci_high']
spec_h[idx] = point['spec_h_value']
spec_h_err[idx] = point['spec_h_err']
spec_h_ci_low[idx] = point['spec_h_ci_low']
spec_h_ci_high[idx] = point['spec_h_ci_high']
tau_int[idx] = point['tau_int']
tau_int_ci_low[idx] = point['tau_int_ci_low']
tau_int_ci_high[idx] = point['tau_int_ci_high']
spec_h_samples[idx, 0] = point['spec_h_value']
entropy_bundle = summarize_entropy_observable(
temperatures=temperatures,
specific_heat_samples=spec_h_samples,
specific_heat_err=spec_h_err,
confidence=DEFAULT_CONFIDENCE_LEVEL,
method='blocking',
)
entropy = np.asarray(entropy_bundle['value'], dtype=float)
entropy_ci_low = np.asarray(entropy_bundle['ci_low'], dtype=float)
entropy_ci_high = np.asarray(entropy_bundle['ci_high'], dtype=float)
undefined_autocorr_flag = (~np.isfinite(tau_int)).astype(np.uint8)
rel_width = (tau_int_ci_high - tau_int_ci_low) / np.maximum(np.abs(tau_int), 1e-12)
tau_interval_unstable_flag = np.asarray(
np.isfinite(rel_width) & (rel_width > 1.0),
dtype=np.uint8,
)
low_effective_sample_flag = np.asarray(
np.isfinite(avg_m_n_eff) & (avg_m_n_eff < 20.0),
dtype=np.uint8,
)
confidence_level = DEFAULT_CONFIDENCE_LEVEL
n_seeds = 1
requested_n_seeds = 1
retained_n_seeds = 1
data_source = (
'Computed a lightweight checkerboard-Metropolis fallback sweep '
f'because {cache_reason}.'
)
exclude_from_thermo_plots = undefined_autocorr_flag > 0
avg_m_plot = np.where(exclude_from_thermo_plots, np.nan, avg_m)
avg_m_ci_low_plot = np.where(exclude_from_thermo_plots, np.nan, avg_m_ci_low)
avg_m_ci_high_plot = np.where(exclude_from_thermo_plots, np.nan, avg_m_ci_high)
avg_e_plot = np.where(exclude_from_thermo_plots, np.nan, avg_e)
avg_e_ci_low_plot = np.where(exclude_from_thermo_plots, np.nan, avg_e_ci_low)
avg_e_ci_high_plot = np.where(exclude_from_thermo_plots, np.nan, avg_e_ci_high)
susc_plot = np.where(exclude_from_thermo_plots, np.nan, susc)
susc_ci_low_plot = np.where(exclude_from_thermo_plots, np.nan, susc_ci_low)
susc_ci_high_plot = np.where(exclude_from_thermo_plots, np.nan, susc_ci_high)
spec_h_plot = np.where(exclude_from_thermo_plots, np.nan, spec_h)
spec_h_ci_low_plot = np.where(exclude_from_thermo_plots, np.nan, spec_h_ci_low)
spec_h_ci_high_plot = np.where(exclude_from_thermo_plots, np.nan, spec_h_ci_high)
entropy_plot = np.where(exclude_from_thermo_plots, np.nan, entropy)
entropy_ci_low_plot = np.where(exclude_from_thermo_plots, np.nan, entropy_ci_low)
entropy_ci_high_plot = np.where(exclude_from_thermo_plots, np.nan, entropy_ci_high)
print(data_source)
print(f'Temperature points: {temperatures.size}')
print(
f'retained seeds = {retained_n_seeds} / requested seeds = {requested_n_seeds}, '
f'confidence level = {confidence_level:.2f}'
)
print(f'Excluded points in thermodynamic plots: {int(np.sum(exclude_from_thermo_plots))}')
Convergence Diagnostics
The two-start protocol evolves a random-start and an ordered-start simulation in parallel and declares convergence once their smoothed magnetization traces overlap. Below \(T_{\mathrm{BKT}}\) the quasi-ordered phase presents a slowly relaxing landscape, so the random-start trace may take longer to settle. Near and above \(T_{\mathrm{BKT}}\) the collapse is generally rapid. The panels below show raw and smoothed traces at three representative temperatures.
[ ]:
diagnostic_temps = [0.5 * T_BKT, 1.0 * T_BKT, 1.5 * T_BKT]
diagnostic_titles = [r'Below $T_{\mathrm{BKT}}$', r'Near $T_{\mathrm{BKT}}$', r'Above $T_{\mathrm{BKT}}$']
diagnostic_steps = 1500
fig, axes = plt.subplots(1, 3, figsize=(14, 4.5), sharey=True)
for ax, t, title in zip(axes, diagnostic_temps, diagnostic_titles):
sim_r = XYSimulation(size=32, temp=t, update='checkerboard', init_state='random', seed=42)
sim_o = XYSimulation(size=32, temp=t, update='checkerboard', init_state='ordered', seed=42)
m_r_list, _ = sim_r.run(n_steps=diagnostic_steps)
m_o_list, _ = sim_o.run(n_steps=diagnostic_steps)
m_r = np.asarray(m_r_list, dtype=float)
m_o = np.asarray(m_o_list, dtype=float)
tau_est = estimate_relaxation_time_two_start(m_r, m_o, smooth_window=60, dwell_window=60)
smooth_r = np.convolve(m_r, np.ones(60)/60, mode='valid')
smooth_o = np.convolve(m_o, np.ones(60)/60, mode='valid')
is_stuck = _detect_quasi_steady_stuck(smooth_r, smooth_o, lattice_size=32)
ax.plot(m_r, color=PALETTE['mag'], alpha=0.3)
ax.plot(m_o, color=PALETTE['eng'], alpha=0.3)
ax.plot(np.arange(59, diagnostic_steps), smooth_r, color=PALETTE['mag'], lw=2, label='Random start')
ax.plot(np.arange(59, diagnostic_steps), smooth_o, color=PALETTE['eng'], lw=2, label='Ordered start')
status_text = 'Converged'
if tau_est < diagnostic_steps:
ax.axvline(tau_est, color='0.2', ls='--', lw=1.5, label='Convergence achieved')
else:
if t < T_BKT and is_stuck:
status_text = 'Stuck (Quasi-steady)'
else:
status_text = 'Not converged'
ax.set_title(f'{title} ($T = {t:.3f}$)\nStatus: {status_text}')
ax.set_xlabel('Monte Carlo steps')
if ax is axes[0]:
ax.set_ylabel(r'Order parameter $|M|$')
ax.grid(alpha=0.25)
axes[0].legend(fontsize=9, loc='lower right')
fig.suptitle('Equilibration Traces: Random vs Ordered Starts (XY Model)', y=1.02)
fig.tight_layout()
plt.show()
Magnetization and Energy
The magnetization per spin \(|M|\) in the XY model behaves differently from the Ising case. True long-range order is forbidden at any finite temperature by the Mermin-Wagner theorem, so \(|M|\) never reaches a sharp plateau at unity. Instead, the finite-size magnetization decays gradually from a value set by the system size at low temperature toward zero above \(T_{\mathrm{BKT}}\). The energy per spin decreases smoothly toward the ground-state limit of \(-2J\) as \(T \to 0\) and increases as thermal fluctuations populate higher-energy spin configurations [1] [2].
[ ]:
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].plot(temperatures, avg_m_plot, color=PALETTE['mag'], lw=2)
axes[0].fill_between(temperatures, avg_m_ci_low_plot, avg_m_ci_high_plot, color=PALETTE['mag'], alpha=0.2)
axes[0].axvline(T_BKT, color='0.3', ls='--', lw=1.2)
axes[0].set_xlabel('Temperature $T$')
axes[0].set_ylabel(r'$|M|$')
axes[0].set_title('Order parameter')
axes[0].grid(alpha=0.25)
axes[1].plot(temperatures, avg_e_plot, color=PALETTE['eng'], lw=2)
axes[1].fill_between(temperatures, avg_e_ci_low_plot, avg_e_ci_high_plot, color=PALETTE['eng'], alpha=0.2)
axes[1].axvline(T_BKT, color='0.3', ls='--', lw=1.2)
axes[1].set_xlabel('Temperature $T$')
axes[1].set_ylabel(r'Energy per spin $E$')
axes[1].set_title('Energy density')
axes[1].grid(alpha=0.25)
fig.suptitle('XY model local-update temperature sweep: mean observables', y=1.02)
fig.tight_layout()
plt.show()
Susceptibility and Specific Heat
In the XY model the susceptibility \(\chi\) and specific heat \(C_v\) both show broad humps near \(T_{\mathrm{BKT}}\) rather than the sharp peaks observed in the Ising model. The BKT transition is driven by vortex unbinding, which is a topological mechanism that does not produce the usual power-law divergences in bulk thermodynamic quantities. The specific heat maximum, in particular, tends to sit slightly above \(T_{\mathrm{BKT}}\) and reflects the energy cost of proliferating free vortices rather than a critical-point singularity [2] [3].
[ ]:
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].plot(temperatures, susc_plot, color=PALETTE['chi'], lw=2)
axes[0].fill_between(temperatures, susc_ci_low_plot, susc_ci_high_plot, color=PALETTE['chi'], alpha=0.2)
axes[0].axvline(T_BKT, color='0.3', ls='--', lw=1.2)
axes[0].set_xlabel('Temperature $T$')
axes[0].set_ylabel(r'Susceptibility $\chi$')
axes[0].set_title('Magnetization fluctuations')
axes[0].grid(alpha=0.25)
axes[1].plot(temperatures, spec_h_plot, color=PALETTE['cv'], lw=2)
axes[1].fill_between(temperatures, spec_h_ci_low_plot, spec_h_ci_high_plot, color=PALETTE['cv'], alpha=0.2)
axes[1].axvline(T_BKT, color='0.3', ls='--', lw=1.2)
axes[1].set_xlabel('Temperature $T$')
axes[1].set_ylabel(r'Specific heat $C_v$')
axes[1].set_title('Energy fluctuations')
axes[1].grid(alpha=0.25)
fig.suptitle('XY model critical fluctuations in the local-update sweep', y=1.02)
fig.tight_layout()
plt.show()
Entropy
The relative entropy is reconstructed by integrating \(C_v(T)/T\) downward from the highest sampled temperature, following the same procedure used for the Ising sweep. The low-temperature phase carries the entropy of algebraically correlated spin waves, while above \(T_{\mathrm{BKT}}\) entropy rises more steeply as free vortices contribute additional configurational degrees of freedom.
[ ]:
fig, ax = plt.subplots(figsize=(7.5, 4.8))
ax.plot(temperatures, entropy_plot, color=PALETTE['entropy'], lw=2)
ax.fill_between(temperatures, entropy_ci_low_plot, entropy_ci_high_plot, color=PALETTE['entropy'], alpha=0.2)
ax.axhline(0.0, color='0.4', ls='--', lw=1.0, label=r'Reference at highest sampled $T$')
ax.axvline(T_BKT, color='0.3', ls='--', lw=1.2)
ax.set_xlabel('Temperature $T$')
ax.set_ylabel(r'Relative entropy $S(T) - S(T_{\max})$')
ax.set_title('Entropy from thermodynamic integration (XY model)')
ax.grid(alpha=0.25)
ax.legend(fontsize=9)
fig.tight_layout()
plt.show()
Critical Slowing Down
The integrated autocorrelation time \(\tau_{\mathrm{int}}\) grows near \(T_{\mathrm{BKT}}\) as large-scale spin-wave and vortex fluctuations become harder to decorrelate with local Metropolis updates. The effective sample size \(N_{\mathrm{eff}} \sim N_{\mathrm{meas}} / (2\tau_{\mathrm{int}})\) drops accordingly. Flagged points indicate temperatures where the autocorrelation estimate is undefined (zero-variance window) or unreliable (wide confidence interval or low \(N_{\mathrm{eff}}\)). The BKT slowdown is gentler than the Ising divergence because the transition is infinite-order, but it remains the dominant numerical cost factor for equilibrium sampling [4] [5].
[ ]:
flagged_tau = (undefined_autocorr_flag > 0) | (tau_interval_unstable_flag > 0)
flagged_neff = low_effective_sample_flag > 0
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
axes[0].plot(temperatures, tau_int, color=PALETTE['tau'], lw=2)
axes[0].fill_between(temperatures, tau_int_ci_low, tau_int_ci_high, color=PALETTE['tau'], alpha=0.2)
axes[0].axvline(T_BKT, color='0.3', ls='--', lw=1.2)
if np.any(flagged_tau):
axes[0].scatter(
temperatures[flagged_tau],
tau_int[flagged_tau],
color='black',
marker='x',
s=45,
label='flagged point',
)
axes[0].set_xlabel('Temperature $T$')
axes[0].set_ylabel(r'Integrated autocorrelation time $\tau_{\mathrm{int}}$')
axes[0].set_title('Equilibrium memory')
axes[0].grid(alpha=0.25)
if np.any(flagged_tau):
axes[0].legend(fontsize=9)
axes[1].plot(temperatures, avg_m_n_eff, color=PALETTE['mag'], lw=2)
axes[1].axvline(T_BKT, color='0.3', ls='--', lw=1.2)
if np.any(flagged_neff):
axes[1].scatter(
temperatures[flagged_neff],
avg_m_n_eff[flagged_neff],
color='black',
marker='x',
s=45,
label=r'low $N_{\mathrm{eff}}$',
)
axes[1].set_xlabel('Temperature $T$')
axes[1].set_ylabel(r'Effective sample size $N_{\mathrm{eff}}$')
axes[1].set_title('Independent information per temperature')
axes[1].grid(alpha=0.25)
if np.any(flagged_neff):
axes[1].legend(fontsize=9)
fig.suptitle('Critical slowing down for XY checkerboard Metropolis data', y=1.02)
fig.tight_layout()
plt.show()
Autocorrelation Function at Three Temperatures
The autocorrelation function \(\rho_M(\tau)\) shows how quickly the Markov chain forgets its current magnetization. Below \(T_{\mathrm{BKT}}\) the decay is moderately slow because spin waves carry long-range correlations. Near \(T_{\mathrm{BKT}}\) the decay slows further as vortex-antivortex pairs on all scales couple to the magnetization dynamics. Above the transition, free vortices scatter correlations effectively and \(\rho_M\) drops quickly [4] [5].
[ ]:
autocorr_temperatures = [0.5 * T_BKT, T_BKT, 1.5 * T_BKT]
autocorr_labels = [r'0.5 $T_{\mathrm{BKT}}$', r'$T_{\mathrm{BKT}}$', r'1.5 $T_{\mathrm{BKT}}$']
autocorr_colors = ['#4878CF', '#D65F5F', '#6ACC65']
fig, ax = plt.subplots(figsize=(7.8, 4.8))
for idx, (temp, label, color) in enumerate(zip(autocorr_temperatures, autocorr_labels, autocorr_colors)):
sim_random = XYSimulation(
size=40,
temp=float(temp),
update='checkerboard',
init_state='random',
seed=5_000 + idx,
)
sim_ordered = XYSimulation(
size=40,
temp=float(temp),
update='checkerboard',
init_state='ordered',
seed=5_000 + idx,
)
convergence_equilibrate(
sim_random,
sim_ordered,
chunk_size=120,
max_steps=9_000,
)
mags, _ = sim_random.run(n_steps=1_600)
mags_arr = np.asarray(mags, dtype=float)
try:
corr, tau_here = calculate_autocorr(time_series=mags_arr, max_lag=160)
except ZeroVarianceAutocorrelationError:
corr = np.array([np.nan])
tau_here = float('nan')
lags = np.arange(corr.size)
ax.plot(lags, corr, lw=2, color=color, label=f'{label} ($\\tau_{{int}} \\approx {tau_here:.1f}$)')
ax.axhline(0.0, color='0.4', ls='--', lw=1.0)
ax.set_xlabel(r'Lag $\tau$')
ax.set_ylabel(r'Autocorrelation $\rho_M(\tau)$')
ax.set_title(r'Checkerboard-Metropolis autocorrelation of $|M|$ (XY model)')
ax.set_ylim(-0.05, 1.05)
ax.grid(alpha=0.25)
ax.legend(fontsize=9)
fig.tight_layout()
plt.show()
Synthesis
The full temperature sweep shows the BKT crossover from several complementary perspectives. The finite-size magnetization decays, the energy curve bends, the susceptibility and specific heat develop broad humps, the entropy slope changes, and the autocorrelation time peaks, all in the vicinity of \(T_{\mathrm{BKT}} \approx 0.893\). These signatures are less singular than in the Ising model because the BKT mechanism is topological (vortex unbinding) rather than a conventional symmetry-breaking transition. The companion BKT Transition notebook examines the vortex density, helicity modulus, and correlation function structure that directly probe this topological content.
Bibliography
[1] N. D. Mermin and H. Wagner, “Absence of Ferromagnetism or Antiferromagnetism in One- or Two-Dimensional Isotropic Heisenberg Models,” Physical Review Letters 17, 1133 (1966). APS
[2] J. M. Kosterlitz and D. J. Thouless, “Ordering, metastability and phase transitions in two-dimensional systems,” Journal of Physics C 6, 1181 (1973). IoP Open Access
[3] J. M. Kosterlitz, “The critical properties of the two-dimensional xy model,” Journal of Physics C 7, 1046 (1974). IoP Open Access
[4] A. D. Sokal, “Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms,” lecture notes (1989), reprinted in Functional Integration: Basics and Applications (1997). Springer Link
[5] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics, Oxford University Press (1999). Lecture notes summary