Analytic continuation

Note

Improvements to analytic continuation in AmpForm are currently being developed in Chew-Mandelstam and Analyticity.

Analytic continuation allows one to handle resonances just below threshold (\(m_0 < m_1 + m_2\) in Eq. (3)). In practice, this entails using a specific function for \(\rho\) in Eq. (1).

Definitions

Three usual choices for \(\rho\) are the following:

Import Python libraries

Hide code cell source

import warnings
from typing import Any

import sympy as sp
from IPython.display import Math
from sympy.printing.latex import LatexPrinter

from ampform.dynamics.phasespace import ChewMandelstamIntegral
from ampform.io import aslatex
from ampform.kinematics.phasespace import BreakupMomentum
from ampform.sympy import unevaluated

warnings.filterwarnings("ignore")

1) Break-up momentum

The sqrt() or ComplexSqrt of BreakupMomentumSquared:

from ampform.dynamics import BreakupMomentumSquared

s, m1, m2, L = sp.symbols("s, m1, m2, L", nonnegative=True)
q_squared = BreakupMomentumSquared(s, m1, m2)
Math(aslatex({q_squared: q_squared.evaluate()}))

\[\begin{split}\displaystyle \begin{aligned} 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}\]

2) ‘Normal’ phase space factor

The ‘normal’ PhaseSpaceFactor (the denominator makes the difference to (1)!):

from ampform.dynamics import PhaseSpaceFactor

rho = PhaseSpaceFactor(s, m1, m2)
Math(aslatex({rho: rho.evaluate()}))

\[\begin{split}\displaystyle \begin{aligned} \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}\]

3) ‘Complex’ phase space factor

A PhaseSpaceFactorComplex that uses ComplexSqrt:

from ampform.dynamics import PhaseSpaceFactorComplex

rho_c = PhaseSpaceFactorComplex(s, m1, m2)
Math(aslatex({rho_c: rho_c.evaluate()}))

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

4) ‘Analytic continuation’ of the phase space factor

The following ‘case-by-case’ analytic continuation for decay products with an equal mass, EqualMassPhaseSpaceFactor:

from ampform.dynamics import EqualMassPhaseSpaceFactor

rho_ac = EqualMassPhaseSpaceFactor(s, m1, m2)
Math(aslatex({rho_ac: rho_ac.evaluate()}))

\[\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 > \left(m_{1} + m_{2}\right)^{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} \\ \end{aligned}\end{split}\]

with

Show code cell source

Hide code cell source

from ampform.dynamics import PhaseSpaceFactorAbs

rho_hat = PhaseSpaceFactorAbs(s, m1, m2)
Math(aslatex({rho_hat: rho_hat.evaluate()}))

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

(Mind the absolute value.)

5) Chew-Mandelstam for \(S\)-waves

A PhaseSpaceFactorSWave that uses chew_mandelstam_s_wave():

from ampform.dynamics import PhaseSpaceFactorSWave

rho_cm = PhaseSpaceFactorSWave(s, m1, m2)
Math(aslatex({rho_cm: rho_cm.evaluate()}))

\[\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} \\ \end{aligned}\end{split}\]

Cut structure

When analytically continued into the complex plane, the breakup momentum and phase space factor exhibit distinct cut structures, depending on the definition of the square root in the numerator. The separated square root definition,

\[\begin{split} \begin{aligned} \rho_\alpha(s) &= \frac{\sqrt{s-(m_{1,\alpha}-m_{2,\alpha})^2}\sqrt{s-(m_{1,\alpha}+m_{2,\alpha})^2}}{s} \\ q_\alpha(s) &= \frac{\sqrt{s-(m_{1,\alpha}-m_{2,\alpha})^2}\sqrt{s-(m_{1,\alpha}+m_{2,\alpha})^2}}{2\sqrt{s}} \,, \end{aligned} \end{split}\]

leads to a cleaner cut structure than

\[\begin{split} \begin{aligned} \rho_\alpha(s) &= \frac{2q_\alpha(s)}{\sqrt{s}} = \frac{\sqrt{(s-(m_{1,\alpha}-m_{2,\alpha})^2) (s-(m_{1,\alpha}+m_{2,\alpha})^2})}{s} \\ q_\alpha(s) &= \frac{\sqrt{(s-(m_{1,\alpha}-m_{2,\alpha})^2)(s-(m_{1,\alpha}+m_{2,\alpha})^2})}{2\sqrt{s}} \,. \end{aligned} \end{split}\]

Here we investigate the cut structure of each of these definitions using AmpForm and SymPy.

Definition of phase space factor with double cut structure

Hide code cell source

@unevaluated
class PhaseSpaceFactorOld(sp.Expr):
    s: Any
    m1: Any
    m2: Any

    def evaluate(self) -> sp.Expr:
        s, m1, m2 = self.args
        q_squared = BreakupMomentumSquared(s, m1, m2)
        return 2 * sp.sqrt(q_squared) / sp.sqrt(s)

    def _latex_repr_(self, printer: LatexPrinter, *args) -> str:
        s = self.args[0]
        s_latex = printer._print(s)
        return Rf"\rho\left({s_latex}\right)"

