Peskin & Schroeder Solutions

Name: Quantum Field Theory — Free Online Course
Availability: InStock

Video solutions to problems from "An Introduction to Quantum Field Theory"

Problems Solved

Chapters Covered

Video Quality

Free

Open Access

About the Textbook

"An Introduction to Quantum Field Theory" by Michael E. Peskin and Daniel V. Schroeder is one of the most widely used graduate-level textbooks in quantum field theory. First published in 1995, it has become the standard reference for QFT courses worldwide.

The textbook covers: relativistic wave equations, Feynman diagrams, QED, renormalization, non-Abelian gauge theories, the Standard Model, and more advanced topics. These video solutions help students work through the challenging problem sets.

Chapter 2: The Klein-Gordon Field

5 problems solved

Ch. 2Problem 2.1 a)

Classical electromagnetism as a field theory

Ch. 2Problem 2.1 b)

Noether's theorem application

Ch. 2Problem 2.2 a)

Complex scalar field

Ch. 2Problem 2.2 b)

Conserved current for complex field

Ch. 2Problem 2.2 c)

Quantization of complex field

Chapter 3: The Dirac Field

3 problems solved

Ch. 3Problem 3.1

Lorentz group representations and spinor transformations

Ch. 3Problem 3.2

Gordon identity derivation

Ch. 3Problem 3.3

Dirac bilinear transformation properties

Chapter 4: Interacting Fields and the S-Matrix

3 problems solved

Ch. 4Problem 4.1

Wick's theorem application to time-ordered products

Ch. 4Problem 4.2

Feynman propagator computation in position and momentum space

Ch. 4Problem 4.3

Scattering amplitude in phi-four theory

Chapter 5: Elementary Processes of QED

3 problems solved

Ch. 5Problem 5.1

e+e- to mu+mu- differential cross section

Ch. 5Problem 5.2

Compton scattering and Klein-Nishina formula

Ch. 5Problem 5.3

Bhabha scattering: e+e- to e+e-

Chapter 6: Radiative Corrections

2 problems solved

Ch. 6Problem 6.1

Vertex correction and anomalous magnetic moment (g-2)

Ch. 6Problem 6.2

Electron self-energy and mass renormalization

Chapter 7: Renormalization

2 problems solved

Ch. 7Problem 7.1

Superficial degree of divergence in QED

Ch. 7Problem 7.2

Running coupling constant and renormalization group

Chapter 9: Functional Methods

2 problems solved

Ch. 9Problem 9.1

Generating functional for the free scalar field

Ch. 9Problem 9.2

Path integral for the quantum harmonic oscillator

Chapter 3: The Dirac Field — Solution Details

The Dirac equation is the cornerstone of relativistic quantum mechanics for spin-1/2 particles. The central equation is:

$$\boxed{(i\gamma^\mu \partial_\mu - m)\psi = 0}$$

where the gamma matrices satisfy the Clifford algebra anticommutation relations:

$$\{\gamma^\mu, \gamma^\nu\} = 2g^{\mu\nu}\mathbb{I}_{4\times 4}$$

Problem 3.1: Lorentz Group Representations

This problem explores how spinors transform under Lorentz transformations. Under an infinitesimal Lorentz transformation $\Lambda^\mu{}_\nu = \delta^\mu{}_\nu + \omega^\mu{}_\nu$, a Dirac spinor transforms as:

$$\psi(x) \to S(\Lambda)\,\psi(\Lambda^{-1}x), \quad S(\Lambda) = \exp\!\Bigl(-\frac{i}{4}\omega_{\mu\nu}\,\sigma^{\mu\nu}\Bigr)$$

where $\sigma^{\mu\nu} = \frac{i}{2}[\gamma^\mu, \gamma^\nu]$ are the generators of the Lorentz group in the spinor representation.

Key result: The $(\frac{1}{2},0) \oplus (0,\frac{1}{2})$ representation of the Lorentz group is reducible into left-handed and right-handed Weyl spinors, connected through parity.

Problem 3.2: Gordon Identity

The Gordon identity decomposes the vector current into a convection current and a spin (magnetic moment) contribution:

$$\boxed{\bar{u}(p')\gamma^\mu u(p) = \bar{u}(p')\left[\frac{(p'+p)^\mu}{2m} + \frac{i\sigma^{\mu\nu}(p'-p)_\nu}{2m}\right]u(p)}$$

