K-matrix

While ampform does not yet provide a generic way to formulate an amplitude model with \(\boldsymbol{K}\)-matrix dynamics, the (experimental) kmatrix module makes it fairly simple to produce a symbolic expression for a parameterized \(\boldsymbol{K}\)-matrix with an arbitrary number of poles and channels and play around with it interactively. For more info on the \(\boldsymbol{K}\)-matrix, see the classic paper by Chung [4], PDG2020, §Resonances, or this instructive presentation [5].

Section Physics summarizes [4], so that the kmatrix module can reference to the equations. It also points out some subtleties and deviations.

Note

The \(\boldsymbol{K}\)-matrix approach was originally worked worked out in [TR-005] K-matrix, [TR-009] Lorentz-invariant K-matrix, and [TR-010] P-vector. Those reports contained a few mistakes, which have been addressed here.

%matplotlib widget
import logging
import re
import warnings

import graphviz
import ipywidgets
import matplotlib.pyplot as plt
import mpl_interactions.ipyplot as iplt
import numpy as np
import sympy as sp
from IPython.display import Image
from matplotlib import cm
from mpl_interactions.controller import Controls

import symplot
from ampform.dynamics import (
    BlattWeisskopfSquared,
    PhaseSpaceFactor,
    kmatrix,
    relativistic_breit_wigner,
)

logging.basicConfig()
logging.getLogger().setLevel(logging.ERROR)

warnings.filterwarnings("ignore")

BlattWeisskopfSquared.max_angular_momentum = 3

Physics

The \(\boldsymbol{K}\)-matrix formalism is used to describe coupled, two-body formation processes of the form \(c_j d_j \to R \to a_i b_i\), with \(i,j\) representing each separate channel and \(R\) an intermediate state by which these channels are coupled.

dot = """
digraph {
    rankdir=LR;
    node [shape=point, width=0];
    edge [arrowhead=none];
    "Na" [shape=none, label="aᵢ"];
    "Nb" [shape=none, label="bᵢ"];
    "Nc" [shape=none, label="cⱼ"];
    "Nd" [shape=none, label="dⱼ"];
    { rank=same "Nc", "Nd" };
    { rank=same "Na", "Nb" };
    "Nc" -> "N0";
    "Nd" -> "N0";
    "N1" -> "Na";
    "N1" -> "Nb";
    "N0" -> "N1" [label="R"];
    "N0" [shape=none, label=""];
    "N1" [shape=none, label=""];
}
"""
graphviz.Source(dot)
../../_images/k-matrix_9_0.svg

A small adaptation allows us to describe a coupled, two-body production process of the form \(R \to a_ib_i\) (see Production processes).

In the following, \(n\) denotes the number of channels and \(n_R\) the number of poles. In the kmatrix module, we use \(0 \leq i,j < n\) and \(1 \leq R \leq n_R\).

Partial wave expansion

In amplitude analysis, the main aim is to express the differential cross section \(\frac{d\sigma}{d\Omega}\), that is, the intensity distribution in each spherical direction \(\Omega=(\phi,\theta)\) as we can observe in experiments. This differential cross section can be expressed in terms of the scattering amplitude \(A\):

(1)\[ \frac{d\sigma}{d\Omega} = \left|A(\Omega)\right|^2. \]

We can now further express \(A\) in terms of partial wave amplitudes by expanding it in terms of its angular momentum components:1

(2)\[ A(\Omega) = \frac{1}{q_i}\sum_L\left(2L+1\right) T_L(s) {D^{*L}_{\lambda\mu}}\left(\phi,\theta,0\right) \]

with \(L\) the total angular momentum of the decay product pair, \(\lambda=\lambda_a-\lambda_b\) and \(\mu=\lambda_c-\lambda_d\) the helicity differences of the final and initial states, \(D\) a Wigner-\(D\) function, and \(T_J\) an operator representing the partial wave amplitude.

The above sketch is just with one channel in mind. The same holds true though for a number of channels \(n\), with the only difference being that the \(T\) operator becomes an \(n\times n\) \(\boldsymbol{T}\)-matrix.

Transition operator

The important point is that we have now expressed \(A\) in terms of an angular part (depending on \(\Omega\)) and a dynamical part \(\boldsymbol{T}\) that depends on the Mandelstam variable \(s\).

The dynamical part \(\boldsymbol{T}\) is usually called the transition operator. It describes the interacting part of the more general scattering operator \(\boldsymbol{S}\), which describes the (complex) amplitude \(\langle f|\boldsymbol{S}|i\rangle\) of an initial state \(|i\rangle\) transitioning to a final state \(|f\rangle\). The scattering operator describes both the non-interacting amplitude and the transition amplitude, so it relates to the transition operator as:

(3)\[ \boldsymbol{S} = \boldsymbol{I} + 2i\boldsymbol{T} \]

with \(\boldsymbol{I}\) the identity operator. Just like in [4], we use a factor 2, while other authors choose \(\boldsymbol{S} = \boldsymbol{I} + i\boldsymbol{T}\). In that case, one would have to multiply Eq. (2) by a factor \(\frac{1}{2}\).

Ensuring unitarity

Knowing the origin of the \(\boldsymbol{T}\)-matrix, there is an important restriction that we need to comply with when we further formulate a parametrization: unitarity. This means that \(\boldsymbol{S}\) should conserve probability, namely \(\boldsymbol{S}^\dagger\boldsymbol{S} = \boldsymbol{I}\). Luckily, there is a trick that makes this easier. If we express \(\boldsymbol{S}\) in terms of an operator \(\boldsymbol{K}\) by applying a Cayley transformation:

(4)\[ \boldsymbol{S} = (\boldsymbol{I} + i\boldsymbol{K})(I - i\boldsymbol{K})^{-1}, \]

unitarity is conserved if \(\boldsymbol{K}\) is real. With some matrix jumbling, we can derive that the \(\boldsymbol{T}\)-matrix can be expressed in terms of \(\boldsymbol{K}\) as follows:

(5)\[ \boldsymbol{T} = \boldsymbol{K} \left(\boldsymbol{I} - i\boldsymbol{K}\right)^{-1} = \left(\boldsymbol{I} - i\boldsymbol{K}\right)^{-1} \boldsymbol{K}. \]

Lorentz-invariance

The description so far did not take Lorentz-invariance into account. For this, we first need to define a two-body phase space matrix \(\boldsymbol{\rho}\):

(6)\[\begin{split} \boldsymbol{\rho} = \begin{pmatrix} \rho_0 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \rho_{n-1} \end{pmatrix}. \end{split}\]

with \(\rho_i\) given by (2) in PhaseSpaceFactor for the final state masses \(m_{a,i}, m_{b,i}\). The Lorentz-invariant amplitude \(\boldsymbol{\hat{T}}\) and corresponding Lorentz-invariant \(\boldsymbol{\hat{K}}\)-matrix can then be computed from \(\boldsymbol{T}\) and \(\boldsymbol{K}\) with:2

