D^+\to\pi^+\pi^+\pi^- with Dalitz-plot decomposition

\(D^+\to\pi^+\pi^+\pi^-\) with Dalitz-plot decomposition#

This report demonstrates a Bose-symmetric isobar model for \(D^+\to\pi^+\pi^+\pi^-\) using QRules for transition generation, AmpForm-DPD for symbolic helicity amplitudes, and TensorWaves for numerical evaluation. The amplitude is checked for Bose symmetry and used to plot an illustrative Dalitz density.

The resonance content follows two published analyses:

CLEO [CLEO Collaboration et al., 2007], Table II: \(\rho(770)\), \(f_2(1270)\), \(f_0(980)\), \(f_0(1370)\), \(f_0(1500)\), and a low-mass \(\sigma\) component, represented here by \(f_0(500)\).
LHCb [The LHCb collaboration et al., 2023], Table 2: also includes \(\omega(782)\), \(\rho(1450)\), and \(\rho(1700)\). Its S wave is determined with a quasi-model-independent parametrization, rather than separate scalar isobars. States tested without significant improvement, such as \(\rho_3(1690)\), are not included here.

This nine-resonance example combines CLEO’s scalar-isobar content with the additional LHCb vector states. It retains the package Breit–Wigner dynamics, masses, widths, and radii. CLEO instead uses a complex \(\sigma\) pole and a Flatté function for \(f_0(980)\); LHCb uses Gounaris–Sakurai functions for the \(\rho\) states. The S wave, especially near the \(K\bar K\) threshold, therefore remains an approximation. The source fit-fraction ratios set the relative component strengths below; the result is a literature-informed illustration, not a reproduction of either fit. No CP violation is assumed.

Generate the decay chains#

Label \(D^+\) by 0 and \((\pi^+_1,\pi^+_2,\pi^-_3)\) by \((1,2,3)\). The charm decay is weak, so imposing strong-interaction conservation on the entire chain would incorrectly forbid it. We allow weak interactions in QRules; the selected neutral resonances also have strong-decay-compatible \(J^P\) for their two-pion daughters.

