Blatt–Weisskopf from Hankel function

Blatt–Weisskopf from Hankel function#

Introduction#

As of AmpForm v0.15, the implementation of BlattWeisskopfSquared contains hard-coded polynomials, see implementation here. The motivation for this can be found in the citations mentioned in its API documentation. However, as noted by @mmikhasenko in ComPWA/ampform#417, the polynomials can be derived from the spherical[1] Hankel functions of the first kind. Von Hippel and Quigg[2] derived a generalization of the centrifugal barrier factor \(F_L\), also called form factor, that was introduced by [Blatt and Weisskopf, 1979], showing that

\[F_\ell^2(z^2) = \frac{1}{z^2\left|h^{(1)}_\ell\left(z\right)\right|^2}\,,\]

where \(h_\ell^{(1)}\) is a Hankel function of the first kind. They also noted that, if \(z\in\mathbb{R}\),

(1)#\[ h_\ell^{(1)}(z) = \left(- i\right)^{\ell+1} \frac{e^{iz}}{z} \sum_{k=0}^\ell \frac{(\ell+k)!}{(\ell-k)! \, k!} \left(\frac{i}{2z}\right)^k. \]

In the following, we call \(F_\ell(z)\) the unnormalized Blatt–Weisskopf form factor. Following Chung and other resources (see e.g. [Chung et al., 1995], p. 415), AmpForm implements a unitless, normalized Blatt–Weisskopf factor \(B_L\), meaning that \(B_L(1)=1\).[3] It can be defined in terms of \(F_L\) as

\[B_L^2(z^2) = \frac{F_L^2(z^2)}{F_L^2(1)} = \frac{\left|h^{(1)}_L(1)\right|^2}{z^2\left|h^{(1)}_L(z)\right|^2}\,.\]

Note

As of writing, AmpForm uses \(z\) as argument in BlattWeisskopfSquared. This means we have to work with a square root and assume that \(z \geq 0\), meaning

(2)#\[ B_L^2(z) = \frac{\left|h^{(1)}_L(1)\right|^2}{z\left|h^{(1)}_L\left(\sqrt{z}\right)\right|^2}\,. \]

Hankel function of the first kind#

Built-in SymPy function#

SymPy offers a Hankel function of the first kind, scipy.special.hankel1.

z = sp.Symbol("z", nonnegative=True, real=True)
ell = sp.Symbol(R"\ell", integer=True, nonnegative=True)
sp.hankel1(ell, z)

\[\displaystyle H^{(1)}_{\ell}\left(z\right)\]

This function is the general[1] Hankel function \(H_\ell\) and the class does not offer algebraic simplifications for specific values or assumptions of \(\ell\) and \(z\).

display(
    sp.hankel1(ell, z).doit(),
    sp.hankel1(ell, 0).doit(),
    sp.hankel1(0, z).doit(),
    sp.hankel1(0, 0).doit(),
)

\[\displaystyle H^{(1)}_{\ell}\left(z\right)\]

\[\displaystyle H^{(1)}_{\ell}\left(0\right)\]

\[\displaystyle H^{(1)}_{0}\left(z\right)\]

\[\displaystyle H^{(1)}_{0}\left(0\right)\]

Custom class definition#

To implement Equation (1) for the spherical Hankel function, we have to define a custom @unevaluated expression class. The following class evaluates to the sum given in Equation (1). We introduce a special sympy.Sum class that does not ‘unfold’ on symbolic input for \(\ell\) if doit() is called (see Nested doit).

@unevaluated
class SphericalHankel1(sp.Expr):
    l: sp.Symbol | int
    z: sp.Symbol | float
    _latex_repr_ = R"h_{{{l}}}^{{(1)}}\left({z}\right)"

    def evaluate(self) -> sp.Expr:
        l, z = self.args
        k = sp.Dummy("k", integer=True, nonnegative=True)
        return (
            (-sp.I) ** (1 + l)
            * (sp.exp(z * sp.I) / z)
            * SymbolicSum(
                sp.factorial(l + k)
                / (sp.factorial(l - k) * sp.factorial(k))
                * (sp.I / (2 * z)) ** k,
                (k, 0, l),
            )
        )