(7)\[\begin{split} \begin{eqnarray} \boldsymbol{T} & = & \sqrt{\boldsymbol{\rho}} \; \boldsymbol{\hat{T}} \sqrt{\boldsymbol{\rho}} \\ \boldsymbol{K} & = & \sqrt{\boldsymbol{\rho}} \; \boldsymbol{\hat{K}} \sqrt{\boldsymbol{\rho}}. \end{eqnarray} \end{split}\]

With these definitions, we can deduce that:

(8)\[ \boldsymbol{\hat{T}} = \boldsymbol{\hat{K}} (\boldsymbol{I} - i\boldsymbol{\rho}\boldsymbol{\hat{K}})^{-1} = (\boldsymbol{I} - i\boldsymbol{\hat{K}}\boldsymbol{\rho})^{-1} \boldsymbol{\hat{K}}. \]

Production processes

As noted in the intro, the \(\boldsymbol{K}\)-matrix describes scattering processes of type \(cd \to ab\). It can however be generalized to describe production processes of type \(R \to ab\). Here, the amplitude is described by a final state \(F\)-vector of size \(n\), so the question is how to express \(F\) in terms of transition matrix \(\boldsymbol{T}\).

dot = """
digraph {
    rankdir=LR;
    node [shape=point, width=0];
    edge [arrowhead=none];
    "a" [shape=none, label="aᵢ"];
    "b" [shape=none, label="bᵢ"];
    "R" [shape=none, label="R"];
    "N" [shape=none, label=""];
    "R" -> "N";
    "N" -> "a";
    "N" -> "b";
    { rank=same "Na", "Nb" };
}
"""
graphviz.Source(dot)
../../_images/k-matrix_20_0.svg

One approach by [6] is to transform \(\boldsymbol{T}\) into \(F\) (and its relativistic form \(\hat{F}\)) through the production amplitude \(P\)-vector:

(9)\[\begin{split} \begin{eqnarray} F & = & \left(\boldsymbol{I}-i\boldsymbol{K}\right)^{-1}P \\ \hat{F} & = & \left(\boldsymbol{I}-i\boldsymbol{\hat{K}\boldsymbol{\rho}}\right)^{-1}\hat{P}, \end{eqnarray} \end{split}\]

where we can compute \(\boldsymbol{\hat{K}}\) from (7):

(10)\[ \hat{\boldsymbol{K}} = \sqrt{\boldsymbol{\rho}^{-1}} \boldsymbol{K} \sqrt{\boldsymbol{\rho}^{-1}}. \]

Another approach by [7] further approximates this by defining a \(Q\)-vector:

(11)\[ Q = \boldsymbol{K}^{-1}P \quad \mathrm{and} \quad \hat{Q} = \boldsymbol{\hat{K}}^{-1}\hat{P} \]

that is taken to be constant (just some ‘fitting’ parameters). The \(F\)-vector can then be expressed as:

(12)\[ F = \boldsymbol{T}Q \quad \mathrm{and} \quad \hat{F} = \boldsymbol{\hat{T}}\hat{Q} \]

Note that for all these vectors, we have:

(13)\[ F=\sqrt{\boldsymbol{\rho}}\hat{F},\quad P=\sqrt{\boldsymbol{\rho}}\hat{P},\quad\mathrm{and}\quad Q=\sqrt{\boldsymbol{\rho}^{-1}}\hat{Q}. \]

Pole parametrization

After all these matrix definitions, the final challenge is to choose a correct parametrization for the elements of \(\boldsymbol{K}\) and \(P\) that accurately describes the resonances we observe.3 There are several choices, but a common one is the following summation over the poles \(R\):4

(14)\[\begin{split} \begin{eqnarray} K_{ij} &=& \sum_R\frac{g_{R,i}g_{R,j}}{m_R^2-s} + c_{ij} \\ \hat{K}_{ij} &=& \sum_R \frac{g_{R,i}(s)g_{R,j}(s)}{\left(m_R^2-s\right)\sqrt{\rho_i\rho_j}} + \hat{c}_{ij} \end{eqnarray} \end{split}\]

with \(c_{ij}, \hat{c}_{ij}\) some optional background characterization and \(g_{R,i}\) the residue functions. The residue functions are often further expressed as:

(15)\[\begin{split} \begin{eqnarray} g_{R,i} &=& \gamma_{R,i}\sqrt{m_R\Gamma^0_{R,i}} \\ g_{R,i}(s) &=& \gamma_{R,i}\sqrt{m_R\Gamma_{R,i}(s)} \end{eqnarray} \end{split}\]

with \(\gamma_{R,i}\) some real constants and \(\Gamma^0_{R,i}\) the partial width of each pole. In the Lorentz-invariant form, the fixed width \(\Gamma^0\) is replaced by an “energy dependent” EnergyDependentWidth \(\Gamma(s)\).5 The width for each pole can be computed as \(\Gamma^0_R = \sum_i\Gamma^0_{R,i}\).

The production vector \(P\) is commonly parameterized as:6

(16)\[\begin{split} \begin{eqnarray} P_i &=& \sum_R \frac{\beta^0_R\,g_{R,i}(s)}{m_R^2-s} \\ \hat{P}_i &=& \sum_R \frac{\beta^0_R\,g_{R,i}(s)}{\left(m_R^2-s\right)\sqrt{\rho_i}} \\ &=& \sum_R \frac{ \beta^0_R\gamma_{R,i}m_R\Gamma^0_R B_{R,i}(q(s)) }{m_R^2-s} \end{eqnarray} \end{split}\]

with \(B_{R,i}(q(s))\) the centrifugal damping factor (see BlattWeisskopfSquared) for channel \(i\) and \(\beta_R^0\) some (generally complex) constants that describe the production information of the decaying state \(R\). Usually, these constants are rescaled just like the residue functions in (15):

(17)\[ \beta^0_R = \beta_R\sqrt{m_R\Gamma^0_R}. \]

Implementation

Non-relativistic K-matrix

A non-relativistic \(\boldsymbol{K}\)-matrix for an arbitrary number of channels and an arbitrary number of poles can be formulated with the NonRelativisticKMatrix.formulate() method:

n_poles = sp.Symbol("n_R", integer=True, positive=True)
k_matrix_nr = kmatrix.NonRelativisticKMatrix.formulate(
    n_poles=n_poles, n_channels=1
)
k_matrix_nr[0, 0]
\[\displaystyle \frac{\sum_{R=1}^{n_{R}} \frac{{\Gamma}_{R,0} {\gamma}_{R,0}^{2} {m}_{R}}{- s + {m}_{R}^{2}}}{- i \sum_{R=1}^{n_{R}} \frac{{\Gamma}_{R,0} {\gamma}_{R,0}^{2} {m}_{R}}{- s + {m}_{R}^{2}} + 1}\]

Notice how the \(\boldsymbol{K}\)-matrix reduces to a relativistic_breit_wigner() in the case of one channel and one pole (but for a residue constant \(\gamma\)):

