$D_s$ decays to $\eta$ and $\eta^\prime$ final states: a phenomenological analysis

We consider the semileptonic and nonleptonic $D_s$ decay modes to final states with $\eta$ and $\eta^\prime$. We use QCD sum rules to determine the $D_s \to \eta$ form factor $f_+^\eta$, and a generalized factorization ansatz to compute nonleptonic decays. We propose a parameterization of possible OZI suppressed contributions producing the $\eta^\prime$ in the final state, compatible with current data; such a scheme can be further constrained improving the precision of the measurement of the $D_s$ decay rates, as expected by the ongoing experiments.


Introduction
The exclusive D s decays to final states containing η and η ′ represent nearly 30% of the total decay rate of the D s meson. Therefore, D s could be a suitable system where to gather information on important aspects of the η − η ′ phenomenology, namely the longstanding issue of the η − η ′ mixing. Moreover, D s can be used to further investigate some unsettled aspects of nonleptonic heavy meson decays, such as the anomalously large η ′ production observed in several heavy meson decay channels. Examples are B − → K − η ′ and D 0 → K 0 η ′ , the measured decay rates of which are substantially larger than what can be expected by naive theoretical calculations.
The current experimental situation concerning the D s transitions to η and η ′ is summarized in Table 1 [1], mainly using the results obtained by the CLEO Collaboration in the past few years [2]. The experimental results are expected to be improved in the Table 1: Experimental rates and branching fractions of semileptonic and nonleptonic D s decays to final states containing η and η ′ .
The results in Table 1 have inspired several considerations. First, it has been proposed that information on the η − η ′ mixing could be obtained just considering the semileptonic decay modes. As a matter of fact, writing the hadronic matrix element governing the transition D + s → ηℓ + ν in terms of form factors: (q = p − p ′ ) and a similar expression for D + s → η ′ ℓ + ν, the ratio B(D + s → η ′ ℓ + ν) B(D + s → ηℓ + ν) could be used to access the η −η ′ mixing angle through the ratios of the form factors f η ′ ± (q 2 )/f η ± (q 2 ) which are related to the η − η ′ mixing scheme [3,4]. In particular, information could be gathered on the mixing scheme in the flavour basis [5,6], which consists in writing the η and η ′ states as combinations of |η q >= 1 √ 2 |uu + dd > and |η s >= |ss >: It has been shown [5] that in this scheme a single angle is essentially required, since |φ s − φ q |/(φ s + φ q ) ≪ 1, a result confirmed by a QCD sum rule calculation [7]. Therefore, one can safely assume φ s ≃ φ q ≃ φ; the most recent estimates of φ give values close to 40 0 [6,8]. In the flavour scheme, the semileptonic form factors relative to D + s → ηℓ + ν and D + s → η ′ ℓ + ν satisfy the relation so that the possibility of a direct comparison with the results for φ obtained from the analyses of other channels involving η − η ′ particles could be envisaged. The situation is particularly simple in the case of semileptonic D + s decays to positrons or antimuons, where essentially only the form factors f η(η ′ ) + (q 2 ) are involved. However, in order to pursue this program, one has to neglect possible contributions to the semileptonic decay amplitude from diagrams where η and η ′ are produced through gluon emission; we shall consider this problem below.
As for nonleptonic decays, naive factorization, using the semileptonic D s → η and D s → η ′ form factors and the Wilson coefficients relevant for the transitions in Table 1, does not allow to predict all the branching fractions of D + s → η (′) π + and D + s → η (′) ρ + [9]. The same conclusion is obtained by analyzing the various decay channels in terms of transition amplitudes related by SU(3) F symmetry to analogous amplitudes for D decays [10], or accounting for some effects of the inelastic final state rescattering [11]. In particular, the prediction for the rate of the decay mode D + s → η ′ ρ + is lower than the experimental measurement by more than a factor of two. This is disappointing: in Cabibbo favoured hadronic D s decays the final state contains a single isospin mode, thus ruling out possible interference effects due to the elastic final state interactions; moreover, the conservation of G-parity does not allow to include inelastic effects of intermediate states consisting of ordinary mesons in the D + s → η ′ ρ + mode [12]. Therefore, a different mechanism must be invoked to explain the enhanced η ′ production. It has been suggested that the enhancement could be due to OZI suppressed diagrams with the η ′ produced by gluons and the cs pair annihilating to a charged W [3]. This mechanism would not affect substantially the η production, since the coupling of the gluons to η is estimated to be smaller than the coupling to η ′ [13]. However, a mechanism of this type, violating the OZI rule, could also affect the semileptonic D + s → η ′ ℓ + ν transition, spoiling the possibility of using the relation (3) to gather information on the angle φ from the semileptonic decay rates. Moreover, these effects could be also present in other systems, namely in D decays, although in such cases the annihilation amplitudes are Cabibbo suppressed.
The effects of the gluon production of the η and η ′ , although plausible, are not included in ordinary analyses since they are difficult to take into account in a quantitative way.
Nevertheless, their investigation is of particular relevance, and we shall try to perform it in a phenomenological way.
In this paper, we compute the form factor f η + (q 2 ) relative to D + s → ηℓ + ν, showing that the result is in agreement with the experimental measurement in Table 1. On the other hand, assuming the standard value of the η − η ′ mixing angle together with the naive factorization, other results in Table 1 are not reproduced. Therefore, we adopt a generalized factorization ansatz, fitting the relevant parameters from the experiment; moreover, we assume that the effect of the process producing the η ′ through the annihilation of the cs pair numerically modifies the D s → η ′ form factors. This enables us to investigate whether the experimental results can be reproduced by this assumption and how the reduction of the experimental uncertainty can be used to test various consequences of our ansatz.
2 QCD sum rule calculation of f η + (q 2 ) Let us first compute the form factor f η + (q 2 ) using a nonperturbative method, such as the QCD sum rule technique [14]. We adopt the usual strategy of considering a three-point function: with J η 5 =siγ 5 s the pseudoscalar quark density probing the strangeness content of the η, J µ =sγ µ c the weak current inducing the c → s transition, and J Ds 5 =ciγ 5 s a quark current having the D s quantum numbers. The momenta P and q are defined as P = p + p ′ and q = p − p ′ , respectively. For the invariant function Π + (p 2 , p ′2 , q 2 ) a double dispersion relation in the variables p 2 , p ′2 can be written down: where possible subtraction terms have been omitted. The spectral function ρ(s 1 , s 2 , q 2 ) contains, for low values of s 1 and s 2 , a double δ−function corresponding to the transition Isolating such a contribution, and neglecting possible subtraction terms which we discuss later on, we can write: (6) In (6) we have assumed that the contribution of higher resonances and continuum of states starts from the effective thresholds s 0 1 and s 0 2 . The hadronic parameter A represents the matrix element: while the projection of the J Ds 5 current on the D s state is given by the matrix element The correlator (4) can be computed in QCD for large Euclidean values of p 2 and p ′2 by an Operator Product Expansion, expanding the T-product in (4) as a sum of a perturbative contribution and non perturbative terms, proportional to vacuum expectation values of quark and gluon gauge invariant operators of increasing dimension, the vacuum condensates. In practice, only the first few condensates numerically contribute, the most important ones being the dimension 3 <ss > and dimension 5 <sgσGs > condensates. The QCD expression for Π + reads: Invoking quark-hadron duality, i.e. assuming that the hadronic and the perturbative QCD spectral densities give the same result when integrated above the thresholds s 0 1 and s 0 2 , we get the sum rule: The variables r and r ′ are defined as r = p 2 − m 2 c and r ′ = p ′ − m 2 s . The domain D is bounded by the curves where ∆ = s 1 − m 2 c + m 2 s , and by the lines s 2 = s 0 2 and s 1 = s 0 1 . Eq. (10) can be improved by applying to both its sides a Borel transform, defined as follows: where F is a generic function of Q 2 . The application of such a procedure to the sum rule amounts to exploiting the relation: with M 2 a Borel parameter. The operation, applied independently to the variables −p 2 and −p ′2 , improves the convergence of the series in the OPE in the r.h.s. of Eq. (9) by factorials in n, and, for suitably chosen values of the Borel parameters, enhances the contribution of the low-lying states in the hadronic representation of the correlator Π + .
Moreover, since the Borel transform of a polynomial vanishes, the procedure allows to get rid of subtraction terms in the dispersion relations, which are polynomials in p 2 or p ′2 . Therefore, a final sum rule can be worked out, keeping only the contribution of the lowest dimensional condensates: In the numerical analysis of (17)  and m s = 140 MeV [16,15]. As for the D s decay constant, we used f Ds = 225 MeV [15], while for the parameter A we adopted the two-point QCD sum rule result A = 0.115 GeV 2 [7]. The obtained sum rule shows stability to the variation of the Borel parameter in the region 2.5 GeV 2 ≤ M 2 1 ≤ 3.5 GeV 2 and 1.6 GeV 2 ≤ M 2 2 ≤ 2.4 GeV 2 , with the thresholds s 0 1 and s 0 2 in the ranges s 0 1 = 5.9 − 6.1 GeV 2 and s 0 2 = 0.9 − 1.1 GeV 2 , respectively. The form factor f η + (q 2 ), obtained in the range of momentum transfer 0 ≤ q 2 ≤ 0.5 GeV 2 , is depicted in Fig. 1; it can be fitted by a linear expression with A = 0.14 GeV −2 and B = 0.50±0.04. This expression is consistent, in the considered range of momentum transfer, with a polar form f η , with the mass of the pole M P ≃ 1.9 GeV .
In the following, we shall consider the form factor f η + (q 2 ) as a theoretical input in a phenomenological analysis of D s transitions.