Show code cell source

Hide code cell source

q = BreakupMomentum(s, m1, m2)
Math(aslatex({q: q.doit()}))

\[\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}} \\ \end{aligned}\end{split}\]

Show code cell source

Hide code cell source

rho = PhaseSpaceFactor(s, m1, m2)
Math(aslatex({rho: rho.doit()}))

\[\begin{split}\displaystyle \begin{aligned} \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}\]

Show code cell source

Hide code cell source

q = BreakupMomentumSquared(s, m1, m2)
Math(aslatex({q: sp.sqrt(q).doit()}))

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

Show code cell source

Hide code cell source

rho = PhaseSpaceFactorOld(s, m1, m2)
Math(aslatex({rho: rho.doit()}))

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

Show code cell source

Hide code cell source

import matplotlib.pyplot as plt
import numpy as np
from matplotlib_inline.backend_inline import set_matplotlib_formats

set_matplotlib_formats("svg")

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

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

rho_func = sp.lambdify(s, sp.I * PhaseSpaceFactor(s, m1, m2).doit().subs(parameters))
rho_old_func = sp.lambdify(
    s, sp.I * PhaseSpaceFactorOld(s, m1, m2).doit().subs(parameters)
)
q_func = sp.lambdify(s, BreakupMomentum(s, m1, m2).doit().subs(parameters))
q_old_func = sp.lambdify(
    s, sp.sqrt(BreakupMomentumSquared(s, m1, m2)).doit().subs(parameters)
)

fig, axes = plt.subplots(figsize=(9, 7), ncols=2, nrows=2, sharey=True)
fig.patch.set_facecolor("none")
ax1, ax2, ax3, ax4 = axes.flatten()
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.set_xlabel(R"$\mathrm{Re}\,s$", labelpad=-12)
    ax.set_xticks([0])
    ax.hlines(0, thr_neg, thr_pos, color="black", lw=1.5, zorder=10)
    ax.scatter([thr_neg, thr_pos], [0, 0], color="black", s=25, zorder=10)
ax1.set_ylabel(R"$\mathrm{Im}\,s$")
ax1.set_ylim(-y_max, +y_max)
ax1.set_yticks([0])

ax3.set_ylabel(R"$\mathrm{Im}\,s$")
ax3.set_ylim(-y_max, +y_max)
ax3.set_yticks([0])
style = dict(
    cmap="coolwarm",
    rasterized=True,
    vmin=-z_max,
    vmax=+z_max,
    zorder=-10,
)
mesh = ax1.pcolormesh(X, Y, rho_func(S).real, **style)
cbar = fig.colorbar(mesh, ax=ax1, pad=0.01)
cbar.ax.set_ylabel(R"$\mathrm{Re}\,i\rho$", labelpad=0, rotation=270)
cbar.ax.set_yticks([-z_max, +z_max])
cbar.ax.set_yticklabels(["$-$", "$+$"])
mesh = ax2.pcolormesh(X, Y, q_func(S).imag, **style)
cbar = fig.colorbar(mesh, ax=ax2, pad=0.01)
cbar.ax.set_ylabel(R"$\mathrm{Im}\,q$", labelpad=0, rotation=270)
cbar.ax.set_yticks([-z_max, +z_max])
cbar.ax.set_yticklabels(["$-$", "$+$"])
mesh = ax3.pcolormesh(X, Y, rho_old_func(S).real, **style)
cbar = fig.colorbar(mesh, ax=ax3, pad=0.01)
cbar.ax.set_ylabel(R"$\mathrm{Re}\,i\rho$", labelpad=0, rotation=270)
cbar.ax.set_yticks([-z_max, +z_max])
cbar.ax.set_yticklabels(["$-$", "$+$"])
mesh = ax4.pcolormesh(X, Y, q_old_func(S).imag, **style)
cbar = fig.colorbar(mesh, ax=ax4, pad=0.01)
cbar.ax.set_ylabel(R"$\mathrm{Im}\,q$", labelpad=0, rotation=270)
cbar.ax.set_yticks([-z_max, +z_max])
cbar.ax.set_yticklabels(["$-$", "$+$"])
ax1.plot(X[0], rho_func(X[0] + ϵi).real, c="darkblue")
ax1.plot(X[0], rho_func(X[0] + ϵi).imag, c="darkgreen")
ax2.plot(X[0], q_func(X[0] + ϵi).real, c="darkblue")
ax2.plot(X[0], q_func(X[0] + ϵi).imag, c="darkgreen")
ax3.plot(X[0], rho_old_func(X[0] + ϵi).real, c="darkblue")
ax3.plot(X[0], rho_old_func(X[0] + ϵi).imag, c="darkgreen")
ax4.plot(X[0], q_old_func(X[0] + ϵi).real, c="darkblue")
ax4.plot(X[0], q_old_func(X[0] + ϵi).imag, c="darkgreen")
ax1.text(
    0.99,
    0.48,
    "real",
    c="darkblue",
    transform=ax1.transAxes,
    ha="right",
    va="top",
)
ax1.text(
    0.99,
    0.97,
    "imag",
    c="darkgreen",
    transform=ax1.transAxes,
    ha="right",
    va="top",
)
ax1.set_title("Phase space factor double sqrt")
ax2.set_title("Break-up momentum factor double sqrt")
ax3.set_title("Phase space factor single sqrt")
ax4.set_title("Break-up momentum factor single sqrt")
fig.tight_layout()
plt.show()