k_matrix_1r = kmatrix.NonRelativisticKMatrix.formulate(n_poles=1, n_channels=1)
k_matrix_1r[0, 0].doit().simplify()
\[\displaystyle - \frac{{\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1}}{s + i {\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1} - {m}_{1}^{2}}\]

Now let’s investigate the effect of using a \(\boldsymbol{K}\)-matrix to describe two poles in one channel and see how it compares with the sum of two Breit-Wigner functions (two ‘resonances’). Two Breit-Wigner ‘poles’ with the same parameters would look like this:

s, m1, m2, Gamma1, Gamma2 = sp.symbols("s, m1, m2, Gamma1, Gamma2")
bw1 = relativistic_breit_wigner(s, m1, Gamma1)
bw2 = relativistic_breit_wigner(s, m2, Gamma2)
bw = bw1 + bw2
bw
\[\displaystyle \frac{\Gamma_{1} m_{1}}{- i \Gamma_{1} m_{1} + m_{1}^{2} - s} + \frac{\Gamma_{2} m_{2}}{- i \Gamma_{2} m_{2} + m_{2}^{2} - s}\]

while a \(\boldsymbol{K}\)-matrix parametrizes the two poles as:

k_matrix_2r = kmatrix.NonRelativisticKMatrix.formulate(n_poles=2, n_channels=1)
k_matrix = k_matrix_2r[0, 0].doit()
# reformulate terms
denominator, nominator = k_matrix.args
term1 = nominator.args[0] * denominator
term2 = nominator.args[1] * denominator
k_matrix = term1 + term2
k_matrix
\[\displaystyle \frac{{\Gamma}_{2,0} {\gamma}_{2,0}^{2} {m}_{2}}{\left(- s + {m}_{2}^{2}\right) \left(- i \left(\frac{{\Gamma}_{2,0} {\gamma}_{2,0}^{2} {m}_{2}}{- s + {m}_{2}^{2}} + \frac{{\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1}}{- s + {m}_{1}^{2}}\right) + 1\right)} + \frac{{\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1}}{\left(- s + {m}_{1}^{2}\right) \left(- i \left(\frac{{\Gamma}_{2,0} {\gamma}_{2,0}^{2} {m}_{2}}{- s + {m}_{2}^{2}} + \frac{{\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1}}{- s + {m}_{1}^{2}}\right) + 1\right)}\]

To simplify things, we can set the residue constants \(\gamma\) to one. Notice how the \(\boldsymbol{K}\)-matrix has introduced some coupling (‘interference’) between the two terms.

def remove_residue_constants(expression):
    expression = symplot.substitute_indexed_symbols(expression)
    residue_constants = filter(
        lambda s: re.match(r"^\\?gamma", s.name),
        expression.free_symbols,
    )
    return expression.xreplace({gamma: 1 for gamma in residue_constants})


display(
    remove_residue_constants(bw),
    remove_residue_constants(k_matrix),
)
\[\displaystyle \frac{\Gamma_{1} m_{1}}{- i \Gamma_{1} m_{1} + m_{1}^{2} - s} + \frac{\Gamma_{2} m_{2}}{- i \Gamma_{2} m_{2} + m_{2}^{2} - s}\]
\[\displaystyle \frac{\Gamma_{1,0} m_{1}}{\left(m_{1}^{2} - s\right) \left(- i \left(\frac{\Gamma_{1,0} m_{1}}{m_{1}^{2} - s} + \frac{\Gamma_{2,0} m_{2}}{m_{2}^{2} - s}\right) + 1\right)} + \frac{\Gamma_{2,0} m_{2}}{\left(m_{2}^{2} - s\right) \left(- i \left(\frac{\Gamma_{1,0} m_{1}}{m_{1}^{2} - s} + \frac{\Gamma_{2,0} m_{2}}{m_{2}^{2} - s}\right) + 1\right)}\]

Now, just like in Inspect model interactively, we use symplot to visualize the difference between the two expressions. The important thing is that the Argand plot on the right shows that the \(\boldsymbol{K}\)-matrix conserves unitarity.

Note that we have to call symplot.substitute_indexed_symbols() to turn the Indexed instances in this Matrix expression into Symbols before calling this function. We also call symplot.rename_symbols() so that the residue \(\gamma\)’s get a name that does not have to be dummified by lambdify().

# Prepare expressions
m = sp.Symbol("m")
k_matrix = symplot.substitute_indexed_symbols(k_matrix)
rename_gammas = lambda s: re.sub(  # noqa: E731
    r"\\([Gg])amma_{([0-9]),0}", r"\1amma\2", s
)
gamma1, gamma2 = sp.symbols("gamma1 gamma2")
bw = symplot.rename_symbols(bw, rename_gammas)
k_matrix = symplot.rename_symbols(k_matrix, rename_gammas)
bw = bw.xreplace({s: m**2})
k_matrix = k_matrix.xreplace({s: m**2})

# Prepare sliders and domain
np_kmatrix, sliders = symplot.prepare_sliders(k_matrix, m)
np_bw = sp.lambdify((m, Gamma1, Gamma2, gamma1, gamma2, m1, m2), bw.doit())

m_min, m_max = 0, 3
domain_1d = np.linspace(m_min, m_max, 200)
domain_argand = np.linspace(m_min - 2, m_max + 2, 1_000)
sliders.set_ranges(
    m1=(0, 3, 100),
    m2=(0, 3, 100),
    Gamma1=(0, 2, 100),
    Gamma2=(0, 2, 100),
    gamma1=(0, 2),
    gamma2=(0, 2),
)
sliders.set_values(
    m1=1.1,
    m2=1.9,
    Gamma1=0.2,
    Gamma2=0.3,
    gamma1=1,
    gamma2=1,
)
if STATIC_WEB_PAGE:
    # Concatenate flipped domain for reverse animation
    domain = np.linspace(1.0, 2.7, 50)
    domain = np.concatenate((domain, np.flip(domain[1:])))
    sliders._sliders["m1"] = domain


def create_argand(func):
    def wrapped(**kwargs):
        values = func(domain_argand, **kwargs)
        argand = np.array([values.real, values.imag])
        return argand.T

    return wrapped


# Create figure
fig, axes = plt.subplots(
    ncols=2,
    figsize=1.2 * np.array((8, 3.8)),
    tight_layout=True,
)
ax_intensity, ax_argand = axes
m_label = "$m_{a+b}$"
ax_intensity.set_xlabel(m_label)
ax_intensity.set_ylabel("$|A|^2$")
ax_argand.set_xlabel("Re($A$)")
ax_argand.set_ylabel("Im($A$)")

# Plot intensity
controls = iplt.plot(
    domain_1d,
    lambda *args, **kwargs: np.abs(np_kmatrix(*args, **kwargs) ** 2),
    label="$K$-matrix",
    **sliders,
    ylim="auto",
    ax=ax_intensity,
)
iplt.plot(
    domain_1d,
    lambda *args, **kwargs: np.abs(np_bw(*args, **kwargs) ** 2),
    label="Breit-Wigner",
    controls=controls,
    ylim="auto",
    ax=ax_intensity,
)
plt.legend(loc="upper right")
iplt.axvline(controls["m1"], c="gray", linestyle="dotted")
iplt.axvline(controls["m2"], c="gray", linestyle="dotted")