resonance	\(J^P\)	mass (MeV)	width (MeV)	\(L_\mathrm{dec}^\mathrm{min}\)	\(L_\mathrm{prod}^\mathrm{min}\)
\(f_{0}(1370) \to \pi^{+} \pi^{-}\)	\(0^+\)	1,350	350	0	0
\(f_{0}(1370) \to \pi^{+} \pi^{-}\)	\(0^+\)	1,350	350	0	0
\(f_{0}(1500) \to \pi^{+} \pi^{-}\)	\(0^+\)	1,522	108	0	0
\(f_{0}(1500) \to \pi^{+} \pi^{-}\)	\(0^+\)	1,522	108	0	0
\(f_{0}(500) \to \pi^{+} \pi^{-}\)	\(0^+\)	600	450	0	0
\(f_{0}(500) \to \pi^{+} \pi^{-}\)	\(0^+\)	600	450	0	0
\(f_{0}(980) \to \pi^{+} \pi^{-}\)	\(0^+\)	990	60	0	0
\(f_{0}(980) \to \pi^{+} \pi^{-}\)	\(0^+\)	990	60	0	0
\(f_{2}(1270) \to \pi^{+} \pi^{-}\)	\(2^+\)	1,275	186	2	2
\(f_{2}(1270) \to \pi^{+} \pi^{-}\)	\(2^+\)	1,275	186	2	2
\(\omega(782) \to \pi^{+} \pi^{-}\)	\(1^-\)	782	8	1	1
\(\omega(782) \to \pi^{+} \pi^{-}\)	\(1^-\)	782	8	1	1
\(\rho(1450)^{0} \to \pi^{+} \pi^{-}\)	\(1^-\)	1,465	400	1	1
\(\rho(1450)^{0} \to \pi^{+} \pi^{-}\)	\(1^-\)	1,465	400	1	1
\(\rho(1700)^{0} \to \pi^{+} \pi^{-}\)	\(1^-\)	1,720	250	1	1
\(\rho(1700)^{0} \to \pi^{+} \pi^{-}\)	\(1^-\)	1,720	250	1	1
\(\rho(770)^{0} \to \pi^{+} \pi^{-}\)	\(1^-\)	775	147	1	1
\(\rho(770)^{0} \to \pi^{+} \pi^{-}\)	\(1^-\)	775	147	1	1

        flowchart LR
    T0_N0["$$D^{+}$$"]
    T0_1["$$1: \pi^{+}$$"]
    T0_2["$$2: \pi^{+}$$"]
    T0_3["$$3: \pi^{-}$$"]
    T0_N1@{ shape: text, label: " " }
    T0_4("$$\begin{array}{ll} f_{0}(500) & \omega(782) \\\ f_{0}(980) & \rho(770)^{0} \\\ f_{2}(1270) & \rho(1450)^{0} \\\ f_{0}(1370) & \rho(1700)^{0} \\\ f_{0}(1500) &  \end{array}$$")
    T0_N0 --- T0_4
    T0_4 --- T0_N1
    T0_N0 --- T0_2
    T0_N1 --- T0_1
    T0_N1 --- T0_3
    T1_N0["$$D^{+}$$"]
    T1_1["$$1: \pi^{+}$$"]
    T1_2["$$2: \pi^{+}$$"]
    T1_3["$$3: \pi^{-}$$"]
    T1_N1@{ shape: text, label: " " }
    T1_4("$$\begin{array}{ll} f_{0}(500) & \omega(782) \\\ f_{0}(980) & \rho(770)^{0} \\\ f_{2}(1270) & \rho(1450)^{0} \\\ f_{0}(1370) & \rho(1700)^{0} \\\ f_{0}(1500) &  \end{array}$$")
    T1_N0 --- T1_4
    T1_4 --- T1_N1
    T1_N0 --- T1_1
    T1_N1 --- T1_2
    T1_N1 --- T1_3

Formulate and display the model#

Define

\[s=\sigma_1=(p_2+p_3)^2,\qquad t=\sigma_2=(p_1+p_3)^2,\qquad u=\sigma_3=(p_1+p_2)^2,\]

The invariant relation, scattering angle, and dynamics definitions below are rendered from the model and package expression objects. All external spins are zero, so the Wigner functions reduce to Legendre polynomials.

The dynamics factor \(X_R\) is the product of the resonance lineshape and the production and decay form factors. \(\hat{B}_J^2\) is AmpForm’s normalized Blatt–Weisskopf function, with \(\hat{B}_J^2(1)=1\); the hat distinguishes it from the unnormalized PDG convention. Each resonance \(R\) has its own decay radius \(R_R\) and the production vertex uses \(R_{D^+}\), all at the package default of \(1\,\mathrm{GeV}^{-1}\). Angular momentum requires \(L_\mathrm{production}=L_\mathrm{decay}=J_R\); canonical-helicity supplies these orbital angular momenta to the dynamics builder.

The dynamics come from AmpForm-DPD’s BreitWignerBuilder with numerator="mass_width", so each propagator is multiplied by \(m_R\Gamma_R\). The decay form factor is evaluated at \(s\) and the production form factor uses the running resonance mass \(\sqrt{s}\) as daughter mass, as the rendered definitions below show. Neither is normalized at the resonance pole. The Breit–Wigner, running width, and form-factor implementations are those of the packages.

The Bose-symmetric amplitude and intensity are

\[\mathcal A(s,t)=\sum_R c_R\left[X_R(s)P_{J_R}(z(s,t))+X_R(t)P_{J_R}(z(t,s))\right],\qquad I(s,t)=|\mathcal A(s,t)|^2,\]