Note

When defining the break-up momentum and the phasespace factor with a single square root in the numerator a more complex cut structure emerges when continuing the functions into the complex \(s\)-plane.

Dispersion integral

To get an analytic phasespace factor for higher angular momenta, the one has to compute the dispersion integral. According to PDG, Rev. Resonances the once-substracted dispersion integral is given by:

\[ \Sigma_a(s+0 i)=\frac{s-s_{\mathrm{thr}_a}}{\pi} \int_{s_{\mathrm{thr}_a}}^{\infty} \frac{\rho_a\left(s^{\prime}\right) n_a^2\left(s^{\prime}\right)}{\left(s^{\prime}-s_{\mathrm{thr}_a}\right)\left(s^{\prime}-s-i 0\right)} \mathrm{d} s^{\prime} \]

integral_expr = ChewMandelstamIntegral(s, m1, m2, L)
integral_expr.doit(deep=False)

\[\displaystyle \frac{s - \left(m_{1} + m_{2}\right)^{2}}{\pi} \int\limits_{\left(m_{1} + m_{2}\right)^{2}}^{\infty} \frac{\mathcal{F}_{L}\left(x, m_{1}, m_{2}\right)^{2} \rho\left(x\right)}{\left(x - \left(m_{1} + m_{2}\right)^{2}\right) \left(- i \epsilon - s + x\right)}\, dx\]

integral_s_wave_func = sp.lambdify(
    [s, m1, m2, integral_expr.epsilon],
    integral_expr.subs(L, 0).doit(),
)
integral_s_wave_func = np.vectorize(integral_s_wave_func)

Same definition for P-wave function

Hide code cell source

integral_p_wave_func = sp.lambdify(
    [s, m1, m2, integral_expr.epsilon],
    integral_expr.subs(L, 1).doit(),
)
integral_p_wave_func = np.vectorize(integral_p_wave_func)

Show code cell source

Hide code cell source

s_values = np.linspace(-0.15, 1.4, num=200)
s_wave_values = integral_s_wave_func(s_values, m1_val, m2_val, epsilon=1e-5)
p_wave_values = integral_p_wave_func(s_values, m1_val, m2_val, epsilon=1e-5)

l_val = [0, 1]
fig, axes = plt.subplots(figsize=(6, 7), nrows=2, sharex=True)
fig.patch.set_facecolor("none")
ax1, ax2 = axes
fig.suptitle(f"Symbolic dispersion integrals for $m_1={m1_val:.2f}, m_2={m2_val:.2f}$")
for ax in axes:
    ax.axhline(0, linewidth=0.5, c="black")
    ax.axvline(thr_pos, linestyle="--", color="red", alpha=0.7)
    ax.patch.set_facecolor("none")
    ax.set_title(f"$L = {l_val}$")
    ax.set_ylabel(R"$16\pi \; \Sigma(s)$")
axes[-1].set_xlabel("$s$ (GeV$^2$)")

ax1.set_title("$S$-wave ($L=0$)")
ax1.plot(s_values, 16 * np.pi * s_wave_values.real, label="real")
ax1.plot(s_values, 16 * np.pi * s_wave_values.imag, label="imaginary")

ax2.set_title("$P$-wave ($L=1$)")
ax2.plot(s_values, 16 * np.pi * p_wave_values.real, label="real")
ax2.plot(s_values, 16 * np.pi * p_wave_values.imag, label="imaginary")

ax1.legend(framealpha=0)
fig.tight_layout()
plt.show()

Interactive visualization

%matplotlib widget

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())
np_breakup_momentum = sp.lambdify(
    (m, m1, m2),
    2 * ComplexSqrt(q_squared.subs(s, m**2).doit()),
)

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

plot_domain = np.linspace(0, 3, 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, 8],
    sharex=True,
    sharey=True,
)
fig.canvas.footer_visible = False
fig.canvas.header_visible = False
fig.canvas.toolbar_visible = False

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

ylim = (-0.1, 1.4)


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


q_math = ComplexSqrt(sp.Symbol("q^2")) / (8 * sp.pi)
ax_q.set_title(f"${sp.latex(q_math)}$")
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)}$")
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)$")
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)$")
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,
)

fig.tight_layout()
plt.legend(loc="upper right")
plt.show()