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 (1st quadrant of the complex plane) and the unphysical Riemann sheets are only defined for negative imaginary part (4th quadrant of the complex plane). For the single-channel case there are two Riemann sheets, one physical and one unphysical.

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.sympy import unevaluated
from IPython.display import Math

warnings.filterwarnings("ignore")

Phase space factor definitions#

Hide code cell source
# @title
from ampform.kinematics.phasespace import Kallen


@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}\]
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#

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())
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)
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()
../_images/66d7b830843c23860236fe49fb49b5c66670b1df6e50b172e2eacd01eaa678be.svg

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#

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",
}
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()