Lecture 17 – Collision theory#

Hide code cell content
import warnings
from collections.abc import Callable
from typing import Any

import numpy as np
import plotly.graph_objects as go
import sympy as sp
from ampform.io import aslatex
from ampform.sympy import unevaluated
from IPython.display import Math
from plotly.colors import DEFAULT_PLOTLY_COLORS
from plotly.subplots import make_subplots

warnings.filterwarnings("ignore")
/home/runner/work/strong2020-salamanca/strong2020-salamanca/docs/lecture17/.venv/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

This notebook is an attempt to recreate the Mathematica notebook provided by Miguel Albaladejo. Another nice tutorial about the complex plane is this Julia notebook by Mikhail Mikhasenko.

Riemann sheets#

Square root example#

There are multiple solutions for \(x\) to the equation \(y^2 = x\) – the fact that we usually take \(y = \sqrt{x}\) to be the solution to this equation is just a matter of convention. It would be more complete to represent the solution as a set of points in the complex plane. In this case, we have the set \(S = \left\{\left(z, w\right)\in\mathbb{C}^2 | w^2=z\right\}\). This is set forms a Riemann surface in \(\mathbb{C}^2\) space.

Hide code cell content
def plot_riemann_surfaces(
    funcs: list[Callable],
    func_unicode: str,
    boundaries: tuple[complex, float] | tuple[complex, complex] = (0, 1),
    resolution: int | tuple[int, int] = 50,
    colorize: bool = True,
    mask: Callable[[np.ndarray, np.ndarray], bool] | None = None,
) -> None:
    X, Y = create_meshgrid(boundaries, resolution)
    Z = X + Y * 1j
    T = [f(Z) for f in funcs]
    if mask is not None:
        the_mask = np.full(Z.shape, fill_value=False)
        for t in T:
            the_mask |= mask(Z, t)
        if np.all(the_mask):
            msg = "All points were masked away"
            raise ValueError(msg)
        X[the_mask] = np.nan
        Y[the_mask] = np.nan
        Z[the_mask] = np.nan
        for t in T:
            t[the_mask] = np.nan

    vmax = max(max(t.imag.max(), t.real.max()) for t in T)
    style = lambda i, t: dict(
        cmin=-vmax,
        cmax=+vmax,
        showscale=colorize,
        colorscale=(
            "RdBu_r"
            if colorize
            else [[0, "rgb(0, 0, 0)"], [1, DEFAULT_PLOTLY_COLORS[i - 1]]]
        ),
        surfacecolor=t.real if colorize else np.ones(shape=t.shape),
    )
    S_im = [
        go.Surface(x=X, y=Y, z=t.imag, **style(i, t), name=f"Sheet {i}")
        for i, t in enumerate(T, 1)
    ]
    S_re = [
        go.Surface(x=X, y=Y, z=t.real, **style(i, t), name=f"Sheet {i}")
        for i, t in enumerate(T, 1)
    ]
    fig = make_subplots(
        cols=2,
        specs=[[{"type": "surface"}, {"type": "surface"}]],
        subplot_titles=(f"Im {func_unicode}", f"Re {func_unicode}"),
    )
    for i in range(len(funcs)):
        fig.add_trace(S_im[i], col=1, row=1)
        fig.add_trace(S_re[i], col=2, row=1)
    fig.update_layout(height=550, width=1_000)
    fig.update_traces(colorbar=dict(title="Re/Im"))
    fig.show()


def create_meshgrid(
    boundaries: tuple[complex, float] | tuple[complex, complex] = (0, 1),
    resolution: int | tuple[int, int] = 50,
) -> tuple[np.ndarray, np.ndarray]:
    if isinstance(resolution, tuple):
        x_res, y_res = resolution
    else:
        x_res, y_res = resolution, resolution
    box_min, box_max = boundaries
    if isinstance(box_max, float | int):
        pos, r_max = box_min, box_max
        R, Θ = np.meshgrid(
            np.linspace(0, r_max, num=x_res),
            np.linspace(-np.pi, +np.pi, num=y_res),
        )
        X = R * np.cos(Θ) + pos
        Y = R * np.sin(Θ) + pos
        return X, Y
    x1 = complex(box_min).real
    x2 = complex(box_max).real
    y1 = complex(box_min).imag
    y2 = complex(box_max).imag
    return np.meshgrid(
        np.linspace(x1, x2, num=x_res),
        np.linspace(y1, y2, num=y_res),
    )


