Chapter 1: Field Quantization and Canonical Formalism

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq

# ===================================
# Time evolution of 1D Klein-Gordon field
# ===================================

def klein_gordon_evolution(phi_init, L=10.0, N=128, m=1.0, T=5.0, dt=0.01):
    """Solve Klein-Gordon equation time evolution using spectral method

    Args:
        phi_init: Initial field configuration
        L: System size
        N: Number of lattice points
        m: Mass
        T: Total time
        dt: Time step

    Returns:
        x, t_array, phi_xt: Spatial coordinates, time array, field spacetime evolution
    """
    x = np.linspace(0, L, N, endpoint=False)
    k = 2 * np.pi * fftfreq(N, L/N)  # Momentum space

    # Dispersion relation: ω(k) = sqrt(k^2 + m^2)
    omega_k = np.sqrt(k**2 + m**2)

    # Initial condition: phi(x,0) and pi(x,0) = ∂_t phi(x,0)
    phi = phi_init.copy()
    pi = np.zeros_like(phi)  # Initial momentum is zero

    # Time evolution array
    n_steps = int(T / dt)
    t_array = np.linspace(0, T, n_steps)
    phi_xt = np.zeros((n_steps, N))
    phi_xt[0] = phi

    for i in range(1, n_steps):
        # Time evolution in Fourier space (split-step method)
        phi_k = fft(phi)
        pi_k = fft(pi)

        # Time evolution operator: exp(-iωt) and exp(iωt)
        phi_k_new = phi_k * np.cos(omega_k * dt) + (pi_k / omega_k) * np.sin(omega_k * dt)
        pi_k_new = pi_k * np.cos(omega_k * dt) - phi_k * omega_k * np.sin(omega_k * dt)

        phi = ifft(phi_k_new).real
        pi = ifft(pi_k_new).real
        phi_xt[i] = phi

    return x, t_array, phi_xt

# Example run: Time evolution of Gaussian wave packet
L, N = 10.0, 128
x = np.linspace(0, L, N, endpoint=False)
phi_init = np.exp(-((x - L/2)**2) / 0.5)  # Gaussian

x, t_array, phi_xt = klein_gordon_evolution(phi_init, m=1.0)

print(f"Number of time steps: {len(t_array)}")
print(f"Energy conservation check: φ(t=0) range [{phi_xt[0].min():.3f}, {phi_xt[0].max():.3f}]")
print(f"                          φ(t=T) range [{phi_xt[-1].min():.3f}, {phi_xt[-1].max():.3f}]")

from sympy import *
from sympy.physics.quantum import *

# ===================================
# Commutation relations of creation-annihilation operators
# ===================================

class AnnihilationOp(Operator):
    """Annihilation operator a"""
    pass

class CreationOp(Operator):
    """Creation operator a†"""
    pass

def commutator(A, B):
    """Commutator [A, B]"""
    return A*B - B*A

# Symbol definition
a = AnnihilationOp('a')
a_dag = CreationOp('a†')

# Verification of canonical commutation relations
print("Canonical commutation relation check:")
print(f"Assume [a, a†] = 1")

# Number operator N = a†a
print("\nProperties of number operator N = a†a:")
print("[a, N] = [a, a†a] = [a, a†]a + a†[a, a] = a")
print("[a†, N] = [a†, a†a] = [a†, a†]a + a†[a†, a] = -a†")

# Action on Fock states
n = Symbol('n', integer=True, positive=True)
print("\nAction on Fock state |n⟩:")
print(f"a |n⟩ = √n |n-1⟩")
print(f"a† |n⟩ = √(n+1) |n+1⟩")
print(f"N |n⟩ = n |n⟩")

import numpy as np
from itertools import combinations_with_replacement

# ===================================
# States and energy calculation in Fock space
# ===================================

def fock_state_energy(k_list, m=1.0):
    """Calculate energy of a Fock state

    Args:
        k_list: List of momenta (each element is a 3D vector)
        m: Particle mass

    Returns:
        Energy eigenvalue
    """
    energy = 0.0
    for k in k_list:
        k_mag = np.linalg.norm(k)
        omega_k = np.sqrt(k_mag**2 + m**2)
        energy += omega_k
    return energy

