FEAT: implement general dispersion integral (#480)

* BREAK: move break-up momentum to `kinematics` module
* DOC: update link to Chung paper

Author: Zeyna777

.cspell.json

@@ -251,6 +251,7 @@
     "hlines",
     "htmlcov",
     "imap",
+    "infty",
     "ipynb",
     "ipython",
     "isort",
@@ -289,6 +290,7 @@
     "topness",
     "unevaluatable",
     "unitarity",
+    "unphysical",
     "venv",
     "weisskopf",
     "yticklabels",

docs/_extend_docstrings.py

@@ -21,18 +21,20 @@
 from sympy.printing.numpy import NumPyPrinter
 
 from ampform.dynamics.phasespace import (
-    BreakupMomentum,
-    BreakupMomentumComplex,
-    BreakupMomentumKallen,
-    BreakupMomentumSplitSqrt,
-    BreakupMomentumSquared,
     PhaseSpaceFactor,
     PhaseSpaceFactorKallen,
     PhaseSpaceFactorSplitSqrt,
     PhaseSpaceFactorSWave,
 )
 from ampform.io import aslatex
 from ampform.kinematics.lorentz import ArraySize, FourMomentumSymbol
+from ampform.kinematics.phasespace import (
+    BreakupMomentum,
+    BreakupMomentumComplex,
+    BreakupMomentumKallen,
+    BreakupMomentumSplitSqrt,
+    BreakupMomentumSquared,
+)
 from ampform.sympy._array_expressions import ArrayMultiplication
 from ampform.sympy._cache import get_readable_hash, make_hashable
 
@@ -204,18 +206,12 @@ def extend_BoostZMatrix() -> None:
     )
 
 
-def extend_PhaseSpaceFactorComplex() -> None:
-    from ampform.dynamics.phasespace import PhaseSpaceFactorComplex
+def extend_ChewMandelstamIntegral() -> None:
+    from ampform.dynamics.phasespace import ChewMandelstamIntegral
 
-    s, m_a, m_b = sp.symbols("s, m_a, m_b")
-    expr = PhaseSpaceFactorComplex(s, m_a, m_b)
+    s, m_a, m_b, ell = sp.symbols(R"s m_a m_b \ell")
+    expr = ChewMandelstamIntegral(s, m_a, m_b, ell)
     _append_latex_doit_definition(expr)
-    _append_to_docstring(
-        PhaseSpaceFactorComplex,
-        """
-    with :math:`q^2(s)` defined as :eq:`BreakupMomentumSquared`.
-    """,
-    )
 
 
 def extend_ComplexSqrt() -> None:
@@ -409,6 +405,20 @@ def extend_PhaseSpaceFactorAbs() -> None:
     )
 
 
+def extend_PhaseSpaceFactorComplex() -> None:
+    from ampform.dynamics.phasespace import PhaseSpaceFactorComplex
+
+    s, m_a, m_b = sp.symbols("s, m_a, m_b")
+    expr = PhaseSpaceFactorComplex(s, m_a, m_b)
+    _append_latex_doit_definition(expr)
+    _append_to_docstring(
+        PhaseSpaceFactorComplex,
+        """
+    with :math:`q^2(s)` defined as :eq:`BreakupMomentumSquared`.
+    """,
+    )
+
+
 def extend_Phi() -> None:
     from ampform.kinematics.angles import Phi
 

docs/bibliography.bib

