Continuation into other Riemann sheets

The scattering amplitude can be analytically continued into the complex plane from the observable lineshape on the real axis. Its analytic structure is inherited from the phase space factor, where each two-body threshold results in a branch point. From this branch point onwards, the branch cut continues to \(+\infty\) along the real axis, giving rise to two solutions called Riemann sheets. Each additional two-body threshold branch point gives rise to another pair of Riemann sheets. The sheet continued from the physical amplitude is called the physical sheet. Continuing this sheet over the branch cut leads us to the unphysical Riemann sheets, which contains the poles that characterize a hadronic state. This notebook aims to calculate the unphysical Riemann sheets for \(L = 1\) using the Chew–Mandelstam dispersion integral.

See also

TR-026 and TR-027.

Import Python libraries

Hide code cell source

import warnings

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from IPython.display import Math
from matplotlib_inline.backend_inline import set_matplotlib_formats

from ampform.dynamics.form_factor import FormFactor
from ampform.dynamics.phasespace import (
    ChewMandelstamIntegral,
    PhaseSpaceFactor,
    chew_mandelstam_s_wave,
)
from ampform.io import aslatex, improve_latex_rendering

improve_latex_rendering()
plt.rc("font", size=14)
set_matplotlib_formats("svg")
warnings.filterwarnings("ignore")

\(T\) matrix definition

According to PDG 2025, §50. Resonances, Eq. (50.41), the scattering amplitude is given by

\[ M = n \left( I - i K \rho n^2 \right)^{-1} K n \,, \]

where \(n\) is a matrix with barrier factor and \(\rho\) is a matrix of phase space factors. Both are diagonal matrices containing the respective factors for each channel. In order to have an amplitude function that is analytic over the complex plane, the factor \(i\rho n^2\) is replaced by the once-subtracted Chew–Mandelstam dispersion integral \(\Sigma^l\).

Matrix symbol definitions

Hide code cell source

class DiagonalMatrix(sp.DiagonalMatrix):
    def _latex(self, printer, *args):
        return printer._print(self.args[0])


n = 2
I = sp.Identity(n)
K = sp.MatrixSymbol("K", n, n)
CM = DiagonalMatrix(sp.MatrixSymbol(R"\Sigma^l", n, n))
FF = DiagonalMatrix(sp.MatrixSymbol(R"n", n, n))
Math(aslatex({CM: CM.as_explicit()}))

\[\begin{split}\displaystyle \begin{aligned} \Sigma^{l} \;&=\; \left[\begin{matrix}{\Sigma^{l}}_{0,0} & 0\\0 & {\Sigma^{l}}_{1,1}\end{matrix}\right] \\ \end{aligned}\end{split}\]

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

\[\displaystyle n \left(\mathbb{I} + K \Sigma^{l}\right)^{-1} K n\]

T1.as_explicit()[0, 0].simplify(doit=False)

\[\displaystyle \frac{\left(\left({K}_{1,1} {\Sigma^{l}}_{1,1} + 1\right) {K}_{0,0} - {K}_{0,1} {K}_{1,0} {\Sigma^{l}}_{1,1}\right) {n}_{0,0}^{2}}{{K}_{0,0} {K}_{1,1} {\Sigma^{l}}_{0,0} {\Sigma^{l}}_{1,1} + {K}_{0,0} {\Sigma^{l}}_{0,0} - {K}_{0,1} {K}_{1,0} {\Sigma^{l}}_{0,0} {\Sigma^{l}}_{1,1} + {K}_{1,1} {\Sigma^{l}}_{1,1} + 1}\]

Symbol definitions for pole parametrization

Hide code cell source

s = sp.Symbol("s")
ma1 = sp.Symbol("m_{a1}")
mb1 = sp.Symbol("m_{b1}")
ma2 = sp.Symbol("m_{a2}")
mb2 = sp.Symbol("m_{b2}")
m0 = sp.Symbol("m0")
g1 = sp.Symbol(R"g^{0}_1")
g2 = sp.Symbol(R"g^{0}_2")

Symbolic dispersion integral

The general form of the Chew–Mandelstam dispersion integral is the following:

s, m1, m2, L = sp.symbols("s, m1, m2, L", nonnegative=True)
integral_expr = ChewMandelstamIntegral(s, m1, m2, L)
integral_expr.doit(deep=False)

