Chapter 4: Critical Phenomena and Scaling Theory

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve # van der Waals equation of state def P_vdw(V, T, a, b, R): """van der Waals pressure""" return R * T / (V - b) - a / V**2 # Parameters for CO₂ R = 8.314e-6 # MPa·m³/(mol·K) a = 0.3658 # MPa·m⁶/mol² b = 4.267e-5 # m³/mol # Critical point T_c = 8 * a / (27 * R * b) P_c = a / (27 * b**2) V_c = 3 * b print("=== van der Waals Critical Point (CO₂) ===") print(f"T_c = {T_c:.2f} K = {T_c - 273.15:.2f} °C") print(f"P_c = {P_c:.4f} MPa = {P_c*10:.2f} bar") print(f"V_c = {V_c*1e6:.2f} cm³/mol") print(f"\nExperimental values (CO₂):") print(f" T_c = 304.13 K = 30.98 °C") print(f" P_c = 7.38 MPa = 73.8 bar") print() # Behavior of physical quantities near the critical point epsilon_range = np.logspace(-3, -0.5, 50) # (T - T_c) / T_c # Computing the order parameter (density difference) def compute_order_parameter(epsilon, a, b, R, V_c, T_c): """Order parameter m = (ρ_L - ρ_G) / ρ_c""" T = T_c * (1 + epsilon) # Determine liquid/gas densities via the Maxwell construction (simplified) # Here we use the mean-field approximation m ~ ε^(1/2) beta_mf = 0.5 # critical exponent of mean-field theory m = epsilon**beta_mf return m # Susceptibility (compressibility) def compute_susceptibility(epsilon, a, b, R, V_c, T_c): """Susceptibility χ ~ ε^(-γ)""" T = T_c * (1 + epsilon) V = V_c # Isothermal compressibility κ_T = -1/V (∂V/∂P)_T dP_dV = -R * T / (V - b)**2 + 2 * a / V**3 kappa_T = -1 / (V * dP_dV) # Normalize by the value near the critical point T_ref = T_c * 1.001 dP_dV_ref = -R * T_ref / (V - b)**2 + 2 * a / V**3 kappa_T_ref = -1 / (V * dP_dV_ref) chi = kappa_T / kappa_T_ref return chi # Computation m_values = [compute_order_parameter(-eps, a, b, R, V_c, T_c) for eps in epsilon_range] chi_values = [compute_susceptibility(eps, a, b, R, V_c, T_c) for eps in epsilon_range] # Visualization fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Order parameter ax1 = axes[0] ax1.loglog(epsilon_range, m_values, 'b-', linewidth=2.5, label='van der Waals') # Theoretical power law beta_mf = 0.5 ax1.loglog(epsilon_range, epsilon_range**beta_mf, 'r--', linewidth=2, label=f'Power law ε^{beta_mf}') ax1.set_xlabel('ε = |T - T_c| / T_c') ax1.set_ylabel('Order parameter m') ax1.set_title('Critical behavior of the order parameter') ax1.legend() ax1.grid(True, alpha=0.3, which='both') # Susceptibility ax2 = axes[1] ax2.loglog(epsilon_range, chi_values, 'g-', linewidth=2.5, label='van der Waals') # Theoretical power law gamma_mf = 1.0 ax2.loglog(epsilon_range, epsilon_range**(-gamma_mf), 'm--', linewidth=2, label=f'Power law ε^{-gamma_mf}') ax2.set_xlabel('ε = |T - T_c| / T_c') ax2.set_ylabel('Susceptibility χ') ax2.set_title('Critical behavior of the susceptibility') ax2.legend() ax2.grid(True, alpha=0.3, which='both') plt.tight_layout() plt.savefig('critical_vdw_exponents.png', dpi=300, bbox_inches='tight') plt.show() # Summary of critical exponents print("=== Critical Exponents of van der Waals (Mean-Field Theory) ===\n") critical_exponents_mf = { 'α (heat capacity)': 0, 'β (order parameter)': 0.5, 'γ (susceptibility)': 1.0, 'δ (critical isotherm)': 3.0, 'ν (correlation length)': 0.5, 'η (correlation function)': 0, } print(f"{'Exponent':<25} {'Mean-field':<15} {'3D Ising (exp.)':<15}") print("-" * 55) for name, value_mf in critical_exponents_mf.items(): # Experimental values for the 3D Ising model ising_3d = { 'α (heat capacity)': 0.110, 'β (order parameter)': 0.326, 'γ (susceptibility)': 1.237, 'δ (critical isotherm)': 4.789, 'ν (correlation length)': 0.630, 'η (correlation function)': 0.036, } value_ising = ising_3d[name] print(f"{name:<25} {value_mf:<15.3f} {value_ising:<15.3f}") print("\nVerification of the scaling law (mean-field theory):") alpha, beta, gamma = 0, 0.5, 1.0 print(f" α + 2β + γ = {alpha} + 2×{beta} + {gamma} = {alpha + 2*beta + gamma}") print(f" (theoretical value: 2)")