@@ -64,7 +64,7 @@ @misc{chungPrimerKmatrixFormalism1995
   author = {Chung, Suh-Urk},
   year = {1995},
   month = mar,
-  url = {http://www.ep.ph.bham.ac.uk/exp/WA102/pics/chung5.ps.gz}
+  url = {https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=88b101a5300736f78293cf10116c32e5d25e3c91}
 }
 
 @techreport{chungSpinFormalismsUpdated2014,

docs/conf.py

@@ -288,6 +288,7 @@ def _get_dataclasses(module):
 linkcheck_anchors = False
 linkcheck_ignore = [
     "http://www.curtismeyer.com",
+    "https://citeseerx.ist.psu.edu/document",
     "https://doi.org/10.1002",  # 403 for onlinelibrary.wiley.com
     "https://doi.org/10.1093",  # 403 for PTEP
     "https://doi.org/10.1103",  # 403 for Phys Rev D
@@ -344,6 +345,7 @@ def _get_dataclasses(module):
     ".rst": "restructuredtext",
 }
 suppress_warnings = [
+    "myst.directive_unknown",
     "myst.domains",
     # skipping unknown output mime type: application/json
     # https://github.com/ComPWA/ampform/runs/8132373732?check_suite_focus=true#step:5:127

docs/usage/dynamics.ipynb

@@ -77,6 +77,7 @@
     "dynamics/custom\n",
     "dynamics/analytic-continuation\n",
     "dynamics/k-matrix\n",
+    "dynamics/riemann-sheets\n",
     "```\n",
     "\n",
     "```{autolink-skip}\n",
@@ -674,7 +675,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.9"
+   "version": "3.13.11"
   }
  },
  "nbformat": 4,

docs/usage/dynamics/analytic-continuation.ipynb

@@ -99,8 +99,9 @@
     "from IPython.display import Markdown, Math\n",
     "from matplotlib_inline.backend_inline import set_matplotlib_formats\n",
     "\n",
-    "from ampform.dynamics.phasespace import BreakupMomentum\n",
+    "from ampform.dynamics.phasespace import ChewMandelstamIntegral\n",
     "from ampform.io import aslatex\n",
+    "from ampform.kinematics.phasespace import BreakupMomentum\n",
     "\n",
     "\n",
     "def display_doit(exprs: list[sp.Expr], deep: bool = False) -> Math:\n",
@@ -365,9 +366,11 @@
     "S = X + 1j * Y\n",
     "ϵi = 1e-7j\n",
     "\n",
-    "parameters = {m1: 0.2, m2: 0.6}\n",
-    "thr_neg = (parameters[m1] - parameters[m2]) ** 2\n",
-    "thr_pos = (parameters[m1] + parameters[m2]) ** 2\n",
+    "m1_val = 0.2\n",
+    "m2_val = 0.6\n",
+    "parameters = {m1: m1_val, m2: m2_val}\n",
+    "thr_neg = (m1_val - m2_val) ** 2\n",
+    "thr_pos = (m1_val + m2_val) ** 2\n",
     "\n",
     "fig, axes = plt.subplots(\n",
     "    dpi=200,\n",
@@ -546,6 +549,114 @@
     "})"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Dispersion integral"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To get an analytic phasespace factor for higher angular momenta, the one has to compute the dispersion integral. According to [PDG 2025, §50. Resonances, Eq. (50.45)](https://pdg.lbl.gov/2025/reviews/rpp2025-rev-resonances.pdf#page=16) the once-substracted dispersion integral is given by:\n",
+    "\n",
+    "$$\n",
+    "\\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}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L = sp.Symbol(\"L\", integer=True, nonnegative=True)\n",
+    "integral_expr = ChewMandelstamIntegral(s, m1, m2, L)\n",
+    "integral_expr.doit(deep=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "integral_s_wave_func = sp.lambdify(\n",
+    "    [s, m1, m2, integral_expr.epsilon],\n",
+    "    integral_expr.subs(L, 0).doit(),\n",
+    ")\n",
+    "integral_s_wave_func = np.vectorize(integral_s_wave_func)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "mystnb": {
+     "code_prompt_show": "Same definition for P-wave function"
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "integral_p_wave_func = sp.lambdify(\n",
+    "    [s, m1, m2, integral_expr.epsilon],\n",
+    "    integral_expr.subs(L, 1).doit(),\n",
+    ")\n",
+    "integral_p_wave_func = np.vectorize(integral_p_wave_func)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "s_values = np.linspace(-0.15, 1.4, num=200)\n",
+    "s_wave_values = integral_s_wave_func(s_values, m1_val, m2_val, epsilon=1e-5)\n",
+    "p_wave_values = integral_p_wave_func(s_values, m1_val, m2_val, epsilon=1e-5)\n",
+    "\n",
+    "l_val = [0, 1]\n",
+    "fig, axes = plt.subplots(figsize=(6, 7), nrows=2, sharex=True)\n",
+    "fig.patch.set_facecolor(\"none\")\n",
+    "ax1, ax2 = axes\n",
+    "fig.suptitle(f\"Symbolic dispersion integrals for $m_1={m1_val:.2f}, m_2={m2_val:.2f}$\")\n",
+    "for ax in axes:\n",
+    "    ax.axhline(0, linewidth=0.5, c=\"black\")\n",
+    "    ax.axvline(thr_pos, linestyle=\"--\", color=\"red\", alpha=0.7)\n",
+    "    ax.patch.set_facecolor(\"none\")\n",
+    "    ax.set_title(f\"$L = {l_val}$\")\n",
+    "    ax.set_ylabel(R\"$16\\pi \\; \\Sigma(s)$\")\n",
+    "axes[-1].set_xlabel(\"$s$ (GeV$^2$)\")\n",
+    "\n",
+    "ax1.set_title(\"$S$-wave ($L=0$)\")\n",
+    "ax1.plot(s_values, 16 * np.pi * s_wave_values.real, label=\"real\")\n",
+    "ax1.plot(s_values, 16 * np.pi * s_wave_values.imag, label=\"imaginary\")\n",
+    "\n",
+    "ax2.set_title(\"$P$-wave ($L=1$)\")\n",
+    "ax2.plot(s_values, 16 * np.pi * p_wave_values.real, label=\"real\")\n",
+    "ax2.plot(s_values, 16 * np.pi * p_wave_values.imag, label=\"imaginary\")\n",
+    "\n",
+    "ax1.legend(framealpha=0)\n",
+    "fig.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -813,7 +924,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.9"
+   "version": "3.13.11"
   }
  },
  "nbformat": 4,