Lecture 17 – Collision theory#
Import Python libraries
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.
Code for plotting Riemann sheets with plotly
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)
Show 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\):
Show 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#
Define special square root for negative values
@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}\]
Show Riemann surface of ⁺√z
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),
)
Show Riemann surface of ⁺√z
Define G and σ expression classes
@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}\]
Show 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}\]
Convert expressions to numerical functions
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)
Show 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#
Define expression classes
@unevaluated(real=False)
class S(sp.Expr):
s: Any
m: Any
Mρ: Any
GV: Any
fπ: 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
Mρ: Any
GV: Any
fπ: 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, Mρ, GV, fπ, g0, sign = self.args
return 1 / (1 / V1(s, m, Mρ, GV, fπ) - G(s, m, g0, sign))
@unevaluated(real=False)
class V1(sp.Expr):
s: Any
m: Any
Mρ: Any
GV: Any
fπ: Any
_latex_repr_ = R"V_1\left({s}\right)"
def evaluate(self) -> sp.Expr:
s, m, Mρ, GV, fπ = self.args
return -(2 * p2(s, m)) / (3 * fπ**2) * (
1 - GV**2 / fπ**2 * 2 * s / (s - Mρ**2)
) - GV**2 / fπ**4 * p2(s, m) * h(Mρ**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, Mρ, GV, fπ = sp.symbols("a m_rho, G_V f_pi")
_exprs = [
S(s, m, Mρ, GV, fπ, g0, sign=-1),
T(s, m, Mρ, GV, fπ, g0, sign=-1),
T(s, m, Mρ, GV, fπ, g0, sign=+1),
V1(s, m, Mρ, GV, fπ),
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}\]
Show code cell source
gv = sp.Symbol("g_v")
substitutions = {
fπ: 87.3,
GV: sp.sqrt(gv**2 * fπ**2) / 2,
gv: 1,
m: 139,
Mρ: 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}\]
Show code cell source
T_exprs = [
T(s, m, Mρ, GV, fπ, g0, sign).doit().xreplace(substitutions).xreplace(substitutions)
for sign in [-1, +1]
]
T_funcs = [sp.lambdify(s, expr) for expr in T_exprs]
Show 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()
Show 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()