import numpy as np import matplotlib.pyplot as plt # van der Waals equation in reduced variables def P_reduced_vdw(V_r, T_r): """Reduced pressure P_r = f(V_r, T_r)""" return 8 * T_r / (3 * V_r - 1) - 3 / V_r**2 # Data for several substances substances = { 'CO₂': {'T_c': 304.13, 'P_c': 7.38, 'V_c': 94.0}, # K, MPa, cm³/mol 'H₂O': {'T_c': 647.1, 'P_c': 22.06, 'V_c': 56.0}, 'N₂': {'T_c': 126.2, 'P_c': 3.39, 'V_c': 90.1}, 'CH₄': {'T_c': 190.6, 'P_c': 4.60, 'V_c': 99.0}, } # Reduced isotherms T_r_values = [0.9, 1.0, 1.1, 1.2, 1.5] colors = ['blue', 'red', 'green', 'orange', 'purple'] fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Plot in reduced variables (universal) ax1 = axes[0] V_r_range = np.linspace(0.4, 5, 200) for T_r, color in zip(T_r_values, colors): P_r_values = [] for V_r in V_r_range: P_r = P_reduced_vdw(V_r, T_r) if P_r > 0: P_r_values.append(P_r) else: P_r_values.append(np.nan) ax1.plot(V_r_range, P_r_values, color=color, linewidth=2, label=f'T_r = {T_r:.1f}') ax1.axhline(1.0, color='gray', linestyle='--', alpha=0.5, label='P_r = 1') ax1.axvline(1.0, color='gray', linestyle='--', alpha=0.5, label='V_r = 1') ax1.plot(1.0, 1.0, 'ko', markersize=10, label='Critical point') ax1.set_xlabel('V_r = V / V_c') ax1.set_ylabel('P_r = P / P_c') ax1.set_title('Law of corresponding states (reduced variables)\nSame curve for all substances') ax1.set_xlim([0, 5]) ax1.set_ylim([0, 3]) ax1.legend(fontsize=9) ax1.grid(True, alpha=0.3) # Plot in real variables (substance-dependent) ax2 = axes[1] # Example at T_r = 1.1 T_r = 1.1 for name, props in list(substances.items())[:3]: # CO₂, H₂O, N₂ T_c = props['T_c'] P_c = props['P_c'] V_c = props['V_c'] T = T_r * T_c V_range = V_r_range * V_c P_values = [] for V_r in V_r_range: P_r = P_reduced_vdw(V_r, T_r) if P_r > 0: P_values.append(P_r * P_c) else: P_values.append(np.nan) ax2.plot(V_range, P_values, linewidth=2, label=f'{name} ({T:.0f} K)') ax2.set_xlabel('Molar volume V (cm³/mol)') ax2.set_ylabel('Pressure P (MPa)') ax2.set_title(f'Isotherms in real variables (T_r = {T_r})\nDifferent for each substance') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('critical_corresponding_states.png', dpi=300, bbox_inches='tight') plt.show() # Compressibility factor Z = PV/(RT) print("=== Verification of the Law of Corresponding States ===\n") print("Compressibility factor at the critical point: Z_c = P_c V_c / (R T_c)") print() print(f"{'Substance':<10} {'Z_c (exp.)':<15} {'Z_c (vdW)':<15}") print("-" * 40) R = 8.314 # J/(mol·K) = 8.314e-6 MPa·m³/(mol·K) Z_c_vdw = 3 / 8 # universal van der Waals value for name, props in substances.items(): T_c = props['T_c'] P_c = props['P_c'] * 1e6 # MPa → Pa V_c = props['V_c'] * 1e-6 # cm³/mol → m³/mol Z_c_exp = P_c * V_c / (R * T_c) print(f"{name:<10} {Z_c_exp:<15.4f} {Z_c_vdw:<15.4f}") print(f"\nvan der Waals theory: Z_c = 3/8 = {Z_c_vdw:.4f} (same for all substances)") print("Experimental values: Z_c ≈ 0.27-0.29 (varies slightly by substance)") print("\nThe law of corresponding states holds approximately (affected by quantum effects and molecular structure)")

