Dalitz plot analysis of the $D^+\to K^-K^+K^+$ decay

TL;DR

This paper presents the first detailed Dalitz plot analysis of the decay D+ to K-K+K+ using LHCb data, revealing the dominant S-wave component and providing insights into K-K+ scattering amplitudes.

Contribution

It introduces the first amplitude analysis of the D+ to K-K+K+ decay, employing isobar and chiral Lagrangian models to study the resonant structure and scattering amplitudes.

Findings

S-wave component is dominant in the decay

Small contribution from phi(1020) meson

Negligible tensor resonance contributions

Abstract

The resonant structure of the doubly Cabibbo-suppressed decay $D^+ \to K^-K^+K^+$ is studied for the first time. The measurement is based on a sample of pp-collision data, collected at a centre-of-mass energy of 8 TeV with the LHCb detector and corresponding to an integrated luminosity of 2 fb$^-1$. The amplitude analysis of this decay is performed with the isobar model and a phenomenological model based on an effective chiral Lagrangian. In both models the S-wave component in the $K^-K^+$ system is dominant, with a small contribution of the $\phi(1020)$ meson and a negligible contribution from tensor resonances. The $K^-K^+$ scattering amplitudes for the considered combinations of spin (0,1) and isospin (0,1) of the two-body system are obtained from the Dalitz plot fit with the phenomenological decay amplitude.

Figures

Tables

Table 1. Table 1: Results from the D + → K − ​ K + ​ K + → superscript 𝐷 superscript 𝐾 superscript 𝐾 superscript 𝐾 {{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}} Dalitz plot fit with the isobar models A, B and C. Magnitudes, | c k | subscript 𝑐 𝑘 |c_{k}| , phases, arg ⁡ ( c k ) subscript 𝑐 𝑘 \arg(c_{k}) (in degrees), and fit fractions (in %) are given with statistical uncertainties only.

		Model A		Model B		Model C
$\varphi \left(1020\right)K^{+}$	Magnitude	1 [fixed]		1 [fixed]		1 [fixed]
	Phase	0 [fixed]		0 [fixed]		0 [fixed]
	Fraction	6.17 $\pm$	0.47	6.40 $\pm$	0.47	6.40 $\pm$	0.48
$f_0\left(980\right)K^{+}$	Magnitude	3.12 $\pm$	0.10	2.64 $\pm$	0.08	2.84 $\pm$	0.13
	Phase	$-58.9$ $\pm$	4.9	$-36.5$ $\pm$	7.6	$-25.9$ $\pm$	8.4
	Fraction	23.7 $\pm$	3.0	17.7 $\pm$	2.1	20.4 $\pm$	1.5
$f_0\left(1370\right)K^{+}$	Magnitude	3.46 $\pm$	0.46	2.33 $\pm$	0.35	–
	Phase	13.1 $\pm$	7.7	42 $\pm$	10	–
	Fraction	25.4 $\pm$	5.0	18.7 $\pm$	1.5	–
	$f_0\left(1370\right)$ mass [$\mathrm{GeV}$ ]	1.422 $\pm$	0.015	1.401 $\pm$	0.009	–
	$f_0\left(1370\right)$ width [$\mathrm{GeV}$ ]	0.324 $\pm$	0.038	0.178 $\pm$	0.031	–
NR	Magnitude	–		8.8 $\pm$	1.3	11.7 $\pm$	1.8
	Phase	–		$-5.5$ $\pm$	6.5	$-39.0$ $\pm$	4.4
	Fraction	–		18.4 $\pm$	5.9	32.7 $\pm$	8.2
$a_0\left(980\right)K^{+}$	Magnitude	–		–		5.9 $\pm$	0.4
	Phase	–		–		48.5 $\pm$	3.0
	Fraction	–		–		53.5 $\pm$	7.4
	S-wave fraction	92 $\pm$	11	91 $\pm$	13	93 $\pm$	12
	Fractions sum	55.4 $\pm$	5.9	61.2 $\pm$	6.4	113 $\pm$	11

Table 2. Table 2: Fit results with model A, given in terms of the magnitudes | c k | subscript 𝑐 𝑘 |c_{k}| , phases, arg ⁡ ( c k ) subscript 𝑐 𝑘 \arg(c_{k}) (in degrees), and fit fractions (in %). For each measurement, the first uncertainty is statistical, the second systematic and the third is a systematic uncertainty due to model.

Component	Magnitude				Phase [deg.]				Fraction (%)
$\varphi \left(1020\right)K^{+}$	1.0 [fixed]				0.0 [fixed]				6.17 $\pm$	0.47  $\pm$	0.19  $\pm$	0.41
$f_0\left(980\right)K^{+}$	3.12 $\pm$	0.10  $\pm$	0.13  $\pm$	0.33	$-58.9$ $\pm$	4.9  $\pm$	2.3  $\pm$	2.0	23.7 $\pm$	3.0  $\pm$	2.1  $\pm$	3.3
$f_0\left(1370\right)K^{+}$	3.46 $\pm$	0.46  $\pm$	0.32  $\pm$	0.73	13.1 $\pm$	7.7  $\pm$	1.6  $\pm$	3.2	25.4 $\pm$	5.0  $\pm$	3.4  $\pm$	3.8
sum									55.4 $\pm$	5.9  $\pm$	0.4  $\pm$	0.6

Table 3. Table 3: Results of the D + → K − ​ K + ​ K + → superscript 𝐷 superscript 𝐾 superscript 𝐾 superscript 𝐾 D^{+}\rightarrow K^{-}K^{+}K^{+} Dalitz plot fit with the Triple-M amplitude.

parameter	value
$F$	${94.3}_{-1.7}^{+2.8}$ $\pm$	1.5$\mathrm{MeV}$
$m_{a_0}$	${947.7}_{-5.0}^{+5.5}$ $\pm$	6.6$\mathrm{MeV}$
$m_{S_o}$	${992.0}_{-7.5}^{+8.5}$ $\pm$	8.6$\mathrm{MeV}$
$m_{S_1}$	${1330.2}_{-6.5}^{+5.9}$ $\pm$	5.1$\mathrm{MeV}$
$m_{\varphi }$	${1019.54}_{-0.10}^{+0.10}$ $\pm$	0.51$\mathrm{MeV}$
$G_{\varphi }$	${0.464}_{-0.009}^{+0.013}$ $\pm$	0.007 0
$c_d$	$-{78.9}_{-2.7}^{+4.2}$ $\pm$	1.9$\mathrm{MeV}$
$c_m$	${106.0}_{-4.6}^{+7.7}$ $\pm$	3.3$\mathrm{MeV}$
${\tilde{c}}_d$	$-{6.15}_{-0.54}^{+0.55}$ $\pm$	0.19$\mathrm{MeV}$
${\tilde{c}}_m$	$-{10.8}_{-1.5}^{+2.0}$ $\pm$	0.4$\mathrm{MeV}$

Table 4. Table 4: Relative fractions (%) of the various components of the Triple-M amplitude. The uncertainties correspond to the combined statistical and systematic uncertainties.

${\mathrm{FF}}_{\mathrm{NR}}$	${\mathrm{FF}}^{00}$	${\mathrm{FF}}^{01}$	${\mathrm{FF}}^{10}$	${\mathrm{FF}}^{11}$	${\mathrm{FF}}_{S-\mathrm{wave}}$
14 $\pm$ 1	29 $\pm$ 1	131 $\pm$ 2	7.1 $\pm$ 0.9	0.26 $\pm$ 0.01	94 $\pm$ 1

Table 5. Table 5: Systematic uncertainties (%) on the results of the isobar model fit. For comparison, the statistical uncertainties are listed in the last column.

parameter	binning	sim. stat.	bkg.	total	model	stat.
$\left|c_{f_0\left(980\right)}\right|$	1.0	1.4	3.8	4.2	11	3.0
${\delta }_{f_0\left(980\right)}$	3.1	1.9	1.4	3.9	3.4	8.2
$\left|c_{f_0\left(1370\right)}\right|$	3.5	3.5	7.8	9.2	21	13
${\delta }_{f_0\left(1370\right)}$	9.3	5.2	4.4	12	24	59
$M_{f_0\left(1370\right)}$	0.1	0.3	0.6	0.7	2.0	1.0
${\Gamma }_{f_0\left(1370\right)}$	3.7	3.0	3.1	5.7	12	12