# Argand plots
iplt.scatter(
    create_argand(np_kmatrix),
    label="$K$-matrix",
    controls=controls,
    parametric=True,
    s=1,
    ax=ax_argand,
)
iplt.scatter(
    create_argand(np_bw),
    label="Breit-Wigner",
    controls=controls,
    parametric=True,
    s=1,
    ax=ax_argand,
)
plt.legend(loc="upper right");
../../_images/k-matrix_41_0.png

Relativistic K-matrix

Relativistic \(\boldsymbol{K}\)-matrices for an arbitrary number of channels and an arbitrary number of poles can be formulated with the RelativisticKMatrix.formulate() method:

L = sp.Symbol("L", integer=True, negative=False)
n_poles = sp.Symbol("n_R", integer=True, positive=True)
rel_k_matrix_nr = kmatrix.RelativisticKMatrix.formulate(
    n_poles=n_poles, n_channels=1, angular_momentum=L
)
rel_k_matrix_nr[0, 0]
\[\displaystyle \frac{\overline{\sqrt{\rho\left(s\right)}} \sqrt{\rho\left(s\right)} \sum_{R=1}^{n_{R}} \frac{\Gamma\left(s\right) {\gamma}_{R,0}^{2} {m}_{R}}{- s + {m}_{R}^{2}}}{- i \rho\left(s\right) \sum_{R=1}^{n_{R}} \frac{\Gamma\left(s\right) {\gamma}_{R,0}^{2} {m}_{R}}{- s + {m}_{R}^{2}} + 1}\]

Again, as in Non-relativistic K-matrix, the \(\boldsymbol{K}\)-matrix reduces to something of a relativistic_breit_wigner(). This time, the width has been replaced by a EnergyDependentWidth and some PhaseSpaceFactors have been inserted that take care of the decay into two decay products:

rel_k_matrix_1r = kmatrix.RelativisticKMatrix.formulate(
    n_poles=1, n_channels=1, angular_momentum=L
)
symplot.partial_doit(rel_k_matrix_1r[0, 0], sp.Sum).simplify(doit=False)
\[\displaystyle - \frac{\overline{\sqrt{\rho\left(s\right)}} \Gamma_{1,0}\left(s\right) {\gamma}_{1,0}^{2} {m}_{1} \sqrt{\rho\left(s\right)}}{s + i \Gamma_{1,0}\left(s\right) {\gamma}_{1,0}^{2} {m}_{1} \rho\left(s\right) - {m}_{1}^{2}}\]

Note that another difference with relativistic_breit_wigner_with_ff() is an additional phase space factor in the denominator. That one disappears in P-vector.

The \(\boldsymbol{K}\)-matrix with two poles becomes (neglecting the \(\sqrt{\rho_0(s)}\)):

rel_k_matrix_2r = kmatrix.RelativisticKMatrix.formulate(
    n_poles=2, n_channels=1, angular_momentum=L
)
rel_k_matrix_2r = symplot.partial_doit(rel_k_matrix_2r[0, 0], sp.Sum)
rel_k_matrix_2r = symplot.substitute_indexed_symbols(rel_k_matrix_2r)
s, m_a, m_b = sp.symbols("s, m_a0, m_b0")
rho = PhaseSpaceFactor(s, m_a, m_b)
rel_k_matrix_2r = rel_k_matrix_2r.xreplace(
    {
        sp.sqrt(rho): 1,
        sp.conjugate(sp.sqrt(rho)): 1,
    }
)
denominator, nominator = rel_k_matrix_2r.args
term1 = nominator.args[0] * denominator
term2 = nominator.args[1] * denominator
rel_k_matrix_2r = term1 + term2
rel_k_matrix_2r
\[\displaystyle \frac{\gamma_{1,0}^{2} m_{1} \Gamma_{1,0}\left(s\right)}{\left(m_{1}^{2} - s\right) \left(- i \left(\frac{\gamma_{1,0}^{2} m_{1} \Gamma_{1,0}\left(s\right)}{m_{1}^{2} - s} + \frac{\gamma_{2,0}^{2} m_{2} \Gamma_{2,0}\left(s\right)}{m_{2}^{2} - s}\right) \rho\left(s\right) + 1\right)} + \frac{\gamma_{2,0}^{2} m_{2} \Gamma_{2,0}\left(s\right)}{\left(m_{2}^{2} - s\right) \left(- i \left(\frac{\gamma_{1,0}^{2} m_{1} \Gamma_{1,0}\left(s\right)}{m_{1}^{2} - s} + \frac{\gamma_{2,0}^{2} m_{2} \Gamma_{2,0}\left(s\right)}{m_{2}^{2} - s}\right) \rho\left(s\right) + 1\right)}\]

This again shows the interference introduced by the \(\boldsymbol{K}\)-matrix, when compared with a sum of two Breit-Wigner functions.

P-vector

For one channel and an arbitrary number of poles \(n_R\), the \(F\)-vector gets the following form:

n_poles = sp.Symbol("n_R", integer=True, positive=True)
kmatrix.NonRelativisticPVector.formulate(n_poles=n_poles, n_channels=1)[0]
\[\displaystyle \frac{\sum_{R=1}^{n_{R}} \frac{{\Gamma}_{R,0} {\beta}_{R} {\gamma}_{R,0} {m}_{R}}{- s + {m}_{R}^{2}}}{- i \sum_{R=1}^{n_{R}} \frac{{\Gamma}_{R,0} {\gamma}_{R,0}^{2} {m}_{R}}{- s + {m}_{R}^{2}} + 1}\]

The RelativisticPVector looks like:

kmatrix.RelativisticPVector.formulate(
    n_poles=n_poles, n_channels=1, angular_momentum=L
)[0]
\[\displaystyle \frac{\overline{\sqrt{\rho\left(s\right)}} \sqrt{\rho\left(s\right)} \sum_{R=1}^{n_{R}} \frac{\sqrt{B_{L}^2\left(q^2\left(s\right)\right)} {\Gamma}_{R,0} {\beta}_{R} {\gamma}_{R,0} {m}_{R}}{- s + {m}_{R}^{2}}}{\overline{\sqrt{\rho\left(s\right)}} - i \sqrt{\rho\left(s\right)} \sum_{R=1}^{n_{R}} \frac{\Gamma\left(s\right) {\gamma}_{R,0}^{2} {m}_{R}}{- s + {m}_{R}^{2}}}\]

As in Non-relativistic K-matrix, if we take \(n_R=1\), the \(F\)-vector reduces to a Breit-Wigner function, but now with an additional factor \(\beta\).