def cut_t(
    cutoff: float | tuple[float, float],
) -> Callable[[np.ndarray, np.ndarray], bool]:
    if isinstance(cutoff, tuple):
        re_cut, im_cut = cutoff
    else:
        re_cut, im_cut = cutoff, cutoff
    return lambda _, t: (np.abs(t.real) > re_cut) | (np.abs(t.imag) > im_cut)
Hide code cell source
plot_riemann_surfaces(
    funcs=[lambda z: -np.sqrt(z), lambda z: +np.sqrt(z)],
    func_unicode="±√z",
)
plot_riemann_surfaces(
    funcs=[
        lambda z: -1 / np.sqrt(z),
        lambda z: +1 / np.sqrt(z),
    ],
    func_unicode="1/±√z",
    mask=cut_t(10),
)

Note also that since \(y = e^{x + 2n \pi i}\) for \(\forall n \in \mathbb{Z}\), we have that \(x = \log(y) + 2n\pi i\):

Hide code cell source
plot_riemann_surfaces(
    funcs=[
        lambda z: np.log(z) - 2j * np.pi,
        lambda z: np.log(z) + 2j * np.pi,
        np.log,
    ],
    func_unicode="log z",
    boundaries=(0, np.e**2),
    mask=cut_t((np.e, np.nan)),
)
Video explainers

Definition of the G(s) functions#

Hide code cell source
@unevaluated(real=False)
class SignedSqrt(sp.Expr):
    z: Any
    _latex_repr_ = R"\sqrt[+]{{{z}}}"

    def evaluate(self) -> sp.Expr:
        z = self.z
        return sp.sqrt(abs(z)) * sp.exp(sp.I * PosArg(z) / 2)


@unevaluated
class PosArg(sp.Expr):
    z: Any
    _latex_repr_ = R"\arg^+\left({z}\right)"

    def evaluate(self) -> sp.Expr:
        z = self.z
        arg = sp.arg(z)
        return sp.Piecewise(
            (arg + 2 * sp.pi, sp.im(z) < 0),
            (arg, True),
        )


z = sp.Symbol("z", complex=True)
Math(aslatex({e: e.evaluate() for e in [SignedSqrt(z), PosArg(z)]}))
\[\begin{split}\displaystyle \begin{array}{rcl} \sqrt[+]{z} &=& e^{\frac{i \arg^+\left(z\right)}{2}} \sqrt{\left|{z}\right|} \\ \arg^+\left(z\right) &=& \begin{cases} \arg{\left(z \right)} + 2 \pi & \text{for}\: \operatorname{im}{\left(z\right)} < 0 \\\arg{\left(z \right)} & \text{otherwise} \end{cases} \\ \end{array}\end{split}\]
Hide code cell source
plot_riemann_surfaces(
    funcs=[sp.lambdify(z, SignedSqrt(z).doit())],
    func_unicode="⁺√z",
    mask=lambda z, _: (np.abs(z.imag) < 1e-5) & (z.real > 0),
    resolution=(30, 301),
)
Hide code cell output
Hide code cell source
@unevaluated(real=False)
class G(sp.Expr):
    s: Any
    m: Any
    g0: Any
    sign: int = +1

    def _latex_repr_(self, printer, *args) -> str:
        s = printer._print(self.args[0])
        sign = self.args[-1]
        number = "I" if sign < 0 else "II"
        return f"G_{{{number}}}({s})"

    def evaluate(self) -> sp.Expr:
        s, m, g0, sign = self.args
        sigma = Sigma(s, m)
        g = (g0 - sigma * sp.log((sigma - 1) / (sigma + 1))) / (16 * sp.pi**2)
        return sp.Piecewise(
            (g, sign < 0),
            (G(s, m, g0, sign=-1) + 2 * sp.I * sigma / (16 * sp.pi), True),
        )


@unevaluated(real=False)
class Sigma(sp.Expr):
    s: Any
    m: Any
    _latex_repr_ = R"\sigma\left({s}\right)"

    def evaluate(self) -> sp.Expr:
        s, m = self.args
        return SignedSqrt(1 - 4 * m**2 / s)


