TR-026

Visualization of the Riemann sheets for the single-channel \(T\) matrix with one resonance pole

This report investigates and reproduces the Riemann sheets shown in Fig. 50.1 and 50.2 of the PDG. The lineshape parametrization is directly derived with the \(K\)-matrix formalism. The transition from the first physical sheet to the second unphysical sheet is derived using analytic continuation.

🚧 ampform#67

Single-channel Riemann sheets

The \(T\) function can be extended into the complex plane. This results in \(2^n\) Riemann sheets for \(n\) channels, each starting at the threshold \(s_{thr}=(m_1+m_2)^{2}\) of the two final state particles, the so-called branching point of the respective channel going along the so-called branch cut along the real axis where the function is not uniquely defined to \(+\infty\). This choice of the direction of the brach cut is most commonly used in particle physics. The physical Riemann sheet is defined for positive imaginary part (1^st quadrant of the complex plane) and the unphysical Riemann sheets are only defined for negative imaginary part (4^th quadrant of the complex plane). For the single-channel case there are two Riemann sheets, one physical and one unphysical.

Import Python libraries Hide code cell content

# @title
from __future__ import annotations

import warnings
from typing import Any

import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go
import sympy as sp
from ampform.io import aslatex
from ampform.kinematics.phasespace import Kallen
from ampform.sympy import unevaluated
from IPython.display import Math

warnings.filterwarnings("ignore")

Phase space factor definitions

Show code cell source Hide code cell source

# @title
@unevaluated(real=False)
class PhaseSpaceFactor(sp.Expr):
    s: Any
    m1: Any
    m2: Any
    _latex_repr_ = R"\rho_{{{m1}, {m2}}}\left({s}\right)"

    def evaluate(self) -> sp.Expr:
        s, m1, m2 = self.args
        return sp.sqrt((s - ((m1 + m2) ** 2)) * (s - (m1 - m2) ** 2) / s**2)


@unevaluated(real=False)
class PhaseSpaceCM(sp.Expr):
    s: Any
    m1: Any
    m2: Any
    _latex_repr_ = R"\rho^\mathrm{{CM}}_{{{m1},{m2}}}\left({s}\right)"

    def evaluate(self) -> sp.Expr:
        s, m1, m2 = self.args
        return -16 * sp.pi * sp.I * ChewMandelstam(s, m1, m2)


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

    def evaluate(self) -> sp.Expr:
        s, m1, m2 = self.args
        q = BreakupMomentum(s, m1, m2)
        return (
            1
            / (16 * sp.pi**2)
            * (
                (2 * q / sp.sqrt(s))
                * sp.log((m1**2 + m2**2 - s + 2 * sp.sqrt(s) * q) / (2 * m1 * m2))
                - (m1**2 - m2**2) * (1 / s - 1 / (m1 + m2) ** 2) * sp.log(m1 / m2)
            )
        )


@unevaluated(real=False)
class BreakupMomentum(sp.Expr):
    s: Any
    m1: Any
    m2: Any
    _latex_repr_ = R"q\left({s}\right)"

    def evaluate(self) -> sp.Expr:
        s, m1, m2 = self.args
        return sp.sqrt(Kallen(s, m1**2, m2**2)) / (2 * sp.sqrt(s))


s, m1, m2 = sp.symbols("s m1 m2")
rho_expr = PhaseSpaceFactor(s, m1, m2)
rho_cm_expr = PhaseSpaceCM(s, m1, m2)
cm_expr = ChewMandelstam(s, m1, m2)
q_expr = BreakupMomentum(s, m1, m2)
kallen = Kallen(*sp.symbols("x:z"))
src = aslatex({
    e: e.doit(deep=False) for e in [rho_expr, rho_cm_expr, cm_expr, q_expr, kallen]
})
Math(src)

\[\begin{split}\displaystyle \begin{array}{rcl} \rho_{m_{1}, m_{2}}\left(s\right) &=& \sqrt{\frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{s^{2}}} \\ \rho^\mathrm{CM}_{m_{1},m_{2}}\left(s\right) &=& - 16 i \pi \Sigma\left(s\right) \\ \Sigma\left(s\right) &=& \frac{- \left(m_{1}^{2} - m_{2}^{2}\right) \left(- \frac{1}{\left(m_{1} + m_{2}\right)^{2}} + \frac{1}{s}\right) \log{\left(\frac{m_{1}}{m_{2}} \right)} + \frac{2 \log{\left(\frac{m_{1}^{2} + m_{2}^{2} + 2 \sqrt{s} q\left(s\right) - s}{2 m_{1} m_{2}} \right)} q\left(s\right)}{\sqrt{s}}}{16 \pi^{2}} \\ q\left(s\right) &=& \frac{\sqrt{\lambda\left(s, m_{1}^{2}, m_{2}^{2}\right)}}{2 \sqrt{s}} \\ \lambda\left(x, y, z\right) &=& x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2} \\ \end{array}\end{split}\]