f_vector_1r = kmatrix.NonRelativisticPVector.formulate(n_poles=1, n_channels=1)
symplot.partial_doit(f_vector_1r[0], sp.Sum)
\[\displaystyle \frac{{\Gamma}_{1,0} {\beta}_{1} {\gamma}_{1,0} {m}_{1}}{\left(1 - \frac{i {\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1}}{- s + {m}_{1}^{2}}\right) \left(- s + {m}_{1}^{2}\right)}\]

And when we neglect the phase space factors \(\sqrt{\rho_0(s)}\), the RelativisticPVector reduces to a relativistic_breit_wigner_with_ff()!

rel_f_vector_1r = kmatrix.RelativisticPVector.formulate(
    n_poles=1, n_channels=1, angular_momentum=L
)
rel_f_vector_1r = symplot.partial_doit(rel_f_vector_1r[0], sp.Sum)
rel_f_vector_1r = symplot.substitute_indexed_symbols(rel_f_vector_1r)
s, m_a, m_b = sp.symbols("s, m_a0, m_b0")
rho = PhaseSpaceFactor(s, m_a, m_b)
rel_f_vector_1r.xreplace(
    {
        sp.sqrt(rho): 1,
        sp.conjugate(sp.sqrt(rho)): 1,
    }
).simplify(doit=False)
\[\displaystyle - \frac{\Gamma_{1,0} \gamma_{1,0} \beta_{1} m_{1} \sqrt{B_{L}^2\left(q^2\left(s\right)\right)}}{i \gamma_{1,0}^{2} m_{1} \Gamma_{1,0}\left(s\right) - m_{1}^{2} + s}\]

Note that the \(F\)-vector approach introduces additional \(\beta\)-coefficients. These can constants can be complex and can introduce phase differences form the production process.

f_vector_2r = kmatrix.NonRelativisticPVector.formulate(n_poles=2, n_channels=1)
f_vector = f_vector_2r[0].doit()
denominator, nominator = f_vector.args
term1 = nominator.args[0] * denominator
term2 = nominator.args[1] * denominator
f_vector = term1 + term2
f_vector
\[\displaystyle \frac{{\Gamma}_{2,0} {\beta}_{2} {\gamma}_{2,0} {m}_{2}}{\left(- s + {m}_{2}^{2}\right) \left(- i \left(\frac{{\Gamma}_{2,0} {\gamma}_{2,0}^{2} {m}_{2}}{- s + {m}_{2}^{2}} + \frac{{\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1}}{- s + {m}_{1}^{2}}\right) + 1\right)} + \frac{{\Gamma}_{1,0} {\beta}_{1} {\gamma}_{1,0} {m}_{1}}{\left(- s + {m}_{1}^{2}\right) \left(- i \left(\frac{{\Gamma}_{2,0} {\gamma}_{2,0}^{2} {m}_{2}}{- s + {m}_{2}^{2}} + \frac{{\Gamma}_{1,0} {\gamma}_{1,0}^{2} {m}_{1}}{- s + {m}_{1}^{2}}\right) + 1\right)}\]

Now again let’s compare the compare this with a sum of two relativistic_breit_wigner()s, now with the two additional \(\beta\)-constants.

beta1, beta2 = sp.symbols("beta1 beta2")
bw_with_phases = beta1 * bw1 + beta2 * bw2
display(
    bw_with_phases,
    remove_residue_constants(f_vector),
)
\[\displaystyle \frac{\Gamma_{1} \beta_{1} m_{1}}{- i \Gamma_{1} m_{1} + m_{1}^{2} - s} + \frac{\Gamma_{2} \beta_{2} m_{2}}{- i \Gamma_{2} m_{2} + m_{2}^{2} - s}\]
\[\displaystyle \frac{\Gamma_{1,0} \beta_{1} m_{1}}{\left(m_{1}^{2} - s\right) \left(- i \left(\frac{\Gamma_{1,0} m_{1}}{m_{1}^{2} - s} + \frac{\Gamma_{2,0} m_{2}}{m_{2}^{2} - s}\right) + 1\right)} + \frac{\Gamma_{2,0} \beta_{2} m_{2}}{\left(m_{2}^{2} - s\right) \left(- i \left(\frac{\Gamma_{1,0} m_{1}}{m_{1}^{2} - s} + \frac{\Gamma_{2,0} m_{2}}{m_{2}^{2} - s}\right) + 1\right)}\]
# Prepare expressions
f_vector = symplot.substitute_indexed_symbols(f_vector)
rename_gammas = lambda s: re.sub(  # noqa: E731
    r"\\([Gg])amma_{([0-9]),0}", r"\1amma\2", s
)
c1, c2, phi1, phi2 = sp.symbols("c1 c2 phi1 phi2", real=True)
bw_with_phases = symplot.rename_symbols(bw_with_phases, rename_gammas)
f_vector = symplot.rename_symbols(f_vector, rename_gammas)
substitutions = {
    s: m**2,
    beta1: c1 * sp.exp(sp.I * phi1),
    beta2: c2 * sp.exp(sp.I * phi2),
}
bw_with_phases = bw_with_phases.xreplace(substitutions)
f_vector = f_vector.xreplace(substitutions)

# Prepare sliders and domain
np_f_vector, sliders = symplot.prepare_sliders(f_vector, m)
np_bw = sp.lambdify(
    (m, Gamma1, Gamma2, c1, c2, gamma1, gamma2, m1, m2, phi1, phi2),
    bw_with_phases.doit(),
)

# Set plot domain
x_min, x_max = 0, 3
y_min, y_max = -0.5, +0.5
plot_domain = np.linspace(x_min, x_max, num=150)
plot_domain_argand = np.linspace(x_min - 2, x_max + 2, num=400)
x_values = np.linspace(x_min, x_max, num=160)
y_values = np.linspace(y_min, y_max, num=80)
X, Y = np.meshgrid(x_values, y_values)
plot_domain_complex = X + Y * 1j

# Set slider values and ranges
sliders.set_ranges(
    c1=(0, 2, 100),
    c2=(0, 2, 100),
    phi1=(0, np.pi, np.pi / 20),
    phi2=(0, np.pi, np.pi / 20),
    m1=(0, 3, 300),
    m2=(0, 3, 300),
    Gamma1=(-2, 2, 100),
    Gamma2=(-2, 2, 100),
    gamma1=(0, 2),
    gamma2=(0, 2),
)
sliders.set_values(
    c1=1,
    c2=1,
    m1=1.4,
    m2=1.7,
    Gamma1=0.2,
    Gamma2=0.3,
    gamma1=1,
    gamma2=1,
)


def create_argand(func):
    def wrapped(**kwargs):
        values = func(plot_domain_argand, **kwargs)
        argand = np.array([values.real, values.imag])
        return argand.T

    return wrapped