import numpy as np import matplotlib.pyplot as plt def ising_1d_exact(T, J, h, k_B=1.0): """Exact solution of the 1D Ising model Args: T: temperature J: exchange interaction h: external magnetic field k_B: Boltzmann constant Returns: m: magnetization chi: susceptibility U: internal energy C: heat capacity """ beta = 1 / (k_B * T) # Magnetization (case h ≠ 0) if abs(h) > 1e-10: # Exact solution for the magnetization exp_beta_J = np.exp(beta * J) exp_beta_h = np.exp(beta * h) exp_minus_beta_h = np.exp(-beta * h) # Contributions to the partition function Z_plus = exp_beta_J * exp_beta_h + np.exp(-beta * J) * exp_minus_beta_h Z_minus = exp_beta_J * exp_minus_beta_h + np.exp(-beta * J) * exp_beta_h m = (Z_plus - Z_minus) / (Z_plus + Z_minus) else: m = 0 # for h = 0, m = 0 at T > 0 # Susceptibility (derivative at h = 0) exp_2beta_J = np.exp(2 * beta * J) chi = beta * exp_2beta_J / (1 + exp_2beta_J) # Internal energy U = -J * np.tanh(beta * J) # Heat capacity C = k_B * (beta * J)**2 / np.cosh(beta * J)**2 return m, chi, U, C # Parameters J = 1.0 k_B = 1.0 h = 0.1 # small external field # Temperature range T_range = np.linspace(0.1, 5.0, 100) m_values = [] chi_values = [] U_values = [] C_values = [] for T in T_range: m, chi, U, C = ising_1d_exact(T, J, h, k_B) m_values.append(m) chi_values.append(chi) U_values.append(U) C_values.append(C) # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Magnetization ax1 = axes[0, 0] ax1.plot(T_range, m_values, 'b-', linewidth=2.5) ax1.set_xlabel('Temperature T/J') ax1.set_ylabel('Magnetization m') ax1.set_title(f'Magnetization of the 1D Ising model (h = {h}J)') ax1.grid(True, alpha=0.3) ax1.axhline(0, color='gray', linestyle='--', alpha=0.5) # Susceptibility ax2 = axes[0, 1] ax2.plot(T_range, chi_values, 'r-', linewidth=2.5) ax2.set_xlabel('Temperature T/J') ax2.set_ylabel('Susceptibility χ') ax2.set_title('Susceptibility of the 1D Ising model (h = 0)') ax2.grid(True, alpha=0.3) # Internal energy ax3 = axes[1, 0] ax3.plot(T_range, U_values, 'g-', linewidth=2.5) ax3.set_xlabel('Temperature T/J') ax3.set_ylabel('Internal energy U/J') ax3.set_title('Internal energy of the 1D Ising model') ax3.grid(True, alpha=0.3) # Heat capacity ax4 = axes[1, 1] ax4.plot(T_range, C_values, 'm-', linewidth=2.5) ax4.set_xlabel('Temperature T/J') ax4.set_ylabel('Heat capacity C/k_B') ax4.set_title('Heat capacity of the 1D Ising model') ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('critical_ising_1d_exact.png', dpi=300, bbox_inches='tight') plt.show() # Display results print("=== 1D Ising Model (Exact Solution) ===\n") print("Key properties:") print(" - No phase transition (T_c = 0)") print(" - For h = 0, m = 0 at T > 0") print(" - Susceptibility diverges as T → 0") print(" - Heat capacity remains finite (no peak)") print() # Low- and high-temperature limits T_low = 0.1 m_low, chi_low, U_low, C_low = ising_1d_exact(T_low, J, 0, k_B) print(f"Low-temperature limit (T = {T_low}J):") print(f" Susceptibility χ = {chi_low:.4f}") print(f" Internal energy U/J = {U_low:.4f} → -1 (fully aligned)") print() T_high = 5.0 m_high, chi_high, U_high, C_high = ising_1d_exact(T_high, J, 0, k_B) print(f"High-temperature limit (T = {T_high}J):") print(f" Susceptibility χ = {chi_high:.4f} → 0") print(f" Internal energy U/J = {U_high:.4f} → 0 (random)")