Table 6. Table 6: Systematic uncertainties (%) on the results of the Triple-M fit. For comparison, the statistcal uncertainties are listed in the last column.

parameter	binning	sim. stat.	PID	bkg.	fit bias	total	stat.
$F$	0.53	0.07	0.09	1.5	0.11	1.6	1.8
$m_{a_0}$	0.54	0.14	-	0.40	0.16	0.70	0.54
$m_{So}$	0.60	0.21	-	0.56	0.21	0.87	0.82
$m_{S1}$	0.16	0.15	-	0.13	0.04	0.26	0.41
$m_{\varphi }$	0.002	0.001	-	-	0.002	0.003	0.005
$G_{\varphi }$	0.86	0.25	0.02	1.2	0.15	1.5	1.9
$c_d$	0.18	0.08	0.09	2.4	0.13	2.4	3.3
$c_m$	0.16	0.11	-	2.7	0.10	2.7	4.7
${\tilde{c}}_d$	0.13	0.15	-	2.6	1.1	3.1	8.8
${\tilde{c}}_m$	0.19	0.11	0.08	2.8	1.9	3.4	13

Table 7. Table 7: Blatt–Weisskopf form factors for angular momentum L = 0 , 1 , 2 𝐿 0 1 2 L=0,1,2 with two distinct formulations.

L	$B_L$	$B_L^{\prime }$
0	1	1
1	$\sqrt{\frac{2z}{1+z}}$	$\sqrt{\frac{1+z_0}{1+z}}$
2	$\sqrt{\frac{13z^2}{1+z}}$	$\sqrt{\frac{{\left(z_0-3\right)}^2+9z_0}{{\left(z-3\right)}^2+9z}}$

Equations

$$\mathcal{L}=\prod ^{N_{\mathrm{cand}}}{f}_{\mathrm{S}}\times S_{\mathrm{PDF}}(s_{12},s_{13})+(1-f_{\mathrm{S}})\times B_{\mathrm{PDF}}(s_{12},s_{13}),$$

$$S_{\mathrm{PDF}}(s_{12},s_{13})=\frac{1}{N_{\mathrm{S}}}{\left|\mathcal{T}(s_{12},s_{13})\right|}^2\epsilon (s_{12},s_{13}),$$

$$N_{\mathrm{S}}=\int \mathrm{d}{s}_{12}\mathrm{d}{s}_{13}{\left|\mathcal{T}(s_{12},s_{13})\right|}^2\epsilon (s_{13},s_{13}).$$

$${\Delta }_i=\frac{(N_{\mathrm{pred}}^i-N_{\mathrm{obs}}^i)}{{\sigma }_i},$$

$$\mathcal{T}(s_{12},s_{13})=c_{\mathrm{NR}}+\sum _k{c}_k{T}_k(s_{12},s_{13}).$$

$$T_k(s_{12},s_{13})=F_D^L(s_{12})F_R^L(s_{12})\times \mathcal{S}(s_{12},s_{13})\times M_R(s_{12})+(2\leftrightarrow 3).$$

$${\mathrm{FF}}_k=\frac{\int \mathrm{d}{s}_{12}\mathrm{d}{s}_{13}|c_k{T}_k(s_{12},s_{13}){|}^2}{\int \mathrm{d}{s}_{12}\mathrm{d}{s}_{13}{\left|{\sum }_i{c}_i{T}_i(s_{12},s_{13})\right|}^2}.$$

$$⟨K^{-}{K}^{+}{K}^{+}|\mathcal{T}|D^{+}⟩=-\left[\frac{G_F}{\sqrt{2}}{\sin }^2{\theta }_C\right]⟨K^{-}{K}^{+}{K}^{+}|A_{\mu }|0⟩⟨0|A_{\mu }|D^{+}⟩,$$

$$\mathcal{T}=T_{\mathrm{NR}}+\sum _{J,I}{T}^{JI},J,I=0,1.$$

$$T_{\mathrm{NR}}$$

$$\mathcal{T}=\left[T^S+T^P+(2\leftrightarrow 3)\right],$$

$$T^S=T_{NR}^S+T^{00}+T^{01}$$

$$T^P=T_{NR}^P+T^{11}+T^{10}.$$

$${\mathrm{FF}}_{\mathrm{NR}}=\frac{\int \mathrm{d}{s}_{12}\mathrm{d}{s}_{13}|T_{\mathrm{NR}}(s_{12},s_{13}){|}^2}{\int \mathrm{d}{s}_{12}\mathrm{d}{s}_{13}{\left|\mathcal{T}(s_{12},s_{13})\right|}^2},{\mathrm{FF}}^{JI}=\frac{\int \mathrm{d}{s}_{12}\mathrm{d}{s}_{13}|T^{JI}(s_{12},s_{13}){|}^2}{\int \mathrm{d}{s}_{12}\mathrm{d}{s}_{13}{\left|\mathcal{T}(s_{12},s_{13})\right|}^2}.$$

$$\mathcal{S}=(-2|p_1||p_3|{)}^L{P}_L(\cos {\theta }_{13}^{R_{12}}),$$

$$M_R(s_{12})=\frac{1}{s_{12}-m_R^2+i{m}_R\Gamma (s_{12})},$$

$$\Gamma (s_{12})={\Gamma }_R{\left(\frac{q}{q_0}\right)}^{2L+1}\frac{m_R}{\sqrt{s_{12}}}{\left(\frac{F_R^L(z)}{F_R^L(z_0)}\right)}^2$$

$$M_R(s_{12})=\frac{1}{s_{12}-m_R^2+i{m}_R({\rho }_{\pi \pi }{g}_{\pi }^2+{\rho }_{KK}{g}_K^2)},$$

$$\begin{array}{rl}{\rho }_{\pi \pi }& =\sqrt{\left(\frac{s_{12}}{4}-m_{\pi }^2\right)}+\sqrt{\left(\frac{s_{12}}{4}-m_{{\pi }^0}^2\right)},\\ {\rho }_{KK}& =\sqrt{\left(\frac{s_{12}}{4}-m_K^2\right)}+\sqrt{\left(\frac{s_{12}}{4}-m_{K^0}^2\right)}.\end{array}$$

$$\mathcal{T}$$

$$C=\left\{\left[\frac{G_F}{\sqrt{2}}{\sin }^2{\theta }_C\right]\frac{2F_D}{F}\frac{m_K^2}{(m_D^2-m_K^2)}\right\},$$

$$T_{\mathrm{NR}}=C\left[(s_{12}-m_K^2)+(s_{13}-m_K^2)\right].$$

$$T^{(1,1)}=-\frac{1}{4}\left[{\bar{\Gamma }}_{KK}^{(1,1)}-{\Gamma }_{c|KK}^{(1,1)}\right](s_{13}-s_{23}),$$

$${\bar{\Gamma }}_{KK}^{(1,1)}=\frac{1}{D_{\rho }(s_{12})}\left[M_{21}^{(1,1)}{\Gamma }_{(0)\pi \pi }^{(1,1)}+\left(1-M_{11}^{(1,1)}\right){\Gamma }_{(0)KK}^{(1,1)}\right],$$

$$T^{(1,0)}=-\frac{1}{4}\left[{\bar{\Gamma }}_{KK}^{(1,0)}-{\Gamma }_{c|KK}^{(1,0)}\right](s_{13}-s_{23}),$$

$${\bar{\Gamma }}_{KK}^{(1,0)}=\frac{1}{D_{\varphi }(s_{12})}{\Gamma }_{(0)KK}^{(1,0)},$$

$$T^{(0,1)}=-\frac{1}{2}\left[{\bar{\Gamma }}_{KK}^{(0,1)}-{\Gamma }_{c|KK}^{(0,1)}\right],$$

