Lorentz-invariant K-matrix

Lorentz-invariant K-matrix#

Physics#

The Lorentz-invariant description \(\boldsymbol{\hat{T}}\) of the \(\boldsymbol{T}\)-matrix is:

(1)#\[ \boldsymbol{T} = \sqrt{\boldsymbol{\rho^\dagger}} \, \boldsymbol{\hat{T}} \sqrt{\boldsymbol{\rho}} \]

with the phase space factor matrix \(\boldsymbol{\rho}\) defined as:

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

and

(3)#\[ \rho_i = \frac{2q_i}{m} = \sqrt{ \left[1-\left(\frac{m_{i,a}+m_{i,b}}{m}\right)^2\right] \left[1-\left(\frac{m_{i,a}-m_{i,b}}{m}\right)^2\right] } \]

This results in a similar transformation for the \(\boldsymbol{K}\)-matrix

(4)#\[ \boldsymbol{K} = \sqrt{\boldsymbol{\rho^\dagger}} \; \boldsymbol{\hat{K}} \sqrt{\boldsymbol{\rho}} \]

with (compare Eq. (5) in TR-005):

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

It’s common to integrate these phase space factors into the parametrization of \(K_{ij}\) as well:

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

Compare this with Eq. (6) in TR-005.

In addition, one often uses an “energy dependent” coupled_width() \(\Gamma_R(m)\) instead of a fixed width \(\Gamma_R\) as done in TR-005.

Implementation#

Wrapping expressions#

To keep a nice rendering, we wrap the expressions for phase_space_factor() and coupled_width() into a class that derives from Expr (see e.g. the implementation of BlattWeisskopfSquared). Note that we need to use partial_doit() to keep these expression symbols after evaluating the Sum.

@implement_doit_method()
class PhaseSpaceFactor(UnevaluatedExpression):
    is_commutative = True

    def __new__(
        cls,
        s: sp.Symbol,
        m_a: sp.Symbol,
        m_b: sp.Symbol,
        i: int,
        **hints,
    ) -> PhaseSpaceFactor:
        return create_expression(cls, s, m_a, m_b, i, **hints)

    def evaluate(self) -> sp.Expr:
        s, m_a, m_b, *_ = self.args
        return phase_space_factor_complex(s, m_a, m_b)

    def _latex(self, printer: LatexPrinter, *args) -> str:
        s = printer._print(self.args[0])
        i = self.args[-1]
        return Rf"\rho_{{{i}}}({s})"


@implement_doit_method()
class CoupledWidth(UnevaluatedExpression):
    is_commutative = True

    def __new__(
        cls,
        s: sp.Symbol,
        mass0: sp.IndexedBase,
        gamma0: sp.IndexedBase,
        m_a: sp.IndexedBase,
        m_b: sp.IndexedBase,
        angular_momentum: int,
        R: int | sp.Symbol,
        i: int,
        **hints,
    ) -> CoupledWidth:
        return create_expression(
            cls, s, mass0, gamma0, m_a, m_b, angular_momentum, R, i, **hints
        )

    def evaluate(self) -> sp.Expr:
        s, mass0, gamma0, m_a, m_b, angular_momentum, R, i = self.args

        def phsp_factor(s, m_a, m_b):
            return PhaseSpaceFactor(s, m_a, m_b, i)

        return coupled_width(
            s,
            mass0[R],
            gamma0[R, i],
            m_a[i],
            m_b[i],
            angular_momentum=angular_momentum,
            meson_radius=1,
            phsp_factor=phsp_factor,
        )

    def _latex(self, printer: LatexPrinter, *args) -> str:
        s = printer._print(self.args[0])
        R = self.args[-2]
        i = self.args[-1]
        return Rf"{{\Gamma_{{{R},{i}}}}}({s})"

And here is what the equations look like:

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

width_expr = CoupledWidth(m**2, M, Gamma, m_a, m_b, 0, R, i)
phsp_expr = PhaseSpaceFactor(m**2, m_a[i], m_b[i], i)

\[\displaystyle \begin{align*} {\Gamma_{R,i}}(m^{2}) = & \frac{B_{0}^2\left(\frac{\left(m^{2} - \left({m_{a}}_{i} - {m_{b}}_{i}\right)^{2}\right) \left(m^{2} - \left({m_{a}}_{i} + {m_{b}}_{i}\right)^{2}\right)}{4 m^{2}}\right) {\Gamma}_{R,i} \rho_{i}(m^{2})}{B_{0}^2\left(\frac{\left(- \left({m_{a}}_{i} - {m_{b}}_{i}\right)^{2} + {m}_{R}^{2}\right) \left(- \left({m_{a}}_{i} + {m_{b}}_{i}\right)^{2} + {m}_{R}^{2}\right)}{4 {m}_{R}^{2}}\right) \rho_{i}({m}_{R}^{2})} \end{align*}\]

\[\displaystyle \begin{align*} \rho_{i}(m^{2}) = & \frac{\sqrt[\mathrm{c}]{\frac{\left(m^{2} - \left({m_{a}}_{i} - {m_{b}}_{i}\right)^{2}\right) \left(m^{2} - \left({m_{a}}_{i} + {m_{b}}_{i}\right)^{2}\right)}{4 m^{2}}}}{8 \pi m} \end{align*}\]

Note

In PhaseSpaceFactor, we used PhaseSpaceFactorComplex instead of PhaseSpaceFactor, meaning that we choose the positive square root when values under the square root are negative. The only reason for doing this is, so that there is output in the figure under Visualization. The choice for which square root to choose has to do with analyticity (see TR-004) and choosing which Riemann sheet to connect to. This issue is ignored in this report.

Generalization#

