Scattering Amplitudes in Quantum Field Theory

Let's implememnt the calculations from Scattering Amplitudes in Quantum Field Theory in Python.

from sympy import Eq, Symbol, sympify, Matrix, sqrt, zeros, Basic, ones, pprint
from sympy.physics import msigma, mgamma

Now let's define the symbols for the energy and mass.

E = Symbol("E", real=True)
m = Symbol("m", real=True)

The Pauli matrices are defined as follows:

sigma1 = msigma(1)
sigma2 = msigma(2)
sigma3 = msigma(3)

$$
\sigma^{\mu} = \begin{pmatrix}
0 & \sigma^{\mu} \\
\sigma^{\mu} & 0
\end{pmatrix}
$$

The gamma matrices are defined as follows:

gamma0 = mgamma(0)
gamma1 = mgamma(1)
gamma2 = mgamma(2)
gamma3 = mgamma(3)
gamma5 = mgamma(5)

These are defined as:

$$
\gamma^{\mu} = \begin{pmatrix}
0 & \sigma^{\mu} \\
\sigma^{\mu} & 0
\end{pmatrix}
$$

The u-spinor is a solution to the Dirac equation for particles with positive energy.

$$
u(p, r) = \begin{pmatrix}
\chi \\
\frac{\vec{\sigma} \cdot \vec{p}}{E + m} \chi
\end{pmatrix}
$$

def u(p, r):
    """
    Calculate the u-spinor for a particle with momentum p and spin r.

    Args:
    p (tuple): Momentum components (p1, p2, p3)
    r (int): Spin state (1 or 2)

    Returns:
    Matrix: The u-spinor
    """
    if r not in [1, 2]:
        raise ValueError("Value of r should lie between 1 and 2")
    p1, p2, p3 = p
    if r == 1:
        ksi = Matrix([[1], [0]])
    else:
        ksi = Matrix([[0], [1]])
    a = (sigma1 * p1 + sigma2 * p2 + sigma3 * p3) / (E + m) * ksi
    if a == 0:
        a = zeros(2, 1)
    return sqrt(E + m) * Matrix([[ksi[0, 0]], [ksi[1, 0]], [a[0, 0]], [a[1, 0]]])

The v-spinor is a solution to the Dirac equation for antiparticles.

$$
v(p, r) = \begin{pmatrix}
\frac{\vec{\sigma} \cdot \vec{p}}{E + m} \chi \\
\chi
\end{pmatrix}
$$

def v(p, r):
    """
    Calculate the v-spinor for an antiparticle with momentum p and spin r.

    Args:
    p (tuple): Momentum components (p1, p2, p3)
    r (int): Spin state (1 or 2)

    Returns:
    Matrix: The v-spinor
    """
    if r not in [1, 2]:
        raise ValueError("Value of r should lie between 1 and 2")
    p1, p2, p3 = p
    if r == 1:
        ksi = Matrix([[1], [0]])
    else:
        ksi = -Matrix([[0], [1]])
    a = (sigma1 * p1 + sigma2 * p2 + sigma3 * p3) / (E + m) * ksi
    if a == 0:
        a = zeros(2, 1)
    return sqrt(E + m) * Matrix([[a[0, 0]], [a[1, 0]], [ksi[0, 0]], [ksi[1, 0]]])

The common notation in QFT is the slash notation of a four-momentum. This is defined as follows:

$$
p\!\!\!/ = \gamma^{\mu} p_{\mu} = \gamma^{\mu} \left( p_0, \vec{p} \right)
$$

def pslash(p):
    """
    Calculate the slash notation of a four-momentum.

    Args:
    p (tuple): Spatial momentum components (p1, p2, p3)

    Returns:
    Matrix: The slash notation of the four-momentum
    """
    p1, p2, p3 = p
    p0 = sqrt(m**2 + p1**2 + p2**2 + p3**2)
    return gamma0 * p0 - gamma1 * p1 - gamma2 * p2 - gamma3 * p3

Then we have the trace of a matrix defined as follows:

def Tr(M):
    return M.trace()

We can define several useful functions to display equations and calculate values.

def eq(lhs, rhs):
    """
    Pretty print an equation.

    Args:
    lhs (str): Left-hand side of the equation
    rhs: Right-hand side of the equation
    """
    pprint(Eq(sympify(lhs), rhs))

Now can get into the scattering amplitude. This process represents electron-positron annihilation followed by pair production.

$$
e^{-} + e^{+} \rightarrow \gamma \rightarrow e^{-} + e^{+}
$$

a = Symbol("a", real=True)
b = Symbol("b", real=True)
c = Symbol("c", real=True)

p = (a, b, c)

# Verify orthogonality of spinors
assert u(p, 1).D * u(p, 2) == Matrix(1, 1, [0])
assert u(p, 2).D * u(p, 1) == Matrix(1, 1, [0])

Now we define the momentum components.

# Define momentum components
p1, p2, p3 = [Symbol(x, real=True) for x in ["p1", "p2", "p3"]]
pp1, pp2, pp3 = [Symbol(x, real=True) for x in ["pp1", "pp2", "pp3"]]
k1, k2, k3 = [Symbol(x, real=True) for x in ["k1", "k2", "k3"]]
kp1, kp2, kp3 = [Symbol(x, real=True) for x in ["kp1", "kp2", "kp3"]]

p = (p1, p2, p3)
pp = (pp1, pp2, pp3)
k = (k1, k2, k3)
kp = (kp1, kp2, kp3)

mu = Symbol("mu")

For the intermediate steps, we need to calculate the following expressions:

$$
e = (p\!\!\!/ + m) (k\!\!\!/ - m)
$$

$$
f = p\!\!\!/ + m
$$

$$
g = p\!\!\!/ - m
$$

e = (pslash(p) + m * ones(4)) * (pslash(k) - m * ones(4))
f = pslash(p) + m * ones(4)
g = pslash(p) - m * ones(4)

eq("Tr(f*g)", Tr(f * g))

Now we can calculate the scattering amplitude.

$$
M = \sum_{\mu=0}^{3} M_{\mu}
$$

M0 = [
    (v(pp, 1).D * mgamma(mu) * u(p, 1)) * (u(k, 1).D * mgamma(mu, True) * v(kp, 1))
    for mu in range(4)
]
M = M0[0] + M0[1] + M0[2] + M0[3]
M = M[0]
if not isinstance(M, Basic):
    raise TypeError("Invalid type of variable")

d = Symbol("d", real=True)  # d=E+m

Finally, this calculation represents the matrix element of the scattering process.

eq("M", M)
print("-" * 40)
M = ((M.subs(E, d - m)).expand() * d**2).expand()

print(latex(M))

$$
M = \left(k_{3} \sqrt{E + m} \overline{\frac{1}{\sqrt{E + m}}} +
\frac{kp_{3} \overline{\sqrt{E + m}}}{\sqrt{E + m}}\right) \left(\frac{p_{3} \overline{\sqrt{E + m}}}{\sqrt{E + m}} + \ pp_{3} \sqrt{E + m} \overline{\frac{1}{\sqrt{E + m}}}\right) + ...
$$