$${\bar{\Gamma }}_{KK}^{(0,1)}=\frac{1}{D_{a_0}(s_{12})}\left[M_{21}^{(0,1)}{\Gamma }_{(0)\pi 8}^{(0,1)}+\left(1-M_{11}^{(0,1)}\right){\Gamma }_{(0)KK}^{(0,1)}\right],$$

$$T^{(0,0)}=-\frac{1}{2}\left[{\bar{\Gamma }}_{KK}^{(0,0)}-{\Gamma }_{c|KK}^{(0,0)}\right],$$

$${\bar{\Gamma }}_{KK}^{(0,0)}=\frac{1}{D_S(s_{12})}\{\left[M_{21}^{(0,0)}\left(1-M_{33}^{(0,0)}\right)+M_{23}^{(0,0)}{M}_{31}^{(0,0)}\right]{\Gamma }_{(0)\pi \pi }^{(0,0)}$$

$$+\left[\left(1-M_{11}^{(0,0)}\right)\left(1-M_{33}^{(0,0)}\right)-M_{13}^{(0,0)}{M}_{31}^{(0,0)}\right]{\Gamma }_{(0)KK}^{(0,0)}$$

$$+\left[M_{23}^{(0,0)}\left(1-M_{11}^{(0,0)}\right)+M_{13}^{(0,0)}{M}_{21}^{(0,0)}\right]{\Gamma }_{(0)88}^{(0,0)}\}.$$

$${\Gamma }_{(0)\pi \pi }^{(1,1)}=C\left\{\left[\frac{\sqrt{2}{G}_V^2}{F^2}\right]\frac{s_{12}^2}{s_{12}^2-m_{\rho }^2}+{\left[-\frac{1}{\sqrt{2}}\right]}_c\right\},$$

$${\Gamma }_{(0)KK}^{(1,1)}=C\left\{\left[\frac{G_V^2}{F^2}\right]\frac{s_{12}^2}{s_{12}^2-m_{\rho }^2}+{\left[-\frac{1}{2}\right]}_c\right\}.$$

Full text

EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

​​​

CERN-EP-2018-336

LHCb-PAPER-2018-039

Dalitz plot analysis of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay

LHCb collaboration†††Authors are listed at the end of this paper.

The resonant structure of the doubly Cabibbo-suppressed decay ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ is studied for the first time. The measurement is based on a sample of $pp$-collision data, collected at a centre-of-mass energy of 8 TeV with the LHCb detector and corresponding to an integrated luminosity of 2$\mbox{\,fb}^{-1}$. The amplitude analysis of this decay is performed with the isobar model and a phenomenological model based on an effective chiral Lagrangian. In both models the S-wave component in the $K^{-}K^{+}$ system is dominant, with a small contribution of the $\phi(1020)$ meson and a negligible contribution from tensor resonances. The $K^{+}K^{-}$ scattering amplitudes for the considered combinations of spin (0,1) and isospin (0,1) of the two-body system are obtained from the Dalitz plot fit with the phenomenological decay amplitude.

Published in JHEP 04 (2019) 063.

© 2024 CERN for the benefit of the LHCb collaboration. CC-BY-4.0 licence.

1 Introduction

The theoretical treatment of weak decays of charm mesons is very challenging. The charm quark is not light enough for the reliable application of chiral perturbation theory, which is successfully applied in predictions of kaon decays. The charm quark is also not heavy enough for the reliable application of the factorisation approach and heavy-quark expansion tools, as used in predictions of properties of b hadrons. The description of charm meson decays relies on approximate symmetries and phenomenological models. For such models, the knowledge of branching fractions and the resonant structures, in the case of multi-body decays, are key inputs. In this paper, the first determination of the resonant structure of the doubly Cabibbo-suppressed decay ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ is presented.111Charge conjugation is implied throughout the paper. The analysis is based on a data sample of $pp$ collisions collected with the LHCb detector, corresponding to an integrated luminosity of 2 $\mbox{\,fb}^{-1}$ at a centre-of-mass energy of 8 TeV. The determination of the resonant structure of this decay is complementary to the recent LHCb measurement of its branching fraction [1], based on the same data set.

The amplitude analysis of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay is performed using two methods. The Dalitz plot is fitted with the isobar model, in which the decay amplitude is a coherent sum of resonant and nonresonant amplitudes [2]. The Dalitz plot is also fitted with a phenomenological model derived from an effective chiral Lagrangian with resonances [3]. This phenomenological model, referred to as the multi-meson model, or Triple-M, includes the effects of coupled channels — $\pi\pi$, $K^{+}K^{-}$, $\pi\eta$, $\eta\eta$ and $\rho\pi$ — in the final state interactions (FSI), in four considered combinations of spin J and isospin I ($J\!=\!0,1$; $I\!=\!0,1$). Given the small phase space of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay and the lack of tensor resonances with significant coupling to $K^{+}K^{-}$, the contribution from D-wave is expected to be suppressed.

An additional motivation for the Dalitz plot analysis of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay is to obtain the $K^{+}K^{-}$ scattering amplitudes. Most information currently available on $\pi\pi$ and $K\pi$ scattering is obtained indirectly from meson-nucleon interactions [4, 5, 6]. In the regime where the momentum transferred to the nucleon is small enough, the interaction is assumed to be dominated by the one-pion-exchange amplitude. The asymptotically free incoming meson interacts with a virtual pion, resulting in what is generally referred as $\pi\pi$ and $K\pi$ scattering data. The resulting $\pi\pi\rightarrow\pi\pi$ and $K\pi\rightarrow K\pi$ phase shifts are affected by ambiguities and large systematic uncertainties. The $\pi\pi\rightarrow K{\kern 1.99997pt\overline{\kern-1.99997ptK}{}}$ scattering was studied both in $\pi p$ and $\pi n$ reactions [7, 8], and in $p\bar{p}$ annihilation at rest [9]. For the $K{\kern 1.99997pt\overline{\kern-1.99997ptK}{}}\rightarrow K{\kern 1.99997pt\overline{\kern-1.99997ptK}{}}$ scattering, no meson-nucleon data exists.

Three-body decays of $D$ mesons into kaons and pions are an interesting alternative for light-meson spectroscopy, as they are complementary to the meson-nucleon reactions. Large data sets from the $B$-factories and LHCb exist for these decays. However, it is necessary to isolate the physics of two-body systems from the rich dynamics of three-body decays, which involve the weak decay of the $c$ quark, the formation of the mesons and their FSI. This is achieved with the Triple-M decay amplitude, in which these three stages are included. The FSI are described in terms of the $K^{+}K^{-}$ scattering amplitudes for the considered spin-isospin combinations, allowing the determination of these amplitudes from a fit to the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ Dalitz plot.

This paper is organised as follows. A brief description of the LHCb detector is presented in Sec. 2. The signal selection is presented in Sec. 3. In Sec. 4, the efficiency determination and background model are discussed. The formalism for the Dalitz plot fit is presented in Sec. 5. In Sec. 6, the results of the fit with the isobar model are presented, whilst the results of the Dalitz plot fit with the Triple-M amplitude are presented in Sec. 7. Systematic uncertainties are discussed in Sec. 8. A summary and the conclusions are presented in Sec. 9.

2 Detector and simulation

The LHCb detector [10, 11] is a single-arm forward spectrometer covering the pseudorapidity range $2<\eta<5$, designed for the study of particles containing $b$ or $c$ quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the $pp$ interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about $4{\mathrm{\,Tm}}$, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum, $p$, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200$\mathrm{\,Ge\kern-1.00006ptV}$.222Natural units with $\hbar=c=1$ are used in this paper. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), is measured with a resolution of (15+29/\mbox{p_{\rm T}}){\,\upmu\rm m}, where $p_{\rm T}$ is the component of the momentum transverse to the beam, in GeV. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [12]. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and pre-shower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multi-wire proportional chambers.