\[\displaystyle \frac{s - \left(m_{1} + m_{2}\right)^{2}}{\pi} \int\limits_{\left(m_{1} + m_{2}\right)^{2}}^{\infty} \frac{\mathcal{F}_{L}\left(x, m_{1}, m_{2}\right)^{2} \rho\left(x\right)}{\left(x - \left(m_{1} + m_{2}\right)^{2}\right) \left(- i \epsilon - s + x\right)}\, dx\]

For S-waves (\(L=0\)), it has an analytic solution, given by

chew_mandelstam_s_wave(s, m1, m2)

\[\displaystyle \frac{\frac{2 q^\mathrm{c}\left(s\right)}{\sqrt{s}} \log{\left(\frac{m_{1}^{2} + m_{2}^{2} + 2 \sqrt{s} q^\mathrm{c}\left(s\right) - s}{2 m_{1} m_{2}} \right)} - \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)}}{\pi}\]

\(T\) matrix parametrization

Show code cell source

Hide code cell source

L_val = 1
T1_expressions = {
    K[0, 0]: (g1 * g2 * m0) / (s - m0**2),
    K[0, 1]: (g1 * g2 * m0) / (s - m0**2),
    K[1, 0]: (g1 * g2 * m0) / (s - m0**2),
    K[1, 1]: (g2 * g2 * m0) / (s - m0**2),
    FF[0, 0]: FormFactor(s, ma1, mb1, angular_momentum=L_val),
    FF[1, 1]: FormFactor(s, ma2, mb2, angular_momentum=L_val),
}
Math(aslatex(T1_expressions))

\[\begin{split}\displaystyle \begin{aligned} {K}_{0,0} \;&=\; \frac{g^{0}_1 g^{0}_2 m_{0}}{- m_{0}^{2} + s} \\ {K}_{0,1} \;&=\; \frac{g^{0}_1 g^{0}_2 m_{0}}{- m_{0}^{2} + s} \\ {K}_{1,0} \;&=\; \frac{g^{0}_1 g^{0}_2 m_{0}}{- m_{0}^{2} + s} \\ {K}_{1,1} \;&=\; \frac{\left(g^{0}_2\right)^{2} m_{0}}{- m_{0}^{2} + s} \\ {n}_{0,0} \;&=\; \mathcal{F}_{1}\left(s, m_{a1}, m_{b1}\right) \\ {n}_{1,1} \;&=\; \mathcal{F}_{1}\left(s, m_{a2}, m_{b2}\right) \\ \end{aligned}\end{split}\]

T1_expr = T1.as_explicit().xreplace(T1_expressions)
T1_expr[0, 0].simplify(doit=False)

\[\displaystyle \frac{g^{0}_1 g^{0}_2 m_{0} \left(- g^{0}_1 g^{0}_2 m_{0} {\Sigma^{l}}_{1,1} + \left(g^{0}_2\right)^{2} m_{0} {\Sigma^{l}}_{1,1} - m_{0}^{2} + s\right) \mathcal{F}_{1}\left(s, m_{a1}, m_{b1}\right)^{2}}{g^{0}_1 \left(g^{0}_2\right)^{2} m_{0}^{2} \left(- g^{0}_1 + g^{0}_2\right) {\Sigma^{l}}_{0,0} {\Sigma^{l}}_{1,1} - g^{0}_2 m_{0} \left(m_{0}^{2} - s\right) \left(g^{0}_1 {\Sigma^{l}}_{0,0} + g^{0}_2 {\Sigma^{l}}_{1,1}\right) + \left(m_{0}^{2} - s\right)^{2}}\]

Sheets II, III, and IV

In the case of two channels, there are four Riemann sheets: one physical and three unphysical ones. The physical sheet is calculated using the analytic solution of the Chew–Mandelstam function. The other sheets are reached by adding the discontinuity across the branch cut:

For higher angular momentum \(l\):

\[ \operatorname{Disc} \varSigma_\ell(s) \;=\; 2i\,\rho(s)\,n_\ell^2(s) \,, \]

for instance,

\[ T^{II}=n\left(\mathbb{I} + K (\varSigma_\ell + 2 i \rho)\right)^{-1} K n \,. \]

Depending on the centre-of-mass energy, different Riemann sheets connect smoothly to the physical one. Therefore, two cases are studied: one where the resonance mass is above the threshold of the second and first channel, and another where the resonance mass is between the threshold of the first and second channel. For the 2 channel one gets:

