Dalitz plot studies of $D^0 \to K^0_S K^+ K^-$ decays in a factorization approach

The $\textit{BABAR}$ Collaboration data of the $D^0 \to K^0_S K^+ K^-$ process are analyzed within a quasi two-body factorization framework. In earlier studies, assuming $D^0$ transitions to two kaons and the transitions between one kaon and two kaons to proceed through the dominant intermediate resonances, we approximated them as being proportional to the kaon form factors. To obtain good fits, one has to multiply the scalar-kaon form factors, derived from unitary relativistic coupled-channel models or in a dispersion relation approach, by phenomenological energy-dependent functions. The final state kaon-kaon interactions in the $S$-, $P$- and $D$- waves are taken into account. All $S$-wave channels are treated in a unitary way. The $K^+K^-$ and $\bar K^0 K^+$ $S$-wave effective mass squared distributions, corrected for phase space, are shown, in a model-independent manner, to be significantly different. Then the $f_0(980)$ resonance must be included at variance with the BABAR analysis. The best fit has 19 free parameters and indicates i) the dominance of annihilation amplitudes, ii) a large dominance of the $f_0(980)$ meson in the near threshold $K^+K^-$ invariant mass distribution and iii) a sizable branching fraction to the $K_S^0 \ [\rho(770)^+ + \rho(1450)^+ + \rho(1700)^+] $ final states. An appendix provides an update of the determination of the isoscalar-scalar meson-meson amplitudes based on an enlarged set of data. A second appendix proposes two alternative fits based on scalar-kaon form factors calculated from the Muskhelishvili-Omn\`es dispersion relation approach. These fits have $\chi^2$ quite close to that of the best fit but they show important contributions from both the $f_0$ and $a_0^0$ mesons and a weaker role of the $\rho^+$ mesons.


I. INTRODUCTION
Measurements of the D 0 -D 0 mixing parameters, through Dalitz-plot time dependent amplitude analyses of the weak process D 0 → K 0 S K + K − , have been performed by the Belle [1] and BABAR collaborations [2]. Such studies could help in the understanding of the origin of mixing and may indicate the presence of new physics contribution. As predicted by the standard model in the charm sector, the violation of the CP symmetry should be small for these D 0 decays. In Refs. [1,2] the description of the D 0 → K 0 S K + K − decay amplitude has been performed using the isobar model developed in [3], extended in [4] and [2,5]. The isobar model has also been applied in the experimental analysis based on the data taken from the BESIII experiment [6,7].
The Cabibbo-Kobayashi-Maskawa (CKM) angle γ (or φ 3 ) has been evaluated from the analyses of the B ± → D 0 K ± , with D 0 → K 0 S π + π − and D 0 → K 0 S K + K − decays [8][9][10][11]. This angle can be also measured using some knowledge on the strong-phase difference between D 0 andD 0 → K 0 S K + K − decay amplitudes obtained by the CLEO Collaboration [12]. This method has been used by the Belle [13], LHCb [14] and BESIII [15] collaborations. and of the D 0 decay mechanism into K 0 S K + K − . The experimental analyses like that of Ref. [2] rely mainly upon the use of the isobar model. For a given reaction, this model has basically two fitted parameters for each part of the decay amplitude. In this approach one can take into account many existing resonances coupled to the interacting pairs of mesons. In Refs. [2,5], the authors introduce explicitly eight resonances a 0 (980) 0 , a 0 (980) + , a 0 (980) − , φ(1020) , f 2 (1270) , f 0 (1370) , a 0 (1450) 0 , a 0 (1450) + . Their analysis rely on 17 free parameters.
However, the decay amplitudes are not unitary and unitarity is not preserved in the three-body decay channels; it is also violated in the two-body subchannels. Furthermore, it is difficult to differentiate the S-wave amplitudes from the nonresonant background terms. Their interferences are noteworthy and then some two-body branching fractions, extracted from the data, could be unreliable. One of the difficulty in the experimental analyses based on the isobar model is the choice of the resonances needed to reach a good agreement with the Dalitz plot data. In Ref. [4] the BABAR collaboration authors have added the scalar a 0 (1450) to their model developed in 2005 [3]. In the recent BESIII analysis [7] the Dalitz plot is described with six resonances: a 0 (980) 0 , a 0 (980) + , φ(1020) , a 2 (1320) + , a 2 (1320) − , a 0 (1450) − .
Extending our previous work on the D 0 → K 0 S π + π − decays [16], we analyze, in the quasi-two-body factorization framework, the D 0 → K 0 S K + K − data provided by the BABAR collaboration [2]. As in our earlier studies, we assume that two of the three final-state mesons form a single state which originates from a quark-antiquark, qq, pair and then apply the factorization procedure to these quasitwo-body final states. Starting from the weak effective Hamiltonian, we derive tree and annihilation (W -meson exchange) amplitudes both being either Cabibbo favored (CF) with c → sdu transition or doubly Cabibbo suppressed (DCS) with c → dus transition.
In the factorization approach, the CF and DCS amplitudes are expressed as superpositions of appropriate effective coefficients and two products of two transition matrix elements. The kaon form factors do not appear explicitly except the isovector ones that enter in only one term of the CF tree amplitude 1 . In all other terms of our amplitudes, one has to evaluate either, for the tree ones, the matrix elements of D 0 transitions to two-kaon states or, for the annihilation ones, the transitions between one kaon and two kaon-states. Similarly to previous studies [17], assuming these transitions to proceed through the dominant intermediate resonances, we have approximated them as being proportional to the isoscalar or isovector kaon form factors.
Taking advantage of the coupling between the ππ and the KK channels and extending the derivation of the unitary isoscalar-scalar pion form factor [18] to that of the kaon one, two of the present authors (L. L. and R. K.) together with two collaborators, have recently studied, in the quasi-twobody QCD factorization approach, the B ± → K + K − K ± decays [19,20]. These S-wave form factors are derived using a unitary relativistic three coupled-channel model including ππ, KK and effective (2π)(2π) interactions together with chiral symmetry constraints. They include the contributions of the f 0 (980) and f 0 (1400) resonances and require the knowledge of the isoscalar-scalar meson-meson amplitudes from the two kaon threshold to energies above the D 0 mass. The parameters of these amplitudes derived in the late nineties by three of us (R. K., L. L. and B. L.) [21,22] have been updated using new precise low energy ππ data together with an enlarged set of data as is shown in The calculation requires also the knowledge of a contribution proportional to the isovector-scalar form factor; it is represented by a function calculated from a unitary model with relativistic twocoupled channel πη and KK equations based on the two-channel model of the a 0 (980) and a 0 (1450) resonances built in [23,24].
their masses, widths and partial decay widths into KK [26]. Consequently, as in Ref. [2], we consider only the f 2 (1270) to represent the D wave. It is an "effective" f 2 (1270) which takes into account both tensor mesons.
In the present approach, a best fit is achieved in which the data are reproduced with amplitudes that are optimized notably by adjusting the isoscalar-and isovector-scalar form factors. It is interesting then to see what could be the outcome of a model based, for instance, on the scalar form factors calculated from the Muskhelishvili-Omnès dispersion relation approach [27][28][29][30]. As in our best fit model (χ 2 /ndf = 1.25), we have to introduce energy dependent phenomenological functions multiplying the scalar form factors to obtain two fits with almost as good χ 2 /ndf . Branching fractions of these two alternative models are compared to those of our best fit model.
The paper is organized as follows. A full derivation of the D 0 → K 0 S K + K − decay amplitude, in the framework of the quasi-two-body factorization approach is given in Sec. II. Based only on the experimental data of the BABAR Collaboration [2] and without any model assumptions, we show in Sec. III that the f 0 (980) contribution, at variance with the BABAR analysis, has to be included in the decay amplitude. Section IV presents the result of our best fit to the Dalitz plot data sample of Ref. [2]. Some discussion and concluding remarks can be found in Sec. V. Appendix A presents the update of the description of the ππ,KK and effective (2π)(2π) S-wave isospin zero scattering amplitudes. Two alternative models for the D 0 → K 0 S K + K − decay amplitude, using kaon scalar-form factors derived from the dispersion relation approach, are presented in Appendix B.
The decay amplitudes for the D 0 → K 0 S K + K − process can be evaluated as matrix elements of the effective weak Hamiltonian [31] where the coefficients V CKM are given in terms of Cabibbo-Kobayashi-Maskawa quark-mixing matrix elements and G F denotes the Fermi coupling constant. The C i (µ) are the Wilson coefficients of the four-quark operators O i (µ) at a renormalization scale µ, chosen to be equal to the c-quark mass m c .
The transition matrix elements that occur in the present work require two specific values of the V CKM coupling matrix elements [16]: where λ is the sine of the Cabibbo angle and where Λ 1 is associated to Cabibbo favored (CF) transitions while Λ 2 is associated to doubly Cabibbo suppressed (DCS) amplitudes. The amplitudes are functions of the Mandelstam invariants where p 0 , p + and p − are the four-momenta of the K 0 S , K + and K − mesons, respectively. Energymomentum conservation implies where p D 0 is the D 0 four-momentum and m D 0 = 1864.83 MeV, m K 0 = 497.611 MeV and m K = 493.677 MeV denote the D 0 , K 0 and charged kaon masses (Ref. [26]).

A. Tree and annihilation CF and DCS amplitudes
The full amplitude is the superposition of two tree CF and DCS amplitudes, T CF (s 0 , s − , s + ) and T DCS (s 0 , s − , s + ) and of two annihilation (i.e., exchange of W meson between the c and u quarks of the D 0 ) CF and DCS amplitudes, A CF (s 0 , s − , s + ) and A DCS (s 0 , s − , s + ). Thus, one writes the full amplitude as Although the three variables s 0 , s − , s + appear as arguments of the amplitudes, all amplitudes depend only on two of them because of the relation (4). Assuming that the factorization approach [31][32][33][34] holds, the diagram illustrated in Fig. 1 is pro- L final states and the diagram in Fig. 2, proportional to Λ 1 a 2 (m c ) with quasi-two-body K 0 [K + K − ] I L with angular momentum L = S, P, D and isospin I = 0, 1 states, yield the tree CF amplitude T CF (s 0 , s − , s + ) which reads, with |0 denoting the vacuum state, The short hand notation of the last line of Eq. (6) implies the L summation 2 . It will be used wherever necessary. In the case of the creation of a K + K − pair we indicate by a subscript q the qq pair from which it originates (here a uu one, as seen in Fig. 2). We shall therefore use the notation [K + K − ] I L,q and/or [K + K − ] I q whenever necessary. We have also introduced the short-hand notation which will be used throughout the text (γ and γ 5 are Dirac's matrices). In deriving Eq. (6) small CP violation effects in K 0 S decays are neglected and we use At leading order in the strong coupling constant α S , the effective QCD factorization coefficients a 1 (m c ) and a 2 (m c ) are expressed as where N C = 3 is the number of colors. Higher order vertex and hard scattering corrections are not discussed in the present work and we introduce effective values for these coefficients. From now on, the simplified notation a 1 ≡ a 1 (m c ) and a 2 ≡ a 2 (m c ) will be used. As in Ref. [16], we take a 1 = 1.1 and a 2 = −0.5.
2 e.g., The corresponding DCS annihilation amplitudes, A DCS (s 0 , s − , s + ), (see Fig. 6), are easily obtained 3 In the amplitudes (11) and (12) we neglect the terms with the quasi two-boby K + (p+) [K 0 L final states. One expects a small contribution of these terms because there exist no strangeness -2 ( from the CF amplitudes in Eq. (11) Let us now review in detail all the amplitudes that will have to be evaluated. We will follow closely the construction detailed in Ref. [16].

B. Explicit tree amplitudes
In the following, starting from the expressions given in Eqs. (6) for the CF amplitudes and in (10) for the DCS ones, we will express the different three-body matrix elements entering in the amplitudes in terms of vertex functions noted G R S,P, Further on, we will need to introduce transition form factors for which, as in Ref. [16], we will assume isospin and charge conjugation symmetry so that the following equations arise: (s), In the above equations F 0 and F 1 denote scalar and vector transition form factors of two pseudoscalar mesons while A 0 's are transition form factors of pseudoscalar and vector mesons.

Scalar amplitudes
Following a derivation similar to that developped in Ref. [16], the isoscalar-scalar CF amplitude associated to the K 0 [K + K − ] I S,u final states can be described by (see Fig. 2) where f K 0 is the K 0 decay constant and the sum over R S runs over the possibly contributing resonances in the isoscalar-scalar channel. It can be seen here that we have approximated the three-body matrix here we can assume that its variation from one resonance to the other is small and we can choose for its value that of the transition to f 0 (980). Unless otherwise specified, by f 0 in F D 0 f 0 0 (m 2 K 0 ) we mean f 0 (980). We may parametrize the sum over R S by introducing the isoscalar-scalar form factor, Γ n * 2 (s 0 ), where n denotes a nonstrange quark pair and which can be built following the method discussed in Ref. [18]. We then apply the following approximation where χ n is a constant complex factor. Hence The real transition form factor, F D 0 f 0 0 (m 2 K 0 ), can be obtained from Ref. [35]. This amplitude has to be associated with the corresponding isoscalar-scalar S,u DCS amplitude (see Fig. 4) approximated by Recombining the two amplitudes (16) and (17), we have We now turn to the three isovector-scalar tree amplitudes. The isovector-scalar amplitude, associated to the a 0 (980) − and a 0 (1450) − resonances can be written as (see Fig. 3) the a − 0 resonances being built from ud pairs and R S [K 0 K − ] 1 |du = 1. In Eq. (19) f K + denotes the charged kaon decay constant. Parametrizing, as above, the sum over R S as we get The function G 1 (s) can be built using the isospin 1 coupled KK and πη channel description of the a 0 (980) and a 0 (1450) performed in Ref. [23]. As previously for the isoscalar-scalar case, we have assumed here that the variation of the D 0 → R S transition form factor from one resonance to the other is small and we choose it to be that of the lowest resonance, i.e., R S ≡ a 0 (980) − , denoted simply In the case ot the isovector-scalar [K 0 K + ] 1 S K − CF amplitude (a 1 term of Eq. (6), related to the a 0 (980) + and a 0 (1450) + resonances, one has 4 T CF where F is the kaon isovector-scalar form factor and denoted also as F second relation of the Eqs. (13). For the transition form factor F D 0 K − 0 (s + ), following Ref. [36], we use the parameterization: where M V = 2.11 GeV. We have then We proceed similarly for the isovector-scalar CF S,u amplitudes (see Figs. 2 and 4). It is given by where the sum over R S runs over the contributing resonances in that channel, i.e., a 0 (980) 0 and where we assume (isospin invariance) that, to describe the uū transition to the isovector-scalar K + K − state, it is the same function G 1 (s) as that introduced in Eqs. (20) for the dū transitions to the isovector-scalar K 0 K − state. In Eq. (26) a 0 0 means a 0 (980) 0 . We may now rewrite In a similar way, the related isovector-scalar DCS amplitude reads Recombining Eqs. (27) and (28) we get

Vector amplitudes
Let us now study the vector tree amplitudes associated to the K 0 [K + K − ] I P,u channel. The isoscalar-vector CF amplitude can be built from the ω and φ resonances (see Fig. 2). It reads where m R P [K + K − ] denotes the mass of the contributing resonances. Now we introduce the following parametrization in terms of the kaon vector form factor F (s 0 ) and, for the same reasons as introduced in the scalar case [Eqs. (14) and (21)] This approximation relies on the fact that the mixing angle of the vector meson nonet is very close to the ideal mixing angle, θ V = 35.3 • , so that the φ resonance amplitude gives an almost nul contribution.
Note that f ω denotes the decay constant for the ω(782) meson. We have then The associated isoscalar-vector K 0 [K + K − ] 0 P,u DCS amplitude is given by Thus, from Eqs. (32) and (33), the isoscalar-vector amplitude reads The isovector-vector T CF amplitude related to the ρ 0 resonances is given by a similar expression to Eq. (30) (see Fig. 2) where Again, parametrizing the sum over the vertex functions by where f ρ is the charged ρ(770) decay constant, Then comes the contribution of the isovector-vector Fig. 4). It goes as so that the total isovector-vector amplitude is The isovector-vector [K 0 K + ] 1 P K − CF amplitude has the expression (see Fig. 1) The sum over the vertex functions As in Ref. [36] we parametrize with, as before in Eq. (23), M V = 2.11 GeV.

Tensor amplitudes
For the isoscalar-tensor amplitude K 0 [K + K − ] 0 D,u amplitude, one can write (see Fig. 2) but it will be dominated by the f 2 (1270) resonance with mass m f 2 ; it will be described by a Breit-Wigner representation. Linking the vertex function to the form factor and using we have where The three-momenta p 2 , p 0 are defined in the [K + K − ] center-of-mass (c.m.) system. One has and The scalar product p 2 · p 0 which enters the function D(p 2 , p 0 ) is given by the relation One has In Eq. (49), because of the large width of the f 2 meson, an energy dependent total width Γ f 2 (s 0 ) has been introduced (see Eqs. (121) and (122) in Ref. [16]) such that where r = 4.0 GeV −1 and q ≡ |p 2 |. In Eq. (49) the D 0 → f 2 effective transition form factor F D 0 f 2 (m 2 K 0 ) will be treated as a free complex parameter. To this amplitude one has to add the isoscalar-tensor K 0 [K + K − ] 0 D DCS amplitude given by (see Fig. 4) The total isoscalar-tensor amplitude then reads . (57) C. Annihilation amplitudes

Scalar amplitudes
The isoscalar-scalar CF annihilation amplitude corresponding to while the isoscalar-scalar DCS amplitude corresponding to where we have used the equality Thus the total isoscalar-scalar annihilation amplitude reads Following the steps in Sec.II B it leads to where χ s is a complex constant and Γ s * 2 (s 0 ) is the kaon strange scalar form factor.
The isovector-scalar annihilation DCS amplitude, associated to the K + [K 0 K − ] 1 S final states containing the a 0 (980) − and a 0 (1450) − , is given by and hence with, reads The corresponding isovector-scalar annihilation CF amplitude associated to the and contains the a 0 (980) + and a 0 (1450) + . With the approximation we reach The [K + K − ] 1 S final states which would contain the a 0 (980) 0 and a 0 (1450) 0 mesons cannot be formed from a ss pair and thus the corresponding S isovector-scalar DCS amplitude is zero.

Vector amplitudes
We now turn to the vector-annihilation amplitudes. The isoscalar-vector CF amplitude correspond- and is associated to the φ mesons. It may be reexpressed as One has to add the associated DCS amplitude corresponding to K 0 [K + K − ] 0 P,s final states (see Fig. 6 The isovector amplitude corresponding to K − [K 0 K + ] 1 P final states, which contains the ρ(770) + , ρ(1450) + and ρ(1700) + mesons, may be written as Similarly, the isovector-DCS amplitude corresponding to It contains the ρ(770) − , ρ(1450) − and ρ(1700) − mesons.

Tensor amplitudes
Finally we present the tensor amplitudes. The two isoscalar CF and DCS amplitudes associated and They contain the f 2 (1270) meson. In the last equation we have used the relation Hence the total isoscalar-tensor amplitude reads
Since the isovector-scalar form factor F the isovector amplitude associated to the a + 0 resonances in the channel [K 0 (24) and (68)] can be expressed as The amplitude associated to the ρ + resonances [Eqs. (42) and (72)] reads The form factor F The isovector amplitude associated to the a − 0 resonances in the channel [K 0 The isovector amplitude associated to the ρ − resonances is given by Eqs. (45) and (74) where we have applied the relation F Finally, the isoscalar-tensor amplitudes related to the f 2 [Eqs. (57) and (78)] can be recombined to give where can be treated as a complex constant parameter fitted to data.
Our study can provide information on the S-wave content of the KK effective mass densities. In the region of low effective masses, near the KK thresholds, one expects dominant contributions of the S-and P -wave amplitudes which simplifies the partial wave analysis of the experimental Dalitz plot distribution. This analysis has been performed by the BABAR Collaboration for the following three decay reactions: [38]. In the K + K − S-wave effective mass distributions both scalar resonances f 0 (980) and a 0 (980) contribute while in the K 0 K + case only the a 0 (980) + resonance is present.
In the analyses of Refs. [2] and [3] the f 0 (980) contribution has not been introduced. A possible argument in favor of this choice has been formulated in Ref. [3], namely the authors have expected that the presence of the f 0 (980) resonance would lead to an excess in the K + K − mass spectrum with respect to K 0 K + . However, based on the limited statistics of 12540 events they have observed that both KK spectra are approximately equal. Below, using a much larger sample of about 80000 signal events [2], we show that the K + K − and KK spectra in D 0 decays into K 0 S K + K − are significantly different for low KK effective masses. Thus, the contribution of the f 0 (980) resonance is required to obtain a good description of the data of Ref. [2].
To proceed further, the definitions of the KK effective S-wave mass distributions corrected for phase space are needed. If we denote by N (s 0 , s + ) the number of events of the D 0 → K 0 K + K − reaction, the corresponding Dalitz-plot density distribution is given by d 2 N (s 0 , s + )/ds 0 ds + . The K 0 S K + effective mass squared distribution corrected for phase space can be then defined as where s 0max and s 0min are the maximum and minimum s 0 values at fixed s + . If we limit ourselves to the low s + values (for example, up to about 1.05 GeV 2 ) then to a good accuracy the above distribution corresponds to the S-wave part of the total decay amplitude related to the isovector-scalara 0 (980) + resonance. The reason is that the dominant P -wave contribution, related to the φ(1020) resonance, is only present in the K + K − decay channel. Moreover, the low-mass K 0 S K + and K + K − distributions are very well separated on the Dalitz plot [2].
For the low effective K + K − masses one has to subtract the P -wave contribution. Following Ref. [3], this can be done by calculating the spherical harmonic moments and with and where θ is the helicity angle of the K 0 S meson defined with respect to the K + direction in the K + K − center-of-mass frame, s +max and s +min being the maximum and minimum s + values at fixed s 0 . The S-wave K + K − effective mass squared distribution corrected for phase space is then defined as For completeness we give below the kinematical relation for the cosine of the helicity angle with |p + | and |p 0 | defined by Eqs. (51) and (52), respectively.
We have performed a simplified partial wave analysis of the BABAR data published in Ref. [2]. As described above, only the S-and P -waves have been included and the effective K 0 S K + and K + K − masses were smaller than 1.05 GeV 2 . The number of signal events of the D 0 →K 0 K + K − decays was 79900±300. Based on the Dalitz plot density distributions corrected for reconstruction efficiency and background, the S-wave K 0 S K + and K + K − effective mass distribution corrected for phase space are calculated using Eqs. (92) and (96).
The comparison of the calculated S-wave K + K − and K 0 K + distributions is shown in Fig. 7. In the left panel a) a clear surplus of theK 0 K + distribution over the K + K − one is seen below 1.02 GeV.
Above m KK = 1.02 GeV the open circles corresponding to K 0 K + spectrum are in majority located below the closed circles (K + K − events), so we observe a crossing of the two distributions. This effect is statistically significant. It was not so clear in 2005 when the first set of the BABAR data was published. But even then, in Fig. 8 of Ref. [3], one can see the same cross-over tendency as in Fig. 7 although the statistics was lower by a factor larger than 6. In the right panel b), one sees that In conclusion, the shape of the K + K − and K 0 K + S-wave effective mass squared distributions, corrected for phase space, is significantly different, so in the phenomenological analysis of the D 0 → K 0 K + K − data one cannot neglect the f 0 (980) contribution in the decay amplitude.

IV. RESULTS AND DISCUSSION
The differential branching fraction or the Dalitz plot density distribution is defined as where M = 7 i=1 M i is the decay amplitude for the process studied and Γ D 0 is the D 0 width. The decay amplitudes M i have been derived in Sec. II. In Table I one can find some constant parameters which appear in these amplitudes. Parameter Value Reference 1561 [16] f ρ 0.209 [34] f φ 0.22 [39] f D 0 0.2067 [16] To make a comparison of experimental data with model predictions the Dalitz diagram has been divided into five regions as shown in Fig. 8. The fit of the model parameters to the experimental data has been performed using the χ 2 tot function defined as a sum of two components: The value of χ 2 i for each cell i has been defined as in Ref. [40]: where N exp i is a number of experimental signal events in the cell i corrected for the reconstruction efficiency and N th i is the theoretical number of events in the same cell 6 . Including the above corrections one gets the total number of experimental events equal to N exp = 80379. The total number of theoretical events is then taken equal to N exp .
The second component of the χ 2 function is given by In our fit the experimental branching ratio for the decay D 0 → K 0 S K + K − has been taken equal to Br exp = 4.45 × 10 −3 and its error ∆Br exp = 0.34 × 10 −3 . These values agree well with recent values of the Particle Data Group [26]. The theoretical branching fraction Br th is obtained from Eq. (98) after integrations of d 2 Br ds + ds 0 over the variables s + and s 0 . The weight factor w in our fit has been set to 100 in order to obtain a good agreement of the theoretical branching fraction with its experimental value. We have performed many fits with our model using different parameter sets. The best fit χ 2 = 1474.4 has been obtained with the nineteen free parameters which are displayed in Table II.
Since the number of degrees of freedom is ndf = 1196 − 19 = 1177, the χ 2 per degree of freedom is The value of the constant |χ n | has been estimated using a relation derived similarly as Eq. (18) in Ref. [41] in which the coupling constants of the f 0 (980) resonance to the K + K − pair are taken into account instead of the f 0 (980) coupling to the ππ system. However, in the present case one has to include two close f 0 (980) poles sitting on the sheets (-+-) and (-++) (for their complex energy positions, E R 1 and E R 2 , see Table IX in Appendix A). One can generalize Eq.(18) from Ref. [41], valid for the pole position of a single resonance, to the case of two close resonances: where g R 1 K + K − and g   The form factors Γ n 2 (s 0 ) and Γ s 2 (s 0 ) have been calculated in a three-channel model of meson-meson interactions (ππ, KK and an effective 2π 2π), introduced in Ref. [18]. These form factors depend not only on the values of the meson-meson parameters listed in Table VIII in Appendix A but also on two other parameters κ and c defined by Eqs. (28) and (39) in Ref. [18], respectively. Their values κ = 2807.3 MeV and c = 0.109 GeV −4 have been fitted to the B ± → K ± K + K − decay data analyzed in Ref. [20].
In Fig. 9 we show the effective energy dependence E of moduli and phases of the KK isoscalarscalar nonstrange Γ n 2 (E) and strange Γ s 2 (E) form factors. The energy E is equal to the square root of s.
In the above model the kaon threshold energy was set equal to the sum of the charged and neutral kaon masses. However, the M 1 amplitude corresponds to the isoscalar K + K − S-wave state with a threshold lower by about 3.9 MeV in comparison with the K + K 0 S threshold energy. In order to take this effect into account in an approximate way, we introduce the variablē where s is a new parameter which is fitted to the data (see Table II). It corresponds to a zero of P (s 0 ). The third order polynomial in the denominator, with the constant b fixed to 0.0154 GeV −6 , is introduced in order to control asymptotically the high energy behavior of the M 1 amplitude. This denominator replaces the denominator (1 + cE 4 ) with E = √ s 0 in Eq. (39) of [18]. A plot of the function (103) used in our fit is shown as the continuous black line denoted R BF in Fig. 16 (a) where it is compared to the corresponding functions used in the two alternative fits MO P 1(P 2) described in Appendix B. This function reduces the moduli of the amplitudes which depend on the isoscalar-scalar form factors.
The masses and widths of the isovector-scalar resonances a 0 (980) and a 0 (1450) are presented in Table III. They have been fixed during the minimization of the χ 2 function. The parameters of the a 0 (1450) on sheet (− −) were taken from Ref. [26]. However, we have studied the influence of the position of the a 0 (980) pole on sheet (− +) in the complex energy domain on the χ 2 minimum curve.
In this way the a 0 (980) mass and width on sheet (− +) have been determined together with an estimation of their errors. The masses and widths of other two associated a 0 poles are also given in Table III.
The coupled channel model of the a 0 (980) and a 0 (1450) resonances described in Ref. [23] has been implemented. There, the separable πη and KK interactions have been used in the calculation of the S-wave isospin one scattering amplitudes. Altogether the model has five parameters: two  range parameters β 1 and β 2 , two channel coupling constants λ 1 and λ 2 , and the interchannel coupling constant λ 12 (here the channel πη is labeled by 1 and the channel KK by 2). The potential parameters are given in Table IV. There exist direct numerical relations between the four parameters describing the positions of the a 0 (980) and a 0 (1450) resonances in the complex energy plane (Table III) and the four potential parameters β 1 , λ 1 , λ 2 and λ 12 at fixed value of the fifth parameter β 2 . These relations are given in Ref. [42].
The function G 1 (s) is introduced to describe a transition from the uū pair to the KK spin zero isospin one state. Two isovector-scalar resonances a 0 (980) and a 0 (1450) can be formed during that transition. Both resonances are also coupled to the πη state. Therefore it is natural to consider three cases for the transition from the uū pair to the KK state. In the first case the KK pair is directly formed from the uū pair. In the second case the KK pair undergoes the elastic rescattering in the final state. In the third case the intermediate πη pair is formed and then the inelastic transition to the KK state takes place. The interaction between the meson-meson pairs is treated in the framework of the separable potential model fully described in Ref. [42] and used to study the properties of the a 0 resonances (Refs. [23], [24]).
Below we briefly derive the dependence of the G 1 (s) function on the meson-meson transition amplitudes. Labelling by 1 the πη channel and by 2 the KK channel, one can express G 1 (s) as a superposition of three terms: where Here r 1 and r 2 are the coupling constants corresponding to the uū transitions to the πη and KK states, respectively. The function W (s) is the third-degree polynomial where p 1 , p 2 and p 3 are the real parameters included in the list of the model free parameters (see Table II). We keep the same p j , j = 1, 2, 3, parameters for both channels. The fitted polynomial W (s 0 ) is plotted as the continuous black line in Fig. 16(b) where it is compared to the polynomial P F (s 0 ) defined by Eq. (B1) and used in the alternative MO P 1(P 2) fits discussed in Appendix B. The introduction of these polynomials improves the quality of the χ 2 fit, in particular in the region II where the density of events is small. The functions g i (k i ), i = 1, 2, are the vertex functions where m i are the channel reduced masses, k i are the channel momenta and β i are the range parameters.
In the πη channel m 1 = m π m η /(m π + m η ), in the KK channel m 2 = m K /2. We take the neutral and in Eq. (107) we have where the energies are defined as E = √ s, E K (p) = p 2 + m 2 K , E π (p) = p 2 + m 2 π and E η (p) =  Table II. The mass and the width of the φ(1020) resonance seen in Table II are in agreement with the corresponding BABAR values of (1019.55 ± 0.02) MeV and (4.60 ± 0.04) MeV, respectively [2]. The obtained width is higher, by about 0.5 MeV, than the averaged value of (4.249 ± 0.013) MeV given by the Particle Data Group in Ref. [26]. This can be explained by a finite experimental energy resolution.
The branching fraction distributions d 2 Br i ds + ds 0 corresponding to the amplitudes M i , i = 1, ..., 7 are obtained if in Eq. (98) the amplitude M is replaced by M i . One can also define the off-diagonal elements d 2 Br ij ds + ds 0 , i = j, If we integrate over s + and s 0 the differential branching fractions d 2 Br ds + ds 0 , d 2 Br i ds + ds 0 and d 2 Br ij ds + ds 0 then we get the corresponding branching fractions Br, Br i , i = 1 to 7 or the off-diagonal elements Br ij , where i = j. The matrix Br ij is symmetric: Br ij = Br ji .
In Table V Table II and some correlations between the parameters in the amplitudes M 1 and M 3 . Then the branching fraction uncertainties have been obtained from the distributions of the 10 000 values of each branching fraction and of their sum. V: Branching fractions (Br) for different quasi-two-body channels in the best fit to the BABAR data [2].

Amplitude channel
Br i (%) Let us notice the particularly large uncertainties of the branching fraction Br 1 = 60.9 +24.4 −10.6 %. This is due to the fact that the amplitude M 1 consists of three components and contains 9 free parameters.
As seen in Table V Table III the mass of the resonance a 0 (980) equal to 979 +3 −2 on sheet −+ is lower than the K 0 S K + threshold mass of about 991.3 MeV. However, due its finite width of 25 +8 −6 MeV, this resonance, together with the second a 0 (980) resonance on sheet −− at (959 − i 34) MeV, can strongly influence the near threshold s + range of the Dalitz plot density distribution. On the other hand, the mass of the a 0 (1450) resonance lies above the upper range of the K 0 S K + effective mass close to 1371 MeV. However, the a 0 (1450) resonances are wide and they can also affect the distribution of the D 0 → K + K − K 0 S events on the Dalitz plot. The contribution of the quasi-two body channel [K 0 S K + ] P K − is related to nonzero couplings of the P -wave resonances ρ(770) + , ρ(1450) + and ρ(1700) + to K 0 S K + . Although the ρ(770) + mass lies below the K 0 S K + threshold its width is sufficiently large to influence the differential density distribution of the Dalitz plot for s + values above the threshold. The ρ(1450) + width is even larger than that of ρ(770) + , so the whole s + range on the Dalitz plot is sensitive to the strength of its coupling to K 0 S K + . The above three ρ resonances, being wide, cannot create a clear structure or a well distinguished band on the Dalitz plot. This could be a reason why they have not been included in the isobar model analyses [2] and [3].
These results can be compared to the results of the experimental analysis which finds a summed branching fraction of 163.4 %, mainly with 71.1 % coming from the a 0 (980) 0 and a 0 (1450) 0 , 44.1 % from the φ(1020) resonance and 45.1 % from the a 0 (980) + and a 0 (1450) + .  (Eq. (112)) for the best fit to the BABAR data [2]. All numbers are in per cent. In Table VI the diagonal branching fraction terms already shown in Table V 82)]. The first isoscalar term is proportional to the conjugated kaon nonstrange scalar form factor Γ * n 2 (s 0 ) and the second one to the conjugated kaon strange scalar form factor Γ * s 2 (s 0 ). The third term is proportional to the function G 1 (s 0 ) describing the transition from the uū pair of quarks into the K + K − pair of mesons in the isospin one and spin zero state. In Table VII the diagonal and the off-diagonal components of the branching fraction related to the M 1 amplitude are given. They are defined in a similar way as the Br ij components in Eq. (112).
From this Table we see that the major contribution close to 60% is related to the  branching fraction is equal only to 1.19 %). Also the branching fraction equal to 4.48%, corresponding to the isovector-scalar amplitude M I=1 1 , is much smaller than that related to the M s,I=0 1 amplitude. The sum of all the off-diagonal components equals to -4.57 %.
A comparison of the results for this best fit model with those for the alternative MO P 1(P 2) ones can be found in the Appendix B. Dalitz plot projections or one-dimensional effective mass squared distributions of events are calculated by a proper integration of the two-dimensional density distributions. They are shown in Fig. 11.
The errors of the experimental signal weighted event number distributions are the statistical ones.
The histograms correspond to the theoretical distributions normalized to the same total number of events.
The distribution in Fig. 11(a) is strongly dominated by the maximum corresponding to the φ(1020) resonance decaying to the K + K − pair. This decay is in the P -wave and leads to a characteristic two-maximum shape of the Dalitz plot distribution as a function of s + -the square of the K + K − effective mass. Since the branching fraction for the channel [K + K − ] P K 0 S is large (45.5 %) the two Dalitz projections in Figs. 11(b) and(c) have a two-maximum character. There are also two other important contributions related to the amplitudes M 1 and M 3 . However, they do not produce any pronounced structures on the Dalitz plot since both are due to the S-wave in the K + K − or in the Finally let us discuss the low effective mass parts of the K + K − and K 0 S K + distributions (m KK < 1.06 GeV). Since the differential branching fraction [Eq. (98)] is proportional to the Dalitz plot density distribution of events d 2 N ds + ds 0 , one can calculate the theoretical one-dimensional distributions of the event numbers using Eqs. (92)-(94) of Sec. III. They are displayed in Fig. 11 as solid histograms. One can see that the BABAR data agree well with the corresponding lines. This agreement enforces the statement about the significant difference between the K + K − and K 0 S K + effective mass distributions which is due to the dominant f 0 (980) resonance contribution to the KK final state interaction amplitude.

V. CONCLUSIONS
A theoretical model of the D 0 → K 0 S K + K − decay amplitude has been constructed within a quasitwo-body QCD factorization approach introducing scalar kaon form factors to describe the S-wave kaon-kaon final state interaction. In doing so, the contribution of isoscalar-scalar f 0 resonance family, viz. f 0 (980), f 0 (1370) and that of the isovector-scalar a 0 one, viz. a 0 0 (980), a ± 0 (980), a 0 0 (1450), a ± 0 (1450) are taken into account. The isospin zero and one kaon-kaon S-wave interactions have been treated in a unitary way using either coupled channel relativistic equations, or a dispersion relation framework.
The P -and D-waves of the final state kaon-kaon interactions have also been taken into account. Independently of any model assumptions, we have shown that the K + K − andK 0 K + S-wave effective mass squared distributions, corrected for phase space, are significantly different. This means that, in the analyses of the D 0 →K 0 K + K − data, one cannot neglect the contribution of the f 0 (980) resonance and retain only the a 0 (980) contribution.
In Appendix A, we have updated the meson-meson S-wave isospin zero scattering amplitudes.
These include the three coupled, ππ,KK and an effective 2π 2π channels. Using the above amplitudes the new kaon nonstrange and strange form factors Γ n 2 (s 0 ) and Γ s 2 (s 0 ) have been calculated following Ref. [20] and introduced in the data analysis. As seen Fig. 14, these form factors are quite similar to those derived using the Muskhelishvili-Omnès dispersion relation approach [28,29].
In the factorization framework, for the D 0 → K 0 S K + K − process one has to evaluate the matrix elements of the D 0 transitions to two-kaons or the transitions between one kaon and two kaons. The knowledge of these transitions requires that of the three-body strong interaction between the D 0 , K 0 S and K ± mesons and that between the K 0 S , K + and K − mesons. Here, to describe these transitions with the two final kaons in S-wave state, we had to go beyond the simple multiplication of the scalar kaon form factors by a complex constant. And to obtain good fits we have multiplied the isoscalar-scalar kaon form factor by a one free parameter energy-dependent function and introduced into the isovector-scalar function an energy-dependent phenomenological polynomial involving three free parameters.
The undetermined free parameters of our seven D 0 → K 0 S K + K − amplitudes are then related to the strength of the isoscalar-scalar kaon form factor, to the function proportional to the isovector-scalar kaon form factor and to the unknown meson to meson transition form factors. They are obtained through a χ 2 minimization to the BABAR Dalitz plot distribution [2]. It should be pointed out that the low density of events in the central region of this Dalitz plot distributions (see Fig. 8) is difficult to reproduce. Using unitary relativistic equations to built the isoscalar-scalar form factor and a function proportional to the isovector-scalar one, we obtain a best fit (denoted R BF ), with a χ 2 /ndf of 1.25 with 19 free parameters to be compared to that of 1.28 for Ref. [2] which uses 17 free parameters.
In Appendix B, we have studied two alternative fits with scalar-kaon form factors derived in the Muskhelishvili-Omnès dispersion relation framework. All other amplitudes are parametrized as in the best fit model. If the scalar form factors are multiplied by energy dependent phenomenological functions, we obtain two good fits, one, denoted MO P 1 with a χ 2 /ndf of 1.32 and 16 free parameters and another one, MO P 2 , with a χ 2 /ndf of 1.31 and 16 free parameters.
Our fits indicate the dominance of the annihilation amplitudes and for the best fit a large dominance of the [K + K − ] S isospin 0 S-wave contribution and a sizable branching fraction to the [K 0 S K + ] P K − final state with the [K 0 S K + ] pair coupled to ρ(770) + , ρ(1450) + and ρ(1700) + . The alternative fits show important contributions from both the f 0 and a 0 0 mesons and a weaker ρ + mesons role. For all our models, the one-dimensional distributions agree well with that of the BABAR data.
One can estimate the strength of the contributions of the different amplitudes by looking at their branching ratio compared to the sum of their branching ratios. As can be seen in Table XII for the best fit model this sum 7 is 149.5 [126.3, 164.1] % (163.4% in Ref. [2]), which points to sizable interference contributions. The kaon-kaon S-wave interactions, related to the f 0 and a 0 0 resonances, gives a large branching of ∼61 [45, 63] % with a large value (for BF, MO P 1 and MO P 2 , see Table XIII) of ∼60 [23,46] %, for the amplitude proportional to the strange isoscalar-scalar form factor (f 0 7 The numbers in brackets are the corresponding values of the MOP 1 and MOP 2 fits, respectively contributions) and smaller branching ∼5 [16,16] % for the amplitude proportional to the isovectorscalar form factor (a 0 0 contributions). Corresponding figures in the isobar BABAR analysis [2] are ∼71 %, dominated by the a 0 (980) 0 and a 0 (1450) 0 with no f 0 (980) and a f 0 (1370) ∼ 2% .
The branching fraction of the isospin 0 P -wave ∼46 [45,45] %, dominated by the φ(1020) resonance, is similar to that found, ∼44 %, in Ref. [2]. The branching of the isovector amplitude associated to the a + 0 resonances is ∼21 [26,40] % to be compared to ∼ 45 % in Ref. [2]. The branching fraction of the amplitude related to the [ρ(770) + + ρ(1450) + + ρ(1700) + ]K 0 S final state, not introduced in Ref. [2], has a value of ∼22 [8,13] %. One could say that, this contribution with no bumps in the Dalitz plot distribution, is in Ref. [2] taken into account by a part of that of the a + 0 . The charmless hadronic B 0 → K + K − K 0 S studied by the Belle [43] and BABAR [44] Collaborations has the same meson final states as the D 0 decay we have been studied here. A quasi-two-body QCD factorization analysis of this B 0 decay process should allow, to constrain, not only the weak interaction observables but also the scalar kaon form factors, the transitions between one kaon and two kaons and to learn about the B 0 transition to two kaons.

Acknowledgements
We are grateful to François Le Diberder who, at the early stage of this work, has helped us in getting access to the BABAR data. We are deeply indebted to Fernando Martinez-Vidal from the BABAR Collaboration who provided us with experimental information for this study. We thank him for many fruitful exchanges. We are also indebted to Bachir Moussallam for very profitable correspondence and Here we present updated results for the meson-meson S-wave isospin zero scattering amplitudes.
They include the following three coupled channels: ππ, channel 1,KK, channel 2 and effective 2π 2π, channel 3. Our previous fits to the meson-meson scattering data were obtained in the late nineties [21,22]. Since, new precise low energy ππ data have appeared [45]. Moreover, as noticed by Bachir Moussallam [46], we used an assumption valid only below the opening of the third channel, namely the phase of the ππ →KK transition amplitude was set equal to the sum of the elastic ππ andKK phaseshifts. The derivation of the kaon isoscalar-scalar form factors Γ n 2 (s 0 ) and Γ s 2 (s 0 ), used in the present analysis for s 0min = 0.98 GeV 2 ≤ s 0 ≤ s 0max = 1.87 GeV 2 , requires the knowledge of the meson-meson amplitudes at energies above s 0max .
Thus, dropping the above mentioned assumption, we have performed a new analysis based on an enlarged set of data. Using the same three coupled-channel separable potential model as developed in Refs. [21] and [22], we fit the following data: The total number of fitted data is then equal to 140. As in Ref. [21], the fitting method is based on the χ 2 function being a sum of five components related to the five data sets enumerated above.
The resulting χ 2 is equal to 135.04 which, for 14 free model parameters, gives the value χ 2 /ndf = 1.07 when divided by ndf = 140 − 14 = 126 degrees of freedom.
The quality of our fit for the ππ phase shifts and inelasticities is illustrated in Fig. 12. As seen in Fig. 13 a good fit is achieved for the moduli and the phases of the T 12 amplitude. All experimental data sets are well reproduced by our phenomenological model. The resulting separable interaction parameters are listed in Table VIII, their notation being identical to that of Ref. [21].  Positions of the S-matrix poles in the complex energy E plane are given in Table IX. The signs of the imaginary parts of the channel complex momenta k i , i=1,2,3 are indicated in order to mark the corresponding pole position on different Riemann sheets. The total width Γ of a given pole equals to twice |ImE|.
As in the case of solution A (see Table 3 of Ref. [22]) one finds 18 poles. The first four (I to IV), lying on the real axis below the ππ threshold, are related to the S-matrix poles in the absence of  four poles (XI to XIV) attributed to the wider resonance f 0 (1400). The four poles (XV to XVIII) located between 1553 and 1584 MeV are responsible for the structure in the phase of the transition ππ → KK amplitude as can be seen in the right panel of Fig. 13 and in Fig. 6 of Ref. [49]; there is a maximum near 1500 MeV, close to the opening of the third channel, followed by a dip at about 1600 MeV. These latter poles could be related to the f 0 (1370) and f 0 (1500) resonances.
In Table X we present values of the moduli of the channel coupling constants calculated for five typical S matrix poles (for their definitions see Eq. (34) of Ref. [22]). The f 0 (500) poles like that with n • VII are mainly coupled to the ππ channel (i = 1). Also the four poles close to Re E = 1450 MeV have a strong coupling only to the ππ channel. The f 0 (980) poles n •s IX and X are preferentially coupled to theKK channel (i = 2) like the four other poles n •s XV to XVIII. This last group of poles has also a substantial coupling to the (2π)(2π) channel (i = 3) in addition to the strong coupling to theKK (i=3) one. All these poles lie above the opening of the third channel taking place at 2m 3 = 1504 MeV. In this appendix we complete our study by describing the results of two fits of the BABAR-Collaboration Dalitz-plot distribution [2] taking, in the amplitudes with final kaon-kaon states in S wave and isospin 0, the scalar KK form factors derived from the Muskhelishvili-Omnès (MO) approach [27]. The same parametrizations as those described in Sec. II are used for all other amplitudes.
In the MO dispersion-relation framework the isoscalar-scalar Γ n,s 2 (s) form factors have been calculated by B. Moussallam [28,29] from the MO equation using the updated ππ-T matrix of the ππ, KK and effective (2π)(2π) coupled-channel model of Ref. [22] (see Appendix A). In Fig. 14 the moduli of these MO form factors are compared to those derived in Sec. IV from a relativistic coupled-channel model.
In Sec. II one has introduced for the form factors Γ n,s 2 (s) complex phenomenological coefficients of proportionality χ n,s and in Sec. IV, to achieve good fits, notably to reproduce the low density of events in the central region of the Dalitz distribution (see Fig. 8, region II), we have been led to multiply them by the energy-dependent phenomenological functions P i (s 0 ) defined below in Eqs. (B2) and (B3). The isovector-scalar F , is plotted in Fig. 15. The isovector amplitudes associated to the isospin-1 a 0 0 and a + 0 resonances can be expressed in terms of this form factor by using, in the Eqs. (82) and (86), the relation (85) with G 1 (0) = χ 1 . The strength χ 1 is real and to obtain good fits, it was necessary to multiply it by the phenomenological where the free parameters c i , i = 1, 2, 3 and s are real.
An estimation of the phenomenological strength parameters χ n,s using Eq. factor with χ 1 |F  Fig. 15), one obtains χ 1 8.2 GeV −1 . As expected, from our study in Sec. III of the near threshold comparison between the K + K − S-wave effective mass projection with that of the K 0 K + , a good fit without the contribution associated with the isospin 0 f 0 resonances (χ n,s ≡ 0) cannot be obtained.
We also find that improved χ 2 are obtained with the δ 12 parameter of the isovector-scalar F form factor equal to 150 • (see Fig. 15). With the N p =16 free parameters displayed in Table XI, we obtain a fit, denoted as MO P 1 , with a total χ 2 of 1559.7 which corresponds to a χ 2 /ndf =1.32, not as good as that found in the best fit model of Sec. IV. In this fit, the phenomenological function multiplying the Γ n(s) 2 (s 0 ) is chosen to be with the zero s of the function P (s 0 ) [Eq. (103)] at 0 GeV 2 . Fixing |χ n | to 35 (GeV −1 ), χ 1 to 15 (GeV −1 ) a slightly better fit, denoted MO P 2 , with N p = 16 and a total χ 2 of 1546.9 (χ 2 /ndf =1.31) is obtained with   The different branching fractions Br i , i = 1 to 7, of these two fits are compared to those of the best   Table XIII show that the isoscalar f 0 and isovector a 0 0 resonance contributions can be quite different. However, taking into account the large uncertainties in the Br 1 values of the MO P 1 , MO P 2 and R BF fits (see Table XII) the total scalar-resonance contribution in the M 1 amplitude is similar. The corresponding results can be qualitatively interpreted from, the expressions of the amplitudes given in Sec. II D, the values of the different parameters given in Tables II,   XI and Table XIII. The branching fractions can also be partly compared to the fit fractions 9 of the BABAR isobar-model experimental analysis [5].
The dominance of the branching fraction associated to the isoscalar-scalar amplitude M s,I=0 1 for the best fit and to a less extent for the MO P 1 one, can be understood as their moduli |M entering the G 1 (s 0 ) function of the best fit.
In the MO P 2 fit (Br(M n,I=0 1 )= 20.66 %) the f 0 contribution in Γ n 2 (s) is enhanced by a larger |χ n | (35 GeV −1 ) and by the function P 1 (s 0 ) while it is more suppressed in the best fit (1.19 %) than in the MO P 1 one (4.75 %).
A striking difference arises for the modulus of the phenomenological transition form factor : it is quite large, 2.22, for the best fit solution R BF as compared to its magnitude for the two other fits where it is smaller than 0.43. This leads to large value of the branching fraction (59.82 %) for the M s,I=0 1 (s 0 ) amplitude arising from D 0 annihilation via W exchange in the R BF best fit solution (see Table XIII). The |χ s | being the same (26 GeV −1 ), the difference between the MO P 1 (22.69 %) and MO P 2 (45.91 %) branching fractions arises from smaller P 1 (s 0 ) enhancement (for s 0 1.3 GeV 2 ) than that of P 2 (s 0 ), as can be seen in Figs The branching fraction Br 3 , which indicates the isovector a + 0 resonances contribution, has values of 25.9, 40.1 and 20.7 % for the MO P 1 , MO P 2 and R BF fits, respectively (see Table XII). The modulus of the transition form factor F represent the corresponding moduli for the alternative MO P 1 , MO P 2 fits, respectively. As in Fig. 15 for the two vertical dashed lines (the s + limits are very close to the s 0 ones) . the MO P 1 , MO P 2 and R BF fits are equal to 7.7, 13.1 and 21.5 %, respectively (see Table XII). The large contribution in the R BF solution is partly due to the magnitude, 9.38, of the modulus of the transition form factor |A K − ρ + 0 (m 2 D 0 )| to be compared to 5.78 and 7.94 for the MO P 1 and MO P 2 fits, respectively. The Br 4 ratio between that of the R BF and those of the MO P 1 and MO P 2 fits is close to the square of the corresponding |A K − ρ + 0 (m 2 D 0 )| ratios. It can be seen that, to improve the χ 2 of the fits, it seems necessary to increase the ρ + resonances contributions. The f 2 (1270) resonance contributions in Br 7 for our three fits follow the evolution of the square of the |P D | parameter in each fit. They are very small and even smaller than in the BABAR analysis [5].