import numpy as np import matplotlib.pyplot as plt def ising_2d_monte_carlo(L, T, J, h, n_steps, n_equilibrate): """Monte Carlo simulation of the 2D Ising model (Metropolis algorithm) Args: L: lattice size (L×L) T: temperature J: exchange interaction h: external magnetic field n_steps: number of Monte Carlo steps n_equilibrate: number of equilibration steps Returns: m_avg: average magnetization E_avg: average energy spin_config: final spin configuration """ # Initial configuration (random) spins = 2 * np.random.randint(0, 2, size=(L, L)) - 1 beta = 1 / T # k_B = 1 # Energy calculation def compute_energy(spins): E = 0 for i in range(L): for j in range(L): s = spins[i, j] # Sum over nearest neighbors (periodic boundary conditions) neighbors = spins[(i+1) % L, j] + spins[i, (j+1) % L] + \ spins[(i-1) % L, j] + spins[i, (j-1) % L] E += -J * s * neighbors - h * s return E / 2 # correct for double counting # Metropolis algorithm m_history = [] E_history = [] for step in range(n_equilibrate + n_steps): # Pick a spin at random i, j = np.random.randint(0, L, size=2) s = spins[i, j] # Nearest-neighbor spins neighbors = spins[(i+1) % L, j] + spins[i, (j+1) % L] + \ spins[(i-1) % L, j] + spins[i, (j-1) % L] # Energy change dE = 2 * s * (J * neighbors + h) # Metropolis acceptance test if dE < 0 or np.random.rand() < np.exp(-beta * dE): spins[i, j] *= -1 # flip the spin # Measurement (after equilibration) if step >= n_equilibrate: m = np.mean(spins) E = compute_energy(spins) m_history.append(m) E_history.append(E) m_avg = np.mean(m_history) E_avg = np.mean(E_history) return m_avg, E_avg, spins # Parameters L = 20 # lattice size J = 1.0 h = 0.0 n_steps = 5000 n_equilibrate = 1000 # Onsager critical temperature (2D Ising, square lattice) T_c_exact = 2 * J / np.log(1 + np.sqrt(2)) # ≈ 2.269 print(f"=== 2D Ising Model (Metropolis Monte Carlo) ===") print(f"Lattice size: {L}×{L}") print(f"Onsager critical temperature: T_c = {T_c_exact:.4f} J") print() # Temperature range T_range = np.linspace(1.0, 4.0, 20) m_values = [] E_values = [] for T in T_range: m, E, _ = ising_2d_monte_carlo(L, T, J, h, n_steps, n_equilibrate) m_values.append(abs(m)) # absolute value of the spontaneous magnetization E_values.append(E) print(f"T = {T:.2f}: m = {m:.4f}, E = {E:.2f}") # Visualization fig, axes = plt.subplots(1, 3, figsize=(18, 5)) # Magnetization ax1 = axes[0] ax1.plot(T_range, m_values, 'bo-', markersize=5, linewidth=2) ax1.axvline(T_c_exact, color='r', linestyle='--', linewidth=2, label=f'T_c = {T_c_exact:.3f}') ax1.set_xlabel('Temperature T/J') ax1.set_ylabel('Magnetization |m|') ax1.set_title('Spontaneous magnetization of the 2D Ising model') ax1.legend() ax1.grid(True, alpha=0.3) # Energy ax2 = axes[1] ax2.plot(T_range, E_values, 'ro-', markersize=5, linewidth=2) ax2.axvline(T_c_exact, color='r', linestyle='--', linewidth=2) ax2.set_xlabel('Temperature T/J') ax2.set_ylabel('Energy per spin E/(NJ)') ax2.set_title('Energy of the 2D Ising model') ax2.grid(True, alpha=0.3) # Visualization of spin configurations ax3 = axes[2] # Configurations at low (T = 1.5) and high (T = 3.5) temperature T_low = 1.5 T_high = 3.5 _, _, spins_low = ising_2d_monte_carlo(L, T_low, J, h, n_steps, n_equilibrate) _, _, spins_high = ising_2d_monte_carlo(L, T_high, J, h, n_steps, n_equilibrate) # Display side by side combined = np.hstack([spins_low, np.ones((L, 2)), spins_high]) im = ax3.imshow(combined, cmap='coolwarm', interpolation='nearest') ax3.set_title(f'Spin configurations\nLeft: T={T_low} < T_c, Right: T={T_high} > T_c') ax3.axis('off') plt.colorbar(im, ax=ax3, fraction=0.046, pad=0.04) plt.tight_layout() plt.savefig('critical_ising_2d_monte_carlo.png', dpi=300, bbox_inches='tight') plt.show() print(f"\nObservations:") print(f" T < T_c: ordered phase (almost all ↑ or ↓)") print(f" T > T_c: disordered phase (random spin configuration)") print(f" T ≈ T_c: large fluctuations (analogous to critical opalescence)")

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit # 2D Ising Monte Carlo data (from the previous code) # Sample densely near the critical point L = 30 J = 1.0 h = 0.0 T_c_exact = 2 * J / np.log(1 + np.sqrt(2)) # Measure the magnetization for T < T_c T_below = np.linspace(T_c_exact * 0.7, T_c_exact * 0.995, 15) m_below = [] print("=== Fitting the Critical Exponents ===\n") print("Collecting data...") for T in T_below: m, _, _ = ising_2d_monte_carlo(L, T, J, h, n_steps=8000, n_equilibrate=2000) m_below.append(abs(m)) # Measure the susceptibility for T > T_c (simplified: magnetization fluctuations) T_above = np.linspace(T_c_exact * 1.005, T_c_exact * 1.5, 15) chi_above = [] for T in T_above: # Repeat measurements to compute fluctuations m_samples = [] for _ in range(5): m, _, _ = ising_2d_monte_carlo(L, T, J, h, n_steps=3000, n_equilibrate=1000) m_samples.append(m) # Susceptibility χ = β ( - ²) beta = 1 / T m_mean = np.mean(m_samples) m2_mean = np.mean(np.array(m_samples)**2) chi = beta * L * L * (m2_mean - m_mean**2) chi_above.append(chi) # Power-law fitting # m ~ ε^β, χ ~ ε^(-γ) epsilon_below = (T_c_exact - T_below) / T_c_exact epsilon_above = (T_above - T_c_exact) / T_c_exact # β exponent (order parameter) def power_law_beta(eps, beta, A): return A * eps**beta params_beta, _ = curve_fit(power_law_beta, epsilon_below, m_below, p0=[0.3, 1.0], bounds=([0, 0], [1, 10])) beta_fit, A_beta = params_beta # γ exponent (susceptibility) def power_law_gamma(eps, gamma, B): return B * eps**(-gamma) params_gamma, _ = curve_fit(power_law_gamma, epsilon_above, chi_above, p0=[1.2, 10.0], bounds=([0, 0], [3, 1000])) gamma_fit, B_gamma = params_gamma print(f"\nFitting results:") print(f" β (order parameter exponent) = {beta_fit:.3f}") print(f" γ (susceptibility exponent) = {gamma_fit:.3f}") print() print(f"2D Ising theoretical values:") print(f" β = 1/8 = 0.125") print(f" γ = 7/4 = 1.750") print() # Visualization fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # β exponent ax1 = axes[0] ax1.loglog(epsilon_below, m_below, 'bo', markersize=8, label='MC data') ax1.loglog(epsilon_below, power_law_beta(epsilon_below, beta_fit, A_beta), 'r-', linewidth=2, label=f'Fit: β = {beta_fit:.3f}') ax1.loglog(epsilon_below, power_law_beta(epsilon_below, 0.125, A_beta), 'g--', linewidth=2, label='Theory: β = 0.125') ax1.set_xlabel('ε = (T_c - T) / T_c') ax1.set_ylabel('Magnetization m') ax1.set_title('Order parameter exponent β') ax1.legend() ax1.grid(True, alpha=0.3, which='both') # γ exponent ax2 = axes[1] ax2.loglog(epsilon_above, chi_above, 'ro', markersize=8, label='MC data') ax2.loglog(epsilon_above, power_law_gamma(epsilon_above, gamma_fit, B_gamma), 'b-', linewidth=2, label=f'Fit: γ = {gamma_fit:.3f}') ax2.loglog(epsilon_above, power_law_gamma(epsilon_above, 1.75, B_gamma), 'g--', linewidth=2, label='Theory: γ = 1.75') ax2.set_xlabel('ε = (T - T_c) / T_c') ax2.set_ylabel('Susceptibility χ') ax2.set_title('Susceptibility exponent γ') ax2.legend() ax2.grid(True, alpha=0.3, which='both') plt.tight_layout() plt.savefig('critical_exponents_fitting.png', dpi=300, bbox_inches='tight') plt.show() print("Notes:") print(" - Finite-size effects cause deviations from the theoretical values") print(" - Larger lattices (L > 100) allow higher-precision measurements") print(" - Power-law behavior is most pronounced very close to the critical point (ε < 0.01)")