\[\begin{split} \begin{eqnarray} \operatorname{Disc}_{\mathrm{I,II}} \varSigma_\ell(s) &=& 2 i\left[\begin{array}{rr}\rho_1 n_{\ell,1}^2(s) & 0 \\ 0 & 0 \end{array}\right], \\ \operatorname{Disc}_{\mathrm{I,III}} \varSigma_\ell(s) &=& 2 i\left[\begin{array}{rr}\rho_1 n_{\ell,1}^2(s) & 0 \\ 0 & \rho_2 n_{\ell,2}^2(s) \end{array}\right], \\ \operatorname{Disc}_{\mathrm{I,IV}} \varSigma_\ell(s) &=& 2 i\left[\begin{array}{rr}0 & 0 \\ 0& \rho_2 n_{\ell,2}^2(s) \end{array}\right]. \end{eqnarray} \end{split}\]

Show code cell source

Hide code cell source

rho = DiagonalMatrix(sp.MatrixSymbol("rho", n, n))
Math(aslatex({rho: rho.as_explicit()}))

\[\begin{split}\displaystyle \begin{aligned} \rho \;&=\; \left[\begin{matrix}{\rho}_{0,0} & 0\\0 & {\rho}_{1,1}\end{matrix}\right] \\ \end{aligned}\end{split}\]

T_unphysical = FF * (I + K * (CM + sp.I * 2 * rho * FF**2)).inv() * K * FF

Check full expression form

Hide code cell content

T_unphysical.as_explicit()[0, 0].simplify(doit=False)

\[\displaystyle \frac{\left(- \left({\Sigma^{l}}_{1,1} + 2 i {n}_{1,1}^{2} {\rho}_{1,1}\right) {K}_{0,1} {K}_{1,0} + \left({K}_{1,1} {\Sigma^{l}}_{1,1} + 2 i {K}_{1,1} {n}_{1,1}^{2} {\rho}_{1,1} + 1\right) {K}_{0,0}\right) {n}_{0,0}^{2}}{{K}_{0,0} {K}_{1,1} {\Sigma^{l}}_{0,0} {\Sigma^{l}}_{1,1} + 2 i {K}_{0,0} {K}_{1,1} {\Sigma^{l}}_{0,0} {n}_{1,1}^{2} {\rho}_{1,1} + 2 i {K}_{0,0} {K}_{1,1} {\Sigma^{l}}_{1,1} {n}_{0,0}^{2} {\rho}_{0,0} - 4 {K}_{0,0} {K}_{1,1} {n}_{0,0}^{2} {n}_{1,1}^{2} {\rho}_{0,0} {\rho}_{1,1} + {K}_{0,0} {\Sigma^{l}}_{0,0} + 2 i {K}_{0,0} {n}_{0,0}^{2} {\rho}_{0,0} - {K}_{0,1} {K}_{1,0} {\Sigma^{l}}_{0,0} {\Sigma^{l}}_{1,1} - 2 i {K}_{0,1} {K}_{1,0} {\Sigma^{l}}_{0,0} {n}_{1,1}^{2} {\rho}_{1,1} - 2 i {K}_{0,1} {K}_{1,0} {\Sigma^{l}}_{1,1} {n}_{0,0}^{2} {\rho}_{0,0} + 4 {K}_{0,1} {K}_{1,0} {n}_{0,0}^{2} {n}_{1,1}^{2} {\rho}_{0,0} {\rho}_{1,1} + {K}_{1,1} {\Sigma^{l}}_{1,1} + 2 i {K}_{1,1} {n}_{1,1}^{2} {\rho}_{1,1} + 1}\]

Substitute sheet expressions

Hide code cell source

T2_expressions = {
    **T1_expressions,
    rho[0, 0]: PhaseSpaceFactor(s, ma1, mb1),
    rho[1, 1]: 0,
}
T3_expressions = {
    **T1_expressions,
    rho[0, 0]: PhaseSpaceFactor(s, ma1, mb1),
    rho[1, 1]: PhaseSpaceFactor(s, ma2, mb2),
}
T4_expressions = {
    **T1_expressions,
    rho[0, 0]: 0,
    rho[1, 1]: PhaseSpaceFactor(s, ma2, mb2),
}
T2_expr = T_unphysical.as_explicit().xreplace(T2_expressions)
T3_expr = T_unphysical.as_explicit().xreplace(T3_expressions)
T4_expr = T_unphysical.as_explicit().xreplace(T4_expressions)