The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. At the hardware trigger stage, events are required to have a muon with high $p_{\rm T}$ or a hadron, photon or electron with high transverse energy in the calorimeters. The software trigger is divided into two parts. The first employs a partial reconstruction of the candidates from the hardware trigger and a cut-based selection. In the second stage, a full event reconstruction is applied and various dedicated algorithms are used in the selection of specific decays. In this analysis, a dedicated algorithm is used to select ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay candidates.

In the simulation, $pp$ collisions are generated using Pythia [13, *Sjostrand:2007gs] with a specific LHCb configuration [15]. Decays of hadronic particles are described by EvtGen [16], in which final-state radiation is generated using Photos [17]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [18, *Agostinelli:2002hh] as described in Ref. [20].

3 Candidate selection

The ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay candidates are selected offline with requirements that exploit the decay topology by combining three charged particles identified as kaons according to particle-identification (PID) criteria. These particles must form a good-quality decay vertex, detached from the PV. The PV is chosen as that with the smallest value of $\chi^{2}_{\rm IP}$, where $\chi^{2}_{\rm IP}$ is defined as the difference in the vertex-fit $\chi^{2}$ of the PV reconstructed with and without the particle under consideration, in this case the ${D}^{+}$ candidate. The selection of candidates is based on the distance between the PV and the ${D}^{+}$ decay vertex (the flight distance); the IP of the ${D}^{+}$ candidate; the angle between the reconstructed ${D}^{+}$ momentum vector and the vector connecting the PV to the decay vertex; the $\chi^{2}$ of the ${D}^{+}$ decay vertex fit; the distance of closest approach between any two final-state tracks; and the momentum, the transverse momentum and the $\chi^{2}_{\rm IP}$ of the ${D}^{+}$ candidate and of its decay products. The invariant mass of the ${D}^{+}$ candidate is required to be within the interval 1820–1920$\mathrm{\,Me\kern-1.00006ptV}$. In order to suppress the contamination from $D^{+}_{s}\rightarrow K^{-}K^{+}\pi^{+}\pi^{0}$ decays, where the neutral pion is not reconstructed and the charged pion is misidentified as a kaon, more stringent PID requirements are applied to the kaon candidates with the same charge.

A boosted decision tree (BDT) multivariate classifier [21, 22] is used to further reduce the combinatorial background. In order to keep the selection efficiency uniform over the Dalitz plot, the BDT uses only the quantities related to the $D^{+}$ candidate described above. The BDT is trained using simulated ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decays for the signal, and data from the invariant-mass intervals 1820–1840$\mathrm{\,Me\kern-1.00006ptV}$ and 1900–1920$\mathrm{\,Me\kern-1.00006ptV}$ for the background. After the application of all selection requirements, approximately 0.5% of the events include more than one signal candidate. All candidates are retained for further analysis.

The invariant-mass spectrum of the selected ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ sample is shown in Fig. 1. To fit the invariant-mass distribution, the signal probability density function (PDF) is modeled by a sum of two Gaussian functions with a common mean and independent widths that are free parameters. The signal model is validated with simulation. The background PDF is parameterised by an exponential function. The fitted PDF is overlaid with the mass distribution in Fig. 1. For the Dalitz plot analysis, only candidates within the range 1861.4–1879.5 $\mathrm{\,Me\kern-1.00006ptV}$ are considered. This interval corresponds to four times the effective mass resolution, and contains 111 thousand candidates, of which $(90.45\pm 0.07)\%$ correspond to signal.

The Dalitz plot of the candidates in the signal region is shown in the left side of Fig 2. The particle ordering is such that the kaon with charge opposite to that of the $D^{+}$ meson is always particle 1, and the same-sign kaons are randomly assigned particles 2 and 3, i.e. $D^{+}\rightarrow K^{-}(p_{1})K^{+}(p_{2})K^{+}(p_{3})$, where $p_{i}$ are the four-momenta. The Dalitz plot is represented in terms of the square of the invariant masses of the two $K^{-}K^{+}$ combinations, $s_{12}\equiv(p_{1}+p_{2})^{2}$ and $s_{13}\equiv(p_{1}+p_{3})^{2}$. Throughout this paper, the symbol $s_{K^{-}K^{+}}$ is used to represent the invariant mass squared of both $K^{-}K^{+}$ combinations. These Lorentz-invariant quantities are computed constraining the invariant mass of the candidate to the known $D^{+}$ mass [23]. An accumulation of candidates is visible at $s_{K^{-}K^{+}}$ $\sim$1.04$\mathrm{\,Ge\kern-1.00006ptV^{2}}$ which corresponds to the $\phi(1020)K^{+}$ component. The difference in the number of candidates in the regions of the Dalitz plot above and below 1.55$\mathrm{\,Ge\kern-1.00006ptV^{2}}$ (regions I and II in the left side of Fig. 2, respectively) is caused by interference between the $\phi(1020)K^{+}$ and S-wave amplitudes. This interference also shifts the position of the peaks of the $s_{K^{-}K^{+}}$ distributions in the two regions. These two effects are better illustrated in the projections of the Dalitz plot shown in the right side of Fig. 2.

4 Efficiency and background model

4.1 Efficiency variation over the Dalitz plot

In the fit to the Dalitz plot distribution, the variation of the total efficiency across the phase space must be taken into account. The total efficiency is determined from a combination of simulation and methods based on data, and includes the geometrical acceptance of the detector and the reconstruction, selection, PID and trigger efficiencies.

The geometrical acceptance, reconstruction and selection efficiencies are obtained from simulation. The PID efficiency of each $D^{+}$ candidate is determined by multiplying the efficiencies for each of the final-state kaons. The PID efficiencies for the kaons are evaluated from calibration samples of $D^{*+}\rightarrow D^{0}(\rightarrow K^{-}\pi^{+})\pi^{+}$ decays [24] and depend on the particle momentum, pseudorapidity and event charged-particle multiplicity. The trigger efficiency is obtained from simulation, with a correction factor determined from data to account for the small mismatch between the performance of the trigger in data and simulation.

The total efficiency distribution is a two-dimensional histogram with $14\times 14$ uniform bins. A two-dimensional cubic spline is used to smooth this distribution to avoid binning discontinuities, yielding the high-resolution histogram ($300\times 300$ uniform bins), shown in Fig. 3. This histograms is used to weight the signal PDF in the Dalitz plot fit. The binning scheme of the efficiency histogram is a source of systematic uncertainty.

4.2 Background model

The background model is built from the inspection of the mass sidebands of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ signal. The Dalitz plots of candidates from both sidebands, 1820–1840$\mathrm{\,Me\kern-1.00006ptV}$ and 1900–1920$\mathrm{\,Me\kern-1.00006ptV}$, are very similar, with a clear peaking structure, corresponding to random $\phi(1020)K^{+}$ combinations over a smooth distribution.

The Dalitz plot variables are computed from the four-momenta determined by a $D^{+}$ mass constrained fit. This constraint implies an unique boundary of the Dalitz plot, regardless of the value of the invariant mass of the three-kaon system. It also improves the mass resolution of signal candidates, but has the effect to distort and shift any structure present in the Dalitz plot of the background candidates in the sidebands. This effect depends strongly on the invariant mass of the three-kaon system and prevents the determination of the background model from a two-dimensional parameterisation of the Dalitz plots from the sidebands. An alternative method is used instead. Each $m(K^{-}K^{+}K^{+})$ sideband is divided into slices of 5$\mathrm{\,Me\kern-1.00006ptV}$. For each slice, the projections onto the $s_{K^{+}K^{-}}$ axis are fitted using a relativistic Breit–Wigner for the $\phi(1020)$ component (with floated mass and width) and a phase-space distribution. The latter serves as a proxy for both the smooth component spread across the Dalitz plot and the projection of the $\phi$ candidates appearing in the other $s_{K^{+}K^{-}}$ combination. The fraction of the $\phi(1020)$ component is nearly constant in both sidebands, indicating that the background composition is independent of $m(K^{-}K^{+}K^{+})$. A linear interpolation is used to obtain the fraction of peaking background in the signal region and is found to be (20.67$\pm$0.28)%.

