ENH: formulate Blatt–Weisskopf with spherical Hankel

.cspell.json

@@ -244,6 +244,7 @@
     "hankel",
     "helicities",
     "helicity",
+    "Hippel",
     "htmlcov",
     "imap",
     "ipynb",

docs/_extend_docstrings.py

@@ -63,12 +63,19 @@ def extend_docstrings() -> None:
 
 
 def extend_BlattWeisskopfSquared() -> None:
-    from ampform.dynamics import BlattWeisskopfSquared
+    from ampform.dynamics.form_factor import BlattWeisskopfSquared, SphericalHankel1
 
-    z = sp.Symbol("z", real=True)
-    L = sp.Symbol("L", integer=True)
+    z = sp.Symbol("z", nonnegative=True, real=True)
+    L = sp.Symbol("L", integer=True, nonnegative=True)
     expr = BlattWeisskopfSquared(z, angular_momentum=L)
-    _append_latex_doit_definition(expr, deep=True, full_width=True)
+    h1lz = SphericalHankel1(L, z)
+    _append_latex_doit_definition(expr)
+    _append_to_docstring(
+        BlattWeisskopfSquared,
+        f"""\n
+    where :math:`{sp.latex(h1lz)}` is defined by :eq:`SphericalHankel1`.
+    """,
+    )
 
 
 def extend_BoostMatrix() -> None:
@@ -472,6 +479,15 @@ def extend_RotationZMatrix() -> None:
     )
 
 
+def extend_SphericalHankel1() -> None:
+    from ampform.dynamics.form_factor import SphericalHankel1
+
+    z = sp.Symbol("z", nonnegative=True, real=True)
+    ell = sp.Symbol(R"\ell", integer=True, nonnegative=True)
+    expr = SphericalHankel1(ell, z)
+    _append_latex_doit_definition(expr)
+
+
 def extend_Theta() -> None:
     from ampform.kinematics.angles import Theta
 

docs/bibliography.bib