AmpForm’s HelicityAmplitudeBuilder includes identical-particle symmetrization automatically; the fix for this channel was released in v0.15.9 (see the accepted answer in discussion #475). This report uses DalitzPlotDecompositionBuilder from AmpForm-DPD instead.

AmpForm-DPD imposes this symmetry itself. permute_equal_final_states generates the decay chains for both \(\pi^+\) assignments, and formulate then ties the couplings of each resonance’s two subsystems together, using its symmetrization module. The cyclic pair ordering of the DPD paper writes the isobar of subsystem 1 as \((\pi^+_2\pi^-)\) and that of subsystem 2 as \((\pi^-\pi^+_1)\), so the exchange delivers the isobar’s decay products in reversed order. Writing \(z(s,t)=\cos\theta_{23}\), that means \(\cos\theta_{31}(s,t)=-z(t,s)\), and \(P_J(-z)=(-1)^J P_J(z)\) leaves a relative factor \((-1)^{J_R}\) between the two chains. The signs rendered below are exactly that factor. Only the couplings of subsystem 1 survive as free parameters, and the complex-amplitude symmetry check further down verifies the construction.

An overall \(1/\sqrt{2!}\) amplitude convention is absorbed into normalization. The plot covers the full labelled Dalitz domain; integration for an absolute identical-particle rate requires the corresponding \(1/2!\) phase-space factor.

s, t, u = sp.symbols("sigma1:4", nonnegative=True)
m0, m1, m2, m3 = sp.symbols("m:4", nonnegative=True)

dynamics_builder = BreitWignerBuilder(numerator="mass_width")
builder = DalitzPlotDecompositionBuilder(decay)
for name in resonance_names:
    builder.dynamics_choices.register_builder(name, dynamics_builder)
model = builder.formulate(use_coefficients=True)
assert len(model.amplitudes) == 2
Math(aslatex(model.amplitudes, terms_per_line=1))

\[\begin{split}\displaystyle \begin{aligned} A^{1}_{0, 0, 0, 0} \;&=\; \Gamma_{\omega(782)} m_{\omega(782)} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{1}; m_{\omega(782)}, \Gamma_{\omega(782)}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right) \mathcal{H}^\mathrm{\omega(782),1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{\rho(1450)^{0}} m_{\rho(1450)^{0}} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{1}; m_{\rho(1450)^{0}}, \Gamma_{\rho(1450)^{0}}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right) \mathcal{H}^\mathrm{\rho(1450)^{0},1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{\rho(1700)^{0}} m_{\rho(1700)^{0}} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{1}; m_{\rho(1700)^{0}}, \Gamma_{\rho(1700)^{0}}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right) \mathcal{H}^\mathrm{\rho(1700)^{0},1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{\rho(770)^{0}} m_{\rho(770)^{0}} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{1}; m_{\rho(770)^{0}}, \Gamma_{\rho(770)^{0}}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{1}, m_{2}, m_{3}\right) \mathcal{H}^\mathrm{\rho(770)^{0},1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{f_{0}(1370)} m_{f_{0}(1370)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{1}; m_{f_{0}(1370)}, \Gamma_{f_{0}(1370)}\right) \mathcal{H}^\mathrm{f_{0}(1370),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{f_{0}(1500)} m_{f_{0}(1500)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{1}; m_{f_{0}(1500)}, \Gamma_{f_{0}(1500)}\right) \mathcal{H}^\mathrm{f_{0}(1500),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{f_{0}(500)} m_{f_{0}(500)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{1}; m_{f_{0}(500)}, \Gamma_{f_{0}(500)}\right) \mathcal{H}^\mathrm{f_{0}(500),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{f_{0}(980)} m_{f_{0}(980)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{1}; m_{f_{0}(980)}, \Gamma_{f_{0}(980)}\right) \mathcal{H}^\mathrm{f_{0}(980),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{23}\right) \\ \;&+\; \Gamma_{f_{2}(1270)} m_{f_{2}(1270)} \mathcal{R}^\mathrm{BW}_{L=2}\left(\sigma_{1}; m_{f_{2}(1270)}, \Gamma_{f_{2}(1270)}\right) \hat{\mathcal{F}}_{2}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right) \hat{\mathcal{F}}_{2}\left(\sigma_{1}, m_{2}, m_{3}\right) \mathcal{H}^\mathrm{f_{2}(1270),1}_{0, 0, 0, 0} d^{2}_{0,0}\left(\theta_{23}\right) \\ A^{2}_{0, 0, 0, 0} \;&=\; - \Gamma_{\omega(782)} m_{\omega(782)} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{2}; m_{\omega(782)}, \Gamma_{\omega(782)}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{2}, m_{1}, m_{3}\right) \mathcal{H}^\mathrm{\omega(782),1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{31}\right) \\ \;&+\; - \Gamma_{\rho(1450)^{0}} m_{\rho(1450)^{0}} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{2}; m_{\rho(1450)^{0}}, \Gamma_{\rho(1450)^{0}}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{2}, m_{1}, m_{3}\right) \mathcal{H}^\mathrm{\rho(1450)^{0},1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{31}\right) \\ \;&+\; - \Gamma_{\rho(1700)^{0}} m_{\rho(1700)^{0}} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{2}; m_{\rho(1700)^{0}}, \Gamma_{\rho(1700)^{0}}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{2}, m_{1}, m_{3}\right) \mathcal{H}^\mathrm{\rho(1700)^{0},1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{31}\right) \\ \;&+\; - \Gamma_{\rho(770)^{0}} m_{\rho(770)^{0}} \mathcal{R}^\mathrm{BW}_{L=1}\left(\sigma_{2}; m_{\rho(770)^{0}}, \Gamma_{\rho(770)^{0}}\right) \hat{\mathcal{F}}_{1}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right) \hat{\mathcal{F}}_{1}\left(\sigma_{2}, m_{1}, m_{3}\right) \mathcal{H}^\mathrm{\rho(770)^{0},1}_{0, 0, 0, 0} d^{1}_{0,0}\left(\theta_{31}\right) \\ \;&+\; \Gamma_{f_{0}(1370)} m_{f_{0}(1370)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{2}; m_{f_{0}(1370)}, \Gamma_{f_{0}(1370)}\right) \mathcal{H}^\mathrm{f_{0}(1370),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{31}\right) \\ \;&+\; \Gamma_{f_{0}(1500)} m_{f_{0}(1500)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{2}; m_{f_{0}(1500)}, \Gamma_{f_{0}(1500)}\right) \mathcal{H}^\mathrm{f_{0}(1500),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{31}\right) \\ \;&+\; \Gamma_{f_{0}(500)} m_{f_{0}(500)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{2}; m_{f_{0}(500)}, \Gamma_{f_{0}(500)}\right) \mathcal{H}^\mathrm{f_{0}(500),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{31}\right) \\ \;&+\; \Gamma_{f_{0}(980)} m_{f_{0}(980)} \mathcal{R}^\mathrm{BW}_{L=0}\left(\sigma_{2}; m_{f_{0}(980)}, \Gamma_{f_{0}(980)}\right) \mathcal{H}^\mathrm{f_{0}(980),1}_{0, 0, 0, 0} d^{0}_{0,0}\left(\theta_{31}\right) \\ \;&+\; \Gamma_{f_{2}(1270)} m_{f_{2}(1270)} \mathcal{R}^\mathrm{BW}_{L=2}\left(\sigma_{2}; m_{f_{2}(1270)}, \Gamma_{f_{2}(1270)}\right) \hat{\mathcal{F}}_{2}\left(m_{0}^{2}, \sqrt{\sigma_{2}}, m_{2}\right) \hat{\mathcal{F}}_{2}\left(\sigma_{2}, m_{1}, m_{3}\right) \mathcal{H}^\mathrm{f_{2}(1270),1}_{0, 0, 0, 0} d^{2}_{0,0}\left(\theta_{31}\right) \\ \end{aligned}\end{split}\]