import numpy as np import matplotlib.pyplot as plt # Critical exponent data critical_exponents = { 'Mean-field theory': {'alpha': 0, 'beta': 0.5, 'gamma': 1.0, 'delta': 3.0, 'nu': 0.5, 'eta': 0}, '2D Ising (theory)': {'alpha': 0, 'beta': 0.125, 'gamma': 1.75, 'delta': 15, 'nu': 1.0, 'eta': 0.25}, '3D Ising (experiment)': {'alpha': 0.110, 'beta': 0.326, 'gamma': 1.237, 'delta': 4.789, 'nu': 0.630, 'eta': 0.036}, '3D Heisenberg': {'alpha': -0.133, 'beta': 0.365, 'gamma': 1.386, 'delta': 4.80, 'nu': 0.705, 'eta': 0.035}, } # Scaling relations scaling_relations = { 'Rushbrooke': lambda exp: exp['alpha'] + 2*exp['beta'] + exp['gamma'], 'Widom': lambda exp: exp['gamma'] / (exp['beta'] * (exp['delta'] - 1)), 'Fisher': lambda exp: exp['gamma'] / (exp['nu'] * (2 - exp['eta'])), } print("=== Verification of the Scaling Relations ===\n") # Verification for each theory/experiment results = {} for name, exps in critical_exponents.items(): results[name] = {} # Rushbrooke: α + 2β + γ = 2 rushbrooke = exps['alpha'] + 2*exps['beta'] + exps['gamma'] results[name]['Rushbrooke'] = rushbrooke # Widom: γ / (β(δ-1)) = 1 widom = exps['gamma'] / (exps['beta'] * (exps['delta'] - 1)) results[name]['Widom'] = widom # Fisher: γ / (ν(2-η)) = 1 fisher = exps['gamma'] / (exps['nu'] * (2 - exps['eta'])) results[name]['Fisher'] = fisher # Display results print(f"{'Theory/Experiment':<20} {'Rushbrooke':<15} {'Widom':<15} {'Fisher':<15}") print(f"{'(Expected)':<20} {'(2.000)':<15} {'(1.000)':<15} {'(1.000)':<15}") print("-" * 65) for name, vals in results.items(): rush = vals['Rushbrooke'] wid = vals['Widom'] fish = vals['Fisher'] print(f"{name:<20} {rush:<15.4f} {wid:<15.4f} {fish:<15.4f}") # Visualize the deviations fig, ax = plt.subplots(figsize=(10, 6)) relations = ['Rushbrooke', 'Widom', 'Fisher'] expected = [2.0, 1.0, 1.0] x = np.arange(len(relations)) width = 0.2 for i, (name, vals) in enumerate(results.items()): values = [vals[r] for r in relations] ax.bar(x + i*width, values, width, label=name, alpha=0.8) # Lines for expected values for i, (rel, exp_val) in enumerate(zip(relations, expected)): ax.axhline(exp_val, color='red', linestyle='--', linewidth=1.5, alpha=0.5) ax.text(i, exp_val + 0.1, f'Expected: {exp_val}', ha='center', fontsize=9) ax.set_xlabel('Scaling relation') ax.set_ylabel('Computed value') ax.set_title('Verification of scaling relations') ax.set_xticks(x + width * 1.5) ax.set_xticklabels(relations) ax.legend() ax.grid(True, axis='y', alpha=0.3) plt.tight_layout() plt.savefig('critical_scaling_relations.png', dpi=300, bbox_inches='tight') plt.show() # Universality classes print("\n=== Universality Classes ===\n") print("Critical exponents do not depend on the details of the system (lattice") print("structure, details of the interactions); they are determined only by:") print(" 1. Spatial dimension d") print(" 2. Number of order parameter components n") print(" 3. Range of interactions (short- or long-range)") print() print("Examples:") print(" - Ising (n=1): lattice gas, binary alloys, ferromagnets") print(" - Heisenberg (n=3): isotropic ferromagnets") print(" - XY (n=2): superfluidity, superconductivity") print() print("This universality underpins the theoretical foundation of critical phenomena (the renormalization group)")