s, g0 = sp.symbols("s g0", complex=True)
m = sp.Symbol("m", real=True, nonnegative=True)
sigma = Sigma(s, m)
G1 = G(s, m, g0, sign=-1)
G2 = G(s, m, g0, sign=+1)
definitions = {e: e.doit(deep=False) for e in [G1, G2, sigma]}
Math(aslatex(definitions))
\[\begin{split}\displaystyle \begin{array}{rcl} G_{I}(s) &=& \frac{g_{0} - \log{\left(\frac{\sigma\left(s\right) - 1}{\sigma\left(s\right) + 1} \right)} \sigma\left(s\right)}{16 \pi^{2}} \\ G_{II}(s) &=& G_{I}(s) + \frac{i \sigma\left(s\right)}{8 \pi} \\ \sigma\left(s\right) &=& \sqrt[+]{- \frac{4 m^{2}}{s} + 1} \\ \end{array}\end{split}\]
Hide code cell source
substitutions = {
    m: 139,
    g0: 3.0,
}
Math(aslatex(substitutions))
\[\begin{split}\displaystyle \begin{array}{rcl} m &=& 139 \\ g_{0} &=& 3.0 \\ \end{array}\end{split}\]
Hide code cell content
G1_expr = G1.doit().xreplace(substitutions)
G2_expr = G2.doit().xreplace(substitutions)
assert G1_expr.free_symbols == {s}
assert G2_expr.free_symbols == {s}
G1_func = sp.lambdify(s, G1_expr)
G2_func = sp.lambdify(s, G2_expr)
Hide code cell source
plot_riemann_surfaces(
    funcs=[
        lambda z: G1_func(z**2),
        lambda z: G2_func(z**2),
    ],
    func_unicode="G(s)",
    boundaries=(240 - 40j, 320 + 40j),
    colorize=False,
    resolution=(50, 401),
    mask=lambda z, _: np.abs(z.imag) == 0,
)

T-matrix definition#

Hide code cell source
@unevaluated(real=False)
class S(sp.Expr):
    s: Any
    m: Any
    : Any
    GV: Any
    : Any
    g0: Any
    sign: int = +1

    def _latex_repr_(self, printer, *args) -> str:
        s = printer._print(self.args[0])
        sign = self.args[-1]
        number = "I" if sign < 0 else "II"
        return f"S_{{{number}}}({s})"

    def evaluate(self) -> sp.Expr:
        s, m, *_ = self.args
        return 1 - 2 * sp.I * Sigma(s, m) / (16 * sp.pi) * T(*self.args)


@unevaluated(real=False)
class T(sp.Expr):
    s: Any
    m: Any
    : Any
    GV: Any
    : Any
    g0: Any
    sign: int = +1

    def _latex_repr_(self, printer, *args) -> str:
        s = printer._print(self.args[0])
        sign = self.args[-1]
        number = "I" if sign < 0 else "II"
        return f"T_{{{number}}}({s})"

    def evaluate(self) -> sp.Expr:
        s, m, , GV, , g0, sign = self.args
        return 1 / (1 / V1(s, m, , GV, ) - G(s, m, g0, sign))


@unevaluated(real=False)
class V1(sp.Expr):
    s: Any
    m: Any
    : Any
    GV: Any
    : Any
    _latex_repr_ = R"V_1\left({s}\right)"

    def evaluate(self) -> sp.Expr:
        s, m, , GV,  = self.args
        return -(2 * p2(s, m)) / (3 * **2) * (
            1 - GV**2 / **2 * 2 * s / (s - **2)
        ) - GV**2 / **4 * p2(s, m) * h(**2 / (2 * p2(s, m)))


@unevaluated
class h(sp.Expr):
    a: Any
    _latex_repr_ = R"h\left({a}\right)"

    def evaluate(self) -> sp.Expr:
        a = self.args[0]
        return -sp.Mul(
            sp.Rational(2, 3),
            (1 + 6 * a + 3 * a**2),
            evaluate=False,
        ) + a * (2 + 3 * a + a**2) * sp.log(1 + 2 / a)


@unevaluated
class p2(sp.Expr):
    s: Any
    m: Any
    _latex_repr_ = R"p^2\left({s}\right)"

    def evaluate(self) -> sp.Expr:
        s, m = self.args
        return s / 4 - m**2