\[\begin{split}\displaystyle \begin{aligned} f_{0}(1370) \;&=\; 1 \\ f_{0}(1500) \;&=\; 1 \\ f_{0}(500) \;&=\; 1 \\ f_{0}(980) \;&=\; 1 \\ f_{2}(1270) \;&=\; 1 \\ \omega(782) \;&=\; -1 \\ \rho(1450)^{0} \;&=\; -1 \\ \rho(1700)^{0} \;&=\; -1 \\ \rho(770)^{0} \;&=\; -1 \\ \end{aligned}\end{split}\]

The kinematic definitions and spin factors are

\[\begin{split}\displaystyle \begin{aligned} \sigma_{3} \;&=\; m_{0}^{2} + m_{1}^{2} + m_{2}^{2} + m_{3}^{2} - \sigma_{1} - \sigma_{2} \\ \lambda\left(x, y, z\right) \;&=\; x^{2} - 2 x y - 2 x z + y^{2} - 2 y z + z^{2} \\ \theta_{23} \;&=\; \operatorname{acos}{\left(\frac{2 \sigma_{1} \left(- m_{1}^{2} - m_{2}^{2} + \sigma_{3}\right) - \left(m_{0}^{2} - m_{1}^{2} - \sigma_{1}\right) \left(m_{2}^{2} - m_{3}^{2} + \sigma_{1}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{1}^{2}, \sigma_{1}\right)} \sqrt{\lambda\left(\sigma_{1}, m_{2}^{2}, m_{3}^{2}\right)}} \right)} \\ \theta_{31} \;&=\; \operatorname{acos}{\left(\frac{2 \sigma_{2} \left(- m_{2}^{2} - m_{3}^{2} + \sigma_{1}\right) - \left(m_{0}^{2} - m_{2}^{2} - \sigma_{2}\right) \left(- m_{1}^{2} + m_{3}^{2} + \sigma_{2}\right)}{\sqrt{\lambda\left(m_{0}^{2}, m_{2}^{2}, \sigma_{2}\right)} \sqrt{\lambda\left(\sigma_{2}, m_{3}^{2}, m_{1}^{2}\right)}} \right)} \\ d^{0}_{0,0}\left(\operatorname{acos}{\left(z \right)}\right) \;&=\; 1 \\ d^{1}_{0,0}\left(\operatorname{acos}{\left(z \right)}\right) \;&=\; z \\ d^{2}_{0,0}\left(\operatorname{acos}{\left(z \right)}\right) \;&=\; \frac{3 z^{2}}{2} - \frac{1}{2} \\ \end{aligned}\end{split}\]