D s transitions to η and η ′
The form factor f η + (q 2 ) computed above allows us to calculate the semileptonic D + s → ηℓ + ν decay rate. It can also be used to analyze the nonleptonic modes D s → ηπ + and D s → ηρ + if the factorization approximation is adopted. This amounts to consider the effective Hamiltonian with (q 1 q 2 ) V −A =q 1 γ µ (1−γ 5 )q 2 and C 1 and C 2 Wilson coefficients, and factorize the V −A currents appearing in it. As for the modes with η ′ , we further need an input on the η − η ′ mixing, and we choose the angle φ in the flavour basis mixing scheme, with the value φ = 39 0 coming from the measurements of φ → η (′) γ [17]. In Table 2 Table 1; also the result Table 2: Computed semileptonic and nonleptonic D s rates and branching fractions. Nonleptonic rates are obtained using naive factorization. The η − η ′ mixing is described in the flavour basis, with mixing angle φ = 39 0 .
, is within the experimental uncertainty quoted in Table 1. On the other hand, as one can infer by comparing the computed decay rates reported in Table 2 with the experimental measurements in Table  1, the calculations of the nonleptonic modes do not fit all the experimental measurements, as already anticipated by previous analyses.
In order to parameterize the deviation from the factorization approximation, as well as the possible role of the η and η ′ gluon production, we adopt a generalized factorization ansatz, consisting in substituting the combination of the Wilson coefficients a 1 = C 1 + C 2 N c with effective scale-independent parameters a ef f 1 in the factorized amplitudes. The coefficients a ef f 1 should be considered as non-universal, process-dependent parameters [18]. However, since in the decay modes D + s → ηπ + , D + s → η ′ π + , and analogously D + s → ηρ + , D + s → η ′ ρ + , the underlying process is the same, we assume only two process-dependent parameters to describe the deviation from naive factorization: a ef f 1,π describing D + s → ηπ + and D + s → η ′ π + , and a ef f 1,ρ describing D + s → ηρ + and D + s → η ′ ρ + .
As for the possible contribution of OZI suppressed diagrams producing η and η ′ , it is essentially related to the matrix elements 0| GG η (′) , where G is the gluon field andG its dual. Several theoretical investigations suggest that 0| GG |η ≪ 0| GG |η ′ [13,7]; therefore, we assume that such annihilation amplitudes mainly affect the D s transitions to η ′ . A simple parameteterizion consists in modifying the values of the parameters A, B in (18), thus without affecting the shape of the form factor f η ′ + . This seems rather reasonable, since the range of momentum transfer in D s → η ′ transitions is rather narrow (q 2 ≃ 0 for , and a linear expansion is a suitable representation of the form factors. Therefore, in the case of η ′ , we phenomenologically represent the D s → η ′ form factor as It is now possible to use the experimental data in Table 1 to fit all the parameters we have introduced, namely a ef f 1,π , a ef f 1,ρ ,Ā andB. From the decays D + s → ηπ + and D + s → ηρ + we find that the values of a ef f 1,π and a ef f 1,ρ are bound in the ranges: to be compared with the value of a 1 obtained from the Wilson coefficients C 1 and As for the decay mode D + s → η ′ π + , it involves f ef f + (q 2 ); however, only the value f ef f + (0) =B is needed in the approximation M π = 0, allowing us to constrainB in the range Moreover, considering the modes D + s → η ′ ℓ + ν and D + s → η ′ ρ + , we find that the relations

Ds
) constrain the pa-rametersĀ andB in selected regions of the (B,Ā) plane. These regions are delimited by two straight lines, from the datum on the nonleptonic D + s → η ′ ρ + decay rate, and by two ellypses corresponding to the measurement of the semileptonic D + s → η ′ ℓ + ν decay rate.
Considering simultaneously the constraints, all the data in Table 1 can be fitted if, in the (B,Ā) plane, overlap regions exist among the area delimited by the ellypses from Eq.(24), the regions delimited by the straight lines from Eq.(23) and the regions between the vertical lines from Eq. (22). At the present level of accuracy of the experimental data in Table 1 such regions indeed exist. They are depicted in figure 2 and denoted as D 1 , Although it is expected and rather plausible, the existence of such overlap regions was not guaranteed a priori; it shows that we have chosen a sensible scheme to parameterize the decays in Table 1. More important, we expect that an improvement in the accuracy of the experimental data on the D s decay rates would sensibly reduce the size of such overlap regions, and presumably, exclude some of them. Noticeably, already at the present level of accuracy some interesting observations can be drawn. Let us consider, for example, the parameters in the regions D 3 and D 2 . In both the cases the experimental branching fraction of the semileptonic decay mode D s → η ′ ℓν is reproduced. However, a prime difference is that in the region D 3 the parametersĀ andB are opposite in sign, while in the region D 2 they have the same sign. This implies that the relation between the D s → η ′ and D s → η form factors in (3) cannot be satisfied by the parameters in the region D 3 . The same conclusion holds for the region D 1 . The opposite signs betweenĀ andB, as it happens in the regions D 1 and D 3 , have an observable consequence in the spectrum of the semileptonic decay D + s → η ′ ℓ + ν: in this case, a zero in the dΓ(D s → η ′ ℓν)/dq 2 distribution should be observed, as depicted in can be done between the two shapes of the semileptonic distribution.
We can reasonably expect that improved data would restrict the allowed regions in the (B,Ā) plane. It could happen that they do not intersect any more, or that intersection regions could be found with restricted extension, allowing a better determination of the effective parameters introduced in our analysis. The calculation of such parameters remains a challenging task, and we do not attempt it in the present paper. However, it is worth outlining the theoretical framework in which the calculation could be carried out.
Concerning the effective coefficients a ef f 1,π and a ef f 1,ρ , which take into account the deviation from the naive factorization in the corresponding decay modes, their theoretical calculation would consist in a precise determination of nonfactorizable contributions. A step in this direction has been recently performed in the case of some two-body nonleptonic B decays, where the meson picking up the B spectator quark is light, exploiting the large value of the beauty quark mass [19]. In this case, it has been observed that nonfactorizable contributions are of order α s or 1/m b , and a QCD factorization formula has been written for the nonleptonic matrix elements in terms of meson light-cone distribution amplitudes. A possible extension of such a procedure to charm requires the development of a reliable method for computing at least the first (process dependent) 1/m Q correction.
A different approach would consist in considering the corrections to the large N c limit, where factorization becomes exact [20]: also in this case, however, next-to-leading 1/N c terms are generally sizeable, and one needs their actual calculation. Therefore, it seems worth attempting to gain information on the effective coefficients from phenomenological analyses, as done, for example, in [18].
As for the OZI suppressed diagrams producing the η ′ through its coupling to the gluons, together with a weak annihilation of D s , a perturbative calculation could be carried out in QCD, in analogy with the calculation of the η ′ production in quarkonium decays [21,22]. The difference, in the present case, is that one has to account also for the gluon emission from a light (strange) quark, and one cannot exploit the fact that all the quarks involved are heavy, which justifies the application of perturbative QCD methods. The calculation, for small values of q 2 , produces an amplitude for D s → η ′ ℓν of the same form as provided by a linear q 2 representation of the D s − η ′ form factor. An important ingredient in this perturbative calculation is the actual value of the two-gluon-η ′ matrix element < g(k 1 )g(k 2 )|η ′ (p) > describing the vertex η ′ gg for off-shell gluons. Such a matrix element is parameterized by a form factor F (k 2 1 , k 2 2 ) whose value at k 2 1 = k 2 2 = 0 is fixed by the QCD anomaly; as for the momentum dependence, various parameterizations have been proposed in the literature, thus providing different values for the effective parameters A andB introduced in our analysis, which in turn could correspond to various solutions for the spectrum shown in fig.3. One might notice some analogies with the analyses which explain the observed enhancement of the η ′ production in B decays through the mechanism of gluon fusion [23].
All such considerations taken into account, we believe that our proposed scheme, where additional contributions are reabsorbed in the parametrization of the D s → η ′ form factor and in a ef f 1,π , a ef f 1,ρ is useful from the phenomenological point of view, as a starting point for the investigation of the underlying dynamics, and could be extended to other cases.
Before concluding, we want to mention a check of consistency. If we consider the decay mode D + s →K 0 K + , which can be related to D + s → ηπ + through SU(3) F symmetry, and describe the D s →K 0 form factor by f η + (q 2 ), together with f K = 0.160 GeV , we can estimate the effective parameter a ef f 1,K . The experimental measurement B(D + s → K 0 K + ) = (3.6 ± 1.1) × 10 −2 produces a ef f 1,K ∈ [0.72, 1.14], i.e. the effective parameter a ef f 1,K displays a significant overlap with the range determined for a ef f 1,π . In different words, from our analysis and assuming SU(3) F , we would be able to predict rather accurately the experimental datum for B(D s →K 0 K + ).

Conclusions
We have presented a phenomenological analysis of the D s decays to final states containing η and η ′ . Since the theoretical investigations based on SU(3) F symmetry, FSI effects and standard η − η ′ mixing failed in simultaneously reproducing the observed branching ratios for all these decays, we have considered a possible role of annihilation diagrams, in which the η ′ is produced through its coupling to gluons. We have proposed a parametrization of those effects in the D s → η ′ form factor. As for D s → η, we used a theoretical calculation of the form factor f η + (q 2 ) which corresponds to a branching fraction for the decay D s → ηℓν in agreement with data. A fit to all the available experimental results, adopting a generalized factorization scheme for nonleptonic decays, is possible; it constrains the parameters in restricted regions that can be discriminated by making dedicated observations, for example looking at the semileptonic spectrum of the D s → η ′ transitions.
An improvement in the precision of the experimental data on D s decays could support this scheme and be helpful in understanding the dynamics of the η and η ′ production in heavy meson decays.