Analytic continuation

Tip

For more information about analytic continuation with symbolic models, see these Technical Reports: Chew-Mandelstam and Analyticity.

Analytic continuation allows one to handle resonances just below threshold (\(m_R < m_1 + m_1\) for a two-body decay \(R \to 12\)). In practice, this entails using a specific function for \(\rho\) in Eq. (1).

Import Python libraries

Hide code cell source

import timeit
import warnings
from textwrap import dedent

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from IPython.display import Markdown, Math
from matplotlib_inline.backend_inline import set_matplotlib_formats

from ampform.dynamics.phasespace import BreakupMomentum
from ampform.io import aslatex


def display_doit(exprs: list[sp.Expr], deep: bool = False) -> Math:
    return Math(aslatex({x: x.doit(deep=deep) for x in exprs}))


set_matplotlib_formats("svg")
warnings.filterwarnings("ignore")

Common definitions

Common choices for phasespace factor \(\rho\), and the break-up momentum \(q\), are the following:

s, m1, m2 = sp.symbols("s m1 m2", nonnegative=True)

Standard phase space factor (single square root)

Hide code cell source

from ampform.dynamics.phasespace import PhaseSpaceFactor

display_doit([PhaseSpaceFactor(s, m1, m2)])

\[\begin{split}\displaystyle \begin{aligned} \rho\left(s\right) \;&=\; \frac{\sqrt{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}}{s} \\ \end{aligned}\end{split}\]

Phase space factor that is defined along the whole real axis

Hide code cell source

from ampform.dynamics.phasespace import PhaseSpaceFactorComplex
from ampform.sympy.math import ComplexSqrt

rho_c = PhaseSpaceFactorComplex(s, m1, m2)
sqrt = ComplexSqrt(sp.Symbol("z"))
src = aslatex({
    rho_c: rho_c.evaluate(),
    sqrt: sqrt.get_definition(),
})
Math(src)

\[\begin{split}\displaystyle \begin{aligned} \rho^\mathrm{c}\left(s\right) \;&=\; \frac{\sqrt[\mathrm{c}]{s - \left(m_{1} - m_{2}\right)^{2}} \sqrt[\mathrm{c}]{s - \left(m_{1} + m_{2}\right)^{2}}}{s} \\ \sqrt[\mathrm{c}]{z} \;&=\; \begin{cases} i \sqrt{- z} & \text{for}\: z < 0 \\\sqrt{z} & \text{otherwise} \end{cases} \\ \end{aligned}\end{split}\]

The physically correct solution is to use a definition for the phase space factor that is analytic over the whole complex domain and complies with the Schwarz reflection principle. For arbitrary angular momenta \(\ell\), this can only be achieved by (numerically) computing a dispersion integral (via the Chew–Mandelstam function \(\Sigma_\ell\)). An algebraic solution for this integral exists only for S-waves (no angular momentum) and is given by:

Analytic continuation (Chew–Mandelstam) for S-waves

Hide code cell source

from ampform.dynamics.phasespace import BreakupMomentumComplex, PhaseSpaceFactorSWave

display_doit([
    PhaseSpaceFactorSWave(s, m1, m2),
    BreakupMomentumComplex(s, m1, m2),
])

\[\begin{split}\displaystyle \begin{aligned} \rho^\mathrm{CM}\left(s\right) \;&=\; - \frac{i \left(\frac{2 q^\mathrm{c}\left(s\right)}{\sqrt{s}} \log{\left(\frac{m_{1}^{2} + m_{2}^{2} + 2 \sqrt{s} q^\mathrm{c}\left(s\right) - s}{2 m_{1} m_{2}} \right)} - \left(m_{1}^{2} - m_{2}^{2}\right) \left(- \frac{1}{\left(m_{1} + m_{2}\right)^{2}} + \frac{1}{s}\right) \log{\left(\frac{m_{1}}{m_{2}} \right)}\right)}{\pi} \\ q^\mathrm{c}\left(s\right) \;&=\; \frac{\sqrt[\mathrm{c}]{s - \left(m_{1} - m_{2}\right)^{2}} \sqrt[\mathrm{c}]{s - \left(m_{1} + m_{2}\right)^{2}}}{2 \sqrt{s}} \\ \end{aligned}\end{split}\]