The tensor-wave lineshape illustrates the package dynamics. Each definition is obtained by evaluating one layer of its expression object; the production and decay form factors are those of the model, and the barrier factors are shown for all three spins.

\[\begin{split}\displaystyle \begin{aligned} X_{f_{2}(1270)}{\left(\sigma_{1} \right)} \;&=\; \Gamma_{f_{2}(1270)} m_{f_{2}(1270)} \mathcal{R}^\mathrm{BW}_{L=2}\left(\sigma_{1}; m_{f_{2}(1270)}, \Gamma_{f_{2}(1270)}\right) \hat{\mathcal{F}}_{2}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right) \hat{\mathcal{F}}_{2}\left(\sigma_{1}, m_{2}, m_{3}\right) \\ \mathcal{R}^\mathrm{BW}_{L=2}\left(\sigma_{1}; m_{f_{2}(1270)}, \Gamma_{f_{2}(1270)}\right) \;&=\; \frac{1}{m_{f_{2}(1270)}^{2} - i m_{f_{2}(1270)} \Gamma_{21270}\left(\sigma_{1}\right) - \sigma_{1}} \\ \Gamma_{21270}\left(\sigma_{1}\right) \;&=\; \frac{\Gamma_{f_{2}(1270)} \hat{\mathcal{F}}_{2}\left(\sigma_{1}, m_{2}, m_{3}\right)^{2} \rho_{1}\left(\sigma_{1}\right)}{\hat{\mathcal{F}}_{2}\left(m_{f_{2}(1270)}^{2}, m_{2}, m_{3}\right)^{2} \rho_{21270}\left(m_{f_{2}(1270)}^{2}\right)} \\ \hat{\mathcal{F}}_{2}\left(m_{0}^{2}, \sqrt{\sigma_{1}}, m_{1}\right) \;&=\; \sqrt{\hat{B}_{2}^2\left(\left(R_{D^{+}}\right)^{2} q^2_{0}\left(m_{0}^{2}\right)\right)} \\ \hat{\mathcal{F}}_{2}\left(\sigma_{1}, m_{2}, m_{3}\right) \;&=\; \sqrt{\hat{B}_{2}^2\left(R_{f_{2}(1270)}^{2} q^2_{1}\left(\sigma_{1}\right)\right)} \\ q^2_{1}\left(\sigma_{1}\right) \;&=\; \frac{\left(\sigma_{1} - \left(m_{2} - m_{3}\right)^{2}\right) \left(\sigma_{1} - \left(m_{2} + m_{3}\right)^{2}\right)}{4 \sigma_{1}} \\ \hat{B}_{0}^2\left(z\right) \;&=\; 1 \\ \hat{B}_{1}^2\left(z\right) \;&=\; \frac{2 z}{z + 1} \\ \hat{B}_{2}^2\left(z\right) \;&=\; \frac{13 z^{2}}{z^{2} + 3 z + 9} \\ \end{aligned}\end{split}\]