a, , GV,  = sp.symbols("a m_rho, G_V f_pi")
_exprs = [
    S(s, m, , GV, , g0, sign=-1),
    T(s, m, , GV, , g0, sign=-1),
    T(s, m, , GV, , g0, sign=+1),
    V1(s, m, , GV, ),
    h(a),
    p2(s, m),
]
Math(aslatex({e: e.doit(deep=False) for e in _exprs}))
\[\begin{split}\displaystyle \begin{array}{rcl} S_{I}(s) &=& - \frac{i \sigma\left(s\right) T_{I}(s)}{8 \pi} + 1 \\ T_{I}(s) &=& \frac{1}{- G_{I}(s) + \frac{1}{V_1\left(s\right)}} \\ T_{II}(s) &=& \frac{1}{- G_{II}(s) + \frac{1}{V_1\left(s\right)}} \\ V_1\left(s\right) &=& - \frac{G_{V}^{2} h\left(\frac{m_{\rho}^{2}}{2 p^2\left(s\right)}\right) p^2\left(s\right)}{f_{\pi}^{4}} - \frac{2 \left(- \frac{2 G_{V}^{2} s}{f_{\pi}^{2} \left(- m_{\rho}^{2} + s\right)} + 1\right) p^2\left(s\right)}{3 f_{\pi}^{2}} \\ h\left(a\right) &=& a \left(a^{2} + 3 a + 2\right) \log{\left(1 + \frac{2}{a} \right)} - \frac{2 \left(3 a^{2} + 6 a + 1\right)}{3} \\ p^2\left(s\right) &=& - m^{2} + \frac{s}{4} \\ \end{array}\end{split}\]
Hide code cell source
gv = sp.Symbol("g_v")
substitutions = {
    : 87.3,
    GV: sp.sqrt(gv**2 * **2) / 2,
    gv: 1,
    m: 139,
    : 770,
    g0: -3,
}
Math(aslatex(substitutions))
\[\begin{split}\displaystyle \begin{array}{rcl} f_{\pi} &=& 87.3 \\ G_{V} &=& \frac{\sqrt{f_{\pi}^{2} g_{v}^{2}}}{2} \\ g_{v} &=& 1 \\ m &=& 139 \\ m_{\rho} &=& 770 \\ g_{0} &=& -3 \\ \end{array}\end{split}\]
Hide code cell source
T_exprs = [
    T(s, m, , GV, , g0, sign).doit().xreplace(substitutions).xreplace(substitutions)
    for sign in [-1, +1]
]
T_funcs = [sp.lambdify(s, expr) for expr in T_exprs]
Hide code cell source
x = np.linspace(500, 1_100, num=200)
y = np.linspace(1e-5, 150, num=100)
X, Yn = np.meshgrid(x, -y)
X, Yp = np.meshgrid(x, +y)
Zn = X + Yn * 1j
Zp = X + Yp * 1j
Tn = T_funcs[1](Zn**2)
Tp = T_funcs[0](Zp**2)

vmax = 100
sty = lambda t: dict(
    cmin=-vmax,
    cmax=+vmax,
    colorscale="RdBu_r",
    surfacecolor=t.imag,
)
Sn = go.Surface(x=X, y=Yn, z=Tn.real, **sty(Tn), name="Unphysical")
Sp = go.Surface(x=X, y=Yp, z=Tp.real, **sty(Tp), name="Physical", colorbar_title="Re T")
y = Yp[0]
z = x + y * 1j
line = go.Scatter3d(
    x=x,
    y=y,
    z=T_funcs[0](z**2).real,
    marker=dict(size=0),
    line=dict(color="darkgreen", width=1),
)
fig = go.Figure(data=[Sn, Sp, line])
fig.update_layout(height=550, width=600)
fig.update_scenes(
    xaxis_title_text="Re s",
    yaxis_title_text="Im s",
    zaxis_title_text="Im T",
    zaxis_range=[-vmax, +vmax],
)
fig.show()
Hide code cell source
sty = lambda t: dict(
    cmin=-vmax,
    cmax=+vmax,
    colorscale="RdBu_r",
    surfacecolor=t.real,
)
Sn = go.Surface(x=X, y=Yn, z=Tn.imag, **sty(Tn), name="Unphysical")
Sp = go.Surface(x=X, y=Yp, z=Tp.imag, **sty(Tp), name="Physical", colorbar_title="Im T")
y = Yp[0]
z = x + y * 1j
line = go.Scatter3d(
    x=x,
    y=y,
    z=T_funcs[0](z**2).imag,
    marker=dict(size=0),
    line=dict(color="darkgreen", width=1),
)
fig = go.Figure(data=[Sn, Sp, line])
fig.update_layout(height=550, width=600)
fig.update_scenes(
    xaxis_title_text="Re s",
    yaxis_title_text="Im s",
    zaxis_title_text="Re T",
    zaxis_range=[-vmax, +vmax],
)
fig.show()