Tree-Level QED Total Cross-Section: Numerical Evaluation and Energy Dependence

1. From Amplitude to Cross-Section: A Quick Review

Let's first recall the relationship between the differential cross-section $\frac{d\sigma }{d\Omega }$ and the squared scattering amplitude $⟨|\mathcal{M}{|}^2⟩$. In the center-of-mass (CM) frame, for a process $1+2\to 3+4$, the formula is:

where $s=E_{\text{cm}}^2$ is the Mandelstam variable, $|{\vec{p}}_i|$ is the magnitude of the momentum of either incoming particle, and $|{\vec{p}}_f|$ is the magnitude of the momentum of either outgoing particle. The brackets $⟨\dots ⟩$ indicate an average over initial spins and sum over final spins.

The total cross-section, ${\sigma }_{\text{tot}}$, is then found by integrating the differential cross-section over the full solid angle:

To see how these pieces come together, let's watch a video that explains the concept of the differential cross-section and how it can be used. While the video uses Compton scattering as its example and employs symbolic integration, the core physical concepts are directly applicable to our task.

2. The Cross-Section for $e^{+}{e}^{-}\to {\mu }^{+}{\mu }^{-}$

Now, let's specialize to our process of interest. We previously calculated the squared amplitude for $e^{+}{e}^{-}\to {\mu }^{+}{\mu }^{-}$ using trace techniques. The final result for the differential cross-section in the CM frame, neglecting the electron mass ($m_e\ll E$) but keeping the muon mass ($M_{\mu }$), is a classic QED result.

You can find the detailed derivation in the resource below. For our purposes, we will take the final formulas for the differential and total cross-sections as our starting point.

From the resource, we have our key formulas. In the CM frame, let $E$ be the energy of the incoming electron (so $s=(2E{)}^2$).

The differential cross-section is:

The analytical total cross-section (after integrating over $\theta$ and $\varphi$) is:

where $E=\sqrt{s}/2$.

Our task is to integrate the first equation numerically and verify that it matches the second equation.

3. Numerical Integration in Python

To perform the numerical integration, we will use the quad function from the scipy.integrate library. It is a robust and standard tool for 1D definite integrals in Python.

The following short video explains the basics of using quad, including how to pass additional parameters to your function—a feature we will need to pass the energy $E$ to our integrand.

Now, let's set up the calculation. The integrand for the total cross-section is the differential cross-section integrated over $\varphi$, which gives a factor of $2\pi$ since there is no $\varphi$ dependence. The remaining integral is over $\theta$:

So, our integrand function in Python will represent the term in the parentheses: $f(\theta ,E)=2\pi \sin \theta \frac{d\sigma }{d\Omega }(\theta ,E)$.

Your Task: Implement the Calculation

Below is a Python code template. Your task is to fill in the missing parts to perform the numerical integration and plot the results.

Physical Constants:

• Fine-structure constant, $\alpha \approx 1/137.036$
• Muon mass, $M_{\mu }\approx 105.66$ MeV/c²
• We'll use natural units where $\hslash =c=1$. A useful conversion factor for cross-sections is $(\hslash c{)}^2\approx 0.389379{\text{GeV}}^2\cdot \text{mb}$ (millibarns).

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

# --- Physical Constants ---
# Using MeV for energy
ALPHA = 1/137.036
MUON_MASS_MEV = 105.66
# (hbar*c)^2 in MeV^2 * mb
HBARC2_MB = 0.389379 * 1e3 

# 1. Define the differential cross section integrand
# This function should calculate d(sigma)/d(theta) = 2*pi*sin(theta)*d(sigma)/d(Omega)
def integrand(theta, E_cm):
    """
    Calculates the integrand for the total cross section of e+e- -> mu+mu-.
    
    Args:
        theta (float): The scattering angle in radians.
        E_cm (float): The center-of-mass energy in MeV.
        
    Returns:
        float: The value of d(sigma)/d(theta) in MeV^-2.
    """
    s = E_cm**2
    E = E_cm / 2.0
    
    # Threshold condition: reaction is only possible if E > M_mu
    if E <= MUON_MASS_MEV:
        return 0.0
        
    # YOUR CODE HERE:
    # Implement the formula for d(sigma)/d(Omega) and combine it
    # with 2*pi*sin(theta) to get the integrand.
    # dsigma_dOmega = ...
    # return 2 * np.pi * np.sin(theta) * dsigma_dOmega
    
    # Placeholder - replace this
    return 0.0

# 2. Define a function to calculate the total cross-section numerically
def calculate_sigma_numeric(E_cm):
    """
    Calculates the total cross section by numerically integrating the integrand.
    
    Args:
        E_cm (float): Center-of-mass energy in MeV.
    
    Returns:
        float: Total cross section in mb.
    """
    # YOUR CODE HERE:
    # Use scipy.integrate.quad to integrate the 'integrand' function
    # from theta = 0 to theta = pi. Pass E_cm as an argument.
    # The result from quad is a tuple (value, error). We only need the value.
    # sigma_MeV_minus_2, _ = quad(...)
    
    # Convert from MeV^-2 to millibarns (mb)
    # sigma_mb = sigma_MeV_minus_2 * HBARC2_MB
    # return sigma_mb
    
    # Placeholder - replace this
    return 0.0

# 3. Define the analytical formula for comparison
def calculate_sigma_analytical(E_cm):
    """
    Calculates the total cross section using the analytical formula.
    
    Args:
        E_cm (float): Center-of-mass energy in MeV.
        
    Returns:
        float: Total cross section in mb.
    """
    s = E_cm**2
    E = E_cm / 2.0

    if E <= MUON_MASS_MEV:
        return 0.0
        
    # YOUR CODE HERE:
    # Implement the analytical formula for sigma_tot (eq. 69 from the text).
    # sigma_MeV_minus_2 = ...
    
    # Convert from MeV^-2 to mb
    # sigma_mb = sigma_MeV_minus_2 * HBARC2_MB
    # return sigma_mb
    
    # Placeholder - replace this
    return 0.0

