Chapter 5: Renormalization Theory and Effective Field Theory

import numpy as np
from scipy.special import gamma

# ===================================
# Integral formula in dimensional regularization
# ===================================

def dimensional_integral(n, Delta, d=4):
    """Dimensional regularization integral

    I_n(Δ) = ∫ d^d k / (2π)^d  1/(k² + Δ)^n

    Args:
        n: power in denominator
        Delta: mass parameter
        d: spacetime dimension
    """
    epsilon = (4 - d) / 2

    # Formula using Γ function
    prefactor = 1 / (4 * np.pi)**(d / 2)
    gamma_factor = gamma(n - d / 2) / gamma(n)
    delta_factor = Delta**(d / 2 - n)

    I_n = prefactor * gamma_factor * delta_factor

    return I_n

def extract_pole(epsilon, m2, mu2=1.0):
    """Separate pole in ε and finite part

    I ~ 1/ε + log(m²/μ²) + O(ε)
    """
    if epsilon < 1e-6:
        pole = 1 / epsilon
        gamma_E = 0.5772156649
        finite = -gamma_E + np.log(4 * np.pi) - np.log(m2 / mu2)
    else:
        # Evaluation at finite ε
        pole = 1 / epsilon
        finite = -np.log(m2 / mu2)

    return pole, finite

# Example of 1-loop integral
m2 = 1.0  # mass squared
mu2 = 1.0  # renormalization scale
d = 3.99  # d = 4 - 2ε, ε = 0.005

I1 = dimensional_integral(1, m2, d)
epsilon = (4 - d) / 2
pole, finite = extract_pole(epsilon, m2, mu2)

print("Integral with dimensional regularization:")
print("=" * 50)
print(f"Spacetime dimension d = {d} (ε = {epsilon})")
print(f"Mass m² = {m2}")
print(f"Renormalization scale μ² = {mu2}")
print(f"\nIntegral value I₁: {I1:.6e}")
print(f"Pole: 1/ε = {pole:.2f}")
print(f"Finite part: {finite:.6f}")

import numpy as np
from scipy.integrate import odeint

# ===================================
# Renormalization group flow of φ⁴ theory
# ===================================

def beta_phi4(lambda_, d=4):
    """β function of φ⁴ theory (1-loop)

    β(λ) = (4-d)λ + 3λ²/(16π²) + O(λ³)
    """
    epsilon = 4 - d
    beta = epsilon * lambda_ + 3 * lambda_**2 / (16 * np.pi**2)

    return beta

def gamma_phi4(lambda_):
    """Field anomalous dimension (1-loop)"""
    gamma = lambda_ / (16 * np.pi**2)

    return gamma

def rg_flow(lambda_, t, d=4):
    """Differential equation for RG flow

    dλ/dt = β(λ), t = log(μ/μ₀)
    """
    return beta_phi4(lambda_, d)

# Numerical solution of RG flow
lambda_0 = 0.1  # initial coupling constant
t_array = np.linspace(0, 10, 100)  # log(μ/μ₀)

# d=4 (critical dimension)
lambda_d4 = odeint(rg_flow, lambda_0, t_array, args=(4,))

# d=3 (renormalizable)
lambda_d3 = odeint(rg_flow, lambda_0, t_array, args=(3,))

print("RG flow of φ⁴ theory:")
print("=" * 60)
print(f"{'log(μ/μ₀)':<15} {'λ(d=4)':<20} {'λ(d=3)':<20}")
print("-" * 60)

for i in [0, 25, 50, 75, 99]:
    print(f"{t_array[i]:<15.2f} {lambda_d4[i][0]:<20.6f} {lambda_d3[i][0]:<20.6f}")

import numpy as np

# ===================================
# Critical behavior of Ising model
# ===================================

def ising_critical_exponents(d):
    """Critical exponents of Ising model (approximate values)

    Args:
        d: spatial dimension
    """
    exponents = {
        2: {'nu': 1.0, 'beta': 0.125, 'gamma': 1.75, 'alpha': 0.0},
        3: {'nu': 0.63, 'beta': 0.325, 'gamma': 1.24, 'alpha': 0.11},
        4: {'nu': 0.5, 'beta': 0.5, 'gamma': 1.0, 'alpha': 0.0},  # mean field
    }
    return exponents.get(d, exponents[3])