class SymbolicSum(sp.Sum):
    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)


@cache
def _get_indices(expr: sp.Sum) -> set[sp.Symbol]:
    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)}

\[\begin{split}\displaystyle \begin{array}{rcl} h_{\ell}^{(1)}\left(z\right) &=& \frac{\left(- i\right)^{\ell + 1} e^{i z} \sum_{k=0}^{\ell} \frac{\left(\frac{i}{2 z}\right)^{k} \left(k + \ell\right)!}{k! \left(- k + \ell\right)!}}{z} \\ \end{array}\end{split}\]

Indeed, the absolute squared value \(\left|h_\ell^{(1)}\right|^2\) results in a clean fraction of polynomials (after some algebraic simplifications).

\[\begin{split}\displaystyle \begin{array}{rcl} \left|{h_{0}^{(1)}\left(z\right)}\right|^{2} &=& \frac{1}{z^{2}} \\ \left|{h_{1}^{(1)}\left(z\right)}\right|^{2} &=& \frac{z^{2} + 1}{z^{4}} \\ \left|{h_{2}^{(1)}\left(z\right)}\right|^{2} &=& \frac{z^{4} + 3 z^{2} + 9}{z^{6}} \\ \end{array}\end{split}\]

\[\begin{split}\displaystyle \begin{array}{rcl} \left|{h_{0}^{(1)}\left(1\right)}\right|^{2} &=& 1 \\ \left|{h_{1}^{(1)}\left(1\right)}\right|^{2} &=& 2 \\ \left|{h_{2}^{(1)}\left(1\right)}\right|^{2} &=& 13 \\ \end{array}\end{split}\]

Normalized Blatt–Weisskopf form factor#

We now have the required expression classes for re-implementing BlattWeisskopfSquared using Equation (2) (with \(z\) as input, instead of \(z^2\)).

@unevaluated
class BlattWeisskopfSquared(sp.Expr):
    L: sp.Symbol | int
    z: sp.Symbol | float
    _latex_repr_ = R"B^2_{{{L}}}\left({z}\right)"

    def evaluate(self) -> sp.Expr:
        L = self.L
        z = sp.Dummy("z", nonnegative=True, real=True)
        expr = (
            sp.Abs(SphericalHankel1(L, 1)) ** 2
            / sp.Abs(SphericalHankel1(L, sp.sqrt(z))) ** 2
            / z
        )
        if not L.free_symbols:
            expr = expr.doit().simplify()
        return expr.xreplace({z: self.z})

Note

An explicit simplify() is required in order to reproduce the polynomial form upon evaluation. To make the simplification as fast as possible, it is done internally within evaluate() with \(z\) as a dummy variable. This is to avoid performing nested simplifications if \(z\) is in itself an expression (see Nested doit).

\[\begin{split}\displaystyle \begin{array}{rcl} B^2_{L}\left(z\right) &=& \frac{\left|{h_{L}^{(1)}\left(1\right)}\right|^{2}}{z \left|{h_{L}^{(1)}\left(\sqrt{z}\right)}\right|^{2}} \\ \end{array}\end{split}\]

Indeed the polynomials are exactly the same as the original BlattWeisskopfSquared!