The fit to the $s_{K^{+}K^{-}}$ projection has the limitation of being less sensitive to the distribution near the $K^{+}K^{-}$ threshold. The inspection of the Dalitz plot sidebands shows that the smooth background component has more candidates at low values of $s_{K^{+}K^{-}}$ and fewer at low values of $s_{K^{+}K^{+}}\equiv(p_{2}+p_{3})^{2}$, indicating that this smooth distribution is not uniform over the phase space. A model for the smooth component of the background is built assuming a sum of two contributions, random $f_{0}(980)K^{+}$ candidates and a constant amplitude, with equal proportions. The relative fractions of these two terms in the smooth component is treated as a source of systematic uncertainty.

A high-resolution normalised histogram ($300\times 300$ uniform bins) is used in the Dalitz plot fit to represent the background PDF, and is shown in Fig. 5. This histogram is produced from a large simulated sample, using a PDF in which the peaking and smooth components are added incoherently with the estimated relative fractions and weighted by the efficiency function.

5 The Dalitz plot fit procedure

The ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decays are studied through an unbinned maximum-likelihood fit to the observed Dalitz plot distribution. The total PDF is constructed as a sum of signal and background components, and the likelihood function is given by

[TABLE]

where $N_{\rm cand}$ is the total number of candidates and $f_{S}$ is the fraction of signal candidates in the sample, as obtained from the $m(K^{-}K^{+}K^{+})$ fit described in Sec. 3. The background PDF, $B_{\rm PDF}(s_{12},s_{13})$, is described in Sec. 4.2.

The normalised signal PDF is written in terms of the total decay amplitude ${\mathcal{T}}(s_{12},s_{13})$,

[TABLE]

where $\varepsilon(s_{12},s_{13})$ is the detection efficiency, described in Sec. 4.1. The normalisation factor, $N_{\rm S}$, is given by

[TABLE]

For any given model, the amplitude ${\mathcal{T}}(s_{12},s_{13})$ depends on a set of parameters that are floated in the fit. The optimum values for these parameters are determined by minimizing the quantity $-2\ln\mathcal{L}$ using the MINUIT package [25].

In order to compare the fit results of a given model to the Dalitz plot distribution in data, a large simulated sample is generated according to the model, including background and efficiency, normalised to the total number of data candidates. Since there are two identical kaons, the folded Dalitz plot is used, represented as $s_{K^{+}K^{-}}^{\rm high}$ versus $s_{K^{+}K^{-}}^{\rm low}$, which are respectively the higher and the lower values among $s_{12}$ and $s_{13}$. The Dalitz plot distribution is divided into 1024 bins with approximately 110 candidates each and the normalised residuals are computed as

[TABLE]

where, for each bin $i$, $N_{\rm pred}^{i}$ is the predicted number of candidates from the model, $N_{\rm obs}^{i}$ is the number of candidates in the data sample, and $\sigma_{i}$ is the statistical uncertainty from data and simulation added in quadrature. The sum of the square values of $\Delta_{i}$ over all bins is the total $\chi^{2}$ and is used as a metric to compare fit results with different models.

6 Dalitz plot analysis with the isobar model

In the isobar model, the decay amplitude is written as a coherent sum of a constant nonresonant (NR) component and intermediate resonant amplitudes,

[TABLE]

Each resonant amplitude, $T_{k}$, is given by a product of Blatt–Weisskopf penetration factors [26], $F_{D}^{L}$ and $F_{R}^{L}$, accounting for the finite size of the $D^{+}$ meson and the resonance, respectively, the spin amplitude, $\mathcal{S}$, accounting for the conservation of angular momentum, and a function, $M_{R}$, describing the resonance lineshape, which is either a relativistic Breit–Wigner (Eq. A.2) or a Flatté lineshape (Eq. A.4). The Zemach formalism [27] is used for the spin amplitude $\mathcal{S}$. Details of each of these factors are given in Appendix A. Since there are two identical kaons in the final state, the resonant amplitudes are Bose-symmetrised,

[TABLE]

The fit parameters are the complex coefficients $c_{\rm NR}=a_{\rm NR}e^{i\delta_{\rm NR}}$ and $c_{k}=a_{k}e^{i\delta_{k}}$. The results are expressed in terms of the magnitude and phase of the complex coefficient for each component, and the corresponding fit fractions. The fit fractions are computed by integrating the squared modulus of the corresponding amplitude over the phase space, and dividing by the integral of the total amplitude squared,

[TABLE]

The sum of fit fractions for all components is, in general, different from 100% due to the presence of interference; it is less than 100% in the case of net constructive interference or higher than 100% otherwise.

6.1 Signal models

For the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay amplitude, contributions from following resonances are possible: the isoscalars $f_{0}(980)$, $f_{0}(1370)$ and $f_{0}(1500)$; the isovectors $a_{0}(980)$ and $a_{0}(1450)$; the vector $\phi(1020)$; the tensor $f_{2}(1270)$. Contributions from resonances with spin greater than one are suppressed due to the small phase space of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay. In the case of the $f_{2}(1270)$ state, a further suppression is expected due to its small branching fraction to $K^{-}K^{+}$, $(4.6\pm 0.4)$% [23]. The relatively narrow $f^{\prime}_{2}(1525)$ state is neglected since it is well beyond the phase space.

Various combinations of the nonresonant and the possible resonant amplitudes are considered. All models studied contain the $\phi(1020)K^{+}$, which is chosen as the reference amplitude, fixing the phase convention and setting the scale for the magnitudes. The models tested differ by the composition of the S-wave. Near the $K^{+}K^{-}$ threshold, both the $a_{0}(980)$ and $f_{0}(980)$ resonances can contribute. Similarly, at higher $K^{+}K^{-}$ invariant mass, contributions from several scalar resonances are possible.

The $\phi(1020)$ mass and width are fixed to the known values [23]; for the $f_{0}(980)$ state, a Flatté lineshape is used, with parameters from the BESII collaboration [28].

6.2 Results

The simplest model that describes the data, referred to as model A, consists of three intermediate components: $\phi(1020)K^{+}$, $f_{0}(980)K^{+}$, and $f_{0}(1370)K^{+}$. As the $f_{0}(1370)$ state has large uncertainties on its mass and width [23], these parameters are allowed to float in the fit. Its contribution can also be interpreted, within the isobar formalism, as an effective representation for the overlap of two or more broad structures at high $K^{-}K^{+}$ invariant mass.

Further addition of scalar states does not improve the fit quality significantly, creates more complex interference effects, and provides a very similar description of the lineshape and phase behaviour of the total S-wave. For example, in model B, a constant nonresonant contribution is added to the resonant amplitudes of model A. The resulting fit quality is essentially unchanged, with the total $\chi^{2}/{\rm ndof}$ being 1.15 and 1.14 for models A and B, respectively. A similar situation occurs in model C, which has the same amplitudes as in model B plus the $a_{0}(980)K^{+}$ component. In this model, the contribution of the $f_{0}(1370)$ is found to be negligible and the value of $\chi^{2}/{\rm ndof}$ is 1.16. Table 1 summarizes the fit results for these three models. The total S-wave fit fraction includes the interference terms between the various S-wave components. In all cases, the total S-wave in the $K^{+}K^{-}$ system is dominant, a notable feature also observed in other three-body $D$ decays with a pair of identical particles in the final state, such as the $D^{+}\rightarrow K^{-}\pi^{+}\pi^{+}$ and $D^{+}_{(s)}\rightarrow\pi^{-}\pi^{+}\pi^{+}$ decays[23]. The contribution from the $f_{2}(1270)K^{+}$ component is also tested and found to be consistent with zero in all models.

Since model A is the simplest model describing all the general features of the observed Dalitz plot distribution, it is chosen as the baseline result for the fit with the isobar model. The projections of the Dalitz plot, with the model A fit result overlaid, are shown in Fig. 6. The green dashed line represents the phase-space distribution, weighted by the efficiency, evidencing the presence of at least one broad, scalar contribution not consistent with a uniform distribution.

