Scalar resonances in the $D^+\to K^-K^+K^+$ decay

We study theoretically the resonant structure of the double Cabibbo suppressed $D^+\to K^- K^+K^+$ decay. We start from an elementary production diagram, considered subleading in previous approaches, which cannot produce a final $K^- K^+$ pair at the tree level but which we show to be able to provide the strength of the decay through final meson-meson state interaction. The different meson-meson elementary productions are related through SU(3) and the final rescattering is implemented from a suitable implementation of unitary extensions of ChPT which generate dynamically the scalar resonances $f_0(980)$ and $a_0(980)$. We obtain a good agreement with recent experimental data from the LHCb collaboration with a minimal freedom in the fit and show the dominance of the $a_0(980)$ contribution close to the threshold of the $K^- K^+$ spectrum.


I. INTRODUCTION
Weak decay of heavy mesons into hadrons has become an important source of information on the hadron-hadron interaction. In particular the decay of D mesons into three pseudoscalars has drawn the attention of different groups with the aim of learning about the meson-meson interaction [1][2][3][4][5][6]. These works deal about the D → Kππ decay, which is Cabibbo favored and is mainly used to learn about the ππ interaction. Other D decays are studied in [7], using the D → π + π + π − and D → π + K + K − reactions to learn about the ππ and KK interaction, in [8,9] interpreting the D + s → π + π 0 η decay, in [10] studying the single Cabibbo suppressed D + → π + π 0 η decay, or in [11] studying the D + s → π + π 0 a 0 (980)(f 0 (980) reactions. The D 0 → K − π + η reaction is also studied in [12] from where information on the a 0 (980) and κ(K * 0 (700)) is obtained.
The reaction that we study here is D + → K − K + K + , which is doubly Cabibbo suppressed, but which can teach us much about the KK interaction, one of the pseudoscalar interaction channels most poorly known. The reaction has been studied by the LHCb collaboration in [13] and analyzed using two methods, the standard one, the isobar model, and then the triple-M model developed in [14]. The isobar model is the standard method used in the LHCb analysis and in most of the experimental collaborations. The full decay amplitude is written in terms of the only two independent variables T (s 12 , s 13 ) = c N R + k c k T k (s 12 , s 13 ), (1) where c N R is a non resonant background term and T k are * Electronic address: luisroca@um.es † Electronic address: oset@ific.uv.es intermediate resonant amplitudes properly parametrized. The parameters in the different terms and the complex weights c k are obtained by performing a best fit to the Dalitz plot data. The method is efficient to extract information on the role played by different resonances, but has its limitations. Using words of Ref. [14] "This approach, albeit largely employed [15], has conceptual limitations.
The outcome of isobar model analyses are resonance parameters such as fit fractions, masses and widths, which are neither directly related to any underlying dynamical theory nor provide clues to the identification of two-body substructures. Thus, the systematic interpretation of the isobar model results is rather difficult." Steps to make different analyses of the data to allow a better matching with theoretical tools used in the study of meson interactions have been done in [16], and tools to use three body dynamics have also been used in the D → Kππ reactions [1,2,[17][18][19]]. Yet, the majority of analyses rely upon the consideration of two body amplitudes having one of the mesons as spectator and this is also used in [14]. The fact that the three body amplitudes can be constructed from on shell two body amplitudes, since off shell parts are shown to cancel with contact terms present in the theory, makes this approach more realistic [20,21]. However, there are other reasons to neglect terms involving explicitly three particles interacting because after the interaction of a pair of mesons in regions where a resonance appears to be important, the resulting invariant mass of one of the particles of the pair with the third one does not have a given value but usually spans a large region of invariant mass, thus diluting the possible contribution of another resonance, which, however, is taken into consideration with the direct interaction of this original pair considering the third particle as spectator.
The work of [14] uses effective Lagrangians to deal with the weak and strong interaction. For the weak interaction the starting point is the diagram of Fig. 1(b), which involves quark pair annihilation (W-annihilation). The diagram of Fig. 1(a), which involves external emission, is considered as a possible mechanism, but it does no provide K − K + K + upon hadronization of the quarks, and final state interaction is needed to produce this state. For this reason it is neglected in the analysis of [14] and the annihilation mechanism of Fig. 1(a) is used as the starting point. The mesons stemming from the double hadronization of this mechanism are allowed to follow final state interaction. The final state interaction of the mesons in [14] is done using Lagrangians of [22] in which chiral perturbation is used including resonances explicitly. Yet, extra unitarization is considered in [14], using techniques of the chiral unitary approach of [23][24][25], which justifies that the parameters that they get from a fit to the data are not the parameters used in [22]. The unitarization is however done with one approximation with respect to the former works, using the K-matrix approach in which the real part of the loop functions is neglected and only the imaginary part is kept. The use of the K-matrix approach has one intrinsic problem, which is that ignoring the real part of the loops prevents that resonances are dynamically generated. This is why they have to be introduced by hand, and the approach cannot tell us about the nature of these resonances. However, resonances like the f 0 (980) and a 0 (980) appear as a consequence of the meson-meson interaction in the chiral unitary approach of [23][24][25] and the consistency of the approach with data gives support to this picture. We should mention that in a recent work [26], the authors or Ref. [14] are already considering the real parts of the loops.
In the present approach to the D + → K − K + K + reaction, we only have three parameters to fit to the data: the global strength, the strength of the φK + production amplitude and a relative phase of the s-wave to p-wave mechanisms. Since the global strength is irrelevant when we compare to events in the data, and the global strength of the φK + is easily determined from the clean peak in the K + K − mass distribution, our approach has basically one degree of freedom to fit all the data. This contrasts with the 10 free parameters that one has in the approach of [14]. The isobar model has even more parameters.
We have here two of the most important differences between the work of [14] and the one we present here. The starting point for us is not the diagram of Fig. 1(b) but the one of Fig. 1(a). The reason is the following: In the classification of weak decay topologies of [27,28] the order of importance is: external emission, internal emission, W-exchange and annihilation. Given the fact that the KK interaction al low energies is driven by the f 0 (980) and a 0 (980) resonances, the final state interaction is very important and does not destroy the order of the primary decays. In the analysis of the BESIII data for D + s → π + π 0 η, [29], the process removing the ρ + η channel, was supposed to proceed via W-annihilation with a rate of an order of magnitude bigger than for usual W-annihilation processes. Yet, the decay was studied in [8] with an approach similar to the one we follow here and it was shown that the process proceeded via internal emission and final state interaction.
Our approach to these weak processes consist in a first identification of the dominant mechanisms at the quark level, then we proceed with hadronization to produce the mesons that appear in a first step, and then consider the final state interaction of these mesons to produce at the end the desired final state. We must take into account that the hadronization, including new qq pairs to produce mesons, results in a reduction factor in the decay amplitude. Then, in the mechanism of Fig. 1(b) used in [14] one has two hadronizations, while in the mechanism that we have of Fig. 1(a) there is only one hadronization. All these reasons justify that we neglect the mechanism of Fig. 1(b) and start the process with the one of Fig. 1(a). The choice of the starting mechanism is not innocuous: different initial meson channels are produced, which upon final state interaction give rise to the three kaons. This means that coupled channels have to be used in the approach, as also done in [14], but the amount f 0 (980) or a 0 (980) production, for instance, depends on the mechanisms assumed and the weights by which the different meson channels appear in the hadronization.
As we mentioned, in the experiment of Ref. [13], two methods of analysis were made, the first one using the isobar model and the other one using the formalism of [14]. The problems using the isobar model which were exposed in [14] are further evidenced in the study done in [13]. Indeed, when using the isobar model three options, A, B, C are used. In model A they include the φK + , f 0 (980)K + and f 0 (1370)K + channels. In model B they add a nonresonant amplitude. In model C they replace the f 0 (1370)K + channel by the a 0 (980)K + one. They note that both models B and C are equally acceptable. This means that this type of analysis cannot tell us about the relevance of the a 0 (980) resonance in this reaction. By the contrary, the analysis done there using the model of [14] shows a dominant role played by the a 0 (980) resonance. Incidentally, we should also mention that by using fully unitarized amplitudes in coupled channels in the approach that we follow, we do not have to worry about background. The amplitudes contain the resonance pole but provide at the same time background terms away from the resonance peaks.
The strategy of our work, which is widely used (see Ref. [30] for a review on this issue) is to establish the dominant mechanisms at the quark level that produce the desired number of mesons after hadronization including qq pairs with the quantum numbers of the vacuum. Then all pairs of mesons are allowed to undergo final state interaction keeping a third particle as spectator. The amplitudes are properly symmetrized to account for the identity of the particles. The final sate interaction requires to use the meson-meson scattering amplitudes, which we take from the prior study with the chiral unitary approach. This method has a minimum input, basically a global strength and some relative strength from the s-wave to p-wave amplitudes. An agreement with data with this minimum input is considered as giving support for the chiral unitary approach to the meson-meson interaction and in particular for the nature of some of the resonances that it provides as dynamically generated from the meson interaction rather than genuine states of particular quark configurations. The good agreement with data that we will show in the present work will then give support to this kind of molecular picture for the f 0 (980) and a 0 (980) resonances, which adds to the support from many other reactions where these states are produced [30].
FIG. 1: Elementary D + → K + P P process at the quark level.
As briefly discussed in the Introduction, the possible elementary quark topologies at tree level for the D + → K − P P process, where P P stands for a pseudoscalar meson pair, are depicted in Fig. 1. Note that P P in Fig. 1(a) cannot be K + K − at the tree level since it cannot be produced fromdd, but the K + K − pair can be produced via final state interaction (FSI) of the pseudoscalar pair. This necessity for FSI was the main reason for neglecting the diagram (a) in Ref. [14], which is the basis of the LHCb experimental analysis in Ref. [13]. However, our position in the present work is opposite and we are going to argue why we expect the (a) diagram to be dominant. First, diagram (b) represents annihilation and, since the D have spin zero, the W annihilation diagram is suppressed by helicity conservation at the light quark vertex [28]. In addition, diagram (b) requires two hadronizations, each of which reduces the width by about one order of magnitude [31]. Furthermore, diagram (a) relies upon external emission, which has the largest strength for weak interaction [28]. Indeed, one quark is operative and the other one remains spectator which implies a one body operator, versus the two body operator required in the annihilation, and is color favored. On the other hand, the FSI interaction necessary to produce final K + K − is actually required at low invariant masses since that region is much influenced by the a 0 (980) and f 0 (980) resonances which are dynamically generated within the chiral unitary approach (UChPT) from the final P P interaction, as explained in the Introduction. We will come back to the implementation of the FSI trough the UChPT amplitudes in the second part of this section, but first we address the calculation of the elementary production at the quark level in Fig. 1(a).
While thed quark remains as spectator, the c quark becomes a d through the emission of a W boson, which eventually creates the us of a K + . Note that this process involves two Cabibbo suppressed weak transitions, (W cd, W us). The finaldd pair then hadronizes into a final pseudoscalar meson pair, which is implemented by producing an extraqq with the 3 P 0 model [32][33][34]. The weight of the different allowed pseudoscalar pairs produced in the hadronization, can be related, up to a global normalization factor, using the following SU (3) arguments: Let |H be the flavor state of the final hadronic part after the quark-antiquark pair is produced in the hadronization: It can be written as where we have defined The strength of SU (3) comes into play when we associate the matrix M to the usual SU (3) matrix containing the pseudoscalar mesons: where we have used ideal mixing between the singlet and octet to give η and η ′ [35]. Then, the matrix element required in Eq. (3) is Note that, as mentioned above, no K + K − pair is possible in the hadronization fromdd and then it must necessarily be produced in the final state interaction from the five possible pseudoscalar pairs, π − π + , π 0 π 0 , ηη, π 0 η and K 0K 0 , as depicted in Fig. 2. Note that we are not including the scalar f 0 (980) and a 0 (980) resonances as explicit degrees of freedom but they arise naturally in the non-linear dynamics involved when implementing unitarity in coupled channels starting from a lowest order tree level meson-meson chiral potential. This effectively accounts for the resummation shown in Fig. 2 and it is the basis of the chiral unitary approach (UChPT). In the scalar sector, there are several different ways to implement these ideas like the Bethe-Salpeter equation [36], the inverse amplitude method [24,37] or the N/D method [23], but all of them provide similar results. Since in the present work we are going to compare with experimental data from the LHCb collaboration which involves K − K + invariant masses up to 1375 MeV, we use the amplitudes from the N/D approach [23] which provides the largest range of predictability among the aforementioned approaches. All theses approaches rely upon mainly one free parameter coming from the regularization, either a cutoff or a subtraction constant, which is determined from a fit to meson-meson scattering data. The N/D method of [23] is able to extend the applicability range up to higher energies by including in the interaction kernel, in addition to the lowest order ChPT amplitudes, the s-channel exchange of scalar resonances in a chiral symmetric invariant way. These tree level resonances constitute an octet with mass around 1.4 GeV and a singlet around 1 GeV which barely changes the dynamical origin of the a 0 (980) and f 0 (980) but improves the amplitude close to 1400 MeV. In any case, all of the UChPT approaches provide similar results in the region around 1 GeV. In Fig. 3 we show some spin 0 and isospin I = 0 and I = 1 meson-meson scattering amplitudes 1 , which will be needed in the present work. The energy range involved in the present work is from the K + K − threshold, 987 MeV, till M D − m K = 1375 MeV. We can clearly see the shapes for the f 0 (980) and a 0 (980), but note that these shapes are far from being just Breit-Wigners, and this is one of the strong points of UChPT: it provides not only the pole structure of the resonances but the actual scattering amplitude. The a 0 (980) actually corresponds to a cusp at the KK threshold.
We can then write the amplitude corresponding to the process in Fig. 2. If we use the label 1 for the K + coming directly from the D + , label 2 for the K − and 3 for the other K + (see Fig. 2), it can be written as where s ij = (p i + p j ) 2 . In Eq. (6) the sum runs over the five P P allowed channels in Eq. (5), C is an arbitrary global normalization factor to be fitted later on to the experimental LHCb data, h i are the numerical coefficients in front of each P P channel in Eq. (5), t i,K + K − stands for the unitarized (P P ) i → K + K − amplitude in s-wave explained above and G i is the loop function for two pseudoscalar mesons regularized with the same subtraction constant used in the evaluation of t i,K + K − . We can theoretically filter the different isospin contributions taking into account that π − π + , π 0 π 0 and ηη contribute only to I = 0, and π 0 η to I = 1. The K 0K 0 pair contributes to both isospins, but taking into account the isospin decomposition of the different KK states it suffices to substitute, in Eq. (6), for I = 0, and for I = 1.
To the I = 0 amplitude in Eq. (6) we must add the contribution from the φ(1020) meson which, being a genuinē qq resonance in p-wave, is no included in the aforementioned meson-meson scattering amplitudes. Since the amplitude is going to span a large invariant mass region, we consider a full relativistic amplitude as in Ref. [12]: where s 13 = M 2 D + 3m 2 K − s 12 − s 23 , D is a complex arbitrary normalization factor to be fitted later on, and we use an energy dependent p-wave φ width where Γ o is the total width of the φ, p(m) the K momentum in the φ rest frame for a φ invariant mass m, and B(m) is the p-wave Blatt-Weisskopf barrier penetration factor [39] given by In Eq. (12), R stands for the range parameter of the φ a typical value of R = 1.5 GeV −1 , although it is not very relevant. Finally, the three body distribution for the D + → K − K + K + decay is given by where M is the total D + → K − K + K + amplitude considered adding Eqs. (6) and (10) In the results section we will also evaluate the K + K + distribution which is given by but using for s 12 in the argument of M(s 23 , s 12 ), s 12 = M 2 D + 3m 2 K − s 13 − s 23 .

III. RESULTS
Our model has three parameters: one for the global normalization C in Eq. (6) and two for the global weight of the φ meson amplitude, complex D in Eq. (10), which we fit to the experimental [13] K + K − distribution, (only panel-(a) in Fig. 4). The other panels in Fig. 4 represent the K + K + distribution and the distributions s high K + K − and s low K + K − where, according to [13], s high K + K − and s low K + K − represent the highest and lowest values among s 12 and s 23 , see Fig. 5. Theoretically, we evaluate the s low K + K − distribution including θ(s 23 − s 12 ) in the integrand of Eq. (15), with θ the step function.
In all the figures the dots represent the LHCb experimental data [13], the solid line our full model, the dotted line the phase-space, the dashed line the I = 0 contribution and the dashed-dotted line the I = 1. The experimental data in Ref. [13] are not corrected for setup acceptance, however phase-space curves weighted by the efficiency for the different plots in Fig. 4 are provided in the experimental paper. Therefore, we have renormalized each experimental datum such that the phase-space agrees with the theoretical three body distribution.
A first observation from Fig. 4 is that our model fits reasonably well the whole spectrum of the K + K − , K + K + distributions (recall that we have only fitted the data of panel a). This is remarkable, given the little freedom in the fit: just the global normalization factor and the relative complex weight of the φ(1020). The rest is given from the non-trivial unitarization model implied in Eq. (6). It is also worth recalling again that there is no K + K − in the elementary production vertex of Fig. 1(a) and thus all the strength is coming from the FSI starting with meson-meson channels other than K + K − .
On the other hand, by looking at the different isospin contributions in Fig. 4, we see that the I = 1 contribution dominates over the I = 0 one. In particular, close to threshold the accumulation of the strength with respect to the phases-space is mainly due to the I = 1 amplitude, i.e. the effect of the a 0 (980). This is more clearly manifest if we look at Fig. 6, where we zoom in the K + K − mass distribution near the threshold, including, in addition to the theoretical curves of Fig. 4(a), the   contribution considering only the φ meson, only the I = 0 without the φ, and the full model removing the φ. This is also more clearly seen if we theoretically remove the phase-space by plotting, in Fig. 7, the different isospin contributions from the FSI terms, i.e. without the φ, to the squared amplitude |M| 2 of Eq. (13) without the 1 ←→ 3 symmetrization, i.e, M ≡ T (s 23 ), as a function of s 23 . Note that if we included the φ meson or the symmetrization, the amplitude would also depend on the s 12 variable. We see that close to threshold, indicated by the vertical dotted line, the strength is essentially dominated by the a 0 contribution. Below threshold the shapes of the a 0 and f 0 are clearly visible in this plot but are not accessible when considering the actual phase-space. In Fig. 6, at low invariant masses we see a pattern for the f 0 (980) and a 0 (980) contributions different than what was found in Ref. [14]. In both approaches a dominance of the a 0 (980) contribution is found. In Ref. [14] a destructive interference between the a 0 and f 0 contributions was reported. We find a different pattern. The addition of the f 0 (980) resonance increases the contribution of the a 0 (980) at low invariant masses but subtracts for invariant masses higher than 1020 MeV. The interference of the I = 0 and I = 1 contributions is made possible in the present work because both contributions appear in the variables s 12 and s 23 .

IV. SUMMARY
We show theoretically that the D + → K − K + K + decay can be understood from the mechanism that accounts for the final state interaction of an initial pseudoscalar pair (see Fig. 1(a)). Indeed K − K + K + in the final state is not possible at the tree level and hence the rescattering is mandatory. At this point we take advantage of the unitary extensions of chiral perturbation theory which generate dynamically the f 0 (980) and a 0 (980) resonances, without the need to include them as explicit degrees of freedom, and provide the full meson-meson scattering amplitudes, not only the resonances. The relative weights of the initial production of the meson-meson pairs are obtained from SU(3) arguments, and then for the unitarization we use the UChPT amplitudes from the N/D method which allows to extend the range of applicability to the whole final K + K − mass spectrum in this decay. With a minimal freedom, just the global normalization and the weight and phase of the φ(1020) contribution, we are able to fit experimental data on K + K − invariant mass distribution. We also show that the dominant contribution close to threshold comes from the I = 1 (hence the a 0 (980)) contribution in clear dominance over the I = 0 one (the f 0 (980)). The remarkable agreement is a step in favor of considering this mechanism as the leading one in this decay, at odds with other considerations as in the experimental analysis [13] based on the work in Ref. [14] which advocates for the dominance of the mechanism in Fig. 1(b), and to reinforce the dynamical origin of the a 0 (980) and f 0 (980).

V. ACKNOWLEDGMENTS
One of us, E. O., wishes to acknowledge useful discussions with Alberto Correa dos Reis and his encourage-ment for us to study this reaction. This work is partly supported by the National Natural Science Foundation of China under Grant Nos. 11505158, 11847217, 11975083 and 11947413. It is also supported by the Academic Improvement Project of Zhengzhou University. This work is partly supported by the Spanish Ministerio de Economia y Competitividad and European FEDER funds under Contracts No. FIS2017-84038-C2-1-P B and No. FIS2017-84038-C2-2-P B. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824093 for the STRONG-2020 project.