We have been gifted a \(K\)-matrix parameterization for two resonances within one channel, because it is International \(K\)-Matrix Day. The parameter values are:
Find pole positions and the residues of the two resonances with the given parameter values.
from typing import Any
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from ampform.io import aslatex
from ampform.kinematics.phasespace import Kallen
from ampform.sympy import unevaluated
from iminuit import Minuit
from IPython.display import Math
from matplotlib import colorsWe define the phase space factor \(\rho^\mathrm{CM}(s)\) using the Chew–Mandelstam function \(\Sigma(s)\) for \(S\) waves (\(L=0\)). For the discontinuity between the sheets, we use the ‘standard’ phase space factor \(\rho(s)\).
@unevaluated(real=False)
class PhaseSpaceFactor(sp.Expr):
s: Any
m1: Any
m2: Any
_latex_repr_ = R"\rho\left({s}\right)"
def evaluate(self) -> sp.Expr:
s, m1, m2 = self.args
return sp.sqrt(s - (m1 + m2) ** 2) * sp.sqrt(s - (m1 - m2) ** 2) / s
@unevaluated(real=False)
class PhaseSpaceFactorCM(sp.Expr):
s: Any
m1: Any
m2: Any
_latex_repr_ = R"\rho^\mathrm{{CM}}\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))
args = sp.symbols("s m_a m_b")
exprs = [
PhaseSpaceFactor(*args),
PhaseSpaceFactorCM(*args),
ChewMandelstam(*args),
BreakupMomentum(*args),
]
Math(aslatex({expr: expr.doit(deep=False) for expr in exprs}))\(\displaystyle \begin{array}{rcl} \rho\left(s\right) &=& \frac{\sqrt{s - \left(m_{a} - m_{b}\right)^{2}} \sqrt{s - \left(m_{a} + m_{b}\right)^{2}}}{s} \\ \rho^\mathrm{CM}\left(s\right) &=& - 16 i \pi \Sigma\left(s\right) \\ \Sigma\left(s\right) &=& \frac{- \left(m_{a}^{2} - m_{b}^{2}\right) \left(- \frac{1}{\left(m_{a} + m_{b}\right)^{2}} + \frac{1}{s}\right) \log{\left(\frac{m_{a}}{m_{b}} \right)} + \frac{2 \log{\left(\frac{m_{a}^{2} + m_{b}^{2} + 2 \sqrt{s} q\left(s\right) - s}{2 m_{a} m_{b}} \right)} q\left(s\right)}{\sqrt{s}}}{16 \pi^{2}} \\ q\left(s\right) &=& \frac{\sqrt{\lambda\left(s, m_{a}^{2}, m_{b}^{2}\right)}}{2 \sqrt{s}} \\ \end{array}\)
Compare PDG 2024, Figure 50.6:
rho_func = sp.lambdify(args, PhaseSpaceFactor(*args).doit())
rho_cm_func = sp.lambdify(args, PhaseSpaceFactorCM(*args).doit())
x = np.linspace(-0.1, 1.3, num=500)
epsilon = 1e-4j
threshold = dict(m_a=0.13, m_b=0.5)
z_rho = 1j * rho_func(x + epsilon, **threshold)
z_rho_cm = 1j * rho_cm_func(x + epsilon, **threshold)
fig, axes = plt.subplots(figsize=(9, 4), ncols=2, sharey=True)
ax1, ax2 = axes
ax1.plot(x, z_rho.real, c="black", label="Real part", ls="--")
ax1.plot(x, z_rho.imag, c="red", label="Imag part")
ax2.plot(x, z_rho_cm.real, c="black", label="Real part", ls="--")
ax2.plot(x, z_rho_cm.imag, c="red", label="Imag part")
ax1.set_title(R"$i\rho(s)$")
ax2.set_title(R"$i\rho^{CM}(s)$")
s_thr = sum(threshold.values()) ** 2
for ax in axes:
ax.axhline(0, alpha=0.5, c="black", lw=0.3)
ax.axvline(s_thr, c="black", lw=0.5)
ax.legend(loc="lower right")
ax.set_xlabel(R"$s+i\epsilon$ [GeV$^2$]")
ax.set_ylim(-1, +1)
fig.tight_layout()
plt.show(fig)
Construction of the two sheets of the \(T\) matrix using a single-channel \(K\) matrix:
n = 1
K = sp.MatrixSymbol("K", n, n)
rho_cm = sp.MatrixSymbol(R"\rho^\mathrm{CM}", n, n)
rho = sp.MatrixSymbol("rho", n, n)
I = sp.Identity(n)T1 = (I - sp.I * K * rho_cm).inv() * K
T1\(\displaystyle \left(\mathbb{I} + - i K \rho^\mathrm{CM}\right)^{-1} K\)
T1.as_explicit()\(\displaystyle \left[\begin{matrix}\frac{{K}_{0,0}}{- i {K}_{0,0} {\rho^\mathrm{CM}}_{0,0} + 1}\end{matrix}\right]\)
Definition of the second Riemann sheet via the discontinuity of the phasespace factor across the branch cut: \[ \begin{align} \mathrm{Disc}_\mathrm{I-II}(s) &= -2i\rho(s) \\ \rho(s) &= \frac{\sqrt{s-(m_a-m_b)^2}\sqrt{s-(m_a+m_b)^2}}{s} \end{align} \]
T2 = (T1.inv() - 2 * sp.I * rho).inv()
T2\(\displaystyle \left(K^{-1} \left(\mathbb{I} + - i K \rho^\mathrm{CM}\right) + - 2 i \rho\right)^{-1}\)
We now parametrize the \(K\) matrix with two poles and define the standard phase space factor (used for the discrepancy between the sheets), as well as a phase space factor that is defined with the Chew–Mandelstam function:
s, g1, m1, g2, m2 = sp.symbols("s g1 m1 g2 m2")
ma, mb = sp.symbols("m_a m_b")definitions = {
K[0, 0]: g1**2 / (m1**2 - s) + g2**2 / (m2**2 - s),
rho_cm[0, 0]: PhaseSpaceFactorCM(s, ma, mb),
rho[0, 0]: PhaseSpaceFactor(s, ma, mb),
}
Math(aslatex(definitions))\(\displaystyle \begin{array}{rcl} {K}_{0,0} &=& \frac{g_{1}^{2}}{m_{1}^{2} - s} + \frac{g_{2}^{2}}{m_{2}^{2} - s} \\ {\rho^\mathrm{CM}}_{0,0} &=& \rho^\mathrm{CM}\left(s\right) \\ {\rho}_{0,0} &=& \rho\left(s\right) \\ \end{array}\)
Definition of the I. sheet:
T1_expr = T1.as_explicit().subs(definitions)[0].simplify(doit=False)
T1_expr\(\displaystyle \frac{g_{1}^{2} \left(- m_{2}^{2} + s\right) + g_{2}^{2} \left(- m_{1}^{2} + s\right)}{- \left(m_{1}^{2} - s\right) \left(m_{2}^{2} - s\right) + i \left(g_{1}^{2} \left(m_{2}^{2} - s\right) + g_{2}^{2} \left(m_{1}^{2} - s\right)\right) \rho^\mathrm{CM}\left(s\right)}\)
Definition of the II. sheet:
T2_expr = T2.as_explicit().subs(definitions)[0].simplify(doit=False)
T2_expr\(\displaystyle \frac{g_{1}^{2} \left(- m_{2}^{2} + s\right) + g_{2}^{2} \left(- m_{1}^{2} + s\right)}{- \left(m_{1}^{2} - s\right) \left(m_{2}^{2} - s\right) + 2 i \left(g_{1}^{2} \left(m_{2}^{2} - s\right) + g_{2}^{2} \left(m_{1}^{2} - s\right)\right) \rho\left(s\right) + i \left(g_{1}^{2} \left(m_{2}^{2} - s\right) + g_{2}^{2} \left(m_{1}^{2} - s\right)\right) \rho^\mathrm{CM}\left(s\right)}\)
For plotting the \(T\) matrix, we need to define some parameter values.
parameters = {
m1: 1.8,
m2: 1.1,
g1: 0.5,
g2: 0.7,
ma: 0.1,
mb: 0.2,
}T1_func = sp.lambdify(s, T1_expr.doit().subs(parameters))
T2_func = sp.lambdify(s, T2_expr.doit().subs(parameters))Just below the real axis (below the cut), the square value of the unphysical sheet, \(T_\mathrm{II}(s-\epsilon i)\), looks like this:
x = np.linspace(0.0, 6, num=1000)
z = T1_func(x + 1e-9j)
s_thr = (parameters[ma] + parameters[mb]) ** 2
plt.plot(x, np.abs(z) ** 2)
plt.xlabel("s")
plt.ylabel(R"$|T(s)|^2$")
plt.yticks([])
plt.ylim(0, None)
plt.axvline(s_thr, c="red", label=R"$(m_a+m_b)^2$", ls="--", lw=1)
plt.legend()
plt.show()
The two sheets are nicely connected along the real axis and two poles are located below the cut.
X, Y = np.meshgrid(
np.linspace(0, 6, num=500),
np.linspace(1e-8, 2, num=500),
)
S = X + 1j * Y
style = dict(rasterized=True, vmin=-1, vmax=1)
plt.pcolormesh(X, +Y, T1_func(S).imag, **style)
plt.pcolormesh(X, -Y, T2_func(S.conj()).imag, **style)
plt.xlabel("Re(s)")
plt.ylabel("Im(s)")
plt.axvline(s_thr, c="red", lw=1, ls="--")
plt.axhline(0, c="black", lw=0.5)
cbar = plt.colorbar()
cbar.ax.set_title("Im(T)")
plt.show()
The pole positions are at the point where the denominator of the unphysical sheet \(T_\mathrm{II}\) goes to zero.
numerator, denominator = sp.fraction(T2_expr)
denominator\(\displaystyle - \left(m_{1}^{2} - s\right) \left(m_{2}^{2} - s\right) + 2 i \left(g_{1}^{2} \left(m_{2}^{2} - s\right) + g_{2}^{2} \left(m_{1}^{2} - s\right)\right) \rho\left(s\right) + i \left(g_{1}^{2} \left(m_{2}^{2} - s\right) + g_{2}^{2} \left(m_{1}^{2} - s\right)\right) \rho^\mathrm{CM}\left(s\right)\)
denom_func = sp.lambdify(s, denominator.doit().subs(parameters))Z = np.abs(denom_func(S.conj()))
plt.pcolormesh(X, -Y, Z, rasterized=True, norm=colors.LogNorm(vmin=0.1, vmax=4))
cbar = plt.colorbar()
cbar.ax.set_title(R"denominator $T(s)$")
plt.xlabel("Re(s)")
plt.ylabel("Im(s)")
plt.show()
We can find these points by minimizing the absolute value of the denominator of \(T_\mathrm{II}\) using gradient-descent.
s_guess1 = 1.8 - 0.26j
s_guess2 = 3.7 - 0.47jdef fit_pole(s_guess: complex) -> Minuit:
minuit2 = Minuit(cost_function, s_guess.real, s_guess.imag)
minuit2.tol = 0.001
return minuit2.migrad()
def cost_function(s_real: float, s_imag: float):
s = s_real + s_imag * 1j
return np.abs(denom_func(s)) ** 2fit_pole1 = fit_pole(s_guess1)
pole1 = complex(*fit_pole1.values)
pole1(1.8152362212906634-0.31507035582340487j)fit_pole2 = fit_pole(s_guess2)
pole2 = complex(*fit_pole2.values)
pole2(3.7266385502204797-0.4849422565013708j)Plot pole positions:
plt.pcolormesh(X, +Y, T1_func(S).imag, **style)
plt.pcolormesh(X, -Y, T2_func(S.conj()).imag, **style)
plt.plot(pole1.real, pole1.imag, "rx", markersize=10)
plt.plot(pole2.real, pole2.imag, "rx", markersize=10)
plt.xlabel("Re(s)")
plt.ylabel("Im(s)")
plt.axhline(0, c="black", lw=0.5)
plt.axvline(s_thr, c="red", ls="--", lw=1)
cbar = plt.colorbar()
cbar.ax.set_title("Im(T)")
plt.show()
Residues can be computed with a closed integral around the pole. A natural choice is to integrate over a circle \(\gamma=s_\mathrm{p}+re^{i\phi}\) around the pole \(s_\mathrm{p}\). Numerically, we can compute this by evaluating \(T_\mathrm{II}\) at \(N\) points around on that circle. The resulting sum can be simplified to: \[ \begin{align} \mathrm{Res}(T, s_\mathrm{p}) &= \frac{1}{2\pi i} \oint_\gamma T(z) \, dz \\ &= \frac{1}{2\pi i} \int_{-\pi}^{+\pi} T(z) \; \underbrace{i r e^{i\phi}}_{dz/d\phi} \; d\phi \quad \text{with} \quad z = s_\mathrm{p} + r e^{i\phi} \\ &= \frac{r}{2\pi} \int_{-\pi}^{+\pi} T(z) \, e^{i\phi} \, d\phi \\ &= \frac{r}{2\pi} \left(\frac{2\pi}{N}\sum_{i=1}^N T(z_i) \, e^{i\phi_i}\right) \\ &= \frac{r}{N} \sum_{i=1}^N T(z_i) \, e^{i\phi_i} \,. \end{align} \]
Define residue computation:
def compute_residue(f, z0, radius=1e-3, n_points=1_000):
phi = np.linspace(-np.pi, np.pi, n_points, endpoint=False)
z = z0 + radius * np.exp(1j * phi)
return radius / n_points * np.sum(f(z) * np.exp(1j * phi))residue1 = compute_residue(T2_func, pole1)
residue1.round(3)np.complex128(-0.303-0.172j)residue2 = compute_residue(T2_func, pole2)
residue2.round(3)np.complex128(-0.508+0.151j)As a sanity check, we can verify whether \(\lim\limits_{z \to z_0} T(z) = \frac{\mathrm{Res}(T, z_0)}{z - z_0}\), for instance by ploting \(T(z)\) over \(z=z_0+re^{i\phi}\) with \(\phi \in [-\pi, +\pi]\).
phi = np.linspace(-np.pi, np.pi, num=1000)
r = 0.1
z = pole1 + r * np.exp(1j * phi)
plt.plot(phi, T2_func(z).imag, ls="--")
plt.plot(phi, (residue1 / (z - pole1)).imag, ls="--")
plt.xlabel(R"$\phi$")
plt.ylabel("Im(T)")
plt.show()
The residue approximates the coupling of the poles as if they were to be considered as a single pole in a Breit–Wigner parameterization.
\[ T\approx \frac{-\mathrm{Res}}{s_\mathrm{pole}-s} \Rightarrow g_\mathrm{pole} \approx \sqrt{\mathrm{Res}} \]
g_pole1 = np.sqrt(residue1)
g_pole2 = np.sqrt(residue2)Define Breit–Wigner function for the poles:
def breit_wigner(s, g, s_pole):
return g**2 / (s_pole - s)Plot line shape of the pole Breit–Wigner functions and the \(T\) function:
x = np.linspace(0, 6, num=500)
z = T1_func(x + 1e-8j)
z_pole1 = breit_wigner(x, g_pole1, pole1)
z_pole2 = breit_wigner(x, g_pole2, pole2)
plt.plot(x, np.abs(z) ** 2, label="T")
plt.plot(x, np.abs(z_pole1) ** 2, label=R"BW$_\mathrm{Pole1}$", ls="--")
plt.plot(x, np.abs(z_pole2) ** 2, label=R"BW$_\mathrm{Pole2}$", ls="--")
plt.plot(
x,
np.abs(z_pole1 + z_pole2) ** 2,
label=R"$\mathrm{BW}_\mathrm{Pole1}+\mathrm{BW}_\mathrm{Pole2}$",
)
plt.axvline(s_thr, c="black", ls="dotted", lw=1)
plt.xlabel("s")
plt.ylabel(R"Intensity")
plt.ylim(0, None)
plt.legend()
plt.show()