Visualize sheets

2D

Parameter values

Hide code cell source

parameters = {
    ma1: 1.0,
    mb1: 1.5,
    ma2: 1.5,
    mb2: 2.0,
    m0: 3.0,
    g1: 1,
    g2: 1,
}
Math(aslatex(parameters))

\[\begin{split}\displaystyle \begin{aligned} m_{a1} \;&=\; 1.0 \\ m_{b1} \;&=\; 1.5 \\ m_{a2} \;&=\; 1.5 \\ m_{b2} \;&=\; 2.0 \\ m_{0} \;&=\; 3.0 \\ g^{0}_1 \;&=\; 1 \\ g^{0}_2 \;&=\; 1 \\ \end{aligned}\end{split}\]

Lambdify expressions

Hide code cell source

symbols = (s, CM[0, 0], CM[1, 1])
T1_func = sp.lambdify(symbols, T1_expr[0, 0].subs(parameters).doit())
T2_func = sp.lambdify(symbols, T2_expr[0, 0].subs(parameters).doit())
T3_func = sp.lambdify(symbols, T3_expr[0, 0].subs(parameters).doit())
T4_func = sp.lambdify(symbols, T4_expr[0, 0].subs(parameters).doit())

integral_p_wave_func = sp.lambdify(
    [s, m1, m2, integral_expr.epsilon],
    integral_expr.subs(L, L_val).doit(),
)

Helper functions

Hide code cell source

def evaluate(func, s_array):
    args = transform(s_array)
    return func(*args)


def transform(s):
    cm1_array = integral_p_wave_func(s, parameters[ma1], parameters[mb1], ε_integral)
    cm2_array = integral_p_wave_func(s, parameters[ma2], parameters[mb2], ε_integral)
    return s, cm1_array, cm2_array


ε_integral = 1e-8

Define meshgrid

Hide code cell source

x = np.linspace(0, 8, num=300)
y = np.linspace(1e-3, 1, num=100)
X, Y = np.meshgrid(x, y)
Zn = X - Y * 1j
Zp = X + Y * 1j

Evaluate numerical integrals

Hide code cell source

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

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

T3n = evaluate(T3_func, Zn**2)
T3p = evaluate(T3_func, Zp**2)

T4n = evaluate(T4_func, Zn**2)
T4p = evaluate(T4_func, Zp**2)

Show code cell source

Hide code cell source

fig, axes = plt.subplots(figsize=(11, 4), ncols=2, sharey=True)
ax1, ax2 = axes.ravel()
fig.subplots_adjust(wspace=0.1)
fig.suptitle(R"$s_{thr1}<s_{thr2}<m_{res}$", fontsize=20)
ax1.set_title("I and II", y=0.88)
ax2.set_title("I and III", y=0.88)
for ax in axes.ravel():
    ax.set_xlabel(R"$\mathrm{Re}\,\sqrt{s}$")
ax1.set_ylabel(R"$\mathrm{Im}\,\sqrt{s}$")

T_max = 2
style = dict(vmin=-T_max, vmax=+T_max, cmap="RdBu", rasterized=True)
mesh = ax1.pcolormesh(X, Y, T1p.imag, **style)
ax1.pcolormesh(X, -Y, T2n.imag, **style)
ax2.pcolormesh(X, +Y, T1p.imag, **style)
ax2.pcolormesh(X, -Y, T3n.imag, **style)

s_thr1 = parameters[ma1] + parameters[mb1]
s_thr2 = parameters[ma2] + parameters[mb2]
linestyle = dict(ls="dotted", lw=1)
for ax in [ax1, ax2]:
    ax.axhline(0, c="black", **linestyle)
    ax.axvline(s_thr1, c="C0", **linestyle, label=R"$\sqrt{s_\mathrm{thr1}}$")
    ax.axvline(s_thr2, c="C1", **linestyle, label=R"$\sqrt{s_\mathrm{thr2}}$")

linestyle_res = dict(c="r", ls="dotted", label=R"$m_\mathrm{res}$")
for ax in [ax1, ax2]:
    ax.axvline(parameters[m0], **linestyle_res)