Analytic continuation for equal masses

Hide code cell source

from ampform.dynamics.phasespace import (
    BreakupMomentumSquared,
    EqualMassPhaseSpaceFactor,
    PhaseSpaceFactorAbs,
)

display_doit([
    EqualMassPhaseSpaceFactor(s, m1, m1),
    PhaseSpaceFactorAbs(s, m1, m2),
    BreakupMomentumSquared(s, m1, m2),
])

\[\begin{split}\displaystyle \begin{aligned} \rho^\mathrm{eq}\left(s\right) \;&=\; \begin{cases} \frac{i \log{\left(\left|{\frac{\hat{\rho}\left(s\right) + 1}{\hat{\rho}\left(s\right) - 1}}\right| \right)} \hat{\rho}\left(s\right)}{\pi} + \hat{\rho}\left(s\right) & \text{for}\: s > 4 m_{1}^{2} \\\frac{2 i \operatorname{atan}{\left(\frac{1}{\hat{\rho}\left(s\right)} \right)} \hat{\rho}\left(s\right)}{\pi} & \text{otherwise} \end{cases} \\ \hat{\rho}\left(s\right) \;&=\; \frac{2 \sqrt{\left|{q^2\left(s\right)}\right|}}{\sqrt{s}} \\ q^2\left(s\right) \;&=\; \frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{4 s} \\ \end{aligned}\end{split}\]

Cut structure

When continued into the complex plane, the standard (non-analytic) breakup momentum and phase space factor exhibit distinct cut structures, depending on the definition of the square root in the numerator. The distinction between these definitions only becomes relevant below threshold, which is not common in single-channel analyses. Only when resonances lie slightly below threshold (above the pseudothreshold \(|m_1-m_2|\)) is it relevant to work with a function that has cleaner cut structure.

The “split” square root definition,

Show code cell source

Hide code cell source

from ampform.dynamics.phasespace import (
    BreakupMomentumSplitSqrt,
    PhaseSpaceFactorSplitSqrt,
)

display_doit([
    BreakupMomentumSplitSqrt(s, m1, m2),
    PhaseSpaceFactorSplitSqrt(s, m1, m2),
])

\[\begin{split}\displaystyle \begin{aligned} q\left(s\right) \;&=\; \frac{\sqrt{s - \left(m_{1} - m_{2}\right)^{2}} \sqrt{s - \left(m_{1} + m_{2}\right)^{2}}}{2 \sqrt{s}} \\ \rho\left(s\right) \;&=\; \frac{\sqrt{s - \left(m_{1} - m_{2}\right)^{2}} \sqrt{s - \left(m_{1} + m_{2}\right)^{2}}}{s} \\ \end{aligned}\end{split}\]

leads to a cleaner cut structure than a definition with a “single” square root (AmpForm’s default implementation),

Show code cell source

Hide code cell source

display_doit([
    BreakupMomentum(s, m1, m2),
    PhaseSpaceFactor(s, m1, m2),
])

\[\begin{split}\displaystyle \begin{aligned} q\left(s\right) \;&=\; \frac{\sqrt{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}}{2 \sqrt{s}} \\ \rho\left(s\right) \;&=\; \frac{\sqrt{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}}{s} \\ \end{aligned}\end{split}\]

or with a Källén function,

Show code cell source

Hide code cell source

from ampform.dynamics.phasespace import BreakupMomentumKallen, PhaseSpaceFactorKallen

display_doit(
    [
        BreakupMomentumKallen(s, m1, m2),
        PhaseSpaceFactorKallen(s, m1, m2),
    ],
    deep=True,
)

\[\begin{split}\displaystyle \begin{aligned} q\left(s\right) \;&=\; \frac{\sqrt{\frac{m_{1}^{4} - 2 m_{1}^{2} m_{2}^{2} - 2 m_{1}^{2} s + m_{2}^{4} - 2 m_{2}^{2} s + s^{2}}{s}}}{2} \\ \rho\left(s\right) \;&=\; \frac{\sqrt{m_{1}^{4} - 2 m_{1}^{2} m_{2}^{2} - 2 m_{1}^{2} s + m_{2}^{4} - 2 m_{2}^{2} s + s^{2}}}{s} \\ \end{aligned}\end{split}\]