T matrix definition with K matrix

The dynamical part of the scattering amplitude is calculated via \(K\) matrix formalism. In this report the single-channel case with one resonance pole is assumed.

# @title
n = 1
I = sp.Identity(n)
K = sp.MatrixSymbol("K", n, n)
CM = sp.MatrixSymbol(R"{\rho_{cm}}", n, n)

T1 = (I + sp.I * K * CM).inv() * K
T1

\[\displaystyle \left(\mathbb{I} + i K {\rho_{cm}}\right)^{-1} K\]

T1_explicit = T1.as_explicit()
T1_explicit[0, 0]

\[\displaystyle \frac{{K}_{0,0}}{i {K}_{0,0} {{\rho_{cm}}}_{0,0} + 1}\]

Show code cell source Hide code cell source

# @title
g0, m0 = sp.symbols(R"g^{0} m0")
k_expr = (g0**2) / (s - m0**2)
definitions_I = {
    K[0, 0]: k_expr,
    CM[0, 0]: PhaseSpaceCM(s, m1, m2),
}
Math(aslatex(definitions_I))

\[\begin{split}\displaystyle \begin{array}{rcl} {K}_{0,0} &=& \frac{\left(g^{0}\right)^{2}}{- m_{0}^{2} + s} \\ {{\rho_{cm}}}_{0,0} &=& \rho^\mathrm{CM}_{m_{1},m_{2}}\left(s\right) \\ \end{array}\end{split}\]

T1_expr = T1_explicit[0, 0].xreplace(definitions_I)

T1_expr.simplify(doit=False)

\[\displaystyle \frac{\left(g^{0}\right)^{2}}{i \left(g^{0}\right)^{2} \rho^\mathrm{CM}_{m_{1},m_{2}}\left(s\right) - m_{0}^{2} + s}\]

Calculation of the second Riemann sheet

Since the \(T\) function is real below the branch cut it can be shown that the discontinuity above and below the threshold reads as:

\[ CM(s+i\epsilon)-CM(s-i\epsilon)= i\rho -(-i\rho) =2i\rho \]

when \(\epsilon\) goes to zero.
Which leads to:

\[ CM^{-1}_{\mathrm{II}}(s-i\epsilon)= Re(CM^{-1}_{\mathrm{I}}(s-i\epsilon))-i\rho+2i\rho \]

For the the Amplitude for the second sheet is defined as:

\[ A^{-1}_{\mathrm{II}}(s)= A^{-1}_{\mathrm{I}}(s)-2i\rho \]

rho = sp.MatrixSymbol("rho", n, n)
T2 = (T1.inv() + 2 * sp.I * rho).inv()
T2

\[\displaystyle \left(K^{-1} \left(\mathbb{I} + i K {\rho_{cm}}\right) + 2 i \rho\right)^{-1}\]

definitions_II = {
    **definitions_I,
    rho[0, 0]: PhaseSpaceFactor(s, m1, m2),
}
T2_explicit = T2.as_explicit()
T2_expr = T2_explicit[0, 0].xreplace(definitions_II)

T2_expr.simplify(doit=False)

\[\displaystyle \frac{\left(g^{0}\right)^{2}}{i \left(g^{0}\right)^{2} \rho^\mathrm{CM}_{m_{1},m_{2}}\left(s\right) + 2 i \left(g^{0}\right)^{2} \rho_{m_{1}, m_{2}}\left(s\right) - m_{0}^{2} + s}\]

Visualization of the 2 dimensional lineshape

Define numerical functions Hide code cell source

symbols = sp.Tuple(s, m1, m2, m0, g0)
T1_func = sp.lambdify(symbols, T1_expr.doit())
T2_func = sp.lambdify(symbols, T2_expr.doit())

Define meshgrid and parameter values Hide code cell source

# @title
epsilon = 1e-5
x = np.linspace(0, 6, num=200)
y = np.linspace(epsilon, 1, num=100)
X, Y = np.meshgrid(x, y)
Zn = X - Y * 1j
Zp = X + Y * 1j

values = {
    m1: 0.9,
    m2: 0.8,
    m0: 3.1,
    g0: 1.5,
}
args = eval(str(symbols[1:].xreplace(values)))