\[\begin{split}\displaystyle \begin{array}{rcl} B^2_{0}\left(z\right) &=& 1 \\ B^2_{1}\left(z\right) &=& \frac{2 z}{z + 1} \\ B^2_{2}\left(z\right) &=& \frac{13 z^{2}}{z^{2} + 3 z + 9} \\ B^2_{3}\left(z\right) &=& \frac{277 z^{3}}{z^{3} + 6 z^{2} + 45 z + 225} \\ B^2_{4}\left(z\right) &=& \frac{12746 z^{4}}{z^{4} + 10 z^{3} + 135 z^{2} + 1575 z + 11025} \\ B^2_{5}\left(z\right) &=& \frac{998881 z^{5}}{z^{5} + 15 z^{4} + 315 z^{3} + 6300 z^{2} + 99225 z + 893025} \\ B^2_{6}\left(z\right) &=& \frac{118394977 z^{6}}{z^{6} + 21 z^{5} + 630 z^{4} + 18900 z^{3} + 496125 z^{2} + 9823275 z + 108056025} \\ B^2_{7}\left(z\right) &=& \frac{19727003738 z^{7}}{z^{7} + 28 z^{6} + 1134 z^{5} + 47250 z^{4} + 1819125 z^{3} + 58939650 z^{2} + 1404728325 z + 18261468225} \\ B^2_{8}\left(z\right) &=& \frac{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} \\ \end{array}\end{split}\]

Nested doit#

Eventually, the barrier factors take \(z=q/q_R\), with \(q\) the break-up momentum and \(q_R\) an impact factor. Here it becomes crucial that only \(\left|h_\ell^{(1)}(z)\right|^2\) is simplified to a polynomial fraction, not \(q\) itself. The break-up momentum does need to unfold though.

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

Symbolic angular momentum#

BlattWeisskopfSquared(L, z=q2 / qR**2).doit(deep=False)

\[\displaystyle \frac{q_{R}^{2} \left|{h_{L}^{(1)}\left(1\right)}\right|^{2}}{\left|{h_{L}^{(1)}\left(\sqrt{\frac{q^2\left(s\right)}{q_{R}^{2}}}\right)}\right|^{2} q^2\left(s\right)}\]

BlattWeisskopfSquared(L, z=q2 / qR**2).doit(deep=True)

\[\displaystyle \frac{q_{R}^{2} s e^{\sqrt[4]{\frac{1}{q_{R}^{4} s^{2}}} \sin{\left(\frac{\operatorname{atan}_{2}{\left(0,\frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{q_{R}^{2} s} \right)}}{2} \right)} \sqrt{\left|{s - \left(m_{1} - m_{2}\right)^{2}}\right|} \sqrt{\left|{s - \left(m_{1} + m_{2}\right)^{2}}\right|}} \left|{\sqrt{\frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{q_{R}^{2} s}}}\right|^{2} \left|{\sum_{k=0}^{L} \frac{\left(\frac{i}{2}\right)^{k} \left(k + L\right)!}{k! \left(- (k - L)\right)!}}\right|^{2}}{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right) \left|{\sum_{k=0}^{L} \frac{\left(\frac{i}{\sqrt{\frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{q_{R}^{2} s}}}\right)^{k} \left(k + L\right)!}{k! \left(- (k - L)\right)!}}\right|^{2}}\]

Numeric angular momentum#

BlattWeisskopfSquared(L=2, z=q2 / qR**2).doit(deep=False)

\[\displaystyle \frac{13 q^2\left(s\right)^{2}}{q_{R}^{4} \left(9 + \frac{3 q^2\left(s\right)}{q_{R}^{2}} + \frac{q^2\left(s\right)^{2}}{q_{R}^{4}}\right)}\]

BlattWeisskopfSquared(L=2, z=q2 / qR**2).doit(deep=True)

\[\displaystyle \frac{13 \left(s - \left(m_{1} - m_{2}\right)^{2}\right)^{2} \left(s - \left(m_{1} + m_{2}\right)^{2}\right)^{2}}{16 q_{R}^{4} s^{2} \left(9 + \frac{3 \left(s - \left(m_{1} - m_{2}\right)^{2}\right) \left(s - \left(m_{1} + m_{2}\right)^{2}\right)}{4 q_{R}^{2} s} + \frac{\left(s - \left(m_{1} - m_{2}\right)^{2}\right)^{2} \left(s - \left(m_{1} + m_{2}\right)^{2}\right)^{2}}{16 q_{R}^{4} s^{2}}\right)}\]