# Create figure
fig, axes = plt.subplots(
    nrows=3,
    ncols=2,
    figsize=(10, 9),
    gridspec_kw=dict(
        width_ratios=[2.5, 1],
        wspace=0.08,
        left=0.05,
        right=0.99,
        top=0.99,
        bottom=0.05,
    ),
)
fig.canvas.toolbar_visible = False
fig.canvas.header_visible = False
fig.canvas.footer_visible = False
(ax_1d, ax_argand), (ax_2d, ax_empty1), (ax_2d_bw, ax_empty2) = axes
ax_empty1.axis("off")
ax_empty2.axis("off")
for ax in axes.flatten():
    ax.set_yticks([])
ax_argand.set_xticks([])
ax_1d.axes.get_xaxis().set_visible(False)
ax_2d.axes.get_xaxis().set_visible(False)
ax_1d.sharex(ax_2d_bw)
ax_2d.sharex(ax_2d_bw)
ax_2d_bw.set_xlabel("Re $m$")
for ax in (ax_2d, ax_2d_bw):
    ax.set_ylabel("Im $m$")
ax_1d.set_ylabel("$|A|^2$")
ax_argand.set_xlabel("Re($A$)")
ax_argand.set_ylabel("Im($A$)")

for ax in (ax_2d, ax_2d_bw):
    ax.axhline(0, linewidth=0.5, c="black", linestyle="dotted")

# 1D intensity plot
controls = Controls(**sliders)
controls = iplt.plot(
    plot_domain,
    lambda *args, **kwargs: np.abs(np_f_vector(*args, **kwargs) ** 2),
    label="$P$-vector",
    controls=controls,
    ylim="auto",
    ax=ax_1d,
)
iplt.plot(
    plot_domain,
    lambda *args, **kwargs: np.abs(np_bw(*args, **kwargs) ** 2),
    label="Breit-Wigner",
    controls=controls,
    ylim="auto",
    ax=ax_1d,
)
ax_1d.legend(loc="upper right")
mass_line_style = dict(
    c="red",
    alpha=0.3,
)
for name in controls.params:
    if not re.match(r"^m[0-9]+$", name):
        continue
    iplt.axvline(controls[name], ax=ax_1d, **mass_line_style)

# Argand plot
iplt.scatter(
    create_argand(np_f_vector),
    label="$P$-vector",
    controls=controls,
    parametric=True,
    s=1,
    ax=ax_argand,
    xlim="auto",
    ylim="auto",
)
iplt.scatter(
    create_argand(np_bw),
    label="Breit-Wigner",
    controls=controls,
    parametric=True,
    s=1,
    ax=ax_argand,
    xlim="auto",
    ylim="auto",
)

# 3D plot
color_mesh = None
color_mesh_bw = None
pole_indicators = []
pole_indicators_bw = []


def plot3(*, z_cutoff, complex_rendering, **kwargs):
    global color_mesh, color_mesh_bw
    Z = np_f_vector(plot_domain_complex, **kwargs)
    Z_bw = np_bw(plot_domain_complex, **kwargs)
    if complex_rendering == "imag":
        projection = np.imag
        ax_title = "Im $A$"
    elif complex_rendering == "real":
        projection = np.real
        ax_title = "Re $A$"
    elif complex_rendering == "abs":
        projection = np.vectorize(lambda z: np.abs(z) ** 2)
        ax_title = "$|A|$"
    else:
        raise NotImplementedError
    ax_2d.set_title(ax_title + " ($P$-vector)")
    ax_2d_bw.set_title(ax_title + " (Breit-Wigner)")

    if color_mesh is None:
        color_mesh = ax_2d.pcolormesh(X, Y, projection(Z), cmap=cm.coolwarm)
    else:
        color_mesh.set_array(projection(Z)[:-1, :-1])
    color_mesh.set_clim(vmin=-z_cutoff, vmax=+z_cutoff)

    if color_mesh_bw is None:
        color_mesh_bw = ax_2d_bw.pcolormesh(
            X, Y, projection(Z_bw), cmap=cm.coolwarm
        )
    else:
        color_mesh_bw.set_array(projection(Z_bw)[:-1, :-1])
    color_mesh_bw.set_clim(vmin=-z_cutoff, vmax=+z_cutoff)

    if pole_indicators:
        for R, (line, text) in enumerate(pole_indicators, 1):
            mass = kwargs[f"m{R}"]
            line.set_xdata(mass)
            text.set_x(mass + (x_max - x_min) * 0.008)
    else:
        for R in range(1, 2 + 1):
            mass = kwargs[f"m{R}"]
            line = ax_2d.axvline(mass, **mass_line_style)
            text = ax_2d.text(
                x=mass + (x_max - x_min) * 0.008,
                y=0.95 * y_min,
                s=f"$m_{R}$",
                c="red",
            )
            pole_indicators.append((line, text))

    if pole_indicators_bw:
        for R, (line, text) in enumerate(pole_indicators_bw, 1):
            mass = kwargs[f"m{R}"]
            line.set_xdata(mass)
            text.set_x(mass + (x_max - x_min) * 0.008)
    else:
        for R in range(1, 2 + 1):
            mass = kwargs[f"m{R}"]
            line = ax_2d_bw.axvline(mass, **mass_line_style)
            text = ax_2d_bw.text(
                x=mass + (x_max - x_min) * 0.008,
                y=0.95 * y_min,
                s=f"$m_{R}$",
                c="red",
            )
            pole_indicators_bw.append((line, text))


# Create switch for imag/real/abs
name = "complex_rendering"
sliders._sliders[name] = ipywidgets.RadioButtons(
    options=["imag", "real", "abs"],
    description=R"\(s\)-plane plot",
)
sliders._arg_to_symbol[name] = name

# Create cut-off slider for z-direction
name = "z_cutoff"
sliders._sliders[name] = ipywidgets.FloatSlider(
    value=2,
    min=0,
    max=+5,
    step=0.1,
    description=R"\(z\)-cutoff",
)
sliders._arg_to_symbol[name] = name

# Create GUI
sliders_copy = dict(sliders)
slider_groups = []
for R in range(1, 2 + 1):
    vertical_slider_group = [
        sliders_copy.pop(f"c{R}"),
        sliders_copy.pop(f"phi{R}"),
        sliders_copy.pop(f"m{R}"),
        sliders_copy.pop(f"Gamma{R}"),
        sliders_copy.pop(f"gamma{R}"),
    ]
    slider_groups.append(vertical_slider_group)
slider_pairs = np.array(slider_groups).T
h_boxes = [ipywidgets.HBox(tuple(pair)) for pair in slider_pairs]
remaining_sliders = sorted(
    sliders_copy.values(),
    key=lambda s: (str(type(s)), s.description),
)
ui = ipywidgets.VBox(h_boxes + remaining_sliders)
output = ipywidgets.interactive_output(plot3, controls=sliders)
display(ui, output)
../../_images/k-matrix_68_0.png

Interactive visualization

All \(\boldsymbol{K}\)-matrices can be inspected interactively for arbitrary poles and channels with the following applet:

if STATIC_WEB_PAGE:
    L = 0


