Chew-Mandelstam

Chew-Mandelstam#

This report is an attempt formulate the Chew-Mandelstam function described in PDG2023, §50.3.3 (Resonances), pp.14–15 with SymPy, so that it can be implemented in AmpForm.

S-wave#

As can be seen in Eq. (50.44) on PDG2023, §Resonances, p.15, the Chew-Mandelstam function \(\Sigma_a\) for a particle \(a\) decaying to particles \(1, 2\) has a simple form for angular momentum \(L=0\) (\(S\)-wave):

(1)#\[\Sigma_a(s) = \frac{1}{16\pi^2} \left[ \frac{2q_a}{\sqrt{s}} \log\frac{m_1^2+m_2^2-s+2\sqrt{s}q_a}{2m_1m_2} - \left(m_1^2-m_2^2\right) \left(\frac{1}{s}-\frac{1}{(m_1+m_2)^2}\right) \log\frac{m_1}{m_2} \right]\]

The only question is how to deal with negative values for the squared break-up momentum \(q_a^2\). Here, we will use AmpForm’s ComplexSqrt:

\[\begin{split}\displaystyle \begin{array}{rcl} \sqrt[\mathrm{c}]{z} &=& \begin{cases} i \sqrt{- z} & \text{for}\: z < 0 \\\sqrt{z} & \text{otherwise} \end{cases} \\ \end{array}\end{split}\]