# --- Main execution and plotting ---
if __name__ == "__main__":
    # Set up a range of center-of-mass energies to scan
    # Start just above the threshold (2 * M_mu)
    E_cm_min = 2 * MUON_MASS_MEV + 1 
    E_cm_max = 1000  # Scan up to 1 GeV
    energies = np.linspace(E_cm_min, E_cm_max, 200)

    # Calculate cross sections
    sigma_numeric_vals = [calculate_sigma_numeric(E) for E in energies]
    sigma_analytical_vals = [calculate_sigma_analytical(E) for E in energies]
    
    # High-energy approximation: sigma ~ 1/s
    sigma_high_E_approx = (4 * np.pi * ALPHA**2) / (3 * energies**2) * HBARC2_MB

    # Plotting
    plt.style.use('dark_background')
    fig, ax = plt.subplots(figsize=(10, 6))
    
    ax.plot(energies, sigma_analytical_vals, 'w-', label='Analytical Formula', linewidth=2)
    ax.plot(energies, sigma_numeric_vals, 'c--', label='Numerical Integration', linewidth=2)
    ax.plot(energies, sigma_high_E_approx, 'r:', label=r'High-Energy Approx ($\propto 1/s$)', linewidth=2)

    ax.set_xlabel('Center-of-Mass Energy $E_{cm}$ [MeV]', fontsize=12)
    ax.set_ylabel(r'Total Cross-Section $\sigma_{tot}$ [mb]', fontsize=12)
    ax.set_title(r'$e^+e^- \to \mu^+\mu^-$ Cross-Section vs. Energy', fontsize=14)
    ax.legend(fontsize=10)
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.grid(True, which='both', linestyle='--', alpha=0.5)
    
    plt.show()

# In integrand(theta, E_cm):
    term1 = (ALPHA**2) / (4 * s)
    term2 = np.sqrt(1 - (MUON_MASS_MEV**2 / E**2))
    angular_term = (1 + (MUON_MASS_MEV**2 / E**2)) + (1 - (MUON_MASS_MEV**2 / E**2)) * np.cos(theta)**2
    dsigma_dOmega = term1 * term2 * angular_term
    return 2 * np.pi * np.sin(theta) * dsigma_dOmega

# In calculate_sigma_numeric(E_cm):
    sigma_MeV_minus_2, _ = quad(integrand, 0, np.pi, args=(E_cm,))
    sigma_mb = sigma_MeV_minus_2 * HBARC2_MB
    return sigma_mb

# In calculate_sigma_analytical(E_cm):
    term1 = (4 * np.pi * ALPHA**2) / (3 * s)
    term2 = np.sqrt(1 - (MUON_MASS_MEV**2 / E**2))
    term3 = 1 + (MUON_MASS_MEV**2 / (2 * E**2))
    sigma_MeV_minus_2 = term1 * term2 * term3
    sigma_mb = sigma_MeV_minus_2 * HBARC2_MB
    return sigma_mb

When you run the completed script, you should see the plot of the total cross-section as a function of energy. The curves from your numerical integration and the analytical formula should lie perfectly on top of each other, validating your code. You will also see how the exact result converges to the simple $1/s$ behavior at high energies. This falling cross-section is a hallmark of point-like scattering mediated by a single virtual particle.

4. Alternative: Monte Carlo Integration

The quad function uses a sophisticated adaptive quadrature method, which is highly efficient for well-behaved 1D integrals. For more complex problems, especially with more final-state particles (leading to higher-dimensional integrals over the "phase space"), Monte Carlo integration becomes an indispensable tool.

The basic idea is simple: to compute ${\int }_a^{b}f(x)dx$, you sample $N$ random points $x_i$ uniformly in $[a,b]$ and approximate the integral as:

The error in this approximation typically decreases as $1/\sqrt{N}$. While less efficient than quadrature for simple 1D cases, its performance advantage becomes significant in higher dimensions.

We will explore Monte Carlo methods in more detail later in the course, especially when we discuss event generators and lattice field theory.

Conclusion

Congratulations on completing this module! You have progressed from the fundamental principles of QED to calculating and verifying a real-world observable.

Key Takeaways:

• The total cross-section ${\sigma }_{\text{tot}}$ is the integral of the differential cross-section $\frac{d\sigma }{d\Omega }$ over the full solid angle.
• We can use standard numerical libraries like SciPy to perform this integration, allowing us to compute physical observables directly from our theoretical formulas for $|\mathcal{M}{|}^2$.
• For the process $e^{+}{e}^{-}\to {\mu }^{+}{\mu }^{-}$, the total cross-section is zero below the threshold energy $E_{\text{cm}}=2M_{\mu }$, rises to a peak, and then falls proportionally to $1/s$ at high energies.
• Comparing numerical results with analytical formulas is a critical step for verifying implementations in computational physics.

Preview of the Next Module:
So far, all our calculations have been at "tree level"—the lowest non-trivial order in the fine-structure constant $\alpha$. But what about higher-order corrections? These correspond to Feynman diagrams with closed loops. In the next module, "Radiative Corrections and Renormalization in QED", we will venture into this fascinating and challenging territory. Our first step will be to identify the ultraviolet (UV) divergences that appear in one-loop QED diagrams, which are infinities that arise from integrating over arbitrarily high momenta inside the loops. We will then learn the powerful techniques of regularization and renormalization used to tame these infinities and make extraordinarily precise physical predictions.