def plot(
    kmatrix_type: kmatrix.TMatrix,
    n_channels: int,
    n_poles: int,
    angular_momentum: sp.Symbol = 0,
    phsp_factor=PhaseSpaceFactor,
    substitute_sqrt_rho: bool = False,
) -> None:
    # Convert to Symbol: symplot cannot handle IndexedBase
    i, j = sp.symbols("i, j", integer=True, negative=False)
    j = i
    expr = kmatrix_type.formulate(
        n_poles=n_poles,
        n_channels=n_channels,
        angular_momentum=angular_momentum,
        phsp_factor=phsp_factor,
    ).doit()[i, j]
    expr = symplot.substitute_indexed_symbols(expr)
    if substitute_sqrt_rho:
        rho_i = lambda i: phsp_factor(  # noqa
            *sp.symbols(f"s m_a{i} m_b{i}")
        ).doit()
        expr = expr.xreplace(
            {sp.sqrt(rho_i(i)): 1 for i in range(n_channels)}
        ).xreplace(
            {sp.conjugate(sp.sqrt(rho_i(i))): 1 for i in range(n_channels)}
        )
    expr = expr.xreplace({s: m**2})
    expr = symplot.substitute_indexed_symbols(expr)
    np_expr, sliders = symplot.prepare_sliders(expr, m)

    # Set plot domain
    x_min, x_max = 1e-3, 3
    y_min, y_max = -0.5, +0.5

    plot_domain = np.linspace(x_min, x_max, num=500, dtype=np.complex)
    x_values = np.linspace(x_min, x_max, num=160)
    y_values = np.linspace(y_min, y_max, num=80)
    X, Y = np.meshgrid(x_values, y_values)
    plot_domain_complex = X + Y * 1j

    # Set slider values and ranges
    m0_values = np.linspace(x_min, x_max, num=n_poles + 2)
    m0_values = m0_values[1:-1]
    max_angular_momentum = BlattWeisskopfSquared.max_angular_momentum
    if max_angular_momentum is None:
        max_angular_momentum = 8
    if "L" in sliders:
        sliders.set_ranges(L=(0, max_angular_momentum))
    sliders.set_ranges(i=(0, n_channels - 1))
    for R in range(1, n_poles + 1):
        for i in range(n_channels):
            if kmatrix_type in {
                kmatrix.RelativisticKMatrix,
                kmatrix.RelativisticPVector,
            }:
                sliders.set_ranges(
                    {
                        f"m{R}": (0, 3, 100),
                        Rf"\Gamma_{{{R},{i}}}": (-2, +2, 100),
                        Rf"\gamma_{{{R},{i}}}": (0, 10, 100),
                        f"m_a{i}": (0, 1, 0.01),
                        f"m_b{i}": (0, 1, 0.01),
                    }
                )
                sliders.set_values(
                    {
                        f"m{R}": m0_values[R - 1],
                        Rf"\Gamma_{{{R},{i}}}": 2.0
                        * (0.4 + R * 0.2 - i * 0.3),
                        Rf"\gamma_{{{R},{i}}}": 0.25 * (10 - R + i),
                        f"m_a{i}": (i + 1) * 0.25,
                        f"m_b{i}": (i + 1) * 0.25,
                    }
                )
            if kmatrix_type in {
                kmatrix.NonRelativisticPVector,
                kmatrix.NonRelativisticKMatrix,
            }:
                sliders.set_ranges(
                    {
                        f"m{R}": (0, 3, 100),
                        Rf"\Gamma_{{{R},{i}}}": (-1, 1, 100),
                        Rf"\gamma_{{{R},{i}}}": (0, 2, 100),
                    }
                )
                sliders.set_values(
                    {
                        f"m{R}": m0_values[R - 1],
                        Rf"\Gamma_{{{R},{i}}}": (R + 1) * 0.1,
                        Rf"\gamma_{{{R},{i}}}": 1 - 0.1 * R + 0.1 * i,
                    }
                )
            if kmatrix_type in {
                kmatrix.NonRelativisticPVector,
                kmatrix.RelativisticPVector,
            }:
                sliders.set_ranges({f"beta{R}": (-1, 1, 100)})
                sliders.set_values({f"beta{R}": 1})

    # Create interactive plots
    fig, axes = plt.subplots(
        nrows=2,
        figsize=(8, 6),
        sharex=True,
        tight_layout=True,
    )
    fig.canvas.toolbar_visible = False
    fig.canvas.header_visible = False
    fig.canvas.footer_visible = False
    if kmatrix_type in {
        kmatrix.NonRelativisticKMatrix,
        kmatrix.RelativisticKMatrix,
    }:
        fig.suptitle(
            Rf"${n_channels} \times {n_channels}$ $K$-matrix"
            f" with {n_poles} resonances"
        )
    elif kmatrix_type in {
        kmatrix.NonRelativisticPVector,
        kmatrix.RelativisticPVector,
    }:
        fig.suptitle(
            f"$P$-vector for {n_channels} channels and {n_poles} resonances"
        )

    for ax in axes:
        ax.set_xlim(x_min, x_max)
    ax_2d, ax_3d = axes
    ax_2d.set_ylabel("$|T|^{2}$")
    ax_2d.set_yticks([])
    ax_3d.set_xlabel("Re $m$")
    ax_3d.set_ylabel("Im $m$")
    ax_3d.set_xticks([])
    ax_3d.set_yticks([])

    ax_3d.axhline(0, linewidth=0.5, c="black", linestyle="dotted")

    # 2D plot
    def plot(channel: int):
        def wrapped(*args, **kwargs) -> sp.Expr:
            kwargs["i"] = channel
            return np.abs(np_expr(*args, **kwargs)) ** 2

        return wrapped

    controls = Controls(**sliders)
    for i in range(n_channels):
        iplt.plot(
            plot_domain,
            plot(i),
            ax=axes[0],
            controls=controls,
            ylim="auto",
            label=f"channel {i}",
        )
    if n_channels > 1:
        axes[0].legend(loc="upper right")
    mass_line_style = dict(
        c="red",
        alpha=0.3,
    )
    for name in controls.params:
        if not re.match(r"^m[0-9]+$", name):
            continue
        iplt.axvline(controls[name], ax=axes[0], **mass_line_style)

    # 3D plot
    color_mesh = None
    pole_indicators = []
    threshold_indicators = []

    def plot3(*, z_cutoff, complex_rendering, **kwargs):
        nonlocal color_mesh
        Z = np_expr(plot_domain_complex, **kwargs)
        if complex_rendering == "imag":
            Z_values = Z.imag
            ax_title = "Im $T$"
        elif complex_rendering == "real":
            Z_values = Z.real
            ax_title = "Re $T$"
        elif complex_rendering == "abs":
            Z_values = np.abs(Z)
            ax_title = "$|T|$"
        else:
            raise NotImplementedError

        if n_channels == 1:
            axes[-1].set_title(ax_title)
        else:
            i = kwargs["i"]
            axes[-1].set_title(f"{ax_title}, channel {i}")

        if color_mesh is None:
            color_mesh = ax_3d.pcolormesh(X, Y, Z_values, cmap=cm.coolwarm)
        else:
            color_mesh.set_array(Z_values[:-1, :-1])
        color_mesh.set_clim(vmin=-z_cutoff, vmax=+z_cutoff)

        if pole_indicators:
            for R, (line, text) in enumerate(pole_indicators, 1):
                mass = kwargs[f"m{R}"]
                line.set_xdata(mass)
                text.set_x(mass + (x_max - x_min) * 0.008)
        else:
            for R in range(1, n_poles + 1):
                mass = kwargs[f"m{R}"]
                line = ax_3d.axvline(mass, **mass_line_style)
                text = ax_3d.text(
                    x=mass + (x_max - x_min) * 0.008,
                    y=0.95 * y_min,
                    s=f"$m_{R}$",
                    c="red",
                )
                pole_indicators.append((line, text))

        if kmatrix_type is kmatrix.RelativisticKMatrix:
            x_offset = (x_max - x_min) * 0.015
            if threshold_indicators:
                for i, (line_thr, line_diff, text_thr, text_diff) in enumerate(
                    threshold_indicators
                ):
                    m_a = kwargs[f"m_a{i}"]
                    m_b = kwargs[f"m_b{i}"]
                    s_thr = m_a + m_b
                    m_diff = abs(m_a - m_b)
                    line_thr.set_xdata(s_thr)
                    line_diff.set_xdata(m_diff)
                    text_thr.set_x(s_thr)
                    text_diff.set_x(m_diff - x_offset)
            else:
                colors = cm.plasma(np.linspace(0, 1, n_channels))
                for i, color in enumerate(colors):
                    m_a = kwargs[f"m_a{i}"]
                    m_b = kwargs[f"m_b{i}"]
                    s_thr = m_a + m_b
                    m_diff = abs(m_a - m_b)
                    line_thr = ax.axvline(s_thr, c=color, linestyle="dotted")
                    line_diff = ax.axvline(m_diff, c=color, linestyle="dashed")
                    text_thr = ax.text(
                        x=s_thr,
                        y=0.95 * y_min,
                        s=f"$m_{{a{i}}}+m_{{b{i}}}$",
                        c=color,
                        rotation=-90,
                    )
                    text_diff = ax.text(
                        x=m_diff - x_offset,
                        y=0.95 * y_min,
                        s=f"$m_{{a{i}}}-m_{{b{i}}}$",
                        c=color,
                        rotation=+90,
                    )
                    threshold_indicators.append(
                        (line_thr, line_diff, text_thr, text_diff)
                    )
            for i, (_, line_diff, _, text_diff) in enumerate(
                threshold_indicators
            ):
                m_a = kwargs[f"m_a{i}"]
                m_b = kwargs[f"m_b{i}"]
                s_thr = m_a + m_b
                m_diff = abs(m_a - m_b)
                if m_diff > x_offset + 0.01 and s_thr - abs(m_diff) > x_offset:
                    line_diff.set_alpha(0.5)
                    text_diff.set_alpha(0.5)
                else:
                    line_diff.set_alpha(0)
                    text_diff.set_alpha(0)

    # Create switch for imag/real/abs
    name = "complex_rendering"
    sliders._sliders[name] = ipywidgets.RadioButtons(
        options=["imag", "real", "abs"],
        description=R"\(s\)-plane plot",
    )
    sliders._arg_to_symbol[name] = name

    # Create cut-off slider for z-direction
    name = "z_cutoff"
    sliders._sliders[name] = ipywidgets.FloatSlider(
        value=1,
        min=+0.01,
        max=+5,
        step=0.01,
        description=R"\(z\)-cutoff",
    )
    sliders._arg_to_symbol[name] = name

    # Link sliders
    if kmatrix_type is kmatrix.RelativisticKMatrix:
        for i in range(n_channels):
            ipywidgets.dlink(
                (sliders[f"m_a{i}"], "value"),
                (sliders[f"m_b{i}"], "value"),
            )

    # Create GUI
    sliders_copy = dict(sliders)
    h_boxes = []
    for R in range(1, n_poles + 1):
        buttons = [sliders_copy.pop(f"m{R}")]
        if n_channels == 1:
            buttons += [
                sliders_copy.pop(sliders.symbol_to_arg[Rf"\Gamma_{{{R},0}}"]),
                sliders_copy.pop(sliders.symbol_to_arg[Rf"\gamma_{{{R},0}}"]),
            ]
        h_box = ipywidgets.HBox(buttons)
        h_boxes.append(h_box)
    remaining_sliders = sorted(
        sliders_copy.values(), key=lambda s: (str(type(s)), s.description)
    )
    if n_channels == 1:
        remaining_sliders.remove(sliders["i"])
    ui = ipywidgets.VBox(h_boxes + remaining_sliders)
    output = ipywidgets.interactive_output(plot3, controls=sliders)
    display(ui, output)