The distribution of the normalised residuals $\Delta_{i}$ over the Dalitz plot is shown in the left plot of Fig. 7, and their distribution is consistent with a normal Gaussian, as shown in the right plot of Fig. 7. In Table 2 the results including the systematic and model uncertainties, as discussed in Sec. 8, are presented.

The squared modulus and phase of the S-wave amplitude from model A are shown in Fig. 8 as a function of the $K^{+}K^{-}$ mass, with total uncertainties represented as bands. For comparison, the corresponding central results for models B and C are overlaid. Although the S-wave composition is different for these models, the total S-wave description is essentially the same, evidencing that the isobar model fails to disentangle the individual contributions. The $f_{0}(1370)$ parameters are found to be m_{0}=1.422\pm 0.015\pm 0.009\pm 0.028$$\mathrm{\,Ge\kern-1.00006ptV} and \Gamma_{0}=0.324\pm 0.038\pm 0.018\pm 0.038$$\mathrm{\,Ge\kern-1.00006ptV}, where the first uncertainties are statistical and the second systematic.

7 Dalitz plot analysis with the Triple-M amplitude

The basic hypothesis of the Triple-M is the dominance of the annihilation diagram shown in Fig. 9(a). The ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay can also proceed via the diagram in Fig. 9(b), but in this case a $K^{+}K^{-}$ pair could only be produced from the $d\bar{d}$ pair through rescattering, since charged kaons have no $d$-valence quark. The same holds for the production of the $\phi(1020)$ meson which is essentially an $s\bar{s}$ state [23].

Assuming the annihilation diagram is the dominant mechanism for the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay, the Triple-M amplitude is a product of two axial-vector currents,

[TABLE]

where $G_{F}$ is the Fermi decay constant, $\theta_{C}$ is the Cabibbo angle and $A^{\mu}$ are the axial currents. The weak vertex is $\langle\,0\,|A_{\mu}|\,D^{+}(P)\rangle=-i\,\sqrt{2}\,f_{D}\,P_{\mu}$, where $P=p_{1}+p_{2}+p_{3}\,$ is the $D^{+}$ four-momentum and $f_{D}$ is the $D^{+}$ decay constant.

In the Triple-M, the three-kaon system can be produced in two ways, as illustrated in the diagrams in Fig. 10. Diagram (a) represents the production of the three kaons directly from the weak vertex, whereas in diagram (b) two of the three kaons result from the decay of a bare intermediate resonance. Final state interactions are introduced in diagrams (c) and (d). The full black circles indicate the unitarised scattering amplitudes, $A_{K^{+}K^{-}}^{JI}$, representing the scattering $ab\rightarrow K^{+}K^{-}$ with the coupled channels $ab=K^{+}K^{-},\ \pi\pi,\ \eta\pi$ and $\eta\eta$ in a well-defined spin (J ) and isospin (I ) state. The nonresonant component corresponds to diagram (a). Due to the existence of two identical kaons, diagrams (b), (c) and (d) are symmetrised. As in the isobar analysis, contributions of D-wave are expected to be very small and are not included.

The Triple-M decay amplitude therefore has five components,

[TABLE]

The free parameters in the Triple-M amplitude are the couplings and masses of the chiral Lagrangian. There are four couplings, $c_{d},\ c_{m},\ \tilde{c}_{d},\ \tilde{c}_{m}$ in the scalar part, contributing to $T^{00}$ and $T^{01}$ terms; two masses, $m_{So},\ m_{S1}$, for the scalar-isoscalar, $T^{00}$ contribution and one, $m_{a_{0}}$, in the scalar-isovector $T^{01}$ components; one coupling, $G_{V}$, for the vector components, $T^{10}$ and $T^{11}$, and one mass, $m_{\phi}$, in the vector-isoscalar component. In the fit to the data, the combination $G_{\phi}\equiv G_{V}\sin\theta_{\omega-\phi}/F$ is used as free parameter, where $\theta_{\omega-\phi}$ is the $\omega-\phi$ mixing angle. The parameter $F$ is the $SU(3)$ pseudoscalar decay constant, common to all components. For convenience, the formulae of the various components of the Triple-M amplitude are reproduced from Ref. [3] in Appendix B.

Equation 9 resembles that of the isobar model, but there are several significant differences. The free parameters in the Triple-M amplitude are real quantities from the chiral Lagrangian. Some of these parameters appear in different spin-isospin components of the model. In the isobar model the free parameters are the complex coefficients $c_{k}$, from which the individual contributions of the resonances are determined. In the Triple-M amplitude, the relative contributions of the various components are fixed by theory. The nonresonant component is usually represented by an empirical constant in fits with the isobar model. In the Triple-M amplitude, it is a function of the Dalitz plot coordinates and is fully determined by chiral symmetry.

7.1 Fit results

The optimum values of the Triple-M parameters are determined by an unbinned maximum-likelihood fit, as described in Sec. 5. The fitted values of the Triple-M parameters are listed in Table 3, with statistical and systematic uncertainties.

The quality of the fit with the Triple-M amplitude is tested with the metric defined in Eq. 4. The value of $\chi^{2}/{\rm ndof}$ is 1.12. The projections of the Dalitz plot onto the $s_{K^{+}K^{-}}$ and the $s_{K^{+}K^{+}}$ axes, as well as the projections onto the highest and lowest invariant masses squared of the two $K^{+}K^{-}$ combinations, $s_{K^{+}K^{-}}^{\rm high}$ and $s_{K^{+}K^{-}}^{\rm low}$, are shown in Fig. 11, with the fit result superimposed. The projections indicate that the model is in good agreement with the data. The distribution of the normalised residuals over the Dalitz plot is shown in the right panel of Fig. 12. The distribution of normalised residuals, shown in the left panel of Fig. 12, is consistent with a normal Gaussian.

7.2 Interpretation

The resonance masses in the Triple-M are introduced in the denominators, D, of Eqs. B.21–B.24, where the functions M are imaginary and proportional to interaction kernels which contain the bare masses of the effective chiral Lagrangian, $m_{a_{0}}$, $m_{S_{o}}$, $m_{S_{1}}$ and $m_{\phi}$. The Triple-M amplitude is derived assuming that only the imaginary part of the two-body propagators in Eqs. B.25–B.28 is relevant. In this approximation, the bare masses coincide with the masses of the physical states and the association $m_{S_{o}}=m_{f_{0}(980)}$ and $m_{S_{1}}=m_{f_{0}(1370)}$ can be made. As in the case of the isobar model, the masses in the Triple-M correspond to the values of $s_{K^{+}K^{-}}$ for which the real part of the denominator D of Eqs. B.21–B.24 vanishes. At these values of $s_{K^{+}K^{-}}$, only the imaginary parts of the denominators remain, corresponding to the model prediction for the widths. The denominators D would be very similar to those from the isobar model if no coupled channel was considered. The inclusion of coupled channels is, therefore, the main difference between the Triple-M and Breit–Wigner denominators, resulting in widths with different dynamical content.

7.2.1 Resonant structure

The nonresonant contribution in the Triple-M is a three-body amplitude predicted by chiral symmetry. It can be projected into the S- and P-waves rewriting Eq. B.3 as

[TABLE]

where $C$ is a constant common to all components of the Triple-M amplitude, and defined in Eq. B.2. The decay amplitude can then be written as the sum of scalar and vector components

[TABLE]

with

[TABLE]

and

[TABLE]

The relative contribution of each individual component of the Triple-M amplitude is determined by integrating the modulus squared of each term in the right-hand side of Eq. 9 over the phase space of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ decay,

[TABLE]

Similarly, the S-wave contribution can be determined by the integral over the phase space of the modulus squared of the $T^{S}$ component, defined in Eq. 12, and divided by the integral of the modulus squared of the decay amplitude $\mathcal{T}$. The results are summarised in Table 4. There is a large destructive interference between the two scalar below-threshold states, $a_{0}(980)$ and $f_{0}(980)$, yielding an S-wave contribution of $(94\pm 1)$%. The large $a_{0}(980)/f_{0}(980)$ interference may be, in part, due to the fact that in the $K^{+}K^{-}$ mass spectrum these two states have very similar lineshapes, since only the tails are visible. This large interference is also observed in the fit with the isobar model C, yielding similar fit fractions for the S-wave component. A more accurate determination of the relative contribution of the $a_{0}(980)$ and $f_{0}(980)$ resonances could be obtained from a simultaneous analysis of the $D^{+}\rightarrow\pi^{+}\pi^{-}\pi^{+}$ and $D^{+}\rightarrow\eta\pi^{+}\pi^{0}$. The contribution of the $\phi(1020)$ resonance, $(7.1\pm 0.5)\%$, is consistent to that observed in the fit with the isobar model.