import numpy as np import matplotlib.pyplot as plt def renormalize_1d_ising(K): """Real-space renormalization group transformation of the 1D Ising model Args: K: coupling constant K = J/(k_B T) Returns: K_prime: renormalized coupling constant """ # Renormalization over two-spin blocks # Partition function: Z = Σ exp(K s_i s_{i+1}) # Block spin S_I = sign(s_{2I} + s_{2I+1}) # Renormalization of the coupling (1D Ising) # K' = (1/2) ln(cosh(2K)) K_prime = 0.5 * np.log(np.cosh(2 * K)) return K_prime # Renormalization group flow K_initial_values = np.linspace(0.01, 2.0, 20) fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Flow diagram ax1 = axes[0] for K_init in K_initial_values: K_flow = [K_init] K = K_init for _ in range(20): K = renormalize_1d_ising(K) K_flow.append(K) ax1.plot(K_flow, 'o-', markersize=3, linewidth=1, alpha=0.7) ax1.set_xlabel('Renormalization step n') ax1.set_ylabel('Coupling constant K') ax1.set_title('RG flow of the 1D Ising model') ax1.grid(True, alpha=0.3) ax1.text(15, 1.5, 'All flows converge to\nK* = 0 (high-temperature\nfixed point)', fontsize=11, bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.6)) # Plot of K vs K' ax2 = axes[1] K_range = np.linspace(0, 2.0, 100) K_prime_range = [renormalize_1d_ising(K) for K in K_range] ax2.plot(K_range, K_prime_range, 'b-', linewidth=2.5, label="K' = f(K)") ax2.plot(K_range, K_range, 'r--', linewidth=2, label="K' = K (fixed point)") # Fixed point K_fixed = 0 # K* = 0 ax2.plot(K_fixed, K_fixed, 'ro', markersize=10, label=f'Fixed point K* = {K_fixed}') ax2.set_xlabel('K = J/(k_B T)') ax2.set_ylabel("K' (after renormalization)") ax2.set_title('Renormalization group transformation K → K\'') ax2.legend() ax2.grid(True, alpha=0.3) ax2.set_xlim([0, 2.0]) ax2.set_ylim([0, 2.0]) plt.tight_layout() plt.savefig('critical_renormalization_group_1d.png', dpi=300, bbox_inches='tight') plt.show() # Analysis print("=== Renormalization Group (Real-Space RG) ===\n") print("RG transformation of the 1D Ising model:") print(" K' = (1/2) ln(cosh(2K))") print() print("Fixed points:") print(" K* = 0 (high-temperature fixed point, T = ∞)") print(" K* = ∞ (low-temperature fixed point, T = 0)") print() print("Conclusions:") print(" - For K < K* (finite temperature), K → 0 (disordered phase)") print(" - No phase transition at finite temperature (T_c = 0)") print() print("For the 2D Ising model:") print(" - A nontrivial fixed point K_c exists") print(" - K_c = ln(1 + √2)/2 ≈ 0.4407") print(" - T_c = 2J/(k_B ln(1 + √2)) ≈ 2.269 J/k_B") print() # Linearization and critical exponents print("=== Computing Critical Exponents (Linearization) ===\n") print("Linearize near the fixed point K*:") print(" δK_{n+1} = λ δK_n") print(" λ = dK'/dK|_{K=K*} (eigenvalue)") print() # Linearization at K* = 0 K_star = 0 dK_prime_dK = 2 * np.tanh(2 * K_star) # derivative at K* = 0 print(f"1D Ising (K* = 0):") print(f" λ = dK'/dK|_{{K*=0}} = {dK_prime_dK:.4f}") print(f" |λ| < 1 → K* = 0 is stable (attractive fixed point)") print() print("Critical exponent ν of the correlation length:") print(" ξ ~ |T - T_c|^(-ν)") print(" ν = ln(b) / ln(λ) (b: length rescaling factor)") print(" 1D Ising: b = 2, λ = 0 → ν = ∞ (no phase transition)")