T1n = T1_func(Zn**2, *args)
T1p = T1_func(Zp**2, *args)

T2n = T2_func(Zn**2, *args)
T2p = T2_func(Zp**2, *args)

Show code cell source Hide code cell source

# @title
%config InlineBackend.figure_formats = ["svg"]

plt.rcParams.update({"font.size": 16})
fig, axes = plt.subplots(figsize=(15, 6), ncols=2, sharey=True)
ax1, ax2 = axes
for ax in axes:
    ax.set_xlabel(R"$\mathrm{Re}\,\sqrt{s}$")
ax1.set_ylabel(R"$\mathrm{Im}\,T$")

ax1.plot(x, T1n[0].imag, label=R"$T_\mathrm{I}(s-0i)$")
ax1.plot(x, T1p[0].imag, label=R"$T_\mathrm{I}(s+0i)$")
ax1.set_title(f"${sp.latex(rho_cm_expr)}$")
ax1.set_title(R"$T_\mathrm{I}$")

ax2.plot(x, T2n[0].imag, label=R"$T_\mathrm{II}(s-0i)$")
ax2.plot(x, T2p[0].imag, label=R"$T_\mathrm{II}(s+0i)$")
ax2.set_title(R"$T_\mathrm{II}$")

for ax in axes:
    ax.legend()

fig.tight_layout()
plt.show()

The Amplitude for the second sheet is only defined for \(s\) positive real part and negative complex part. It inherits the analytic structure of the phasespace factor \(\rho\) (the branch cut starting form zero and from \(s=s_{thr}\) on the real axis). So it is only defined up to the closest branch cut which is in this case the cut at \(s=s_{thr}\).

Visualization of the Riemann sheets

Show code cell source Hide code cell source

# @title
def sty(sheet_name: str) -> dict:
    sheet_color = sheet_colors[sheet_name]
    n_lines = 12
    return dict(
        cmin=-vmax,
        cmax=+vmax,
        colorscale=[[0, "rgb(0, 0, 0)"], [1, sheet_color]],
        contours=dict(
            x=dict(
                show=True,
                start=x.min(),
                end=x.max(),
                size=(x.max() - x.min()) / n_lines,
                color="black",
                width=1,
            ),
            y=dict(
                show=True,
                start=-y.max(),
                end=+y.max(),
                size=(y.max() - y.min()) / (n_lines // 2),
                color="black",
                width=1,
            ),
        ),
        name=sheet_name,
        opacity=0.4,
        showscale=False,
    )


vmax = 1.6
project = np.imag
sheet_colors = {
    "Physical (T1)": "blue",
    "Unphysical (T2)": "red",
}

Show code cell source Hide code cell source

# @title
Sp = go.Surface(x=X, y=Y, z=-T1p.imag, **sty("Physical (T1)"))
Sn = go.Surface(x=X, y=-Y, z=-T2n.imag, **sty("Unphysical (T2)"))
Sp.name = "Physical sheet I"

s_thr = values[m1] + values[m2]
threshold_filter = x >= s_thr
lineshape = go.Scatter3d(
    x=x[threshold_filter],
    y=np.zeros(threshold_filter.shape),
    z=project(-T1p[0])[threshold_filter],
    line=dict(color="yellow", width=10),
    mode="lines",
    name="Lineshape",
)
point = go.Scatter3d(
    x=[s_thr],
    y=[0],
    z=[0],
    mode="markers",
    marker=dict(color="black", size=6),
    name="Branch point",
)

fig = go.Figure(data=[Sn, Sp, lineshape, point])
fig.update_layout(
    height=550,
    margin=dict(l=0, r=0, t=30, b=0),
    showlegend=True,
    legend=dict(
        orientation="v",
        xanchor="left",
        yanchor="top",
        x=0.05,
        y=0.95,
        font=dict(size=24),
    ),
    title_text="Im(T) with Chew-Mandelstam phase space factor",
    title_font=dict(size=28),
    title=dict(y=0.989),
)
fig.update_scenes(
    camera_center=dict(z=-0.2),
    xaxis_title_text="Re √s",
    yaxis_title_text="Im √s",
    zaxis_title_text="Im T(s)",
    zaxis_range=[-vmax, +vmax],
)

fig.show()

The lineshape, the part that is observed within the experiment, is given as the intersection of the Riemann sheets with real plane. Also note that the second Riemann sheets transitions smoothly into the first one.

Attention

Not that the second Riemann sheet also inherits the singularity at \(s=0\), as it is derived from the common phasespace factor.