ax1.legend(fontsize=11, loc="lower left")
cbar = fig.colorbar(mesh, ax=axes)
cbar.ax.set_ylabel(R"$\mathrm{Im} T(s)$")
fig.show()

Lineshape

Show code cell source

Hide code cell source

x = np.linspace(2.2, 8, num=100)
ε_offset = 1e-2
cn = x - ε_offset * 1j
cp = x + ε_offset * 1j
fig, axes = plt.subplots(figsize=(12, 6), ncols=4, sharey=True)
ax1, ax2, ax3, ax4 = axes.ravel()

ax1.plot(x, evaluate(T1_func, cp**2).imag, label=R"$T_\mathrm{I}(s+0i)$")
ax1.plot(x, evaluate(T1_func, cn**2).imag, label=R"$T_\mathrm{I}(s-0i)$")
ax1.set_title(R"$T_\mathrm{I}$")

ax2.plot(x, evaluate(T2_func, cp**2).imag, label=R"$T_\mathrm{II}(s+0i)$")
ax2.plot(x, evaluate(T2_func, cn**2).imag, label=R"$T_\mathrm{II}(s-0i)$")
ax2.set_title(R"$T_\mathrm{II}$")

ax3.plot(x, evaluate(T3_func, cp**2).imag, label=R"$T_\mathrm{III}(s+0i)$")
ax3.plot(x, evaluate(T3_func, cn**2).imag, label=R"$T_\mathrm{III}(s-0i)$")
ax3.set_title(R"$T_\mathrm{III}$")

ax4.plot(x, evaluate(T4_func, cp**2).imag, label=R"$T_\mathrm{IV}(s+0i)$")
ax4.plot(x, evaluate(T4_func, cn**2).imag, label=R"$T_\mathrm{IV}(s-0i)$")
ax4.set_title(R"$T_\mathrm{IV}$")

for ax in axes:
    ax.axvline(parameters[m0], **linestyle_res)
    ax.axvline(s_thr1, c="C0", **linestyle, label=R"$\sqrt{s_\mathrm{thr1}}$")
    ax.axvline(s_thr2, c="C1", **linestyle, label=R"$\sqrt{s_\mathrm{thr2}}$")
    ax.legend()
    ax.set_xlabel(R"$\mathrm{Re}\,\sqrt{s}$")
    ax.set_ylim(-2, +2)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
ax1.set_ylabel(R"$\mathrm{Im}\,T(s)$ (a.u.)")

fig.tight_layout()
plt.show()

Caution

The integral does not appear to be numerically stable (see ComPWA/ampform#487), as a significant imaginary offset from the real axis (\(\epsilon = 10^{-2}i\)) is required to get rid of the instability fluctuations.

Connection (\(Im(T)\) v.s. \(Im(\sqrt{s})\) )

Depending on the energy, a different unphysical Riemann sheet connects to the physical one. The structures on the connecting sheet have the most influence on the physical line shape.

Show code cell source

Hide code cell source

y = np.linspace(1e-8, 1, num=100)
for s_value in (3.1, 3.7):
    Jn = s_value - y * 1j
    Jp = s_value + y * 1j
    fig, ax = plt.subplots(figsize=(9, 3))
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.plot(+y, evaluate(T1_func, Jp**2).imag, label=R"$T_\mathrm{I}$", c="C0")
    ax.plot(-y, evaluate(T1_func, Jn**2).imag, c="C0")
    ax.plot(+y, evaluate(T2_func, Jp**2).imag, label=R"$T_\mathrm{II}$", c="C1")
    ax.plot(-y, evaluate(T2_func, Jn**2).imag, c="C1")
    ax.plot(+y, evaluate(T3_func, Jp**2).imag, label=R"$T_\mathrm{III}$", c="C2")
    ax.plot(-y, evaluate(T3_func, Jn**2).imag, c="C2")
    ax.plot(+y, evaluate(T4_func, Jp**2).imag, label=R"$T_\mathrm{IV}$", c="C3")
    ax.plot(-y, evaluate(T4_func, Jn**2).imag, c="C3")
    ax.set_title(Rf"Re $\sqrt{{s}}={s_value:g}$")
    ax.set_xlabel(R"$\mathrm{Im}\,\sqrt{s}$")
    ax.set_ylabel(R"$\mathrm{Im}\,T(s)$ (a.u.)")
    ax.legend()
    fig.tight_layout()
    plt.show()