# Retain the package coefficient names; set their numerical values below.
couplings = {
    symbol: value
    for symbol, value in model.parameter_defaults.items()
    if isinstance(symbol, sp.Indexed)
}
fixed_parameters = {
    symbol: value
    for symbol, value in model.parameter_defaults.items()
    if symbol not in couplings
}


def substitute_kinematics(expression):
    expression = expression.doit().xreplace(model.variables).doit()
    expression = expression.subs(u, model.invariants[u])
    return expression.xreplace(fixed_parameters)


amplitude = substitute_kinematics(sum(model.amplitudes.values()))
intensity = substitute_kinematics(model.full_expression)
intensity_function = create_parametrized_function(intensity, couplings, backend="numpy")
amplitude_function = create_parametrized_function(amplitude, couplings, backend="numpy")

M = decay.states[0].mass
mpi = decay.states[1].mass
sigma_sum = M**2 + 3 * mpi**2

coefficient_symbols = {
    name: next(symbol for symbol in couplings if resonance.latex in str(symbol))
    for name, resonance in resonances.items()
}
assert len(set(coefficient_symbols.values())) == len(resonance_names)
component_functions = {
    name: create_function(amplitude.diff(symbol), backend="numpy")
    for name, symbol in coefficient_symbols.items()
}

Physical Dalitz domain and numerical checks#

For a scalar parent, three-body phase space is constant in \(ds\,dt\) inside the kinematic boundary. Thus \(I(s,t)\) is proportional to the Dalitz density without an additional momentum weight. The bounds below are rendered from the same symbolic expressions used for numerical evaluation.

We evaluate only interior points, where scattering angles are defined. The checks below verify Bose symmetry of the complex amplitude and finite nonnegative intensity. An intensity-only symmetry check would miss an antisymmetric amplitude, which is what the exchange signs prevent. The model’s own intensity expression is also compared with the squared modulus of the summed subsystem amplitudes.

Bose symmetry and intensity checks passed at 2000 physical points.

\[\begin{split}\displaystyle \begin{aligned} s_{min} \;&=\; 4 m_{1}^{2} \\ s_{max} \;&=\; \left(m_{0} - m_{1}\right)^{2} \\ t_{-}{\left(\sigma_{1} \right)} \;&=\; \frac{m_{0}^{2}}{2} + \frac{3 m_{1}^{2}}{2} - \frac{\sigma_{1}}{2} - \frac{\sqrt{\lambda\left(m_{0}^{2}, \sigma_{1}, m_{1}^{2}\right) \lambda\left(\sigma_{1}, m_{1}^{2}, m_{1}^{2}\right)}}{2 \sigma_{1}} \\ t_{+}{\left(\sigma_{1} \right)} \;&=\; \frac{m_{0}^{2}}{2} + \frac{3 m_{1}^{2}}{2} - \frac{\sigma_{1}}{2} + \frac{\sqrt{\lambda\left(m_{0}^{2}, \sigma_{1}, m_{1}^{2}\right) \lambda\left(\sigma_{1}, m_{1}^{2}, m_{1}^{2}\right)}}{2 \sigma_{1}} \\ \end{aligned}\end{split}\]