The implementation is quite similar to that of TR-005, with the only difference being additional \(\boldsymbol{\rho}\)-matrix and the insertion of coupled width. Don’t forget to convert back to \(\boldsymbol{T}\) from \(\boldsymbol{\hat{T}}\) with Eq. (1).

def Kij_relativistic(
    m: sp.Symbol,
    M: sp.IndexedBase,
    Gamma: sp.IndexedBase,
    gamma: sp.IndexedBase,
    i: int,
    j: int,
    n_resonances: int | sp.Symbol,
    angular_momentum: int | sp.Symbol = 0,
) -> sp.Expr:
    def residue_function(i):
        return gamma[R, i] * sp.sqrt(
            M[R] * CoupledWidth(m**2, M, Gamma, m_a, m_b, angular_momentum, R, i)
        )

    g_i = residue_function(i)
    g_j = residue_function(j)
    parametrization = (g_i * g_j) / (M[R] ** 2 - m**2)
    return sp.Sum(parametrization, (R, 0, n_resonances - 1))


def relativistic_k_matrix(
    n_resonances: int,
    n_channels: int,
    angular_momentum: int | sp.Symbol = 0,
) -> sp.Matrix:
    # Define symbols
    m = sp.Symbol("m", real=True)
    M = sp.IndexedBase("m", shape=(n_resonances,))
    Gamma = sp.IndexedBase("Gamma", shape=(n_resonances, n_channels))
    gamma = sp.IndexedBase("gamma", shape=(n_resonances, n_channels))
    m_a = sp.IndexedBase("m_a", shape=(n_channels,))
    m_b = sp.IndexedBase("m_b", shape=(n_channels,))
    # Define phase space matrix
    sqrt_rho = sp.zeros(n_channels, n_channels)
    sqrt_rho_dagger = sp.zeros(n_channels, n_channels)
    for i in range(n_channels):
        rho = PhaseSpaceFactor(m**2, m_a[i], m_b[i], i)
        sqrt_rho[i, i] = sp.sqrt(rho)
        sqrt_rho_dagger[i, i] = 1 / sp.conjugate(sqrt_rho[i, i])
    # Define K-matrix and T-matrix
    K = create_symbol_matrix("K", n_channels)
    T_hat = K * (sp.eye(n_channels) - sp.I * rho * K).inv()
    T = sqrt_rho_dagger * T_hat * sqrt_rho
    # Substitute elements
    return T.subs({
        K[i, j]: Kij_relativistic(
            m=m,
            M=M,
            Gamma=Gamma,
            gamma=gamma,
            i=i,
            j=j,
            n_resonances=n_resonances,
            angular_momentum=angular_momentum,
        )
        for i in range(n_channels)
        for j in range(n_channels)
    })


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

Single channel, one resonance (compare relativistic_breit_wigner_with_ff()):

expr = relativistic_k_matrix(n_resonances=1, n_channels=1)[0, 0]
Math(
    sp.multiline_latex(
        lhs=expr,
        rhs=symplot.partial_doit(expr, sp.Sum).simplify(doit=False),
    )
)

\[\displaystyle \begin{align*} \frac{\sqrt{\rho_{0}(m^{2})} \sum_{R=0}^{0} \frac{{\Gamma_{R,0}}(m^{2}) {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}}}{\left(- i \rho_{0}(m^{2}) \sum_{R=0}^{0} \frac{{\Gamma_{R,0}}(m^{2}) {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} + 1\right) \overline{\sqrt{\rho_{0}(m^{2})}}} = &- \frac{{\Gamma_{0,0}}(m^{2}) {\gamma}_{0,0}^{2} {m}_{0} \sqrt{\rho_{0}(m^{2})}}{\left(m^{2} + i {\Gamma_{0,0}}(m^{2}) {\gamma}_{0,0}^{2} {m}_{0} \rho_{0}(m^{2}) - {m}_{0}^{2}\right) \overline{\sqrt{\rho_{0}(m^{2})}}} \end{align*}\]

Two channels, one resonance (‘Flatté’):

expr = relativistic_k_matrix(n_resonances=1, n_channels=2)[0, 0]
symplot.partial_doit(expr, sp.Sum).simplify(doit=False)

\[\displaystyle \frac{{\Gamma_{0,0}}(m^{2}) {\gamma}_{0,0}^{2} {m}_{0} \sqrt{\rho_{0}(m^{2})}}{\left(- m^{2} - i {\Gamma_{0,0}}(m^{2}) {\gamma}_{0,0}^{2} {m}_{0} \rho_{1}(m^{2}) - i {\Gamma_{0,1}}(m^{2}) {\gamma}_{0,1}^{2} {m}_{0} \rho_{1}(m^{2}) + {m}_{0}^{2}\right) \overline{\sqrt{\rho_{0}(m^{2})}}}\]

Single channel, \(n_R\) resonances:

relativistic_k_matrix(n_resonances, n_channels=1)[0, 0]

\[\displaystyle \frac{\sqrt{\rho_{0}(m^{2})} \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma_{R,0}}(m^{2}) {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}}}{\left(- i \rho_{0}(m^{2}) \sum_{R=0}^{n_{R} - 1} \frac{{\Gamma_{R,0}}(m^{2}) {\gamma}_{R,0}^{2} {m}_{R}}{- m^{2} + {m}_{R}^{2}} + 1\right) \overline{\sqrt{\rho_{0}(m^{2})}}}\]

Two channels, \(n_R\) resonances:

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

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

Visualization#

plot_relativistic_k_matrix(
    n_resonances=2,
    n_channels=1,
    angular_momentum=L,
    title="Relativistic $K$-matrix, single channel",
)