Derivation outline: Start from the Dirac equation in momentum space,$(\not\!{p} - m)u(p) = 0$, and use the identity$\gamma^\mu \not\!{p} = p^\mu + i\sigma^{\mu\nu}p_\nu$ with the corresponding equation for $\bar{u}(p')$. Add the two expressions and divide by $2m$.

The Gordon identity is essential for extracting the anomalous magnetic moment in Chapter 6.

Problem 3.3: Dirac Bilinear Transformation Properties

The 16 independent Dirac bilinears $\bar{\psi}\,\Gamma\,\psi$ transform under the Lorentz group as follows:

Bilinear	Components	Transforms as
$\bar{\psi}\psi$	1	Scalar
$\bar{\psi}\gamma^\mu\psi$	4	Vector
$\bar{\psi}\sigma^{\mu\nu}\psi$	6	Antisymmetric tensor
$\bar{\psi}\gamma^\mu\gamma^5\psi$	4	Axial vector
$\bar{\psi}\gamma^5\psi$	1	Pseudoscalar

These 16 bilinears form a complete basis for $4\times 4$ matrices, reflecting the decomposition $1 + 4 + 6 + 4 + 1 = 16$.

Python: Dirac Gamma Matrix Algebra Verification

Numerically verifies the Clifford algebra anticommutation relations$\{\gamma^\mu, \gamma^\nu\} = 2g^{\mu\nu}$ and computes traces of gamma matrix products.

Dirac Gamma Matrix Algebra

Python

Verifies the Clifford algebra {gamma^mu, gamma^nu} = 2g^{mu,nu} and trace identities

script.py79 lines

import numpy as np

# ──────────────────────────────────────────────
# Dirac Gamma Matrices (Dirac/Pauli representation)
# Verify Clifford algebra: {gamma^mu, gamma^nu} = 2 g^{mu nu}
# ──────────────────────────────────────────────

I2 = np.eye(2, dtype=complex)
zero2 = np.zeros((2, 2), dtype=complex)

# Pauli matrices
sigma1 = np.array([[0, 1], [1, 0]], dtype=complex)
sigma2 = np.array([[0, -1j], [1j, 0]], dtype=complex)
sigma3 = np.array([[1, 0], [0, -1]], dtype=complex)
pauli = [sigma1, sigma2, sigma3]

# Gamma matrices in Dirac representation
gamma0 = np.block([[I2, zero2], [zero2, -I2]])
gamma1 = np.block([[zero2, sigma1], [-sigma1, zero2]])
gamma2 = np.block([[zero2, sigma2], [-sigma2, zero2]])
gamma3 = np.block([[zero2, sigma3], [-sigma3, zero2]])
gamma = [gamma0, gamma1, gamma2, gamma3]

# gamma^5 = i gamma^0 gamma^1 gamma^2 gamma^3
gamma5 = 1j * gamma0 @ gamma1 @ gamma2 @ gamma3

# Metric tensor (mostly minus)
g = np.diag([1.0, -1.0, -1.0, -1.0])

print("=" * 55)
print("  Clifford Algebra Verification")
print("  {gamma^mu, gamma^nu} = 2 g^{mu,nu} I_4")
print("=" * 55)

all_ok = True
for mu in range(4):
    for nu in range(4):
        anticomm = gamma[mu] @ gamma[nu] + gamma[nu] @ gamma[mu]
        expected = 2.0 * g[mu, nu] * np.eye(4)
        match = np.allclose(anticomm, expected)
        if not match:
            all_ok = False
        if mu <= nu:
            status = "PASS" if match else "FAIL"
            print(f"  {{gamma^{mu}, gamma^{nu}}} = {2*g[mu,nu]:+.0f} I_4  [{status}]")

print(f"\nAll anticommutation relations: {'PASSED' if all_ok else 'FAILED'}")

print("\n" + "=" * 55)
print("  Useful Trace Identities")
print("=" * 55)

I4 = np.eye(4, dtype=complex)
print(f"  Tr(I_4)           = {np.trace(I4).real:.0f}  (expected 4)")
print(f"  Tr(gamma^mu)      = {np.trace(gamma0).real:.0f}  (expected 0)")
print(f"  Tr(gamma^5)       = {np.trace(gamma5).real:.0f}  (expected 0)")

tr_g0g1 = np.trace(gamma[0] @ gamma[1])
print(f"  Tr(gamma^0 gamma^1) = {tr_g0g1.real:.0f}  (expected 0)")

tr_g0g0 = np.trace(gamma[0] @ gamma[0])
print(f"  Tr(gamma^0 gamma^0) = {tr_g0g0.real:.0f}  (expected 4)")

tr_4gamma = np.trace(gamma[0] @ gamma[1] @ gamma[2] @ gamma[3])
print(f"  Tr(gamma^0 gamma^1 gamma^2 gamma^3) = {tr_4gamma.real:+.0f}i  (expected 0)")

# Tr(gamma^mu gamma^nu) = 4 g^{mu,nu}
print("\n  Tr(gamma^mu gamma^nu) = 4 g^{mu,nu}:")
for mu in range(4):
    for nu in range(4):
        tr_val = np.trace(gamma[mu] @ gamma[nu]).real
        if mu <= nu and abs(tr_val) > 1e-10:
            print(f"    Tr(gamma^{mu} gamma^{nu}) = {tr_val:+.0f}  (4 g^{{{mu},{nu}}} = {4*g[mu,nu]:+.0f})")

print("\n  (gamma^5)^2 = I?", "YES" if np.allclose(gamma5 @ gamma5, I4) else "NO")
print("  {gamma^5, gamma^mu} = 0?", "YES" if all(
    np.allclose(gamma5 @ gamma[mu] + gamma[mu] @ gamma5, np.zeros((4,4)))
    for mu in range(4)) else "NO")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter 4: Interacting Fields — Solution Details

Problem 4.1: Wick's Theorem

Wick's theorem relates time-ordered products to normal-ordered products plus all possible contractions. For a product of $n$ fields:

$$T\{\phi(x_1)\cdots\phi(x_n)\} = N\{\phi(x_1)\cdots\phi(x_n) + \text{all contractions}\}$$

Each contraction of two fields yields the Feynman propagator:

$$\overline{\phi(x)\phi(y)} = D_F(x-y) = \int \frac{d^4k}{(2\pi)^4}\frac{i}{k^2 - m^2 + i\epsilon}e^{-ik\cdot(x-y)}$$

The theorem is applied by systematically enumerating all ways to pair up fields. For four fields, there are 3 distinct complete contractions, corresponding to the three $s$-, $t$-, and $u$-channel diagrams.

Problem 4.2: Feynman Propagator

The Feynman propagator for a massive scalar field can be computed directly. In position space, for a massive field:

$$D_F(x) = \begin{cases} -\frac{m}{4\pi^2\sqrt{-x^2}}K_1(m\sqrt{-x^2}) & \text{for } x^2 < 0 \text{ (spacelike)} \\ \frac{im}{4\pi^2\sqrt{x^2}}K_1(-im\sqrt{x^2}) & \text{for } x^2 > 0 \text{ (timelike)} \end{cases}$$

where $K_1$ is the modified Bessel function. Key properties:

Falls off exponentially for spacelike separations: $D_F \sim e^{-m|x|}$
Satisfies $(\Box + m^2)D_F(x) = -i\delta^{(4)}(x)$
The $i\epsilon$ prescription selects positive-frequency solutions propagating forward in time

Problem 4.3: $\phi^4$ Scattering Amplitude

For $\phi^4$ theory with Lagrangian$\mathcal{L} = \frac{1}{2}(\partial\phi)^2 - \frac{1}{2}m^2\phi^2 - \frac{\lambda}{4!}\phi^4$, the leading-order $2\to 2$ scattering amplitude is simply:

$$i\mathcal{M} = -i\lambda$$

At next-to-leading order ($\mathcal{O}(\lambda^2)$), the three one-loop diagrams contribute:

$$i\mathcal{M}^{(1)} = \frac{(-i\lambda)^2}{2}\left[\int\frac{d^4k}{(2\pi)^4}\frac{i}{k^2-m^2}\frac{i}{(k+p)^2-m^2}\right]_{s,t,u\text{-channels}}$$

The Mandelstam variables are $s = (p_1+p_2)^2$, $t = (p_1-p_3)^2$,$u = (p_1-p_4)^2$, satisfying $s + t + u = 4m^2$.

Fortran: Numerical One-Loop Feynman Parameter Integral

Evaluates the one-loop bubble integral using Feynman parametrization and Gaussian quadrature. The integral $I(p^2) = \int_0^1 dx\, \ln\!\bigl[m^2 - x(1-x)p^2\bigr]$ appears in the $\phi^4$ one-loop correction.

One-Loop Feynman Parameter Integral

Fortran

Numerically evaluates the bubble integral I(p^2) for phi^4 theory using Gaussian quadrature

one_loop_integral.f9064 lines

program one_loop_integral
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: N_GAUSS = 100
  real(dp), parameter :: PI = 3.14159265358979323846_dp
  real(dp) :: m2, p2, x, dx, integral, arg
  real(dp) :: s_vals(7), result
  integer :: i, j

m2 = 1.0_dp  ! mass squared in natural units

write(*,'(A)') '================================================='
  write(*,'(A)') '  One-Loop Bubble Integral (Feynman parameter)'
  write(*,'(A)') '  I(p^2) = int_0^1 dx ln[m^2 - x(1-x) p^2]'
  write(*,'(A)') '  m^2 = 1.0'
  write(*,'(A)') '================================================='
  write(*,'(A)') ''
  write(*,'(A)') '    p^2/m^2     I(p^2)/(16 pi^2)    Threshold?'
  write(*,'(A)') '   ---------   -----------------   -----------'

! Evaluate for various p^2 values (below and above threshold 4m^2)
  s_vals(1) = 0.0_dp
  s_vals(2) = 1.0_dp
  s_vals(3) = 2.0_dp
  s_vals(4) = 3.0_dp
  s_vals(5) = 3.9_dp
  s_vals(6) = 4.0_dp   ! threshold
  s_vals(7) = 5.0_dp

do j = 1, 7
    p2 = s_vals(j) * m2
    integral = 0.0_dp
    dx = 1.0_dp / real(N_GAUSS, dp)

do i = 1, N_GAUSS
      x = (real(i, dp) - 0.5_dp) * dx
      arg = m2 - x * (1.0_dp - x) * p2

if (arg > 0.0_dp) then
        integral = integral + log(arg) * dx
      else
        ! Above threshold: pick up imaginary part
        integral = integral + log(abs(arg)) * dx
      end if
    end do

result = integral / (16.0_dp * PI**2)

if (s_vals(j) < 3.99_dp) then
      write(*,'(F10.2, F20.8, A)') s_vals(j), result, '    below'
    else if (s_vals(j) < 4.01_dp) then
      write(*,'(F10.2, F20.8, A)') s_vals(j), result, '    AT THRESHOLD'
    else
      write(*,'(F10.2, F20.8, A)') s_vals(j), result, '    ABOVE (has Im part)'
    end if
  end do

write(*,'(A)') ''
  write(*,'(A)') '  Note: At p^2 = 4m^2 the argument vanishes at x=1/2,'
  write(*,'(A)') '  producing a branch cut (the two-particle threshold).'
  write(*,'(A)') '  Above threshold, the amplitude acquires an imaginary part'
  write(*,'(A)') '  related to the total cross section via the optical theorem.'

end program one_loop_integral

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Chapter 5: Elementary QED Processes — Solution Details

Problem 5.1: $e^+e^- \to \mu^+\mu^-$ Cross Section

The unpolarized differential cross section for $e^+e^- \to \mu^+\mu^-$in the center-of-mass frame, at leading order in QED:

$$\boxed{\frac{d\sigma}{d\Omega} = \frac{\alpha^2}{4s}(1 + \cos^2\theta)}$$

This is valid in the high-energy limit $s \gg m_\mu^2$. The total cross section integrates to:

$$\sigma_{\text{total}} = \frac{4\pi\alpha^2}{3s}$$

The $(1 + \cos^2\theta)$ angular distribution is characteristic of spin-1/2 pair production via a vector (spin-1) intermediate state.

Problem 5.2: Compton Scattering (Klein-Nishina)

The Klein-Nishina formula for unpolarized Compton scattering$e^-\gamma \to e^-\gamma$ in the lab frame:

$$\frac{d\sigma}{d\Omega} = \frac{\alpha^2}{2m_e^2}\left(\frac{\omega'}{\omega}\right)^2\left[\frac{\omega'}{\omega} + \frac{\omega}{\omega'} - \sin^2\theta\right]$$

where $\omega' = \omega/[1 + (\omega/m_e)(1 - \cos\theta)]$ is the scattered photon energy. In the low-energy limit $\omega \ll m_e$, this reduces to the Thomson cross section $\sigma_T = 8\pi\alpha^2/(3m_e^2)$.

Problem 5.3: Bhabha Scattering

Bhabha scattering ($e^+e^- \to e^+e^-$) receives contributions from both the $s$-channel (annihilation) and $t$-channel (scattering) diagrams. The unpolarized differential cross section is:

$$\frac{d\sigma}{d\Omega} = \frac{\alpha^2}{2s}\left[\frac{1+\cos^4(\theta/2)}{\sin^4(\theta/2)} - \frac{2\cos^4(\theta/2)}{\sin^2(\theta/2)} + \frac{1+\cos^2\theta}{2}\right]$$

The first term is the $t$-channel contribution (divergent as $\theta \to 0$), the last is the $s$-channel, and the middle term is the interference.

Python: Differential Cross Section for $e^+e^- \to \mu^+\mu^-$

Plots the differential cross section as a function of $\cos\theta$at several center-of-mass energies, including finite muon mass effects.

e+e- to mu+mu- Cross Section

Python

Differential and total cross sections with finite muon mass effects at various energies

script.py75 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────
# e+e- -> mu+mu- differential cross section
# Including finite muon mass effects
# ──────────────────────────────────────────────

alpha = 1.0 / 137.036    # fine structure constant
m_mu = 0.10566            # muon mass in GeV
hbarc2 = 0.3894e6         # (hbar c)^2 in GeV^2 pb (for unit conversion)

cos_th = np.linspace(-0.99, 0.99, 500)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# --- Left: dsigma/dOmega vs cos(theta) at different sqrt(s) ---
energies = [0.5, 1.0, 5.0, 10.0, 91.2]  # GeV (last = Z pole, tree-level only)
colors_plot = ['#f472b6', '#c084fc', '#60a5fa', '#34d399', '#fbbf24']

for E, col in zip(energies, colors_plot):
    s = E**2
    if s < (2 * m_mu)**2:
        continue
    beta = np.sqrt(1.0 - 4.0 * m_mu**2 / s)
    # Full formula with muon mass
    dsig = (alpha**2 / (4.0 * s)) * beta * (
        1.0 + (1.0 - beta**2) + (1.0 - (1.0 - beta**2)) * cos_th**2
    )
    dsig_pb = dsig * hbarc2  # convert to pb/sr
    axes[0].plot(cos_th, dsig_pb, color=col, linewidth=2,
                 label=f'$\\sqrt{{s}}$ = {E} GeV')

axes[0].set_xlabel('cos $\\\\theta$', color='#94a3b8', fontsize=13)
axes[0].set_ylabel('d$\\sigma$/d$\\Omega$ [pb/sr]', color='#94a3b8', fontsize=13)
axes[0].set_title('$e^+e^- \\\\to \\mu^+\\mu^-$ Differential Cross Section',
                   color='white', fontsize=14)
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')
axes[0].set_yscale('log')

# --- Right: Total cross section vs sqrt(s) ---
sqrts = np.linspace(2 * m_mu + 0.001, 50.0, 500)
s_arr = sqrts**2
beta_arr = np.sqrt(1.0 - 4.0 * m_mu**2 / s_arr)
sigma_tot = (4.0 * np.pi * alpha**2 / (3.0 * s_arr)) * beta_arr * (
    1.0 + (1.0 - beta_arr**2) / 2.0
)
sigma_tot_pb = sigma_tot * hbarc2

axes[1].plot(sqrts, sigma_tot_pb, color='#f472b6', linewidth=2.5)
axes[1].axvline(x=2*m_mu, color='#fbbf24', linestyle='--', alpha=0.6, label=f'Threshold 2$m_\\mu$ = {2*m_mu:.3f} GeV')
axes[1].set_xlabel('$\\sqrt{s}$ [GeV]', color='#94a3b8', fontsize=13)
axes[1].set_ylabel('$\\sigma_{\\mathrm{tot}}$ [pb]', color='#94a3b8', fontsize=13)
axes[1].set_title('Total Cross Section vs Energy', color='white', fontsize=14)
axes[1].set_yscale('log')
axes[1].set_xscale('log')
axes[1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight',
            facecolor='#0a0a1a', edgecolor='none')
print("Plot saved: e+e- -> mu+mu- cross sections")
print(f"  sigma_tot at sqrt(s)=10 GeV: {(4*np.pi*alpha**2/(3*100))*hbarc2:.2f} pb")
print(f"  sigma_tot at sqrt(s)=91 GeV (Z pole, tree QED only): {(4*np.pi*alpha**2/(3*91.2**2))*hbarc2:.4f} pb")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter 6: Radiative Corrections — Solution Details

Problem 6.1: Vertex Correction and Anomalous Magnetic Moment

The one-loop vertex correction modifies the QED vertex$-ie\gamma^\mu$ to $-ie\Gamma^\mu(p',p)$, where:

$$\Gamma^\mu(p',p) = \gamma^\mu F_1(q^2) + \frac{i\sigma^{\mu\nu}q_\nu}{2m}F_2(q^2)$$

The form factors satisfy $F_1(0) = 1$ (charge renormalization) and$F_2(0) = a_e$ gives the anomalous magnetic moment. Schwinger's famous result at one loop:

$$\boxed{a_e = \frac{g-2}{2} = \frac{\alpha}{2\pi} \approx 0.00116}$$

The calculation uses the Feynman parameter integral:

$$F_2(0) = \frac{\alpha}{\pi}\int_0^1 dx\int_0^{1-x}dy\,\frac{2m^2 z}{m^2(1-z)^2} = \frac{\alpha}{2\pi}$$

where $z = 1 - x - y$. This result, first computed by Julian Schwinger in 1948, was a triumph of renormalized QED.

Problem 6.2: Electron Self-Energy

The one-loop electron self-energy $-i\Sigma(\not\!{p})$ is given by:

$$\Sigma(\not\!{p}) = -\frac{\alpha}{4\pi}\int_0^1 dx\,(4m - 2x\not\!{p})\left[\frac{2}{\epsilon} - \gamma_E + \ln\frac{4\pi\mu^2}{\Delta}\right]$$

where $\Delta = -x(1-x)p^2 + xm_\gamma^2 + (1-x)m^2$ and we use dimensional regularization with $d = 4 - \epsilon$. The physical mass is determined by:

$$m_{\text{phys}} = m_0 + \Sigma(m_{\text{phys}})$$

The self-energy contains both UV divergences (requiring renormalization) and an IR divergence from the photon mass $m_\gamma \to 0$, which cancels against soft bremsstrahlung.

Key Results: g-2 and Lamb Shift

Anomalous magnetic moment: Current QED prediction to 5th order:$a_e = \frac{\alpha}{2\pi} - 0.328\,\frac{\alpha^2}{\pi^2} + 1.181\,\frac{\alpha^3}{\pi^3} - 1.912\,\frac{\alpha^4}{\pi^4} + \cdots$

Lamb shift: The $2S_{1/2} - 2P_{1/2}$ splitting in hydrogen arises from the electron self-energy and vacuum polarization:$\Delta E_{\text{Lamb}} \approx \frac{4\alpha^5 m_e}{3\pi n^3}\left[\ln\frac{m_e}{2\bar{E}} + \frac{3}{8} + \cdots\right] \approx 1057\text{ MHz}$

Fortran: Vertex Correction Feynman Parameter Integration

Numerically evaluates the Feynman parameter integral for the one-loop vertex correction form factor $F_2(q^2)$ and verifies $F_2(0) = \alpha/(2\pi)$.

QED Vertex Correction Form Factor

Fortran

Computes F_2(q^2) via Feynman parameter integration, verifying Schwinger's alpha/(2*pi) result

vertex_correction.f9082 lines

program vertex_correction
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: PI = 3.14159265358979323846_dp
  real(dp), parameter :: alpha = 1.0_dp / 137.036_dp
  real(dp) :: m, q2, x, y, z, dx, dy, denom, integrand
  real(dp) :: F2, F2_exact
  integer :: ix, iy, nx, ny
  real(dp) :: q2_vals(6), results(6)
  integer :: iq

m = 1.0_dp    ! electron mass (natural units)
  nx = 500      ! integration points
  ny = 500
  dx = 1.0_dp / real(nx, dp)

write(*,'(A)') '======================================================='
  write(*,'(A)') '  QED Vertex Correction: Form Factor F_2(q^2)'
  write(*,'(A)') '  One-loop Feynman parameter integration'
  write(*,'(A)') '  F_2(0) should equal alpha/(2*pi) = Schwinger result'
  write(*,'(A)') '======================================================='
  write(*,'(A)') ''

F2_exact = alpha / (2.0_dp * PI)
  write(*,'(A,F12.8)') '  Schwinger exact: alpha/(2*pi) = ', F2_exact
  write(*,'(A)') ''

q2_vals(1) = 0.0_dp
  q2_vals(2) = -0.01_dp
  q2_vals(3) = -0.1_dp
  q2_vals(4) = -0.5_dp
  q2_vals(5) = -1.0_dp
  q2_vals(6) = -2.0_dp

write(*,'(A)') '    q^2/m^2       F_2(q^2)        F_2/F_2(0)'
  write(*,'(A)') '   ---------    ------------     -----------'

do iq = 1, 6
    q2 = q2_vals(iq) * m**2
    F2 = 0.0_dp

do ix = 1, nx
      x = (real(ix, dp) - 0.5_dp) * dx
      ! y ranges from 0 to 1-x
      if (1.0_dp - x < 1.0e-12_dp) cycle
      dy = (1.0_dp - x) / real(ny, dp)

do iy = 1, ny
        y = (real(iy, dp) - 0.5_dp) * dy
        z = 1.0_dp - x - y
        if (z < 0.0_dp .or. z > 1.0_dp) cycle

! Denominator of Feynman parameter integral
        denom = m**2 * (1.0_dp - z)**2 - q2 * x * y
        if (denom < 1.0e-15_dp) cycle

! F_2 integrand: 2m^2 z / denom
        integrand = 2.0_dp * m**2 * z / denom
        F2 = F2 + integrand * dx * dy
      end do
    end do

F2 = (alpha / PI) * F2
    results(iq) = F2

if (iq == 1) then
      write(*,'(F10.3, F16.8, A)') q2_vals(iq), F2, '      (Schwinger)'
    else
      write(*,'(F10.3, F16.8, F14.6)') q2_vals(iq), F2, F2 / results(1)
    end if
  end do

write(*,'(A)') ''
  write(*,'(A,F12.8)') '  Numerical  F_2(0) = ', results(1)
  write(*,'(A,F12.8)') '  Exact      F_2(0) = ', F2_exact
  write(*,'(A,F12.6,A)') '  Relative error    = ', &
    abs(results(1) - F2_exact) / F2_exact * 100.0_dp, ' %'
  write(*,'(A)') ''
  write(*,'(A)') '  F_2(q^2) increases for spacelike q^2 < 0'
  write(*,'(A)') '  (the electron appears "larger" at shorter distances)'

end program vertex_correction

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Chapter 7: Renormalization — Solution Details

Problem 7.1: Superficial Degree of Divergence

For a Feynman diagram in QED with $E_e$ external electron lines,$E_\gamma$ external photon lines, and $V$ vertices, the superficial degree of divergence is:

$$\boxed{D = 4 - \frac{3}{2}E_e - E_\gamma}$$

This remarkable result depends only on the external lines, not the number of loops! The superficially divergent diagrams in QED are:

Diagram	$E_e$	$E_\gamma$	$D$	Type
Vacuum polarization	0	2	2	Quadratic (but gauge: log)
Electron self-energy	2	0	1	Linear (but Dirac: log)
Vertex correction	2	1	0	Logarithmic

Ward identity ($Z_1 = Z_2$) ensures that the vertex and self-energy divergences are related, reducing QED renormalization to just two independent counterterms.

Problem 7.2: Running Coupling Constant

The QED coupling constant "runs" with energy scale due to vacuum polarization. The one-loop beta function gives:

$$\boxed{\alpha(q^2) = \frac{\alpha(0)}{1 - \frac{\alpha(0)}{3\pi}\ln\frac{q^2}{m_e^2}}}$$

The beta function for QED is:

$$\beta(\alpha) = \mu\frac{\partial\alpha}{\partial\mu} = \frac{2\alpha^2}{3\pi} + \mathcal{O}(\alpha^3)$$

Since $\beta > 0$ in QED, the coupling increases at short distances. At$q = M_Z \approx 91$ GeV, the effective coupling is$\alpha(M_Z) \approx 1/128$ rather than $\alpha(0) \approx 1/137$. This is in contrast to QCD where $\beta < 0$ (asymptotic freedom).

Python: Running QED Coupling Constant

Plots the running of $\alpha(q^2)$ from low energies to 100 GeV, including the effect of heavier lepton and quark thresholds.

Running QED Coupling alpha(q)

Python

Plots 1/alpha(q) and alpha(q) from 1 MeV to 100 GeV with fermion threshold effects

script.py109 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────
# Running QED coupling alpha(q) with fermion thresholds
# ──────────────────────────────────────────────

alpha_0 = 1.0 / 137.036
m_e = 0.000511   # GeV
m_mu = 0.10566   # GeV
m_tau = 1.777    # GeV
# Light quarks (constituent-like threshold masses)
m_u = 0.003      # GeV
m_d = 0.005      # GeV
m_s = 0.100      # GeV
m_c = 1.27       # GeV (charm)
m_b = 4.18       # GeV (bottom)

# Fermion charges squared and masses
fermions = [
    ('e',      1.0, m_e),
    ('mu',     1.0, m_mu),
    ('tau',    1.0, m_tau),
    ('u',      (2.0/3.0)**2 * 3, m_u),     # Nc * Q^2
    ('d',      (1.0/3.0)**2 * 3, m_d),
    ('s',      (1.0/3.0)**2 * 3, m_s),
    ('c',      (2.0/3.0)**2 * 3, m_c),
    ('b',      (1.0/3.0)**2 * 3, m_b),
]

def alpha_running(q, alpha0=alpha_0):
    """One-loop running alpha with fermion thresholds."""
    if q < 1e-6:
        return alpha0
    q2 = q**2
    delta = 0.0
    for name, Q2Nc, mf in fermions:
        if q > 2 * mf:  # above threshold
            ratio = q2 / mf**2
            if ratio > 1:
                delta += Q2Nc * np.log(ratio)
    correction = 1.0 - (alpha0 / (3.0 * np.pi)) * delta
    if correction <= 0.05:
        return alpha0 / 0.05  # cap to avoid Landau pole
    return alpha0 / correction

q_vals = np.logspace(-3, 2, 2000)  # 0.001 to 100 GeV
alpha_vals = np.array([alpha_running(q) for q in q_vals])
inv_alpha = 1.0 / alpha_vals

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# --- Left: 1/alpha(q) vs q ---
axes[0].plot(q_vals, inv_alpha, color='#fb923c', linewidth=2.5)
axes[0].axhline(y=137.036, color='#64748b', linestyle='--', alpha=0.5, linewidth=1)
axes[0].text(0.003, 137.5, '1/137.036', color='#94a3b8', fontsize=10)

# Mark thresholds
thresholds = [('$2m_e$', 2*m_e, '#f472b6'), ('$2m_\\mu$', 2*m_mu, '#c084fc'),
              ('$2m_c$', 2*m_c, '#60a5fa'), ('$2m_b$', 2*m_b, '#34d399'),
              ('$2m_\\\\tau$', 2*m_tau, '#fbbf24')]
for label, th, col in thresholds:
    if th > 0.001 and th < 100:
        axes[0].axvline(x=th, color=col, linestyle=':', alpha=0.5)
        axes[0].text(th*1.1, 128, label, color=col, fontsize=8, rotation=90)

axes[0].axhline(y=128.0, color='#fb923c', linestyle='--', alpha=0.3)
axes[0].text(0.003, 128.5, '1/128 (at $M_Z$)', color='#94a3b8', fontsize=9)

axes[0].set_xscale('log')
axes[0].set_xlabel('q [GeV]', color='#94a3b8', fontsize=13)
axes[0].set_ylabel('$1/\\alpha(q)$', color='#94a3b8', fontsize=13)
axes[0].set_title('Running QED Coupling: $1/\\alpha(q)$', color='white', fontsize=14)
axes[0].set_ylim(120, 140)
axes[0].tick_params(colors='#94a3b8')

# --- Right: alpha(q) vs q ---
axes[1].plot(q_vals, alpha_vals, color='#fb923c', linewidth=2.5)
axes[1].set_xscale('log')
axes[1].set_xlabel('q [GeV]', color='#94a3b8', fontsize=13)
axes[1].set_ylabel('$\\alpha(q)$', color='#94a3b8', fontsize=13)
axes[1].set_title('QED Fine Structure Constant', color='white', fontsize=14)
axes[1].tick_params(colors='#94a3b8')

# Annotate key values
alpha_mz = alpha_running(91.2)
axes[1].annotate(f'$\\alpha(M_Z)$ = 1/{1/alpha_mz:.1f}',
                 xy=(91.2, alpha_mz), xytext=(5, alpha_mz*1.3),
                 arrowprops=dict(arrowstyle='->', color='#fbbf24'),
                 color='#fbbf24', fontsize=11)

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight',
            facecolor='#0a0a1a', edgecolor='none')

print("Running QED coupling constant plotted")
print(f"  alpha(0)     = 1/{1/alpha_0:.3f}")
print(f"  alpha(1 GeV) = 1/{1/alpha_running(1.0):.1f}")
print(f"  alpha(M_Z)   = 1/{1/alpha_running(91.2):.1f}")
print(f"  alpha(100 GeV) = 1/{1/alpha_running(100):.1f}")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Chapter 9: Functional Methods — Solution Details

Problem 9.1: Generating Functional for the Free Scalar Field

The generating functional for the free scalar field with source $J(x)$ is:

$$Z_0[J] = \int \mathcal{D}\phi\, \exp\!\left[i\int d^4x\,\left(\frac{1}{2}(\partial\phi)^2 - \frac{1}{2}m^2\phi^2 + J\phi\right)\right]$$

Completing the square in $\phi$ and performing the Gaussian functional integral gives:

$$\boxed{Z_0[J] = Z_0[0]\,\exp\!\left[-\frac{1}{2}\int d^4x\,d^4y\, J(x)\,D_F(x-y)\,J(y)\right]}$$

The $n$-point correlation functions are obtained by functional differentiation:

$$\langle 0|T\{\phi(x_1)\cdots\phi(x_n)\}|0\rangle = \frac{1}{Z_0[0]}\left.\frac{\delta^n Z_0[J]}{i^n\,\delta J(x_1)\cdots\delta J(x_n)}\right|_{J=0}$$

For the free theory, only the 2-point function is nonzero (all higher connected correlators vanish), reproducing the Feynman propagator $D_F(x-y)$.

Problem 9.2: Path Integral for the Harmonic Oscillator

The quantum mechanical propagator for the harmonic oscillator can be computed exactly via the path integral:

$$K(x_f, T; x_i, 0) = \int_{x(0)=x_i}^{x(T)=x_f}\mathcal{D}x(t)\,\exp\!\left[\frac{i}{\hbar}\int_0^T dt\,\frac{m}{2}(\dot{x}^2 - \omega^2 x^2)\right]$$

The result is the famous Mehler kernel:

$$\boxed{K(x_f, T; x_i, 0) = \sqrt{\frac{m\omega}{2\pi i\hbar\sin\omega T}}\,\exp\!\left[\frac{im\omega}{2\hbar\sin\omega T}\bigl((x_i^2+x_f^2)\cos\omega T - 2x_i x_f\bigr)\right]}$$

The path integral is evaluated by separating $x(t) = x_{\text{cl}}(t) + \eta(t)$into the classical path and fluctuations. The classical action contributes the exponential, and the Gaussian integral over $\eta$ gives the prefactor.

Python: Monte Carlo Path Integral for Quantum Harmonic Oscillator

Simulates the quantum harmonic oscillator ground state wavefunction using Euclidean (imaginary time) lattice Monte Carlo. The path integral in Euclidean time gives access to the ground state probability distribution.

Monte Carlo Path Integral Simulation

Python

Euclidean lattice Monte Carlo for the quantum harmonic oscillator ground state

script.py150 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────
# Monte Carlo Path Integral: Quantum Harmonic Oscillator
# Euclidean lattice simulation
# ──────────────────────────────────────────────
np.random.seed(42)

# Parameters (natural units: hbar = m = omega = 1)
N_lat = 40          # number of time slices
a = 0.5             # lattice spacing
N_therm = 2000      # thermalization sweeps
N_meas = 8000       # measurement sweeps
N_skip = 5          # skip between measurements
delta = 1.2         # Metropolis step size

# Euclidean action for harmonic oscillator
# S_E = sum_i [ (x_{i+1} - x_i)^2 / (2a) + a * x_i^2 / 2 ]
def euclidean_action(x, a_lat):
    kinetic = np.sum((np.roll(x, -1) - x)**2) / (2.0 * a_lat)
    potential = a_lat * np.sum(x**2) / 2.0
    return kinetic + potential

# Initialize path
x = np.zeros(N_lat)

# Metropolis update
accept_count = 0
total_count = 0
measurements = []

print("=" * 50)
print("  Monte Carlo Path Integral Simulation")
print("  Quantum Harmonic Oscillator (m=omega=hbar=1)")
print(f"  Lattice: {N_lat} sites, spacing a = {a}")
print("=" * 50)

# Thermalization
for sweep in range(N_therm):
    for i in range(N_lat):
        x_old = x[i]
        S_old_local = (x[i] - x[(i-1) % N_lat])**2 / (2*a) + \
                      (x[(i+1) % N_lat] - x[i])**2 / (2*a) + \
                      a * x[i]**2 / 2.0
        x[i] += delta * (2.0 * np.random.random() - 1.0)
        S_new_local = (x[i] - x[(i-1) % N_lat])**2 / (2*a) + \
                      (x[(i+1) % N_lat] - x[i])**2 / (2*a) + \
                      a * x[i]**2 / 2.0
        dS = S_new_local - S_old_local
        if dS > 0 and np.random.random() > np.exp(-dS):
            x[i] = x_old

print(f"  Thermalization: {N_therm} sweeps done")

# Measurement phase
x_samples = []
E_samples = []

for sweep in range(N_meas * N_skip):
    for i in range(N_lat):
        x_old = x[i]
        S_old_local = (x[i] - x[(i-1) % N_lat])**2 / (2*a) + \
                      (x[(i+1) % N_lat] - x[i])**2 / (2*a) + \
                      a * x[i]**2 / 2.0
        x[i] += delta * (2.0 * np.random.random() - 1.0)
        S_new_local = (x[i] - x[(i-1) % N_lat])**2 / (2*a) + \
                      (x[(i+1) % N_lat] - x[i])**2 / (2*a) + \
                      a * x[i]**2 / 2.0
        dS = S_new_local - S_old_local
        if dS > 0 and np.random.random() > np.exp(-dS):
            x[i] = x_old
        else:
            accept_count += 1
        total_count += 1

if sweep % N_skip == 0:
        x_samples.extend(x.tolist())
        # Measure <x^2> for energy
        E_samples.append(np.mean(x**2))

x_all = np.array(x_samples)
accept_rate = accept_count / total_count

print(f"  Measurements: {N_meas} configurations")
print(f"  Acceptance rate: {accept_rate:.3f}")
print(f"  Total position samples: {len(x_all)}")

# Ground state energy from <x^2>
x2_mean = np.mean(E_samples)
x2_err = np.std(E_samples) / np.sqrt(len(E_samples))
# For SHO: <x^2>_0 = hbar/(2*m*omega) = 0.5 in our units
# E_0 = omega/2 = 0.5
print(f"\n  <x^2> = {x2_mean:.4f} +/- {x2_err:.4f}  (exact: 0.5000)")
print(f"  E_0 estimate = {0.5/x2_mean * 0.5:.4f}  (exact: 0.5000)")

# --- Plots ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
fig.patch.set_facecolor('#0a0a1a')
for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# Left: Histogram vs exact ground state
x_plot = np.linspace(-4, 4, 500)
psi0_sq = np.exp(-x_plot**2) / np.sqrt(np.pi)  # |psi_0(x)|^2

axes[0].hist(x_all, bins=80, range=(-4, 4), density=True,
             color='#22d3ee', alpha=0.6, label='MC samples')
axes[0].plot(x_plot, psi0_sq, color='#fbbf24', linewidth=2.5,
             label='$|\\psi_0(x)|^2$ exact')
axes[0].set_xlabel('x', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('Probability density', color='#94a3b8', fontsize=12)
axes[0].set_title('Ground State Wavefunction', color='white', fontsize=13)
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')

# Middle: Sample path
t_lattice = np.arange(N_lat) * a
axes[1].plot(t_lattice, x, color='#a78bfa', linewidth=1.5, marker='o', markersize=3)
axes[1].axhline(y=0, color='#64748b', linewidth=0.5)
axes[1].set_xlabel('Euclidean time $\\\\tau$', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('x($\\\\tau$)', color='#94a3b8', fontsize=12)
axes[1].set_title('Sample Lattice Path', color='white', fontsize=13)
axes[1].tick_params(colors='#94a3b8')

# Right: Autocorrelation of x^2
E_arr = np.array(E_samples)
E_mean = np.mean(E_arr)
E_var = np.var(E_arr)
max_lag = min(200, len(E_arr) // 2)
autocorr = np.array([np.mean((E_arr[:len(E_arr)-t] - E_mean) *
                              (E_arr[t:] - E_mean)) / E_var
                      for t in range(max_lag)])

axes[2].plot(range(max_lag), autocorr, color='#34d399', linewidth=1.5)
axes[2].axhline(y=0, color='#64748b', linewidth=0.5)
axes[2].set_xlabel('Lag (sweeps)', color='#94a3b8', fontsize=12)
axes[2].set_ylabel('Autocorrelation', color='#94a3b8', fontsize=12)
axes[2].set_title('$\\langle x^2 \\rangle$ Autocorrelation', color='white', fontsize=13)
axes[2].tick_params(colors='#94a3b8')

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight',
            facecolor='#0a0a1a', edgecolor='none')
print("\nPlot saved: Monte Carlo path integral results")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Related Resources

QFT Course

Complete quantum field theory course with 8 parts covering classical fields to the Standard Model.

Griffiths QM Solutions

87 video solutions to problems from Griffiths' Introduction to Quantum Mechanics.

Peskin & Schroeder Solutions

About the Textbook

Chapter 2: The Klein-Gordon Field

Chapter 3: The Dirac Field

Chapter 4: Interacting Fields and the S-Matrix

Chapter 5: Elementary Processes of QED

Chapter 6: Radiative Corrections

Chapter 7: Renormalization

Chapter 9: Functional Methods

Chapter 3: The Dirac Field — Solution Details

Problem 3.1: Lorentz Group Representations

Problem 3.2: Gordon Identity

Problem 3.3: Dirac Bilinear Transformation Properties

Python: Dirac Gamma Matrix Algebra Verification

Dirac Gamma Matrix Algebra

Chapter 4: Interacting Fields — Solution Details

Problem 4.1: Wick's Theorem

Problem 4.2: Feynman Propagator

Problem 4.3: $\phi^4$ Scattering Amplitude

Fortran: Numerical One-Loop Feynman Parameter Integral

One-Loop Feynman Parameter Integral

Chapter 5: Elementary QED Processes — Solution Details

Problem 5.1: $e^+e^- \to \mu^+\mu^-$ Cross Section

Problem 5.2: Compton Scattering (Klein-Nishina)

Problem 5.3: Bhabha Scattering

Python: Differential Cross Section for $e^+e^- \to \mu^+\mu^-$

e+e- to mu+mu- Cross Section

Chapter 6: Radiative Corrections — Solution Details

Problem 6.1: Vertex Correction and Anomalous Magnetic Moment

Problem 6.2: Electron Self-Energy

Key Results: g-2 and Lamb Shift

Fortran: Vertex Correction Feynman Parameter Integration

QED Vertex Correction Form Factor

Chapter 7: Renormalization — Solution Details

Problem 7.1: Superficial Degree of Divergence

Problem 7.2: Running Coupling Constant

Python: Running QED Coupling Constant

Running QED Coupling alpha(q)

Chapter 9: Functional Methods — Solution Details

Problem 9.1: Generating Functional for the Free Scalar Field

Problem 9.2: Path Integral for the Harmonic Oscillator

Python: Monte Carlo Path Integral for Quantum Harmonic Oscillator

Monte Carlo Path Integral Simulation

Recommended QFT Lecture Videos

David Tong: QFT Lectures (Cambridge)

Tobias Osborne: Advanced QFT (Hannover)

Perimeter Institute: QFT for Cosmologists

Leonard Susskind: Particle Physics (Stanford)

Related Resources

QFT Course

Griffiths QM Solutions