def generate_fock_states(k_modes, max_particles=3):
    """Generate many-particle Fock states from allowed momentum modes

    Args:
        k_modes: List of possible momentum modes
        max_particles: Maximum number of particles

    Returns:
        fock_states: List of Fock states (each state is a tuple of momenta)
    """
    fock_states = []
    for n in range(max_particles + 1):
        for state in combinations_with_replacement(range(len(k_modes)), n):
            k_list = [k_modes[i] for i in state]
            fock_states.append(k_list)
    return fock_states

# Example for 1D system: k = 0, ±π/L
L = 5.0
k_modes = [
    np.array([0.0]),
    np.array([np.pi / L]),
    np.array([-np.pi / L])
]

fock_states = generate_fock_states(k_modes, max_particles=2)

print("Low-energy states in Fock space:")
print("-" * 50)

for i, state in enumerate(fock_states[:8]):
    n_particles = len(state)
    energy = fock_state_energy(state, m=1.0)

    if n_particles == 0:
        label = "|0⟩ (vacuum)"
    else:
        k_values = [f"k={k[0]:.3f}" for k in state]
        label = f"|{', '.join(k_values)}⟩"

    print(f"{i+1}. {label:<30} E = {energy:.4f}")

import numpy as np

# ===================================
# Simulate anticommutation relations of Fermi operators
# (finite-dimensional approximation with matrix representation)
# ===================================

def fermi_operators(n_states):
    """Construct creation-annihilation operators for n independent Fermi modes

    Fock space dimension: 2^n (each mode has 2 states: occupied/unoccupied)

    Args:
        n_states: Number of Fermi modes

    Returns:
        c: List of annihilation operators (each element is 2^n × 2^n matrix)
        c_dag: List of creation operators
    """
    dim = 2**n_states  # Fock space dimension
    c = []
    c_dag = []

    for i in range(n_states):
        # Annihilation operator for i-th mode
        op = np.zeros((dim, dim), dtype=complex)

        for state in range(dim):
            if (state >> i) & 1:  # i-th mode is occupied
                new_state = state ^ (1 << i)  # Flip i-th bit

                # Jordan-Wigner sign: parity of occupation to the left
                sign = (-1)**bin(state & ((1 << i) - 1)).count('1')

                op[new_state, state] = sign

        c.append(op)
        c_dag.append(op.conj().T)

    return c, c_dag

def anticommutator(A, B):
    """Anticommutator {A, B} = AB + BA"""
    return A @ B + B @ A

# Example with 3 Fermi modes
n_states = 3
c, c_dag = fermi_operators(n_states)

print("Verification of Fermi operator anticommutation relations:")
print("=" * 50)

# Verification of {c_i, c_j†} = δ_ij
print("\n1. {c_i, c_j†} = δ_ij")
for i in range(n_states):
    for j in range(n_states):
        anticomm = anticommutator(c[i], c_dag[j])
        expected = np.eye(2**n_states) if i == j else np.zeros((2**n_states, 2**n_states))
        is_correct = np.allclose(anticomm, expected)
        print(f"   {{c_{i}, c†_{j}}} = δ_{i}{j}: {is_correct}")

# Verification of {c_i, c_j} = 0
print("\n2. {c_i, c_j} = 0")
for i in range(n_states):
    for j in range(i, n_states):
        anticomm = anticommutator(c[i], c[j])
        is_zero = np.allclose(anticomm, 0)
        print(f"   {{c_{i}, c_{j}}} = 0: {is_zero}")

# Pauli exclusion principle: (c†)^2 = 0
print("\n3. Pauli exclusion principle: (c†_i)^2 = 0")
for i in range(n_states):
    square = c_dag[i] @ c_dag[i]
    is_zero = np.allclose(square, 0)
    print(f"   (c†_{i})^2 = 0: {is_zero}")

import numpy as np
from itertools import combinations

# ===================================
# Numerical verification of Wick's theorem (harmonic oscillator example)
# ===================================

def harmonic_operators(n_max):
    """Creation-annihilation operators in Fock space of harmonic oscillator

    Args:
        n_max: Maximum occupation number (restrict Fock space to |0⟩, |1⟩, ..., |n_max⟩)

    Returns:
        a: Annihilation operator (matrix)
        a_dag: Creation operator (matrix)
    """
    dim = n_max + 1
    a = np.zeros((dim, dim))

    for n in range(1, dim):
        a[n-1, n] = np.sqrt(n)

    a_dag = a.T

    return a, a_dag