Show code cell source

Hide code cell source

x_min, x_max = -0.1, +1.5
y_max = 0.8
z_max = 0.5
x = np.linspace(x_min, x_max, num=500)
X, Y = np.meshgrid(x, np.linspace(-y_max, +y_max, num=300))
S = X + 1j * Y
ϵi = 1e-7j

parameters = {m1: 0.2, m2: 0.6}
thr_neg = (parameters[m1] - parameters[m2]) ** 2
thr_pos = (parameters[m1] + parameters[m2]) ** 2

fig, axes = plt.subplots(
    dpi=200,
    figsize=(12, 5),
    ncols=3,
    nrows=2,
    sharex=True,
    sharey=True,
)
fig.subplots_adjust(bottom=0, hspace=0.02, left=0, right=1, top=1, wspace=0.02)
fig.patch.set_facecolor("none")
for ax in axes.flatten():
    ax.patch.set_facecolor("none")
    ax.spines["bottom"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_visible(False)
    ax.hlines(0, thr_neg, thr_pos, color="black", lw=1.5, zorder=10)
    ax.scatter([thr_neg, thr_pos], [0, 0], c="black", s=25, zorder=10)
for ax in axes[-1]:
    ax.set_xlabel(R"$\mathrm{Re}\,s$", labelpad=-12)
    ax.set_xticks([0])
for ax in axes[:, 0]:
    ax.set_ylabel(R"$\mathrm{Im}\,s$")
    ax.set_ylim(-y_max, +y_max)
    ax.set_yticks([0])

style = dict(
    cmap="coolwarm",
    rasterized=True,
    vmin=-z_max,
    vmax=+z_max,
    zorder=-10,
)
for ax, expr_class, title in [
    (axes[0, 0], BreakupMomentumSplitSqrt, '$q(s)$, "split" square root'),
    (axes[0, 1], BreakupMomentum, "$q(s)$, single square root"),
    (axes[0, 2], BreakupMomentumKallen, "$q(s)$ with Källén function"),
    (axes[1, 0], PhaseSpaceFactorSplitSqrt, R'$\rho(s)$, "split" square root'),
    (axes[1, 1], PhaseSpaceFactor, R"$\rho(s)$, single square root"),
    (axes[1, 2], PhaseSpaceFactorKallen, R"$\rho(s)$ with Källén function"),
]:
    expr = expr_class(s, m1, m2)
    func = sp.lambdify(s, expr.doit().subs(parameters))
    mesh = ax.pcolormesh(X, Y, func(S).imag, **style)
    ax.plot(x, func(x + ϵi).real, c="darkblue")
    ax.plot(x, func(x + ϵi).imag, c="darkgreen")
    ax.set_title(title, y=0.9)
cbar = fig.colorbar(mesh, ax=axes, pad=0.01)
cbar.ax.set_ylabel("imag", labelpad=0, rotation=270)
cbar.ax.set_yticks([-z_max, +z_max])
cbar.ax.set_yticklabels(["$-$", "$+$"])
ax = axes[0, 0]
kwargs = dict(transform=ax.transAxes, ha="right", va="top")
ax.text(0.99, 0.85, "real", c="darkblue", **kwargs)
ax.text(0.99, 0.48, "imag", c="darkgreen", **kwargs)
plt.show()

Numerical precision and performance

The numerical precision of functions decreases as the number of operations (nodes) in the function tree increases. As such, BreakupMomentumSplitSqrt and PhaseSpaceFactorSplitSqrt have a slightly lower numerical precision and is slower, even if its cut structure is cleaner. For this reason, AmpForm uses the definition with the single square root by default.

The implementation of the breakup momentum and the phasespace factor that uses the Källén function and decreases numerical precision by a factor \(10\).

Function for rendering differences with algebraic values

Hide code cell source

def render_floating_point_differences(exprs: dict[str, sp.Expr]) -> None:
    m1_val, m2_val = parameters.values()
    n_eval = 1_000
    fig, ax = plt.subplots(figsize=(7, 3))
    fig.patch.set_facecolor("none")
    ax.patch.set_facecolor("none")
    ax.spines["bottom"].set_position("zero")
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_visible(False)
    ax.axvline(thr_pos, c="red", label="Threshold")
    y_max = 0
    src = dedent(f"""
    | Name | Operations | Mean abs. difference | Standard deviation | {n_eval} evaluations (ms) |
    |:-----|-----------:|---------------------:|-------------------:|------:|
    """)
    for name, expr in exprs.items():
        n_ops = sp.count_ops(expr.doit())
        func = sp.lambdify(args, expr.doit())
        exact = [expr.doit().subs(parameters).subs(s, x).n() for x in s_array]
        numerical = func(s_array, m1_val, m2_val)
        diff = numerical - np.array(exact, dtype=float)
        diff_abs = np.abs(diff)
        mean = sp.sympify(diff_abs.mean()).n(3)
        std = sp.sympify(diff_abs.std()).n(3)
        run = lambda: func(test_array, m1_val, m2_val)  # noqa:B023,E731
        ms = 1e3 * np.array(timeit.repeat(run, number=n_eval, repeat=10))
        time = Rf"${ms.mean():.2f} \pm {ms.std():.2f}$"
        src += f"| {name} | {n_ops:,d} | ${sp.latex(mean)}$ | ${sp.latex(std)}$ | {time} |\n"
        y_max = max(y_max, diff_abs[3:].max())
        ax.fill_between(s_array, diff_abs, alpha=0.7, label=name, step="mid")
    ax.set_ylim(0, y_max)
    ax.legend(
        bbox_to_anchor=(1, 1.1),
        framealpha=0,
        loc="upper right",
    )
    ax.set_xlabel("$s$")
    ax.set_ylabel("Abs. diff. with algebraic value")
    fig.tight_layout()
    plt.show()
    display(Markdown(src))


args = s, m1, m2
s_array = np.linspace(thr_pos + 1e-3, x_max, num=50)
test_array = np.linspace(0, 10, num=1000)

Phase space factor

Hide code cell source

render_floating_point_differences({
    "Källén function": PhaseSpaceFactorKallen(*args),
    "Split square root": PhaseSpaceFactorSplitSqrt(*args),
    "Single square root": PhaseSpaceFactor(*args),
})

Name	Operations	Mean abs. difference	Standard deviation	1000 evaluations (ms)
Källén function	21	\(2.03 \cdot 10^{-16}\)	\(5.54 \cdot 10^{-16}\)	\(13.63 \pm 0.25\)
Split square root	12	\(4.07 \cdot 10^{-17}\)	\(5.05 \cdot 10^{-17}\)	\(11.23 \pm 0.13\)
Single square root	10	\(1.57 \cdot 10^{-17}\)	\(3.51 \cdot 10^{-17}\)	\(9.57 \pm 0.21\)

Breakup momentum

Hide code cell source

render_floating_point_differences({
    "Källén function": BreakupMomentumKallen(*args),
    "Split square root": BreakupMomentumSplitSqrt(*args),
    "Single square root": BreakupMomentum(*args),
})

Name	Operations	Mean abs. difference	Standard deviation	1000 evaluations (ms)
Källén function	22	\(9.41 \cdot 10^{-17}\)	\(2.21 \cdot 10^{-16}\)	\(14.70 \pm 0.23\)
Split square root	15	\(2.37 \cdot 10^{-17}\)	\(2.93 \cdot 10^{-17}\)	\(13.93 \pm 0.15\)
Single square root	13	\(1.75 \cdot 10^{-17}\)	\(2.4 \cdot 10^{-17}\)	\(12.35 \pm 0.20\)

Behavior along the real axis

%matplotlib widget

Lambdify to numerical functions

Hide code cell source

import symplot
from ampform.sympy.math import ComplexSqrt

m = sp.Symbol("m", nonnegative=True)
rho_c = PhaseSpaceFactorComplex(m**2, m1, m2)
rho_cm = PhaseSpaceFactorSWave(m**2, m1, m2)
rho_ac = EqualMassPhaseSpaceFactor(m**2, m1, m2)
np_rho_c, sliders = symplot.prepare_sliders(plot_symbol=m, expression=rho_c.doit())
np_rho_ac = sp.lambdify((m, m1, m2), rho_ac.doit())
np_rho_cm = sp.lambdify((m, m1, m2), rho_cm.doit())
q = BreakupMomentumComplex(m**2, m1, m2)
np_breakup_momentum = sp.lambdify((m, m1, m2), 2 * q.doit())

Run this notebook on Binder or locally on Jupyter Lab to interactively modify the parameters.

Define sliders and plot range

Hide code cell source

plot_domain = np.linspace(0, 3, num=500)
sliders.set_ranges(
    m1=(0, 2, 200),
    m2=(0, 2, 200),
)
sliders.set_values(
    m1=0.3,
    m2=0.75,
)

Show code cell source

Hide code cell source

import mpl_interactions.ipyplot as iplt

import symplot

fig, axes = plt.subplots(
    ncols=2,
    nrows=2,
    figsize=[8, 5],
    sharex=True,
    sharey=True,
)
fig.subplots_adjust(bottom=0.08, hspace=0.1, left=0.01, right=0.99, top=1, wspace=0.05)
fig.canvas.footer_visible = False
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False
fig.patch.set_facecolor("none")
for ax in axes.flatten():
    ax.patch.set_facecolor("none")
    ax.spines["bottom"].set_position("zero")
    ax.spines["left"].set_position("zero")
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_visible(False)

(ax_q, ax_rho), (ax_rho_ac, ax_rho_cm) = axes
for ax in axes[-1]:
    ax.set_xlabel("$m$")
for ax in axes.flatten():
    ax.set_yticks([])
    ax.set_yticks([])


def func_imag(func, *args, **kwargs):
    return lambda *args, **kwargs: func(*args, **kwargs).imag


def func_real(func, *args, **kwargs):
    return lambda *args, **kwargs: func(*args, **kwargs).real


ylim = (-0.1, 1.4)
q_math = ComplexSqrt(sp.Symbol("q^2")) / (8 * sp.pi)
ax_q.set_title(f"${sp.latex(q_math)}$", y=0.85)
controls = iplt.plot(
    plot_domain,
    func_real(np_breakup_momentum),
    label="real",
    **sliders,
    ylim=ylim,
    ax=ax_q,
    alpha=0.7,
)
iplt.plot(
    plot_domain,
    func_imag(np_breakup_momentum),
    label="imaginary",
    controls=controls,
    ylim=ylim,
    ax=ax_q,
    alpha=0.7,
)

ax_rho.set_title(f"${sp.latex(rho_c)}$", y=0.85)
iplt.plot(
    plot_domain,
    func_real(np_rho_c),
    label="real",
    controls=controls,
    ylim=ylim,
    ax=ax_rho,
    alpha=0.7,
)
iplt.plot(
    plot_domain,
    func_imag(np_rho_c),
    label="imaginary",
    controls=controls,
    ylim=ylim,
    ax=ax_rho,
    alpha=0.7,
)

ax_rho_ac.set_title(R"equal mass $\rho^\mathrm{eq}(m^2)$", y=0.85)
iplt.plot(
    plot_domain,
    func_real(np_rho_ac),
    label="real",
    controls=controls,
    ylim=ylim,
    ax=ax_rho_ac,
    alpha=0.7,
)
iplt.plot(
    plot_domain,
    func_imag(np_rho_ac),
    label="imaginary",
    controls=controls,
    ylim=ylim,
    ax=ax_rho_ac,
    alpha=0.7,
)

ax_rho_cm.set_title(R"Chew-Mandelstam $\rho^\mathrm{CM}(m^2)$", y=0.85)
iplt.plot(
    plot_domain,
    func_real(np_rho_cm),
    label="real",
    controls=controls,
    ylim=ylim,
    ax=ax_rho_cm,
    alpha=0.7,
)
iplt.plot(
    plot_domain,
    func_imag(np_rho_cm),
    label="imaginary",
    controls=controls,
    ylim=ylim,
    ax=ax_rho_cm,
    alpha=0.7,
)

ax_rho.legend(loc="upper right")
plt.show()