7.2.2 Decay and scattering amplitudes

The phases of the S-wave amplitude, $T^{S}$, and the $K^{+}K^{-}\rightarrow K^{+}K^{-}$ scattering amplitudes, $A^{0I}_{K^{+}K^{-}}$, for the two allowed isospin states, are shown in Fig. 13 as a function of the $K^{+}K^{-}$ invariant mass. The bands correspond to the statistical and systematic uncertainties added in quadrature. The kink in the phase of $T^{S}$ at $m(K^{+}K^{-})\sim 1.25$ GeV is due to the opening of the $\eta\eta$ channel. The curves of Fig. 13 illustrate the difference between decay and scattering amplitudes. The latter, which depends on spin and isospin, is a substructure of the former, which depends only on spin. The expressions of the various scattering amplitudes, derived in Ref. [3], are reproduced in Appendix C.

The physics of two-body scattering is encompassed by the phase shifts and inelasticities. These quantities are obtained from the scattering amplitudes, following the procedure described in Ref. [3]. The phase shifts, $\delta^{JI}_{K^{+}K^{-}}$, and inelasticities $\eta^{JI}_{K^{+}K^{-}}$, are displayed in Fig. 14 for $J\!=\!0$ and $I\!=\!0,1$.

The interpretation of the phase shifts for $K^{+}K^{-}$ scattering is not as straightforward as in the case of elastic scattering, since for both isospin states, the $\pi\pi\rightarrow K^{+}K^{-}$ and $\pi\eta\rightarrow K^{+}K^{-}$ channels are already open at the $K^{+}K^{-}$ threshold. An interesting feature of the results displayed is that the phase variation of $\delta^{00}_{K^{+}K^{-}}$ is monotonic and spans over more than $180^{\circ}$, with a fast variation starting at $m(K^{+}K^{-})\sim 1.4$ GeV, close to the value of $m_{S_{1}}$ and typical of a resonance at high $K^{+}K^{-}$ mass. A fast variation of the phases is observed near threshold for both $\delta^{00}_{K^{+}K^{-}}$ and $\delta^{01}_{K^{+}K^{-}}$, indicating the contribution from the resonances below threshold.

The $\eta\eta$ channel contributes to $T^{00}$ but not to $T^{01}$ and its effect is visible in the bottom left plot of Fig. 14 as a kink at m(K^{+}K^{-})\sim 1.1$$\mathrm{\,Ge\kern-1.00006ptV}. As elastic scattering corresponds to $\eta^{JI}=1$, one sees that the isoscalar component becomes significantly more inelastic after the mass of the second scalar resonance.

8 Systematic uncertainties

Sources of systematic uncertainties associated to the background model, to the efficiency correction and to possible biases in the fitting procedure are common to the fits with the isobar model and the Triple-M. They are summarized in Tables 5 and 6, respectively. There is an additional source of systematic uncertainties on the results of the fit with the isobar model due to the uncertainties on the parameters defining the $f_{0}(980)$ lineshape, which are fixed in the fit. This additional uncertainty, quoted separately from the experimental uncertainties, is estimated by repeating the fit varying the parameters $g_{\pi}$, $g_{K}$ and $m_{0}$ of Eq. A.4 by one standard deviation, one at a time, and taking the largest deviation as the systematic uncertainty. The radii of the Blatt–Weisskopf form factors are also fixed in the fit. However, they impacts only the $\phi(1020)K^{+}$ amplitude. Fits with alternative values of these parameters are performed. The tested values of the radii are 4 and 6 GeV*-1*, for $F_{D}^{L}$, and 1 and 3 GeV*-1*, for $F_{R}^{L}$. Since no significant deviation from the baseline fit is observed, no systematic uncertainty is assigned.

Two types of systematic uncertainties due to the background are investigated. First, the background level is varied according to the uncertainty from the fit to the $K^{-}K^{+}K^{+}$ invariant mass. The data is fitted changing the fraction of the background by $\pm 1\sigma$. No significant change in the fit parameters is found and no systematic uncertainty is assigned. Uncertainties due to the background modelling are also investigated. The background model is built from inspection of the sidebands of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ signal. It is a combination of a peaking structure and a smooth component. The smooth component corresponds to 80% of the background and is modelled by a sum of a constant term and an $f_{0}(980)K^{+}$ contribution, in equal proportions. A systematic uncertainty due to the modelling of the background is assigned by varying the relative fractions of these two components, fitting the data with these alternative background models and taking the largest variation as systematic uncertainty.

Systematic uncertainties are assigned to small biases in the fit using ensembles of 500 simulated samples. Two sets of samples are generated using the Triple-M amplitude and the isobar model, both with the fitted values of the parameters. In the simulations the signal PDFs are weighted by the efficiency function and the background component is included. Each simulated sample is fitted independently, resulting in distributions of fitted values of the parameters and their respective uncertainties. For each parameter, the mean of the distribution of fitted values is compared to the input. The difference is assigned as the systematic uncertainty due to the fit bias. A small bias is observed in the fit with the Triple-M amplitude, whilst no bias is observed in the fit with the isobar model.

The systematic uncertainty associated to the efficiency variation across the Dalitz plot includes the effect of the uncertainties on the PID efficiency and the hardware trigger correction factors, the effect of the finite size of the simulated sample, and the effect of the binning scheme of the efficiency histogram prior to the two-dimensional spline smoothing. The uncertainties on the PID efficiency are due to the finite size of the calibration samples and imply small systematic uncertainties compared to the other sources of systematics, in the fit with the Triple-M amplitude, and negligible uncertainties in the fit with the isobar model. The uncertainty due to the hardware trigger correction factors is found to be negligible. The effect of the finite size of the simulated sample is assessed by generating a set of alternative histograms from the selection efficiency histogram, prior to the hardware trigger correction and the PID efficiency weighting. The content of each bin of the selection efficiency histogram is varied according to a Poisson distribution. For each of these alternative histograms, an efficiency map is produced and used to fit the data. For each parameter, the root mean square of the distribution of fitted values is assigned as a systematic uncertainty. The systematic uncertainty due to the binning scheme of the efficiency map is accessed by varying the number of bins of the final efficiency histogram. The histograms with alternative binnings are fitted by the two-dimensional cubic spline. The data is fitted with these alternative efficiency maps and the largest variation of each parameter is assigned as systematic uncertainty.

9 Summary and conclusions

In this paper, the first Dalitz plot analysis of the doubly Cabibbo-suppressed decay ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ is performed. The two goals of the analysis are the determinations of the resonant structure of the decay and the $K^{+}K^{-}$ scattering amplitudes. The resonant structure is studied with two different approaches. In the fit with the isobar model, several variations of the decay amplitude are tested. The Dalitz plot analysis is also performed with the Triple-M [3], which is a model derived from a chiral effective Lagrangian. The Triple-M amplitude has a nonresonant component plus the minimal $SU(3)$ content corresponding to four states, the $\phi(1020)$, the $a_{0}(980)$ and two isoscalar states, identified with the $f_{0}(980)$ and $f_{0}(1370)$ resonances. A good description of the data is achieved with both approaches.