def normal_order(ops, n_max):
    """Rearrange product of operators into normal order

    Args:
        ops: List of operators ('a' or 'a_dag')
        n_max: Maximum occupation number of Fock space

    Returns:
        Normal-ordered product of operators (matrix)
    """
    a, a_dag = harmonic_operators(n_max)

    # Creation operators to the left, annihilation operators to the right
    creation_ops = [a_dag for op in ops if op == 'a_dag']
    annihilation_ops = [a for op in ops if op == 'a']

    result = np.eye(n_max + 1)
    for op in creation_ops + annihilation_ops:
        result = result @ op

    return result

def compute_contraction(op1, op2, n_max):
    """Calculate contraction of two operators"""
    a, a_dag = harmonic_operators(n_max)

    if op1 == 'a' and op2 == 'a_dag':
        return a @ a_dag - a_dag @ a  # [a, a†]
    else:
        return np.zeros((n_max + 1, n_max + 1))

# Verify Wick's theorem: expand a a† a a†
n_max = 5
a, a_dag = harmonic_operators(n_max)

# Left side: a a† a a†
lhs = a @ a_dag @ a @ a_dag

# Right side: expansion by Wick's theorem
# :a a† a a†: + Contraction(a,a†) :a a†: + Contraction(a,a†) :a† a: + contraction product

# Normal-ordered product: :a a† a a†: = a†^2 a^2
normal = a_dag @ a_dag @ a @ a

# Contraction calculation
contraction_1 = compute_contraction('a', 'a_dag', n_max) @ (a @ a_dag)
contraction_2 = compute_contraction('a', 'a_dag', n_max) @ (a_dag @ a)
contraction_both = compute_contraction('a', 'a_dag', n_max) @ compute_contraction('a', 'a_dag', n_max)

rhs = normal + contraction_1 + contraction_2 + contraction_both

print("Verification of Wick's theorem: a a† a a†")
print("=" * 50)
print(f"Maximum difference between direct calculation and Wick expansion: {np.max(np.abs(lhs - rhs)):.10f}")
print(f"\nVacuum expectation value ⟨0|a a† a a†|0⟩:")
print(f"  Direct calculation: {lhs[0, 0]:.4f}")
print(f"  Wick's theorem: {rhs[0, 0]:.4f}")

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Phonon dispersion relation for 1D atomic chain
# ===================================

def phonon_dispersion_1d(k, K, M, a):
    """Phonon dispersion relation for 1D atomic chain

    Args:
        k: Wave number (can be array)
        K: Spring constant
        M: Atomic mass
        a: Lattice constant

    Returns:
        omega: Angular frequency
    """
    return 2 * np.sqrt(K / M) * np.abs(np.sin(k * a / 2))

def phonon_dos_1d(omega, K, M, a, n_points=1000):
    """Density of states for 1D phonons

    Args:
        omega: Angular frequency (array)
        K, M, a: System parameters
        n_points: Number of integration points

    Returns:
        dos: Density of states g(ω)
    """
    k_max = np.pi / a
    k = np.linspace(-k_max, k_max, n_points)
    omega_k = phonon_dispersion_1d(k, K, M, a)

    dos = np.zeros_like(omega)
    dk = k[1] - k[0]

    for i, om in enumerate(omega):
        # Approximate δ(ω - ω(k)) with Gaussian of small width
        delta_width = 0.01 * (omega[-1] - omega[0])
        delta_approx = np.exp(-((omega_k - om)**2) / (2 * delta_width**2))
        delta_approx /= (np.sqrt(2 * np.pi) * delta_width)

        dos[i] = np.sum(delta_approx) * dk / (2 * np.pi)

    return dos

# Parameter setting (silicon crystal assumption)
K = 50.0  # N/m
M = 28.0855 * 1.66e-27  # Si atom mass (kg)
a = 5.43e-10  # Lattice constant (m)

# Wave number range
k = np.linspace(-np.pi/a, np.pi/a, 200)
omega = phonon_dispersion_1d(k, K, M, a)

# Density of states in frequency range
omega_range = np.linspace(0, np.max(omega), 100)
dos = phonon_dos_1d(omega_range, K, M, a)

print("Phonon properties of 1D atomic chain:")
print("=" * 50)
print(f"Maximum phonon frequency: {np.max(omega)/(2*np.pi)*1e-12:.2f} THz")
print(f"Sound velocity (long-wavelength limit): {2*np.sqrt(K/M)*a:.2f} m/s")
print(f"Zero-point energy (per mode): {0.5*1.055e-34*np.max(omega)*1e3:.2e} meV")

