Chapter 4: Feynman Diagram Techniques

import numpy as np

# ===================================
# Scattering amplitude by Feynman rules
# ===================================

def propagator(p, m, epsilon=1e-3):
    """Propagator for scalar field"""
    p2 = np.dot(p, p)  # p^2 (Minkowski metric)
    return 1j / (p2 - m**2 + 1j * epsilon)

def tree_amplitude_22(lambda_):
    """Tree amplitude for 2→2 scattering (s-channel)"""
    return -1j * lambda_

def tree_amplitude_s_channel(s, lambda_, m):
    """Amplitude with s-channel exchange"""
    # s-channel: p1 + p2 → (exchanged particle) → p3 + p4
    # Amplitude = (-iλ) × i/(s - m²) × (-iλ)
    return -lambda_**2 / (s - m**2)

def symmetry_factor(diagram_type):
    """Calculate symmetry factor"""
    factors = {
        '4pt_contact': 1,        # 4-point contact
        'tadpole': 2,            # Tadpole
        'self_energy': 2,       # Self-energy
        'vacuum_bubble': 8,     # Vacuum bubble (2-loop)
    }
    return factors.get(diagram_type, 1)

# Parameters
lambda_ = 0.1
m = 1.0
s = 10.0  # Mandelstam s

M_tree = tree_amplitude_22(lambda_)
M_s_channel = tree_amplitude_s_channel(s, lambda_, m)

print("Amplitude Calculation by Feynman Rules:")
print("=" * 50)
print(f"4-point contact: M = {M_tree:.6f}")
print(f"s-channel exchange: M = {M_s_channel:.6f}")
print(f"\nSymmetry factor examples:")
print(f"  4-point contact: 1/{symmetry_factor('4pt_contact')}")
print(f"  Tadpole: 1/{symmetry_factor('tadpole')}")

import numpy as np
from scipy.special import gamma

# ===================================
# 1-loop integral with dimensional regularization
# ===================================

def one_loop_integral_dim_reg(m, d=4, mu=1.0):
    """1-loop integral with dimensional regularization

    I = ∫ d^d k / [(2π)^d (k² - m² + iε)]

    Args:
        m: mass
        d: spacetime dimension (d = 4 - 2ε)
        mu: renormalization scale
    """
    epsilon = (4 - d) / 2

    # Result of dimensional regularization (MS-bar scheme)
    if epsilon < 1e-6:
        # Pole at ε → 0
        gamma_E = 0.5772156649  # Euler constant
        result = -1j * m**2 / (16 * np.pi**2) * (
            1 / epsilon - gamma_E + np.log(4 * np.pi) - np.log(m**2 / mu**2)
        )
    else:
        # Evaluation at finite ε
        result = -1j * m**2 / (16 * np.pi**2) * (m / mu)**(−2 * epsilon)

    return result

def self_energy_phi4(p2, lambda_, m, mu=1.0):
    """1-loop self-energy in φ⁴ theory"""
    I_loop = one_loop_integral_dim_reg(m, d=3.99, mu=mu)
    Sigma = lambda_ / 2 * I_loop

    return Sigma

# Parameters
lambda_ = 0.1
m = 1.0
mu = 1.0  # Renormalization scale
p2 = 0.0  # External momentum

Sigma = self_energy_phi4(p2, lambda_, m, mu)

print("1-Loop Self-Energy (Dimensional Regularization):")
print("=" * 50)
print(f"Coupling constant λ: {lambda_}")
print(f"Mass m: {m}")
print(f"Renormalization scale μ: {mu}")
print(f"\nSelf-energy Σ(p²=0): {Sigma:.6f}")
print(f"Real part (mass shift): {Sigma.real:.6e}")
print(f"Imaginary part (decay width): {Sigma.imag:.6e}")

import numpy as np

# ===================================
# Feynman amplitude for Møller scattering
# ===================================

def moller_amplitude(s, t, e, m_e):
    """Invariant amplitude for Møller scattering (spin-averaged)

    e⁻(p1) + e⁻(p2) → e⁻(p3) + e⁻(p4)

    Args:
        s, t: Mandelstam variables
        e: electric charge
        m_e: electron mass
    """
    # Contributions from t-channel and u-channel
    u = 4 * m_e**2 - s - t

    # |M|² after spin averaging (simplified version)
    M2_t = (2 * e**2 / t)**2 * (s**2 + u**2)
    M2_u = (2 * e**2 / u)**2 * (s**2 + t**2)
    M2_int = -(4 * e**4 / (t * u)) * s**2  # Interference term

    M2_avg = (M2_t + M2_u + M2_int) / 4  # Spin average (1/4)

    return M2_avg

def differential_cross_section_moller(s, theta, alpha=1/137, m_e=0.511):
    """Differential cross section dσ/dΩ for Møller scattering"""
    # Kinematics
    p_cm = np.sqrt(s / 4 - m_e**2)
    t = -2 * p_cm**2 * (1 - np.cos(theta))

    e = np.sqrt(4 * np.pi * alpha)
    M2 = moller_amplitude(s, t, e, m_e)

    # Cross section (correction to Mott formula)
    dsigma_dOmega = M2 / (64 * np.pi**2 * s)

    return dsigma_dOmega

# Parameters (low-energy electron scattering)
E_cm = 2.0  # MeV
m_e = 0.511  # MeV
s = E_cm**2
theta = np.pi / 2  # 90 degree scattering

dsigma = differential_cross_section_moller(s, theta, m_e=m_e)

print("Møller Scattering (e⁻e⁻ → e⁻e⁻):")
print("=" * 50)
print(f"Center of mass energy: {E_cm} MeV")
print(f"Scattering angle: {np.degrees(theta):.0f}°")
print(f"Differential cross section dσ/dΩ: {dsigma:.6e} MeV⁻²")
print(f"(barn units: {dsigma * 0.3894:.6e} mb/sr)")

import numpy as np
from scipy.integrate import quad

# ===================================
# Charge renormalization by vacuum polarization
# ===================================

def vacuum_polarization(q2, m_e, mu=1.0):
    """Vacuum polarization function Π(q²)"""

    def integrand(x):
        return x * (1 - x) * np.log(m_e**2 - x * (1 - x) * q2 + 1e-10)

    integral, _ = quad(integrand, 0, 1)

    alpha = 1 / 137
    Pi = -alpha / (3 * np.pi) * (integral - np.log(mu**2) / 6)

    return Pi

def running_coupling(q2, m_e, alpha_0=1/137):
    """Running coupling constant α(q²) = α₀ / (1 - Π(q²))"""
    Pi = vacuum_polarization(q2, m_e)
    alpha_q = alpha_0 / (1 - Pi)

    return alpha_q

# Energy scales
m_e = 0.511  # MeV
Q_values = [1, 10, 100, 1000]  # MeV

print("Running Fine Structure Constant α(Q²):")
print("=" * 60)
print(f"{'Q (MeV)':<15} {'α(Q²)':<20} {'Δα/α (%)':<20}")
print("-" * 60)

alpha_0 = 1 / 137

for Q in Q_values:
    q2 = Q**2
    alpha_Q = running_coupling(q2, m_e, alpha_0)
    delta_alpha = (alpha_Q - alpha_0) / alpha_0 * 100

    print(f"{Q:<15} {alpha_Q:<20.8f} {delta_alpha:<20.6f}")