def verify_scaling_relations(exponents, d):
    """Verification of scaling relations"""
    nu, beta, gamma, alpha = (exponents['nu'], exponents['beta'],
                              exponents['gamma'], exponents['alpha'])

    # Rushbrooke inequality: α + 2β + γ = 2
    rushbrooke = alpha + 2 * beta + gamma

    # Hyperscaling: νd = 2 - α
    hyperscaling_lhs = nu * d
    hyperscaling_rhs = 2 - alpha

    return rushbrooke, hyperscaling_lhs, hyperscaling_rhs

# Verification for each dimension
dimensions = [2, 3, 4]

print("Critical exponents of Ising model:")
print("=" * 70)

for d in dimensions:
    exp = ising_critical_exponents(d)
    rush, hyp_l, hyp_r = verify_scaling_relations(exp, d)

    print(f"\nd = {d}:")
    print(f"  ν = {exp['nu']:.3f}, β = {exp['beta']:.3f}, "
          f"γ = {exp['gamma']:.3f}, α = {exp['alpha']:.3f}")
    print(f"  Rushbrooke: α + 2β + γ = {rush:.3f} (theoretical value: 2)")
    print(f"  Hyperscaling: νd = {hyp_l:.3f}, 2-α = {hyp_r:.3f}")

import numpy as np

# ===================================
# Fermi theory and electroweak unified theory
# ===================================

def fermi_coupling_from_mw(M_W, g_w):
    """Derive Fermi coupling constant from W boson mass

    G_F = g²/(8M_W²)
    """
    G_F = g_w**2 / (8 * M_W**2)
    return G_F

def effective_vs_full_theory(E, M_W, g_w):
    """Comparison of effective theory and full theory

    Low energy (E << M_W): Fermi theory
    High energy (E ~ M_W): Full electroweak theory
    """
    G_F = fermi_coupling_from_mw(M_W, g_w)

    # Cross section in Fermi theory (E << M_W)
    sigma_fermi = G_F**2 * E**2

    # Cross section in full theory (suppression by propagator)
    sigma_full = (g_w**4 / (E**2 + M_W**2)**2) * E**2

    validity = E / M_W  # validity parameter of effective theory

    return sigma_fermi, sigma_full, validity

# Parameters
M_W = 80.4  # GeV (W boson mass)
g_w = 0.65  # weak coupling constant
G_F = 1.166e-5  # GeV^-2 (experimental value)

energies = [1, 10, 50, 100]  # GeV

print("Comparison of effective theory and full theory:")
print("=" * 70)
print(f"W boson mass: {M_W} GeV")
print(f"Fermi constant G_F: {G_F:.3e} GeV^-2")
print(f"\n{'E (GeV)':<12} {'σ_Fermi':<18} {'σ_full':<18} {'E/M_W':<12}")
print("-" * 70)

for E in energies:
    sigma_f, sigma_full, val = effective_vs_full_theory(E, M_W, g_w)

    print(f"{E:<12} {sigma_f:<18.3e} {sigma_full:<18.3e} {val:<12.4f}")

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Landau-Ginzburg free energy
# ===================================

def landau_free_energy(M, r, u):
    """Landau free energy (uniform field)

    F(M) = r/2 M² + u/4 M⁴
    """
    return r / 2 * M**2 + u / 4 * M**4

def equilibrium_magnetization(r, u):
    """Calculate equilibrium magnetization"""
    if r >= 0:
        # Paramagnetic phase
        return 0.0
    else:
        # Ferromagnetic phase
        return np.sqrt(-r / u)

def susceptibility(r, u, M_eq):
    """Susceptibility χ = ∂M/∂H"""
    if r >= 0:
        # χ ~ 1/r (Curie-Weiss)
        return 1 / r
    else:
        # χ ~ 1/(-2r)
        return 1 / (-2 * r)

# Parameters
u = 1.0
r_values = np.linspace(-2.0, 2.0, 100)
M_eq = [equilibrium_magnetization(r, u) for r in r_values]

# Near critical temperature
T_range = r_values  # r ∝ (T - T_c)

print("Phase transition by Landau theory:")
print("=" * 60)
print(f"{'r (T-Tc)':<15} {'M_eq':<15} {'χ':<15}")
print("-" * 60)

for r in [-1.0, -0.5, 0.5, 1.0]:
    M = equilibrium_magnetization(r, u)
    chi = susceptibility(r, u, M) if r != 0 else np.inf

    print(f"{r:<15.2f} {M:<15.6f} {chi:<15.6f}")

import numpy as np

# ===================================
# Ferroelectric phase transition of BaTiO₃ (Landau theory)
# ===================================

