K-matrix

K-matrix#

This report investigates how to implement \(K\)-matrix dynamics with SymPy. We here describe only the version that is not Lorentz-invariant, because it is simplest and allows us to check whether the case \(n_R=1, n=1\) (single resonance, single channel) reduces to a Breit-Wigner function. We followed the physics as described by PDG2020, §Resonances and [Chung et al., 1995, Peters, 2004, Meyer, 2008]. For the Lorentz-invariant version, see TR-009.

A brief overview of the origin of the \(\boldsymbol{K}\)-matrix is given first. This overview follows [Chung et al., 1995], but skips over quite a few details, as this is only an attempt to provide some context of what is going on.

Physics#

The \(\boldsymbol{K}\)-matrix formalism is used to describe coupled, two-body scattering processes of the type \(c_id_i \to R \to a_ib_i\), with \(i\) representing each separate channel and \(R\) a number of resonances that these channels have in common.

https://user-images.githubusercontent.com/29308176/164994485-fc4843c3-856b-4853-857a-679e258cf7c8.svg

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\) as:

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

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

(2)#\[ A(\Omega) = \frac{1}{2q_i}\sum_J\left(2J+1\right) T^J(s) {D^J_{\lambda\mu}}^*\left(\phi,\theta,0\right) \]

with \(\lambda=\lambda_a-\lambda_b\) and \(\mu=\lambda_c-\lambda_d\) the helicity differences of the final and initial states \(ab,cd\).

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

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. The reason is that it describes the interacting part of the 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:[1]

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

with \(\boldsymbol{I}\) the identity operator. With this in mind, there is an important restriction that the \(T\)-operator needs to comply with: unitarity. This means that \(\boldsymbol{S}\) should conserve probability, namely \(\boldsymbol{S}^\dagger\boldsymbol{S} = \boldsymbol{I}\).

K-matrix formalism#

Now there is a trick to ensure unitarity of \(\boldsymbol{S}\). We can 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 \(K\) is real. Finally, 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} \]

Resonances#

The challenge is now to choose a correct parametrization for the elements of \(\boldsymbol{K}\) so that it correctly describes the resonances we observe. There are several choices, but a common one is the following summation over the resonances \(R\):

(6)#\[ K_{ij} = \sum_R\frac{g_{R,i}^*g_{R,j}}{m_R^2-m^2} \]

with \(g_{R,i}\) the residue functions that can be further expressed as

(7)#\[ g_{R,i}=\gamma_{R,i}\sqrt{m_R\Gamma_R} \]

Implementation#

The challenge is to generate a correct parametrization for an arbitrary number of coupled channels \(n\) and an arbitrary number of resonances \(n_R\). Our approach is to construct an \(n \times n\) sympy.Matrix with Symbols as its elements. We then use substitute these Symbols with certain parametrizations using subs(). In order to generate symbols for \(n_R\) resonances and \(n\) channels, we use indexed symbols.

This approach is less elegant and (theoretically) slower than using MatrixSymbols. That approach is explored in TR-007.

It would be nice to use a Symbol to represent the number of channels \(n\) and specify its value later.

n_channels = sp.Symbol("n", integer=True, positive=True)

Unfortunately, this does not work well in the Matrix class. We therefore set variables \(n\) to a specific int value and define some other Symbols for the rest of the implementation.[2] The value we choose in this example is n_channels=1, because we want to see if this reproduces a Breit-Wigner function.[3]

n_channels = 1
i, j, R, n_resonances = sp.symbols("i j R n_R", integer=True, negative=False)
m = sp.Symbol("m", real=True)
M = sp.IndexedBase("m", shape=(n_resonances,))
Gamma = sp.IndexedBase("Gamma", shape=(n_resonances,))
gamma = sp.IndexedBase("gamma", shape=(n_resonances, n_channels))

The parametrization of \(K_{ij}\) from Eq. (6) can be expressed as follows:

def Kij(
    m: sp.Symbol,
    M: sp.IndexedBase,
    Gamma: sp.IndexedBase,
    gamma: sp.IndexedBase,
    i: int,
    j: int,
    n_resonances: int | sp.Symbol,
) -> sp.Expr:
    g_i = gamma[R, i] * sp.sqrt(M[R] * Gamma[R])
    g_j = gamma[R, j] * sp.sqrt(M[R] * Gamma[R])
    parametrization = (g_i * g_j) / (M[R] ** 2 - m**2)
    return sp.Sum(parametrization, (R, 0, n_resonances - 1))

\[\displaystyle K_{ij} = \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,i} {\gamma}_{R,j} {m}_{R}}{- m^{2} + {m}_{R}^{2}}\]

We now define the \(\boldsymbol{K}\)-matrix in terms of a Matrix with IndexedBase instances as elements that can serve as Symbols. These Symbols will be substituted with the parametrization later. We could of course have inserted the parametrization directly, but this slows down matrix multiplication in the following steps.

K_symbol = sp.IndexedBase("K", shape=(n_channels, n_channels))
K = sp.Matrix([[K_symbol[i, j] for j in range(n_channels)] for i in range(n_channels)])
display(K_symbol[i, j], K)