The resonant structure of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$ is largely dominated by the S-wave, with a approximately 7% contribution from the $\phi(1020)K^{+}$ component. The dominance of the S-wave contribution is also observed in other three-body $D^{+}_{(s)}$ decays with a pair of identical particles in the final state, such as the $D^{+}\rightarrow K^{-}\pi^{+}\pi^{+}$ and $D^{+}_{(s)}\rightarrow\pi^{-}\pi^{+}\pi^{+}$ decays[23]. The possibility of determining the individual components of the S-wave, however, is limited by the lack of structures in the Dalitz plot, other than that from the $\phi(1020)$ resonance, and by the fact that the $f_{0}(980)$ and $a_{0}(980)$ mesons poles lie below the $K^{+}K^{-}$ threshold. In all the models tested, large interference between the various S-wave components is observed. In the fit with isobar model, different combinations of scalar resonances and nonresonant amplitudes yield fits of same quality and a very similar S-wave amplitude. In the fit with the Triple-M, a large $a_{0}(980)$ contribution is observed, with a large destructive interference with the $f_{0}(980)$ component that yields an S-wave fraction of about 94%. The separation between the $f_{0}(980)$ and $a_{0}(980)$ contributions could better achieved with a simultaneous analyses of the ${{D}^{+}}\!\rightarrow{{K}^{-}}{{K}^{+}}{{K}^{+}}$, $D^{+}\rightarrow\pi^{-}\pi^{+}\pi^{+}$ and $D^{+}\rightarrow\eta\pi^{+}\pi^{0}$ decays.

Predicitions for the $K^{+}K^{-}\rightarrow K^{+}K^{-}$ scattering amplitudes are obtained from the Dalitz plot fit using the Triple-M amplitude. This is possible because the model incorporates explicitely coupled channels and isospin degrees of freedom. In this respect, the chiral Lagrangian approach represents an advance towards the description of the hadronic part of weak decays of $D$ mesons in a more fundamental basis.

Acknowledgements

We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); MOST and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG (Germany); INFN (Italy); NWO (Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania); MSHE (Russia); MinECo (Spain); SNSF and SER (Switzerland); NASU (Ukraine); STFC (United Kingdom); NSF (USA). We acknowledge the computing resources that are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN (Italy), SURF (Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and Yandex LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), PL-GRID (Poland) and OSC (USA). We are indebted to the communities behind the multiple open-source software packages on which we depend. Individual groups or members have received support from AvH Foundation (Germany); EPLANET, Marie Sklodowska-Curie Actions and ERC (European Union); ANR, Labex P2IO and OCEVU, and Région Auvergne-Rhône-Alpes (France); Key Research Program of Frontier Sciences of CAS, CAS PIFI, and the Thousand Talents Program (China); RFBR, RSF and Yandex LLC (Russia); GVA, XuntaGal and GENCAT (Spain); the Royal Society and the Leverhulme Trust (United Kingdom); Laboratory Directed Research and Development program of LANL (USA).

Appendix

Appendix A Decay amplitudes in the isobar model

Intermediate decay amplitudes within the Isobar model are given by Eq. 6. Each factor appearing in that equation is presented here.

The form factors $F_{D}^{L}$ and $F_{R}^{L}$, for the $D^{+}$ and the resonance decay, respectively, are parameterized by the Blatt–Weisskopf penetration factors [26], and depend on $L$, the orbital angular momentum involved in the transition. Since both the initial state (the $D^{+}$ meson) and the final state (three kaons) have spin 0, $L$ is equal to the spin of the resonance. In the rest frame of a resonance formed by particles 1 and 2 , R12, $q$ is the modulus of the momentum of particle 1 or 2 (the decay momentum), $q_{0}$ is the decay momentum when $s_{12}=m_{R}^{2}$ ($m_{R}$ being the nominal resonance mass), and $d$ is a measure of the effective radius of the decaying meson, fixed in this work to $5.0\mathrm{\,Me\kern-1.00006ptV}^{-1}$ for the $D$ meson and $1.5\mathrm{\,Me\kern-1.00006ptV}^{-1}$ for the resonance. Defining $z=(qd)^{2}$ and $z_{0}=(q_{0}d)^{2}$, the Blatt–Weisskopf barrier factors are usually written with two different formulations, $B_{L}$ and $B_{L}^{\prime}$ [23], given in Table 7. The $B_{L}^{\prime}$ formulation is used in this analysis, consistent with the energy dependent width given below in Eq. A.3, with the momenta in $F_{D}^{L}$ and $F_{R}^{L}$ computed in the rest frame of the respective decaying particle.

The function $\mathcal{S}(\theta^{R12}_{13})$ describes the angular distribution of the decay particles, with $\theta^{R_{12}}_{13}=\theta^{R_{12}}_{13}(s_{12},s_{13})$ being the angle between particles 1 and 3 momenta measured in the rest frame R12. The Zemach formalism [27] is used for the angular distribution

[TABLE]

where $P_{L}$ is the Legendre polynomial of order $L$. For vector and tensor resonances, this term introduces nodes in the Dalitz plot in regions where the helicity angle is either 90*∘* or 270*∘*.

The relativistic Breit–Wigner function [29] is used as the dynamical function,

[TABLE]

where $m_{R}$ is the mass of the resonance and $\Gamma(s_{12})$ is the mass-dependent width,

[TABLE]

with $\Gamma_{R}$ being the nominal resonance width.

In the case of the $f_{0}(980)$ resonance, the relativistic Breit–Wigner is replaced by the Flatté formula[30]

[TABLE]

where $g_{\pi}$ and $g_{K}$ are dimensionless coupling constants to the $K{\kern 1.99997pt\overline{\kern-1.99997ptK}{}}$ and $\pi\pi$ channels, respectively, and $\rho_{\pi\pi}$ and $\rho_{KK}$ are the corresponding phase-space factors,

[TABLE]

All the above formulation holds equally for the resonances in the system composed by particles 1 and 3, with $s_{12}\rightarrow s_{13}$ and $\theta^{R_{12}}_{13}\rightarrow\theta^{R_{31}}_{12}$ (angular functions convention with cyclic permutation $(12)3\rightarrow(31)2$).

Appendix B The Triple-M Decay amplitude

All formulae presented in this appendix are reproduced from Ref. [3] for convenience. The Triple-M decay amplitude for the $D^{+}\rightarrow K^{-}\,K^{+}\,K^{+}$ decay is given by

[TABLE]

where $T_{\rm NR}$ and the $T^{(J,I)}$ are the nonresonant and resonant contributions, respectively. All components are proportional to the kaon mass squared, $m_{K}^{2}$, included in the common factor

[TABLE]

where $F_{D}$ is the $D^{+}$ decay constant, $F$ is the $SU(3)$ pseudoscalar decay constant, $G_{F}$ is the Fermi decay constant and $\theta_{C}$ is the Cabibbo angle. The nonresonant contribution is a three-body amplitude, and therefore is not Bose-symmetrised. It is written as a real polynomial,

[TABLE]

The amplitudes $T^{(J,I)}$ are

[TABLE]

[TABLE]

The various functions $\Gamma_{(0)ab}^{(J,I)}$ correspond to diagrams ($a$) and ($b$) of Fig. 10, and represent the tree-level production of particles $abK^{+}$ from the weak vertex. The functions $\bar{\Gamma}_{KK}^{(J,I)}$ represent the full decay vertex, from which the decay amplitude is obtained after subtracting the contribution of the contact terms $\Gamma_{c|KK}^{(J,I)}$ to avoid double counting. Their explicit form of the $\Gamma_{(0)ab}^{(J,I)}$ functions are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

In the above equations, $m_{\pi}$ and $m_{D}$ are the $\pi^{+}$ and the $D^{+}$ masses, respectively, and $\theta$ is the $\omega-\phi$ mixing angle. The subscripts $8$ refer to the member of the $SU(3)$ octet with the quantum numbers of the $\eta$. The denominators in Eqs. B.5, B.7, B.9 and B.11 are the model prediction for the resonance line shapes:

[TABLE]

The functions $M_{ij}^{(J,I)}$ read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The imaginary propagators $\bar{\Omega}$ are given by

[TABLE]

The functions ${\cal{K}}_{ab|cd}^{(J,I)}$ are the scattering kernels,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix C Scattering amplitudes

The $K^{+}K^{-}$ scattering amplitudes are written in terms of the denominators $D^{(J,I)}$ as

[TABLE]