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):
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
:
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\):
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
.
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.
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:
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\):
For \(\rho_a\), we use AmpForm’s PhaseSpaceFactor
:
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
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:
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.
SymPy expressions#
In the following, we attempt to implement Equation (2) using SymPy integral.
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!