@@ -173,19 +173,24 @@ @incollection{ParticleDataGroup:2020ssz
   url = {https://pdg.lbl.gov/2010/reviews/rpp2010-rev-dalitz-analysis-formalism.pdf}
 }
 
-@phdthesis{pychyGekoppeltePartialwellenanalyseAnnihilationen2016,
-  title = {{Gekoppelte Partialwellenanalyse von 𝑝̅𝑝-Annihilationen im Fluge in die Endzustände 𝐾⁺𝐾⁻𝜋⁰, 𝜋⁰𝜋⁰𝜂 und 𝜋⁰𝜂𝜂}},
-  author = {Pychy, Julian},
-  year = {2016},
-  month = apr,
-  url = {https://hss-opus.ub.ruhr-uni-bochum.de/opus4/frontdoor/deliver/index/docId/4943/file/diss.pdf},
-  school = {Ruhr-Universität Bochum}
-}
-
 @misc{Richman:1984gh,
   title = {An {{Experimenter}}'s {{Guide}} to the {{Helicity Formalism}}},
   author = {Richman, Jeffrey D.},
   year = {1984},
   month = jun,
   url = {https://inspirehep.net/literature/202987}
 }
+
+@article{VonHippel:1972fg,
+  title = {Centrifugal-{{Barrier Effects}} in {{Resonance Partial Decay Widths}}, {{Shapes}}, and {{Production Amplitudes}}},
+  author = {{von Hippel}, Frank and Quigg, C.},
+  year = {1972},
+  month = feb,
+  journal = {Physical Review D},
+  volume = {5},
+  number = {3},
+  pages = {624--638},
+  issn = {0556-2821},
+  doi = {10.1103/PhysRevD.5.624},
+  url = {https://link.aps.org/doi/10.1103/PhysRevD.5.624}
+}

docs/usage/dynamics.ipynb

@@ -148,7 +148,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "AmpForm uses Blatt-Weisskopf functions $B_L$ as _barrier factors_ (also called _form factors_, see {class}`.BlattWeisskopfSquared`):"
+    "AmpForm uses Blatt-Weisskopf functions $B_L$ as _barrier factors_ (also called _form factors_, see {class}`.BlattWeisskopfSquared` and **[TR-029](https://compwa.github.io/report/029)**):"
    ]
   },
   {
@@ -161,8 +161,8 @@
    "source": [
     "from ampform.dynamics import BlattWeisskopfSquared\n",
     "\n",
-    "L = sp.Symbol(\"L\", integer=True)\n",
-    "z = sp.Symbol(\"z\", real=True)\n",
+    "L = sp.Symbol(\"L\", integer=True, nonnegative=True)\n",
+    "z = sp.Symbol(\"z\", nonnegative=True, real=True)\n",
     "ff2 = BlattWeisskopfSquared(z, L)"
    ]
   },
@@ -183,9 +183,15 @@
    },
    "outputs": [],
    "source": [
+    "from ampform.dynamics.form_factor import SphericalHankel1\n",
     "from ampform.io import aslatex\n",
     "\n",
-    "Math(aslatex({ff2: ff2.doit()}))"
+    "ell = sp.Symbol(R\"\\ell\", integer=True, nonnegative=True)\n",
+    "exprs = [\n",
+    "    BlattWeisskopfSquared(z, L),\n",
+    "    SphericalHankel1(ell, z),\n",
+    "]\n",
+    "Math(aslatex({e: e.doit(deep=False) for e in exprs}))"
    ]
   },
   {
@@ -205,7 +211,7 @@
    "source": [
     "from ampform.dynamics import BreakupMomentumSquared\n",
     "\n",
-    "m, m_a, m_b, d = sp.symbols(\"m, m_a, m_b, d\")\n",
+    "m, m_a, m_b, d = sp.symbols(\"m, m_a, m_b, d\", nonnegative=True)\n",
     "s = m**2\n",
     "q_squared = BreakupMomentumSquared(s, m_a, m_b)\n",
     "ff2 = BlattWeisskopfSquared(q_squared * d**2, angular_momentum=L)"
@@ -243,7 +249,7 @@
    "source": [
     "plot_domain = np.linspace(0.01, 4, 500)\n",
     "sliders.set_ranges(\n",
-    "    L=(0, 8),\n",
+    "    L=(0, 10),\n",
     "    m_a=(0, 2, 200),\n",
     "    m_b=(0, 2, 200),\n",
     "    d=(0, 5),\n",
@@ -359,7 +365,7 @@
    "source": [
     "from ampform.dynamics import relativistic_breit_wigner\n",
     "\n",
-    "m, m0, w0 = sp.symbols(\"m, m0, Gamma0\")\n",
+    "m, m0, w0 = sp.symbols(\"m, m0, Gamma0\", nonnegative=True)\n",
     "rel_bw = relativistic_breit_wigner(s=m**2, mass0=m0, gamma0=w0)\n",
     "rel_bw"
    ]
@@ -537,7 +543,7 @@
     "sliders.set_ranges(\n",
     "    m0=(0, 4, 200),\n",
     "    Gamma0=(0, 1, 100),\n",
-    "    L=(0, 8),\n",
+    "    L=(0, 10),\n",
     "    m_a=(0, 2, 200),\n",
     "    m_b=(0, 2, 200),\n",
     "    d=(0, 5),\n",

src/ampform/dynamics/__init__.py

@@ -3,14 +3,14 @@
 .. seealso:: :doc:`/usage/dynamics` and :doc:`/usage/dynamics/analytic-continuation`
 """
 
-# cspell:ignore asner mhash
 from __future__ import annotations
 
-from typing import TYPE_CHECKING, Any, ClassVar
+from typing import TYPE_CHECKING, Any
 
 import sympy as sp
 
 # pyright: reportUnusedImport=false
+from ampform.dynamics.form_factor import BlattWeisskopfSquared
 from ampform.dynamics.phasespace import (
     BreakupMomentumSquared,
     EqualMassPhaseSpaceFactor,  # noqa: F401
@@ -27,114 +27,9 @@
     from sympy.printing.latex import LatexPrinter
 
 
-@unevaluated
-class BlattWeisskopfSquared(sp.Expr):
-    # cspell:ignore pychyGekoppeltePartialwellenanalyseAnnihilationen
-    r"""Blatt-Weisskopf function :math:`B_L^2(z)`, up to :math:`L \leq 8`.
-
-    Args:
-        z: Argument of the Blatt-Weisskopf function :math:`B_L^2(z)`. A usual
-            choice is :math:`z = (d q)^2` with :math:`d` the impact parameter and
-            :math:`q` the breakup-momentum (see `.BreakupMomentumSquared`).
-
-        angular_momentum: Angular momentum :math:`L` of the decaying particle.
-
-    Note that equal powers of :math:`z` appear in the nominator and the denominator,
-    while some sources have nominator :math:`1`, instead of :math:`z^L`. Compare for
-    instance Equation (50.27) in :pdg-review:`2021; Resonances; p.9`.
-
-    Each of these cases for :math:`L` has been taken from
-    :cite:`pychyGekoppeltePartialwellenanalyseAnnihilationen2016`, p.59,
-    :cite:`Chung:1995dx`, p.415, and :cite:`Chung:1995dx`. For a good overview of where
-    to use these Blatt-Weisskopf functions, see :cite:`ParticleDataGroup:2020ssz`.
-
-    See also :ref:`usage/dynamics:Form factor`.
-    """
-
-    z: Any
-    angular_momentum: Any
-    _latex_repr_ = R"B_{{{angular_momentum}}}^2\left({z}\right)"
-
-    max_angular_momentum: ClassVar[int | None] = None
-    """Limit the maximum allowed angular momentum :math:`L`.
-
-    This improves performance when :math:`L` is a `~sympy.core.symbol.Symbol` and you
-    are note interested in higher angular momenta.
-    """
-
-    def evaluate(self) -> sp.Expr:
-        z: sp.Expr = self.args[0]  # type: ignore[assignment]
-        angular_momentum: sp.Expr = self.args[1]  # type: ignore[assignment]
-        cases: dict[int, sp.Expr] = {
-            0: sp.S.One,
-            1: 2 * z / (z + 1),
-            2: 13 * z**2 / ((z - 3) * (z - 3) + 9 * z),
-            3: 277 * z**3 / (z * (z - 15) * (z - 15) + 9 * (2 * z - 5) * (2 * z - 5)),
-            4: (
-                12746
-                * z**4
-                / (
-                    (z**2 - 45 * z + 105) * (z**2 - 45 * z + 105)
-                    + 25 * z * (2 * z - 21) * (2 * z - 21)
-                )
-            ),
-            5: (
-                998881
-                * z**5
-                / (z**5 + 15 * z**4 + 315 * z**3 + 6300 * z**2 + 99225 * z + 893025)
-            ),
-            6: (
-                118394977
-                * z**6
-                / (
-                    z**6
-                    + 21 * z**5
-                    + 630 * z**4
-                    + 18900 * z**3
-                    + 496125 * z**2
-                    + 9823275 * z
-                    + 108056025
-                )
-            ),
-            7: (
-                19727003738
-                * z**7
-                / (
-                    z**7
-                    + 28 * z**6
-                    + 1134 * z**5
-                    + 47250 * z**4
-                    + 1819125 * z**3
-                    + 58939650 * z**2
-                    + 1404728325 * z
-                    + 18261468225
-                )
-            ),
-            8: (
-                4392846440677
-                * z**8
-                / (
-                    z**8
-                    + 36 * z**7
-                    + 1890 * z**6
-                    + 103950 * z**5
-                    + 5457375 * z**4
-                    + 255405150 * z**3
-                    + 9833098275 * z**2
-                    + 273922023375 * z
-                    + 4108830350625
-                )
-            ),
-        }
-        return sp.Piecewise(*[
-            (expression, sp.Eq(angular_momentum, value))
-            for value, expression in cases.items()
-            if self.max_angular_momentum is None or value <= self.max_angular_momentum
-        ])
-
-
 @unevaluated
 class EnergyDependentWidth(sp.Expr):
+    # cspell:ignore asner
     r"""Mass-dependent width, coupled to the pole position of the resonance.
 
     See Equation (50.28) in :pdg-review:`2021; Resonances; p.9` and

src/ampform/dynamics/form_factor.py

@@ -0,0 +1,102 @@
+"""Implementations of the form factor, or barrier factor."""
+
+from __future__ import annotations
+
+from functools import lru_cache
+from typing import Any
+
+import sympy as sp
+
+from ampform.sympy import unevaluated
+
+
+@unevaluated
+class BlattWeisskopfSquared(sp.Expr):
+    r"""Normalized Blatt-Weisskopf function :math:`B_L^2(z)`, with :math:`B_L^2(1)=1`.
+
+    Args:
+        z: Argument of the Blatt-Weisskopf function :math:`B_L^2(z)`. A usual
+            choice is :math:`z = (d q)^2` with :math:`d` the impact parameter and
+            :math:`q` the breakup-momentum (see `.BreakupMomentumSquared`).
+
+        angular_momentum: Angular momentum :math:`L` of the decaying particle.
+
+    Note that equal powers of :math:`z` appear in the nominator and the denominator,
+    while some sources define an *non-normalized* form factor :math:`F_L` with :math:`1`
+    in the nominator, instead of :math:`z^L`. See for instance Equation (50.27) in
+    :pdg-review:`2021; Resonances; p.9`. We normalize the form factor such that
+    :math:`B_L^2(1)=1` and that :math:`B_L^2` is unitless no matter what :math:`z` is.
+
+    .. seealso:: :ref:`usage/dynamics:Form factor`, :doc:`TR-029<compwa:report/029>`,
+      and :cite:`chungFormulasAngularMomentumBarrier2015`.
+
+    With this, the implementation becomes
+    """
+
+    z: Any
+    angular_momentum: Any
+    _latex_repr_ = R"B_{{{angular_momentum}}}^2\left({z}\right)"
+
+    def evaluate(self) -> sp.Expr:
+        ell = self.angular_momentum
+        z = sp.Dummy("z", nonnegative=True, real=True)
+        expr = (
+            sp.Abs(SphericalHankel1(ell, 1)) ** 2
+            / sp.Abs(SphericalHankel1(ell, sp.sqrt(z))) ** 2
+            / z
+        )
+        if not ell.free_symbols:
+            expr = expr.doit().simplify()
+        return expr.xreplace({z: self.z})
+
+
+@unevaluated
+class SphericalHankel1(sp.Expr):
+    r"""Spherical Hankel function of the first kind for real-valued :math:`z`.
+
+    See :cite:`VonHippel:1972fg`, Equation (A12), and :doc:`TR-029<compwa:report/029>`
+    for more info. `This page
+    <https://mathworld.wolfram.com/SphericalHankelFunctionoftheFirstKind.html>`_
+    explains the difference with the *general* Hankel function of the first kind,
+    :math:`H_\ell^{(1)}`.
+
+    This expression class assumes that :math:`z` is real and evaluates to the following
+    series:
+    """
+
+    l: Any  # noqa: E741
+    z: Any
+    _latex_repr_ = R"h_{{{l}}}^{{(1)}}\left({z}\right)"
+
+    def evaluate(self) -> sp.Expr:
+        l, z = self.args  # noqa: E741
+        k = sp.Dummy("k", integer=True, nonnegative=True)
+        return (
+            (-sp.I) ** (1 + l)  # type:ignore[operator]
+            * (sp.exp(z * sp.I) / z)
+            * _SymbolicSum(
+                sp.factorial(l + k)
+                / (sp.factorial(l - k) * sp.factorial(k))
+                * (sp.I / (2 * z)) ** k,  # type:ignore[operator]
+                (k, 0, l),
+            )
+        )
+
+
+class _SymbolicSum(sp.Sum):
+    """See [TR-029](https://compwa.github.io/report/029.html) for why this class is needed."""
+
+    def doit(self, deep: bool = True, **kwargs) -> sp.Expr:
+        if _get_indices(self):
+            expression = self.args[0]
+            indices = self.args[1:]
+            return _SymbolicSum(expression.doit(deep=deep, **kwargs), *indices)
+        return super().doit(deep=deep, **kwargs)
+
+
+@lru_cache(maxsize=None)
+def _get_indices(expr: sp.Sum) -> set[sp.Basic]:
+    free_symbols = set()
+    for index in expr.args[1:]:
+        free_symbols.update(index.free_symbols)
+    return {s for s in free_symbols if not isinstance(s, sp.Dummy)}