from collections.abc import Callable
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 colorsSee the Finding poles notebook for more explanation on the phase space factors and the procedure for finding a pole and computing a residue in the unphysical sheet.
@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) * sp.sqrt(s - (m1 - m2) ** 2) / s
@unevaluated(real=False)
class PhaseSpaceFactorCM(sp.Expr):
s: Any
m1: Any
m2: Any
_latex_repr_ = R"\rho_{{{m1},{m2}}}^\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))
class DiagonalMatrix(sp.DiagonalMatrix):
def _latex(self, printer, *args): # noqa: ARG002
return printer._print(self.args[0])
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_{m_{a},m_{b}}\left(s\right) &=& \frac{\sqrt{s - \left(m_{a} - m_{b}\right)^{2}} \sqrt{s - \left(m_{a} + m_{b}\right)^{2}}}{s} \\ \rho_{m_{a},m_{b}}^\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}\)
n = 2
K = sp.MatrixSymbol("K", n, n)
CM = DiagonalMatrix(sp.MatrixSymbol(R"\rho^\mathrm{CM}", n, n))
rho = DiagonalMatrix(sp.MatrixSymbol("rho", n, n))
I = sp.Identity(n)T1 = (I - sp.I * K * CM).inv() * K
T1\(\displaystyle \left(\mathbb{I} + - i K \rho^\mathrm{CM}\right)^{-1} K\)
s, g1, g2, m1 = sp.symbols("s g1 g2 m1")
ma, mb, mc, md = sp.symbols("m_a m_b m_c m_d")substitutions = {
K[0, 0]: g1**2 / (m1**2 - s),
K[1, 1]: g2**2 / (m1**2 - s),
K[0, 1]: g2 * g1 / (m1**2 - s),
K[1, 0]: g2 * g1 / (m1**2 - s),
CM[0, 0]: PhaseSpaceFactorCM(s, ma, mb),
CM[1, 1]: PhaseSpaceFactorCM(s, mc, md),
}
Math(aslatex(substitutions))\(\displaystyle \begin{array}{rcl} {K}_{0,0} &=& \frac{g_{1}^{2}}{m_{1}^{2} - s} \\ {K}_{1,1} &=& \frac{g_{2}^{2}}{m_{1}^{2} - s} \\ {K}_{0,1} &=& \frac{g_{1} g_{2}}{m_{1}^{2} - s} \\ {K}_{1,0} &=& \frac{g_{1} g_{2}}{m_{1}^{2} - s} \\ {\rho^\mathrm{CM}}_{0,0} &=& \rho_{m_{a},m_{b}}^\mathrm{CM}\left(s\right) \\ {\rho^\mathrm{CM}}_{1,1} &=& \rho_{m_{c},m_{d}}^\mathrm{CM}\left(s\right) \\ \end{array}\)
Indeed, the expression for sheet I looks like a Flatté function (multichannel Breit–Wigner)!
T1.as_explicit().subs(substitutions).simplify(doit=False)[0, 0]\(\displaystyle - \frac{g_{1}^{2}}{i g_{1}^{2} \rho_{m_{a},m_{b}}^\mathrm{CM}\left(s\right) + i g_{2}^{2} \rho_{m_{c},m_{d}}^\mathrm{CM}\left(s\right) - m_{1}^{2} + s}\)
T1_expr = T1.as_explicit().subs(substitutions)
T1_expr_00 = T1_expr[0, 0]
T1_expr_11 = T1_expr[1, 1]The second and third sheet are again calculated through the discontinuity of the \(T\) matrix across the branch cut:
\[ \begin{align} \operatorname{Disc}_{\mathrm{I,II}} T^{-1} &= 2 i\left[\begin{array}{rr}\rho_1 & 0 \\ 0 & 0 \end{array}\right] \\ \operatorname{Disc}_{\mathrm{I,III}} T^{-1} &= 2 i\left[\begin{array}{rr}\rho_1 & 0 \\ 0 & \rho_2 \end{array}\right] \end{align} \]
Depending on the centre-of-mass energy, different Riemann sheets connect smoothly to the physical one. In our case the resonance mass will be above the threshold for the first and the second channel. Therefore the resonance is expected to be located on the third sheet.
Ti = (T1.inv() - 2 * sp.I * rho).inv()
Ti\(\displaystyle \left(- 2 i \rho + K^{-1} \left(\mathbb{I} + - i K \rho^\mathrm{CM}\right)\right)^{-1}\)
rho_subs_II = {
**substitutions,
rho[0, 0]: PhaseSpaceFactor(s, ma, mb),
rho[1, 1]: 0,
}
rho_subs_III = {
**substitutions,
rho[0, 0]: PhaseSpaceFactor(s, ma, mb),
rho[1, 1]: PhaseSpaceFactor(s, mc, md),
}T2_expr = Ti.as_explicit().subs(rho_subs_II)
T3_expr = Ti.as_explicit().subs(rho_subs_III)parameters = {
m1: 0.8,
g1: 0.4,
g2: 0.3,
ma: 0.1,
mb: 0.15,
mc: 0.25,
md: 0.4,
}def lambdify_matrix_expression(matrix_expr: sp.Matrix) -> np.ndarray:
return np.array([
[sp.lambdify(s, matrix_expr[i, j].doit().subs(parameters)) for i in range(n)]
for j in range(n)
])
T1_func = lambdify_matrix_expression(T1_expr)
T2_func = lambdify_matrix_expression(T2_expr)
T3_func = lambdify_matrix_expression(T3_expr)thr1 = (parameters[ma] + parameters[mb]) ** 2
thr2 = (parameters[mc] + parameters[md]) ** 2Y_max = 1
X, Y = np.meshgrid(
np.linspace(0, 1.5, num=100),
np.linspace(1e-8, Y_max, num=100),
)
S = X + 1j * Y
style = dict(cmap=plt.cm.coolwarm, rasterized=True, vmin=-1, vmax=1)
text_kwargs = dict(c="black", fontsize=16)
fig, axes = plt.subplots(figsize=(8, 6), ncols=2, nrows=2, sharex=True, sharey=True)
ax1, ax2, ax3, ax4 = axes.flatten()
ax1.pcolormesh(X, Y, T1_func[0, 0](S).imag, **style)
ax1.pcolormesh(X, -Y, T3_func[0, 0](S.conj()).imag, **style)
ax1.text(0.32, 0.88, R"$T_\mathrm{I}^{11}$", transform=ax1.transAxes, **text_kwargs) # noqa: RUF027
ax1.text(0.32, 0.06, R"$T_\mathrm{III}^{11}$", transform=ax1.transAxes, **text_kwargs)
ax2.pcolormesh(X, +Y, T1_func[0, 1](S).imag, **style)
ax2.pcolormesh(X, -Y, T3_func[0, 1](S.conj()).imag, **style)
ax2.text(0.32, 0.88, R"$T_\mathrm{I}^{12}$", transform=ax2.transAxes, **text_kwargs) # noqa: RUF027
ax2.text(0.32, 0.06, R"$T_\mathrm{III}^{12}$", transform=ax2.transAxes, **text_kwargs)
ax3.pcolormesh(X, +Y, T1_func[1, 0](S).imag, **style)
ax3.pcolormesh(X, -Y, T2_func[1, 0](S.conj()).imag, **style)
ax3.text(0.32, 0.88, R"$T_\mathrm{I}^{21}$", transform=ax3.transAxes, **text_kwargs) # noqa: RUF027
ax3.text(0.32, 0.06, R"$T_\mathrm{II}^{21}$", transform=ax3.transAxes, **text_kwargs)
ax4.pcolormesh(X, +Y, T1_func[1, 1](S).imag, **style)
ax4.pcolormesh(X, -Y, T2_func[1, 1](S.conj()).imag, **style)
ax4.text(0.32, 0.88, R"$T_\mathrm{I}^{22}$", transform=ax4.transAxes, **text_kwargs) # noqa: RUF027
ax4.text(0.32, 0.06, R"$T_\mathrm{II}^{22}$", transform=ax4.transAxes, **text_kwargs)
for ax in axes.flatten():
ax.vlines(thr1, 0, Y_max, color="C0", ls="--")
ax.vlines(thr2, 0, Y_max, color="C1", ls="--")
for ax in axes[1]:
ax.set_xlabel("Re(s)")
for ax in axes[:, 0]:
ax.set_ylabel("Im(s)")
fig.tight_layout()
plt.show()
def get_denominator_functions(matrix_expr: sp.Matrix, parameters: dict) -> np.ndarray:
return np.array([
[lambdify_denominator(matrix_expr[i, j], parameters) for i in range(n)]
for j in range(n)
])
def lambdify_denominator(expr: sp.Expr, parameters: dict) -> Callable:
_, denominator = sp.fraction(expr)
return sp.lambdify(s, denominator.doit().subs(parameters))
T2_denom_func = get_denominator_functions(T2_expr, parameters)
T3_denom_func = get_denominator_functions(T3_expr, parameters)Z = np.abs(T3_denom_func[0, 0](S.conj()))
plt.title(R"Absolute value of denominator of $T_\mathrm{III}^{1, 1}$")
plt.pcolormesh(X, -Y, Z, rasterized=True, norm=colors.LogNorm(vmin=0.1, vmax=4))
plt.colorbar()
plt.xlabel("Re(s)")
plt.ylabel("Im(s)")
plt.show()
def fit_pole(func, s_guess: complex) -> complex:
def cost_function(s_real: float, s_imag: float) -> float:
s = s_real + s_imag * 1j
return np.abs(func(s)) ** 2
minuit2 = Minuit(cost_function, s_guess.real, s_guess.imag)
minuit2.tol = 0.001
fit_result = minuit2.migrad()
return complex(*fit_result.values)poles_positions_II = np.array([
[
fit_pole(T2_denom_func[0, 0], s_guess=0.8 - 0.3j),
fit_pole(T2_denom_func[0, 1], s_guess=0.8 - 0.3j),
],
[
fit_pole(T2_denom_func[1, 0], s_guess=0.8 - 0.3j),
fit_pole(T2_denom_func[1, 1], s_guess=0.8 - 0.3j),
],
])
poles_positions_II.round(3)array([[0.873-0.102j, 0.873-0.102j],
[0.873-0.102j, 0.873-0.102j]])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))residues_II = np.array([
[
compute_residue(T2_func[0, 0], z0=poles_positions_II[0, 0]),
compute_residue(T2_func[0, 1], z0=poles_positions_II[0, 1]),
],
[
compute_residue(T2_func[1, 0], z0=poles_positions_II[1, 0]),
compute_residue(T2_func[1, 1], z0=poles_positions_II[1, 1]),
],
])
residues_II.round(3)array([[-0.179-0.008j, -0.134-0.006j],
[-0.134-0.006j, -0.101-0.005j]])Interstingly, the residue matrix
\[ R_{ij} = \frac{1}{2\pi i} = \oint T^\mathrm{II}_{ij}(s) \, \mathrm{d}s \]
has determinant \(0\) and is always of rank \(1\), that is, it has only one non-zero eigenvalue.
np.testing.assert_almost_equal(np.linalg.det(residues_II), 0)eigenvalues = np.linalg.eigvals(residues_II)
eigenvalues.round(5)array([-0.27993-0.01324j, 0. -0.j ])We also have the relation
\[ g^T g = \lambda_1 \,, \]
where \(\lambda_1\) is the first eigenvalue on the eigenvalue matrix.