plot(
    kmatrix.RelativisticKMatrix,
    n_poles=2,
    n_channels=1,
    angular_momentum=L,
    phsp_factor=PhaseSpaceFactor,
    substitute_sqrt_rho=False,
)
../../_images/k-matrix_75_0.png

1

Further subtleties arise when taking spin into account, especially for sequential decays. This is where spin formalisms come in.

2

An unpublished primer on the \(\boldsymbol{K}\)-matrix by Chung [8] uses a conjugate transpose of \(\boldsymbol{\rho}\), e.g. \(\boldsymbol{T} = \sqrt{\boldsymbol{\rho^\dagger}} \; \boldsymbol{\hat{T}} \sqrt{\boldsymbol{\rho}}\). This should not matter above threshold, where the phase space factor is real, but could have effects below threshold. This is where things become tricky: on the one hand, we need to ensure that \(\boldsymbol{K}\) remains real (unitarity) and on the other, we need to take keep track of which imaginary square root we choose (Riemann sheet). The latter is often called the requirement of analyticity. This is currently being explored in [TR-003] Chew-Mandelstam and [TR-004] Analyticity (WIP).

3

See PDG2020, §Resonances, p.5 and Section 49.1 for a discussion about what poles and resonances are. See also the intro to Section 5 in [4].

4

Eqs. (51) and (52) in [8] take a complex conjugate of one of the residue functions and one of the phase space factors.

5(1,2)

Unlike Eq. (77) in [4], AmpForm defines EnergyDependentWidth as in PDG2020, §Resonances, p.6, Eq. (49.21). The difference is that the phase space factor denoted by \(\rho_i\) in Eq. (77) in [4] is divided by the phase space factor at the pole position \(m_R\). So in AmpForm, the choice is \(\rho_i \to \frac{\rho_i(s)}{\rho_i(m_R)}\).

6

Just as with 5, we have smuggled a bit in the last equation in order to be able to reproduce Eq. (49.19) in PDG2020, §Resonances, p.6 in the case \(n=1,n_R=1\), on which relativistic_breit_wigner_with_ff() is based.