\[\begin{split}\displaystyle \begin{array}{rcl} \Sigma\left(s\right) &=& \frac{1}{16 \pi^{2}} \left(\frac{2 q_a\left(s\right)}{\sqrt{s}} \log{\left(\frac{m_{1}^{2} + m_{2}^{2} + 2 \sqrt{s} q_a\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) \\ q_a\left(s\right) &=& \sqrt[\mathrm{c}]{q^2\left(s\right)} \\ \end{array}\end{split}\]

It should be noted that this equation is not well-defined along the real axis, that is, for \(\mathrm{Im}(s) = 0\). For this reason, we split \(s\) into a real part \(s'\) with a small imaginary offset (the PDG indicates this with \(s+0i\)). We parametrized this imaginary offset with \(\epsilon\), and for the interactive plot, we do so with a power of \(10\):

\[\displaystyle s \to s^{\prime} + 10^{- \epsilon} i\]

We are now ready to use express the symbolic expressions above as a numerical function. For comparison, we will plot the Chew-Mandelstam function for \(S\)-waves next to AmpForm’s PhaseSpaceFactorComplex.

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

symbols = (s_prime, m1, m2, epsilon)
cm_func = sp.lambdify(symbols, cm_expr.doit().subs(s, s_plus))
rho_func = sp.lambdify(symbols, rho_expr.doit().subs(s, s_plus))

As starting values for the interactive plot, we assume \(\pi\eta\) scattering (just like in the PDG section) and use their masses as values for \(m_1\) and \(m_1\), respectively.

Show code cell source Hide code cell source

fig, axes = plt.subplots(figsize=(11, 4.5), ncols=2, tight_layout=True)
ax1, ax2 = axes
for ax in axes:
    ax.axhline(0, color="black", linewidth=0.5)

real_style = {"label": "Real part", "c": "black", "linestyle": "dashed"}
imag_style = {"label": "Imag part", "c": "red"}
threshold_style = {"label": R"$s_\mathrm{thr}$", "c": "grey", "linewidth": 0.5}

xlim = (plot_domain.min(), plot_domain.max())
ylim = (-1, +1)
y_factor = 16 * np.pi
controls = iplt.axvline(
    lambda *args, **kwargs: (kwargs["m1"].value + kwargs["m2"].value) ** 2,
    **sliders,
    ax=ax1,
    **threshold_style,
)
iplt.axvline(
    lambda *args, **kwargs: (kwargs["m1"].value + kwargs["m2"].value) ** 2,
    controls=controls,
    ax=ax2,
    **threshold_style,
)
to_values = lambda kwargs: {k: v.value for k, v in kwargs.items()}
iplt.plot(
    plot_domain,
    lambda *args, **kwargs: (y_factor * 1j * rho_func(*args, **to_values(kwargs))).real,
    controls=controls,
    ylim=ylim,
    xlim=xlim,
    alpha=0.7,
    ax=ax1,
    **real_style,
)
iplt.plot(
    plot_domain,
    lambda *args, **kwargs: (y_factor * 1j * rho_func(*args, **to_values(kwargs))).imag,
    controls=controls,
    ylim=ylim,
    xlim=xlim,
    alpha=0.7,
    ax=ax1,
    **imag_style,
)

iplt.plot(
    plot_domain,
    lambda *args, **kwargs: y_factor * cm_func(*args, **to_values(kwargs)).real,
    controls=controls,
    ylim=ylim,
    xlim=xlim,
    alpha=0.7,
    ax=ax2,
    **real_style,
)
iplt.plot(
    plot_domain,
    lambda *args, **kwargs: y_factor * cm_func(*args, **to_values(kwargs)).imag,
    controls=controls,
    ylim=ylim,
    xlim=xlim,
    alpha=0.7,
    ax=ax2,
    **imag_style,
)

for ax in axes:
    ax.legend(loc="lower right")
    ax.set_xticks(np.arange(0, 1.21, 0.3))
    ax.set_yticks(np.arange(-1, 1.1, 0.5))
    ax.set_xlim()
    ax.set_xlabel("$s$ (GeV$^2$)")

ax1.set_ylabel(R"$16\pi \; i\rho(s)$")
ax2.set_ylabel(R"$16\pi \; \Sigma(s)$")
ax1.set_title(R"Complex phase space factor $\rho$")
ax2.set_title("Chew-Mandelstam $S$-wave ($L=0$)")

fig.savefig("chew-mandelstam-s-wave.svg")
plt.show()

https://github.com/user-attachments/assets/3454e711-8b7e-4ee4-90a6-d6c30bb9d3f8

Tip

Compare the plots above with Figure 50.6 on PDG2023, §Resonances, p.16.

General dispersion integral#

For higher angular momenta, the PDG notes that one has to compute the dispersion integral given by Eq. (50.44) on PDG2023, §Resonances, p.15:

(2)#\[ \Sigma_a(s+0i) = \frac{s-s_{\mathrm{thr}_a}}{\pi} \int^\infty_{s_{\mathrm{thr}_a}} \frac{ \rho_a(s')n_a^2(s') }{ (s' - s_{\mathrm{thr}_a})(s'-s-i0) } \mathop{}\!\mathrm{d}s' \]

Equation (1) is the analytic solution for \(L=0\).

From Equations (50.33–34) on PDG2023, §Resonances, p.12, it can be deduced that the function \(n_a^2\) is the same as AmpForm’s BlattWeisskopfSquared (note that this function is normalized, whereas the PDG’s \(F_j\) function has \(1\) in the nominator). For this reason, we simply use BlattWeisskopfSquared for the definition of \(n_a^2\):

\[\begin{split}\displaystyle \begin{array}{rcl} n_a^2\left(s\right) &=& B_{L}^2\left(\frac{q^2\left(s\right)}{q_{0}^{2}}\right) \\ B_{L}^2\left(z\right) &=& \begin{cases} \frac{2 z}{z + 1} & \text{for}\: L = 1 \\\frac{13 z^{2}}{9 z + \left(z - 3\right)^{2}} & \text{for}\: L = 2 \\\frac{277 z^{3}}{z \left(z - 15\right)^{2} + \left(2 z - 5\right) \left(18 z - 45\right)} & \text{for}\: L = 3 \\\frac{12746 z^{4}}{25 z \left(2 z - 21\right)^{2} + \left(z^{2} - 45 z + 105\right)^{2}} & \text{for}\: L = 4 \\\frac{998881 z^{5}}{z^{5} + 15 z^{4} + 315 z^{3} + 6300 z^{2} + 99225 z + 893025} & \text{for}\: L = 5 \\\frac{118394977 z^{6}}{z^{6} + 21 z^{5} + 630 z^{4} + 18900 z^{3} + 496125 z^{2} + 9823275 z + 108056025} & \text{for}\: L = 6 \\\frac{19727003738 z^{7}}{z^{7} + 28 z^{6} + 1134 z^{5} + 47250 z^{4} + 1819125 z^{3} + 58939650 z^{2} + 1404728325 z + 18261468225} & \text{for}\: L = 7 \\\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} & \text{for}\: L = 8 \end{cases} \\ \end{array}\end{split}\]

For \(\rho_a\), we use AmpForm’s PhaseSpaceFactor:

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

The symbolic integrand is then formulated as:

s_thr = (m1 + m2) ** 2
integrand = (PhaseSpaceFactor(s_prime, m1, m2) * FormFactor(s_prime, m1, m2, L, q0)) / (
    (s_prime - s_thr) * (s_prime - s - epsilon * sp.I)
)
integrand

\[\displaystyle \frac{n_a^2\left(s^{\prime}\right) \rho\left(s^{\prime}\right)}{\left(s^{\prime} - \left(m_{1} + m_{2}\right)^{2}\right) \left(- i \epsilon - s + s^{\prime}\right)}\]

Next, we lambdify() this integrand to a numpy expression so that we can integrate it efficiently:

integrand_func = sp.lambdify(
    args=(s_prime, s, L, epsilon, m1, m2, q0),
    expr=integrand.doit(),
    modules="numpy",
)

Note

Integrals can be expressed symbolically with SymPy, with some caveats. See SymPy integral.

As discussed in TR-016, scipy.integrate.quad() cannot integrate over complex-valued functions, but scipy.integrate.quad_vec() can. For comparison, we now compute this integral for a few values of \(L>0\):

s_domain = np.linspace(s_min, s_max, num=50)
max_L = 3
l_values = list(range(1, max_L + 1))

It is handy to store the numerical results of each dispersion integral in a dict with \(L\) as keys:

s_thr_val = float(s_thr.subs({m1: m1_val, m2: m2_val}))
integral_values = {
    l_val: quad_vec(
        lambda x: integrand_func(
            x,
            s=s_domain,
            L=l_val,
            epsilon=1e-3,
            m1=m1_val,
            m2=m2_val,
            q0=1.0,
        ),
        a=s_thr_val,
        b=np.inf,
    )[0]
    for l_val in tqdm(l_values, desc="Evaluating integrals")
}

Finally, as can be seen from Eq. (2), the resulting values from the integral have to be shifted with a factor \(\frac{s-s_{\mathrm{thr}_a}}{\pi}\) to get \(\Sigma_a\). We also scale the values with \(16\pi\) so that it can be compared with the plot generated in S-wave.

sigma = {
    l_val: (s_domain - s_thr_val) / np.pi * integral_values[l_val] for l_val in l_values
}
sigma_scaled = {l_val: 16 * np.pi * sigma[l_val] for l_val in l_values}

SymPy expressions#

In the following, we attempt to implement Equation (2) using SymPy integral.

\[\displaystyle \frac{s - \left(m_{1} + m_{2}\right)^{2}}{\pi} \int\limits_{\left(m_{1} + m_{2}\right)^{2}}^{\infty} \frac{B_{L}^2\left(q^2\left(s^{\prime}\right)\right) \rho\left(s^{\prime}\right)}{\left(s^{\prime} - \left(m_{1} + m_{2}\right)^{2}\right) \left(- i \epsilon - s + s^{\prime}\right)}\, ds^{\prime}\]

Warning

We have to keep track of the integration variable (\(s'\) in Equation (2)), so that we don’t run into trouble if we use lambdify() with common sub-expressions. The problem is that the integration variable should not be extracted as a common sub-expression, otherwise the lambdified scipy.integrate.quad_vec() expression cannot handle vectorized input.

To keep the function under the integral simple, we substitute angular momentum \(L\) with a definite value before we lambdify:

UnevaluatableIntegral.abs_tolerance = 1e-4
UnevaluatableIntegral.rel_tolerance = 1e-4
integral_s_wave_func = sp.lambdify(
    [s, m1, m2, epsilon],
    integral_expr.subs(L, 0).doit(),
    # integration symbol should not be extracted as common sub-expression!
    cse=partial(sp.cse, ignore=[s_prime], list=False),
)
integral_s_wave_func = np.vectorize(integral_s_wave_func)

integral_p_wave_func = sp.lambdify(
    [s, m1, m2, epsilon],
    integral_expr.subs(L, 1).doit(),
    cse=partial(sp.cse, ignore=[s_prime], list=False),
)
integral_p_wave_func = np.vectorize(integral_p_wave_func)

def _lambdifygenerated(s, m1, m2, epsilon):
    x0 = pi ** (-1.0)
    x1 = (m1 + m2) ** 2
    x2 = -x1
    return (
        x0
        * (s + x2)
        * quad_vec(
            lambda _Dummy_46: (1 / 16)
            * x0
            * sqrt((_Dummy_46 + x2) * (_Dummy_46 - (m1 - m2) ** 2) / _Dummy_46)
            / (sqrt(_Dummy_46) * (_Dummy_46 + x2) * (_Dummy_46 - 1j * epsilon - s)),
            x1,
            PINF,
            epsabs=0.0001,
            epsrel=0.0001,
            limit=50,
        )[0]
    )

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

CPU times: user 751 ms, sys: 201 μs, total: 751 ms
Wall time: 752 ms
CPU times: user 266 ms, sys: 0 ns, total: 266 ms
Wall time: 266 ms

Note that the dispersion integral for \(L=0\) indeed reproduces the same shape as in S-wave!