import numpy as np

# ===================================
# Magnon dispersion in the Heisenberg model
# ===================================

def magnon_dispersion(k, J, S, a):
    """Magnon dispersion in 1D Heisenberg ferromagnet

    Hamiltonian: H = -J Σ S_i · S_{i+1}

    Args:
        k: Wave number
        J: Exchange interaction constant (J > 0 for ferromagnetism)
        S: Spin quantum number
        a: Lattice constant

    Returns:
        omega: Magnon excitation energy
    """
    return 2 * J * S * (1 - np.cos(k * a))

def magnon_energy_gap(J, S, d, B_ext=0.0):
    """Magnon energy gap including anisotropy and external magnetic field

    Args:
        J: Exchange interaction constant
        S: Spin quantum number
        d: Anisotropy constant
        B_ext: External magnetic field

    Returns:
        gap: Energy gap
    """
    g_factor = 2.0
    mu_B = 9.274e-24  # Bohr magneton (J/T)

    return d * S + g_factor * mu_B * B_ext

# Parameters (iron example)
J = 1.0e-20  # J (Joules)
S = 1.0  # Spin quantum number
a = 2.87e-10  # Lattice constant (m)

# Wave number
k = np.linspace(0, 2*np.pi/a, 100)
omega = magnon_dispersion(k, J, S, a)

# Physical quantity calculation
k_small = 1e8  # Small wave number (1/m)
omega_k_small = magnon_dispersion(k_small, J, S, a)
spin_wave_stiffness = omega_k_small / k_small**2

print("Magnons in Heisenberg ferromagnet:")
print("=" * 50)
print(f"Maximum excitation energy: {np.max(omega)*6.242e18:.2f} eV")
print(f"Spin wave stiffness: {spin_wave_stiffness:.2e} J·m^2")
print(f"Long-wavelength limit energy: E(k) ≈ D k^2, D = {spin_wave_stiffness:.2e}")

import numpy as np
from scipy.integrate import quad

# ===================================
# Bose distribution and thermal properties
# ===================================

def bose_einstein(omega, T):
    """Bose-Einstein distribution function

    Args:
        omega: Energy (angular frequency)
        T: Temperature (K)

    Returns:
        n(ω, T): Average occupation number
    """
    k_B = 1.381e-23  # Boltzmann constant (J/K)
    hbar = 1.055e-34  # Planck constant (J·s)

    if T == 0:
        return 0.0

    x = hbar * omega / (k_B * T)
    if x > 50:  # Prevent overflow
        return 0.0

    return 1.0 / (np.exp(x) - 1)

def thermal_energy(omega_k_func, T, k_range, dim=1):
    """Thermal energy of phonons/magnons

    Args:
        omega_k_func: Function of dispersion relation ω(k)
        T: Temperature (K)
        k_range: (k_min, k_max)
        dim: Dimension

    Returns:
        E: Total thermal energy
    """
    k_B = 1.381e-23
    hbar = 1.055e-34

    def integrand(k):
        omega = omega_k_func(k)
        n_BE = bose_einstein(omega, T)
        return hbar * omega * n_BE

    if dim == 1:
        result, _ = quad(integrand, k_range[0], k_range[1])
        return result / (2 * np.pi)
    else:
        raise NotImplementedError("Only 1D implemented")

# Phonon parameters
K, M, a = 50.0, 28.0855 * 1.66e-27, 5.43e-10
omega_phonon = lambda k: 2 * np.sqrt(K / M) * np.abs(np.sin(k * a / 2))

# Magnon parameters
J, S = 1.0e-20, 1.0
omega_magnon = lambda k: 2 * J * S * (1 - np.cos(k * a)) / 1.055e-34

# Temperature range
temperatures = [10, 50, 100, 300]  # K

print("Thermal energy comparison of phonons and magnons:")
print("=" * 60)
print(f"{'T (K)':<10} {'Phonon (J/m)':<20} {'Magnon (J/m)':<20}")
print("-" * 60)

for T in temperatures:
    E_phonon = thermal_energy(omega_phonon, T, (0, np.pi/a))
    E_magnon = thermal_energy(omega_magnon, T, (0, np.pi/a))

    print(f"{T:<10} {E_phonon:<20.3e} {E_magnon:<20.3e}")