Couplings informed by measured component strengths#

A coefficient of one has a different meaning for each unnormalized lineshape. Copying published magnitudes would also mix different angular and lineshape normalizations. Instead, let \(a_R(s,t)\) be the Bose-symmetrized amplitude with coefficient one, and define

\[ N_R=\int_{\mathcal D}|a_R(s,t)|^2\,ds\,dt,\qquad r_R=\frac{\mathrm{FF}_R^{\mathrm{source}}}{\mathrm{FF}_{\rho(770)}^{\mathrm{source}}}. \]

With \(c_{\rho(770)}=1\), choose

\[c_R=\sqrt{r_R\frac{N_{\rho(770)}}{N_R}}\,e^{i\phi_R}.\]

This preserves ratios of diagonal component integrals, independently of the overall amplitude normalization. CLEO [CLEO Collaboration et al., 2007] supplies the six original components; LHCb [The LHCb collaboration et al., 2023] supplies the three additional vectors, using its own \(\rho(770)\) fraction as their reference. These are separate fits, not a combined experimental result. Their central phase values are adopted only as illustrative starting values in the DPD convention: an exact phase transfer would require matching each source’s angular signs and dynamical phase conventions. No uncertainty propagation or refit is performed.

The integrals use equal-area cells on the physical Dalitz domain. Comparing 500 and 1000 bins per axis checks numerical stability, including the narrow \(\omega\) band. The displayed model fractions are computed with the coherent denominator,

\[ \mathrm{FF}^{\mathrm{model}}_R= \frac{|c_R|^2N_R}{\int_{\mathcal D}|\sum_k c_k a_k|^2\,ds\,dt}. \]

They need not equal the source fractions or sum to one because of interference and the combined resonance content.

Normalize the resonance contributions and set literature-informed couplings

Hide code cell source

def integration_grid(n_bins):
    bin_edges = np.linspace(lower, upper, n_bins + 1)
    bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2
    grid_s, grid_t = np.meshgrid(bin_centers, bin_centers)
    t_low, t_high = dalitz_limits(grid_s)
    mask = (grid_t > t_low) & (grid_t < t_high)
    data = {"sigma1": grid_s[mask], "sigma2": grid_t[mask]}
    cell_area = ((upper - lower) / n_bins) ** 2
    return bin_edges, grid_s, grid_t, mask, data, cell_area


def component_integrals(data, cell_area):
    return {
        name: float(np.sum(np.abs(function(data)) ** 2) * cell_area)
        for name, function in component_functions.items()
    }


*_, coarse_data, coarse_area = integration_grid(500)
coarse_norms = component_integrals(coarse_data, coarse_area)
edges, s_grid, t_grid, physical, grid_data, cell_area = integration_grid(1000)
norms = component_integrals(grid_data, cell_area)
assert all(np.isfinite(value) and value > 0 for value in norms.values())
normalization_change = max(abs(coarse_norms[name] / norms[name] - 1) for name in norms)
assert normalization_change < 0.01

for name, (fraction, phase, source) in literature.items():
    ratio = fraction / reference_fractions[source]
    magnitude = np.sqrt(ratio * norms["rho(770)0"] / norms[name])
    couplings[coefficient_symbols[name]] = magnitude * np.exp(1j * np.deg2rad(phase))
parameters = {str(symbol): value for symbol, value in couplings.items()}
intensity_function.update_parameters(parameters)
amplitude_function.update_parameters(parameters)

