Analytic continuation

Show code cell content

Hide code cell content

import logging
import warnings

import ampform
import matplotlib.pyplot as plt
import qrules
from IPython.display import Markdown, Math

from tensorwaves.data import (
    SympyDataTransformer,
    TFPhaseSpaceGenerator,
    TFUniformRealNumberGenerator,
)
from tensorwaves.function.sympy import create_parametrized_function

logging.getLogger("tensorwaves.data").setLevel(logging.ERROR)  # hide progress bars
warnings.filterwarnings("ignore")

Sometimes qrules finds resonances that lie outside phase space, because resonances can ‘leak’ into phase space. The example below is simplified and the two selected resonances are a bit unusual, but it serves to illustrate how to handle these sub-threshold resonances.

reaction = qrules.generate_transitions(
    initial_state="D0",
    final_state=["K-", "K+", "K0"],
    allowed_intermediate_particles=["a(0)(980)0", "a(2)(1320)0"],
    formalism="canonical-helicity",
)
src = qrules.io.asmermaid(
    reaction,
    collapse_graphs=True,
    render_node=False,
    markdown=True,
)
Markdown(src)

        flowchart LR
    T0_0["$$0: K^{-}$$"]
    T0_1["$$1: K^{+}$$"]
    T0_2["$$2: K^{0}$$"]
    T0_N0["$$D^{0}$$"]
    T0_N1@{ shape: text, label: " " }
    T0_3("$$\begin{gathered} a_{0}(980)^{0} \\\ a_{2}(1320)^{0} \end{gathered}$$")
    T0_N0 --- T0_3
    T0_3 --- T0_N1
    T0_N0 --- T0_2
    T0_N1 --- T0_0
    T0_N1 --- T0_1
    

Of the two resonances, \(a_0(980)\) lies just below threshold ― its mass is smaller than the masses of the two decay products combined:

Show code cell source

Hide code cell source

pdg = qrules.load_pdg()
a_meson = pdg["a(0)(980)0"]
k_plus = pdg["K+"]
k_minus = pdg["K+"]
two_k_mass = round(k_minus.mass + k_plus.mass, 3)
Math(
    Rf"""\begin{{array}}{{c}}
    m_{{{a_meson.latex}}} = {a_meson.mass} \pm {a_meson.width}\;\mathrm{{GeV}} \\
    m_{{{k_plus.latex}}} + m_{{{k_minus.latex}}} = {two_k_mass}\;\mathrm{{GeV}}
    \end{{array}}"""
)

\[\begin{split}\displaystyle \begin{array}{c} m_{a_{0}(980)^{0}} = 0.98 \pm 0.075\;\mathrm{GeV} \\ m_{K^{+}} + m_{K^{+}} = 0.987\;\mathrm{GeV} \end{array}\end{split}\]

To correctly describe the dynamics for this resonance, we should make use of analytic continuation. As opposed to Step 1: Formulate model, we now construct a RelativisticBreitWignerBuilder where we set its phase space factor to PhaseSpaceFactorSWave:

from ampform.dynamics import PhaseSpaceFactorSWave
from ampform.dynamics.builder import (
    RelativisticBreitWignerBuilder,
    create_non_dynamic_with_ff,
    create_relativistic_breit_wigner_with_ff,
)

model_builder = ampform.get_builder(reaction)
analytic_breit_wigner_builder = RelativisticBreitWignerBuilder(
    form_factor=True,
    energy_dependent_width=True,
    phsp_factor=PhaseSpaceFactorSWave,  # ty:ignore[invalid-argument-type]
)
model_builder.dynamics.assign("D0", create_non_dynamic_with_ff)
model_builder.dynamics.assign("a(0)(980)0", analytic_breit_wigner_builder)
model_builder.dynamics.assign("a(2)(1320)0", create_relativistic_breit_wigner_with_ff)
model = model_builder.formulate()

Note that we have reduced the coupling for the sub-threshold resonance \(a_0(980)^0\), so that it doesn’t dominate over the \(a_2(1320)^0\) resonance.

The effect can be seen once we generate data. Despite the fact that the resonance lies outside phase space, there is still a contribution to the intensity:

for symbol in model.parameter_defaults:
    if not str(symbol).startswith("C") or "a_{0}(980)^{0}" not in str(symbol):
        continue
    model.parameter_defaults[symbol] = 0.1

intensity = create_parametrized_function(
    expression=model.expression.doit(),
    parameters=model.parameter_defaults,
    backend="jax",
)
helicity_transformer = SympyDataTransformer.from_sympy(
    model.kinematic_variables, backend="numpy"
)
rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=reaction.initial_state[-1].mass,
    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},
)
phsp_momenta = phsp_generator.generate(100_000, rng)

Show code cell source

Hide code cell source

%config InlineBackend.figure_formats = ['svg']

import numpy as np

phsp = helicity_transformer(phsp_momenta)
intensities = np.array(intensity(phsp))

resonances = sorted(
    reaction.get_intermediate_particles(),
    key=lambda p: p.mass,
)
evenly_spaced_interval = np.linspace(0, 1, len(resonances))
colors = [plt.cm.rainbow(x) for x in evenly_spaced_interval]

fig, ax = plt.subplots()
ax.set_xlabel("$m_{02}$ [GeV]")
for p, color in zip(resonances, colors, strict=True):
    ax.axvline(x=p.mass, linestyle="dotted", label=p.name, color=color)
ax.set_yticks([])
ax.hist(phsp["m_01"], bins=100, alpha=0.5, weights=intensities)
ax.legend()
plt.show()