\[\displaystyle {K}_{i,j}\]

\[\displaystyle \left[\begin{matrix}{K}_{0,0}\end{matrix}\right]\]

The \(\boldsymbol{T}\)-matrix can now be computed from Eq. (5):

T = K * (sp.eye(n_channels) - sp.I * K).inv()
T

\[\displaystyle \left[\begin{matrix}\frac{{K}_{0,0}}{- i {K}_{0,0} + 1}\end{matrix}\right]\]

Next, we need to substitute the elements \(K_{i,j}\) with the parametrization we defined above:

T_subs = T.subs({
    K[i, j]: Kij(m, M, Gamma, gamma, i, j, n_resonances)
    for i in range(n_channels)
    for j in range(n_channels)
})
T_subs

\[\displaystyle \left[\begin{matrix}\frac{\sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}}}{- i \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} + 1}\end{matrix}\right]\]

Warning

It is important to perform doit() after subs(), otherwise the Sum cannot be evaluated and there will be no warning of a failed substitution.

Now indeed, when taking \(n_R=1\), the resulting element from the \(\boldsymbol{T}\)-matrix looks like a Breit-Wigner function (compare relativistic_breit_wigner())!

\[\displaystyle - \frac{{\Gamma}_{0} {\gamma}_{0,0}^{2} {m}_{0}}{m^{2} + i {\Gamma}_{0} {\gamma}_{0,0}^{2} {m}_{0} - {m}_{0}^{2}}\]

Generalization#

The above procedure has been condensed into a function that can handle an arbitrary number of resonances and an arbitrary number of channels.

def create_symbol_matrix(name: str, n: int) -> sp.Matrix:
    symbol = sp.IndexedBase("K", shape=(n, n))
    return sp.Matrix([[symbol[i, j] for j in range(n)] for i in range(n)])


def k_matrix(n_resonances: int, n_channels: int) -> sp.Matrix:
    # Define symbols
    m = sp.Symbol("m", real=True)
    M = sp.IndexedBase("m", shape=(n_resonances,))
    Gamma = sp.IndexedBase("Gamma", shape=(n_resonances,))
    gamma = sp.IndexedBase("gamma", shape=(n_resonances, n_channels))
    # Define K-matrix and T-matrix
    K = create_symbol_matrix("K", n_channels)
    T = K * (sp.eye(n_channels) - sp.I * K).inv()
    # Substitute elements
    return T.subs({
        K[i, j]: Kij(m, M, Gamma, gamma, i, j, n_resonances)
        for i in range(n_channels)
        for j in range(n_channels)
    })

Single channel, single resonance:

k_matrix(n_resonances=1, n_channels=1)[0, 0].doit().simplify()

\[\displaystyle - \frac{{\Gamma}_{0} {\gamma}_{0,0}^{2} {m}_{0}}{m^{2} + i {\Gamma}_{0} {\gamma}_{0,0}^{2} {m}_{0} - {m}_{0}^{2}}\]

Single channel, \(n_R\) resonances

k_matrix(n_resonances=sp.Symbol("n_R"), n_channels=1)[0, 0]

\[\displaystyle \frac{\sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}}}{- i \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} + 1}\]

Two channels, one resonance (Flatté function):

k_matrix(n_resonances=1, n_channels=2)[0, 0].doit().simplify()

\[\displaystyle - \frac{{\Gamma}_{0} {\gamma}_{0,0}^{2} {m}_{0}}{m^{2} + i {\Gamma}_{0} {\gamma}_{0,0}^{2} {m}_{0} + i {\Gamma}_{0} {\gamma}_{0,1}^{2} {m}_{0} - {m}_{0}^{2}}\]

Two channels, \(n_R\) resonances:

expr = k_matrix(n_resonances=sp.Symbol("n_R"), n_channels=2)[0, 0]
Math(sp.multiline_latex("", expr))

\[\begin{split}\displaystyle \begin{align*} \mathtt{\text{}} = & \frac{\left(i \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,1}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} - 1\right) \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}}}{\left(\sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}}\right) \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,1}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} + i \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} + i \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,1}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} - \left(\sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0} {\gamma}_{R,1} {m}_{R}}{- m^{2} + {m}_{R}^{2}}\right)^{2} - 1} \\ & + \frac{i \left(\sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0} {\gamma}_{R,1} {m}_{R}}{- m^{2} + {m}_{R}^{2}}\right)^{2}}{- \left(\sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}}\right) \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,1}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} - i \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} - i \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,1}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} + \left(\sum_{R=0}^{n_{R} - 1} \frac{{\Gamma}_{R} {\gamma}_{R,0} {\gamma}_{R,1} {m}_{R}}{- m^{2} + {m}_{R}^{2}}\right)^{2} + 1} \end{align*}\end{split}\]

Visualization#

Now, let’s use matplotlib, mpl_interactions, and symplot to visualize the \(\boldsymbol{K}\)-matrix for arbitrary \(n\) and \(n_R\).

plot_k_matrix(n_resonances=3, n_channels=1)

record

plot_k_matrix(n_resonances=2, n_channels=2)