# Check the numerical intensity against the sum of the individual complex waves.
component_sum = sum(
    couplings[coefficient_symbols[name]] * function(test_data)
    for name, function in component_functions.items()
)
np.testing.assert_allclose(
    amplitude_function(test_data), component_sum, rtol=1e-11, atol=1e-11
)
np.testing.assert_allclose(
    component_sum,
    amplitude_function({"sigma1": t_test, "sigma2": s_test}),
    rtol=1e-11,
    atol=1e-11,
)
coherent_integral = float(np.sum(intensity_function(grid_data)) * cell_area)
assert np.isfinite(coherent_integral)
assert coherent_integral > 0
model_fractions = {
    name: abs(couplings[coefficient_symbols[name]]) ** 2 * norm / coherent_integral
    for name, norm in norms.items()
}
for name, (fraction, _, source) in literature.items():
    np.testing.assert_allclose(
        model_fractions[name] / model_fractions["rho(770)0"],
        fraction / reference_fractions[source],
        rtol=1e-12,
    )

rows = [
    r"| Resonance | Source | Source FF [%] | Starting phase [deg] | $\lvert c_R\rvert$ | Model FF [%] |",
    "|---|---|---:|---:|---:|---:|",
]
for name, (fraction, phase, source) in literature.items():
    coefficient = couplings[coefficient_symbols[name]]
    rows.append(
        f"| ${resonances[name].latex}$ | {source} | {fraction:g} | {phase:g} | "
        f"{abs(coefficient):.4g} | {100 * model_fractions[name]:.3f} |"
    )
rows.extend([
    f"\nLargest change in component integrals on grid refinement: {100 * normalization_change:.3f}%.",
    f"\nSum of model fit fractions: {100 * sum(model_fractions.values()):.2f}%.",
])
Markdown("\n".join(rows))

Resonance	Source	Source FF [%]	Starting phase [deg]	\(\lvert c_R\rvert\)	Model FF [%]
\(\rho(770)^{0}\)	CLEO	20	0	1	20.817
\(f_{0}(500)\)	CLEO	41.8	-3	0.2335	43.508
\(f_{0}(980)\)	CLEO	4.1	12	0.156	4.268
\(f_{2}(1270)\)	CLEO	18.2	-123	3.46	18.944
\(f_{0}(1370)\)	CLEO	2.6	-21	0.0515	2.706
\(f_{0}(1500)\)	CLEO	3.4	-44	0.1218	3.539
\(\omega(782)\)	LHCb	0.103	-103.3	0.2534	0.082
\(\rho(1450)^{0}\)	LHCb	5.4	47	0.2673	4.324
\(\rho(1700)^{0}\)	LHCb	5.7	-65.7	0.5732	4.564

Largest change in component integrals on grid refinement: 0.311%.

Sum of model fit fractions: 102.75%.

Dalitz plot#

The 1000 × 1000 regular grid avoids Monte Carlo fluctuations and resolves the narrow \(\omega(782)\) contribution. The color is the coherent intensity divided by its maximum on this grid, on a logarithmic scale; it is not a normalized probability density or a measured event count. The two axes distinguish the otherwise identical positive pions. Reflection about \(s=t\) must leave the density unchanged.

../_images/6f1953ee56e22618ed347f131475e3f71f92d59dec551f33bb9da084fd8179cf.svg

The \(\rho(770)\), \(f_0(980)\), and \(f_2(1270)\) bands now have strengths informed by measured component ratios. The heavier states contribute broad structures and interference; their presence does not imply separate visible bands. The small \(\omega(782)\) contribution is concentrated near the \(\rho(770)\) band. The intensity remains symmetric under exchange of the positive pions.

One-dimensional mass projections#

The projections integrate the coherent AmpForm-DPD intensity over the physical Dalitz grid. Each has unit area in its pair-mass-squared variable. Bose symmetry makes the \(s\) and \(t\) projections identical; \(u\) describes the same-charge pion pair. Projection bins align with the integration grid to avoid binning artifacts.

../_images/0fd6c2c4690a788ca0659bc9c73d572fc1643d05267e9f318bba29aa325a47a2.svg