def landau_free_energy_ferro(P, a, b, c, E=0):
    """Landau free energy for ferroelectrics

    F(P) = a/2 P² + b/4 P⁴ + c/6 P⁶ - EP

    Args:
        P: polarization
        a: quadratic coefficient (temperature dependent)
        b, c: higher order coefficients
        E: external electric field
    """
    return a / 2 * P**2 + b / 4 * P**4 + c / 6 * P**6 - E * P

def dielectric_constant(a, b, P_eq):
    """Dielectric constant ε ~ χ"""
    if a > 0:
        # Paraelectric phase (Curie-Weiss law)
        epsilon = 1 / a
    else:
        # Ferroelectric phase
        epsilon = 1 / (a + 3 * b * P_eq**2)

    return epsilon

# Parameters for BaTiO₃ (simplified)
T_c = 393  # K (Curie temperature)
alpha_0 = 0.01  # temperature coefficient
b = 1.0
c = 0.1

temperatures = [300, 350, 400, 450]  # K

print("Ferroelectric phase transition of BaTiO₃:")
print("=" * 60)
print(f"Curie temperature: {T_c} K")
print(f"\n{'T (K)':<12} {'a(T)':<15} {'P_eq':<15} {'ε':<15}")
print("-" * 60)

for T in temperatures:
    a_T = alpha_0 * (T - T_c)  # a ∝ (T - T_c)

    # Equilibrium polarization
    if a_T < 0 and b > 0:
        P_eq = np.sqrt(-a_T / b)
    else:
        P_eq = 0.0

    epsilon = dielectric_constant(a_T, b, P_eq) if a_T != 0 else np.inf

    print(f"{T:<12} {a_T:<15.4f} {P_eq:<15.6f} {epsilon:<15.6f}")

import numpy as np

# ===================================
# Cahn-Hilliard equation (spinodal decomposition)
# ===================================

def cahn_hilliard_growth_rate(k, r, kappa):
    """Linear growth rate of CH equation

    ∂c/∂t = M ∇²(δF/δc)
    ω(k) = -M k² (r + κ k²)

    Args:
        k: wave number
        r: free energy coefficient (spinodal for r < 0)
        kappa: gradient energy coefficient
    """
    M = 1.0  # mobility
    omega = -M * k**2 * (r + kappa * k**2)

    return omega

def fastest_growing_mode(r, kappa):
    """Fastest growing mode

    k_m = sqrt(-r / (2κ))
    """
    if r >= 0:
        return 0.0, 0.0

    k_m = np.sqrt(-r / (2 * kappa))
    omega_m = cahn_hilliard_growth_rate(k_m, r, kappa)

    return k_m, omega_m

# Parameters (spinodal decomposition in alloy)
r = -1.0  # spinodal region
kappa = 1.0

k_array = np.linspace(0.01, 2.0, 100)
omega_array = [cahn_hilliard_growth_rate(k, r, kappa) for k in k_array]

k_m, omega_m = fastest_growing_mode(r, kappa)

print("Dynamics of spinodal decomposition:")
print("=" * 50)
print(f"Free energy coefficient r: {r}")
print(f"Gradient coefficient κ: {kappa}")
print(f"\nFastest growing wave number k_m: {k_m:.6f}")
print(f"Growth rate ω(k_m): {omega_m:.6f}")
print(f"Characteristic length scale λ_m: {2*np.pi/k_m:.6f}")

import numpy as np

# ===================================
# Relaxation time near critical point
# ===================================

def relaxation_time(T, T_c, tau_0=1.0, z=2, nu=0.63):
    """Critical slowing down

    τ ~ ξ^z ~ |T - T_c|^{-zν}

    Args:
        z: dynamic critical exponent
        nu: correlation length exponent
    """
    t_reduced = np.abs(T - T_c) / T_c

    if t_reduced < 1e-10:
        return 1e10  # divergence

    tau = tau_0 * t_reduced**(-z * nu)

    return tau

# Ferromagnetic transition of iron
T_c = 1043  # K
tau_0 = 1e-12  # s
z = 2  # dynamic exponent (Model B)
nu = 0.63  # Ising universality class

temperatures = [T_c + dT for dT in [1, 10, 50, 100]]

print("Critical slowing down:")
print("=" * 60)
print(f"Curie temperature T_c: {T_c} K")
print(f"Dynamic exponent z: {z}")
print(f"Correlation length exponent ν: {nu}")
print(f"\n{'T (K)':<15} {'ΔT (K)':<15} {'τ (s)':<20}")
print("-" * 60)

for T in temperatures:
    dT = T - T_c
    tau = relaxation_time(T, T_c, tau_0, z, nu)

    print(f"{T:<15.1f} {dT:<15.1f} {tau:<20.6e}")