Direct CP asymmetry in $D\to \pi^-\pi^+$ and $D\to K^-K^+$ in QCD-based approach

We present the first QCD-based calculation of hadronic matrix elements with penguin topology determining direct CP-violating asymmetries in $D^0\to \pi^-\pi^+$ and $D^0\to K^- K^+$ nonleptonic decays. The method is based on the QCD light-cone sum rules and does not rely on any model-inspired amplitude decomposition, instead leaning heavily on quark-hadron duality. We provide a Standard Model estimate of the direct CP-violating asymmetries in both pion and kaon modes and their difference and comment on further improvements of the presented computation.


I. INTRODUCTION
Despite years of intense experimental efforts, CP-violation has never been unambiguously observed in the decays of up-type quarks. In the Standard Model (SM) this fact can be explained by the suppression of all CP-violating amplitudes resulting from the smallness of relevant Cabbibo-Kobayashi-Maskawa (CKM) matrix elements. To make matters worst, accurate predictions of up-type CP-violating observables are hard to obtain, since the necessary hadronic matrix elements are dominated by long-distance contributions. In order to calculate these matrix elements one needs to employ a QCD-based method that deals with nonperturbative effects in a model-independent manner. In this letter we shall calculate CP-violating observables in exclusive singly Cabibbo-suppressed (SCS) decays of D-mesons using a variant of light-cone QCD sum rules (LCSRs).
Observables that are sensitive to CP -violation are most often written in terms of asym- formed from the partial rates of a D-meson decay to a final state f and of its CP-conjugated counterpart. Depending on the initial state, the asymmetry in Eq. (1) could be a function of time, if D 0 D 0 -mixing is taken into account.
The measured time-integrated asymmetry contains a direct component (see, e.g., [2]), which will be the main focus of this paper. Direct CP asymmetry occurs when the absolute values of the D → f decay amplitude, which we denote by A f ≡ A(D → f ), and of the corresponding CP-conjugated amplitude Af ≡ A(D → f ) are different. This can be realized if the decay amplitude A f can be separated into at least two different parts, where φ 1 = φ 2 are the weak phases (odd under CP ), and δ 1 = δ 2 are the strong phases (even under CP ). The CP-violating asymmetry is then given by The amplitude pattern of Eq.
A model-independent computation of the amplitude ratio and the strong phase difference in Eq.
(3) is a daunting task in charm physics (for reviews see e.g., [3]). In SM, the sizes of direct CP-violating asymmetries for any final state are always proportional to a small combination of the CKM factors V * cb V ub , which ensures that a dir CP (f ) is small even if the maximal strong phase difference is assumed. Hence, if a larger value of CP-violating asymmetry is observed, it would provide a "smoking gun" signal of new physics in the charm sector of quark-flavor transitions. However, due to a very uncertain hadronic input, the available estimates of a dir CP (f ) obtained with various degree of model dependence are mostly qualitative, predicting this asymmetry in the ballpark of a per mille.
Direct and indirect components of CP -violating asymmetries in the decays of neutral D-mesons can be separated with a careful time-dependent analysis. Since the indirect components to a large extent are independent of the final state, it is advantageous to consider differences of direct CP-violating asymmetries ∆a dir CP , This difference is especially interesting if CP asymmetries in the subtracted amplitudes are predicted to have opposite signs. This is in fact realized in SM for f 1 = K − K + and f 2 = π − π + final states.
Earlier experimental results seemed to indicate a somewhat large asymmetry (4) for these final states, with values reaching the order of −1.0%. If confirmed, this would have indicated a possible new physics contribution to flavor-changing neutral currents (FCNC) in charm sector [4] or a previously unaccounted SM contributions [5].
Both asymmetries imply a very small effect, consistent with zero within current experimental uncertainties. Note again that in the SM opposite signs are expected for the asymmetries in the K − K + and π − π + channels.
New results with smaller experimental uncertainty are expected from the Run II LHCb data, as well as from the Belle II experiment. A model-independent calculation with controlled theoretical uncertainties of direct CP-violating asymmetries in SM is, therefore, compellingly needed. However, this task necessitates a calculation of hadronic matrix elements with strongly interacting and energetic two-meson final states, which is a big challenge even for the most advanced lattice QCD methods. In this situation even an order-of-magnitude QCD-based estimate of the hadronic input should become useful to reliably constrain the expected SM contribution to the CP-violating asymmetries.
The aim of this letter is to estimate the hadronic matrix elements relevant for the direct CP-violating asymmetries in D 0 → P − P + decays (P = π, K), employing a computational method which combines the light-cone operator-product expansion (OPE) in QCD and hadronic dispersion relations. More specifically, we use the approach developed in [11] for B → ππ nonleptonic decays, in particular its application to the penguin-topology matrix elements [12]. We define matrix elements of "penguin topology" as those of the weak effective 4-quark operator containing a quark-antiquark pair not present in the valence-quark content of the final P − P + state.
In what follows, we identify the hadronic matrix elements with penguin topologies which are needed to estimate the direct CP asymmetry in D 0 → π − π + and D 0 → K − K + decays.
We then briefly describe the method of Refs. [11,12], adapting it to the D 0 → P − P + decays. The main result is the QCD LCSR for the D 0 → P − P + hadronic matrix elements with the penguin topology. The calculation, which takes into account SU (3) F -violating O(m s ) effects, is valid at large invariant mass P 2 of the π − π + and/or K − K + final state.
Applying quark-hadron duality, the result is analytically continued to the physical point P 2 = m 2 D . The rest of this letter contains numerical analysis and our estimate of the direct CP-asymmetry in D 0 → P − P + decays, followed by a concluding discussion.
The singly Cabibbo suppressed decays of charmed mesons are driven by an effective where the products of CKM matrix elements λ q are defined as Unitarity of the CKM matrix 1 implies that Since the goal of this paper is to capture the dominant contributions to the decay amplitudes in the SM, we shall only take into account the effective current-current operators, where Γ µ = γ µ (1 − γ 5 ), and q = d, s. We shall neglect the penguin operators Q i=3,...,6,8g with small Wilson coefficients.
Furthermore, we introduce a compact notation for the linear combination of the operators (12) with their Wilson coefficients and the Fermi constant, The dominant contribution to the two-body D 0 → P − P + nonleptonic decay amplitudes is given by the hadronic matrix elements of O q , Applying the CKM unitarity relation of Eq. (11) to Eqs. (14) and (15), and subsequently adding and subtracting a term λ b π − π + |O s |D 0 to right-hand side of Eq. (14), we arrange the decay amplitudes in the following form, where we introduce compact notations for the ratios of hadronic matrix elements and and denote by δ π(K) the difference between the strong phases of the amplitudes P s(d) ππ(KK) and A ππ (KK) . In what follows, we do not attempt to calculate the amplitudes A ππ and A KK , having in mind their complicated form in terms of hadronic matrix elements with various quark topologies. Instead, as follows from Eq. (16), these amplitudes can be extracted to a reasonable precision from the measured partial widths of D 0 → π − π + and D 0 → K − K + decays, neglecting small parts of the amplitudes proportional to λ b .
It is instructive to discuss the flavor SU (3) F -symmetry limit of the decay amplitudes in Eq. (16). In this limit separate hadronic matrix elements with pions and kaons in the final state are equal: A KK = A ππ , r K = r π and δ K = δ π . Still, with λ b = 0, the decay amplitudes in Eq. (16) differ from each other by O(λ b ) terms. This difference can be easily understood in the U -spin symmetry limit, in which the initial D 0 -state is a U -singlet. The effective Hamiltonian of Eq. (9) transforms as a combination of a U -triplet and U -singlet, the latter being proportional to λ b , so that there are two U -spin invariant amplitudes contributing to The representation in Eq. (16) has the advantage that the parts of the amplitudes proportional to λ b /λ s generate direct CP -violating asymmetries, due to the weak phase contained in this combination of CKM parameters. The asymmetries will vanish if either the ratio r π,K → 0 or the strong phase δ π,K → 0. Hence, computation of the ratios r π,K and the phases δ π,K will result in the prediction of direct CP-violating asymmetries of Eq. (4).
It is important to note that we will not be using the expansion of the decay amplitudes in flavor topology diagrams, which is frequently employed in the analysis of two-body nonleptonic decays of charmed meson, starting from the earlier papers [15] (for a more recent analysis see, e.g, [16]). In that approach, flavor symmetries and experimental data are used to fit the "topological amplitudes". In our calculation, such expansion is unnecessary, first of all because we are only calculating the "penguin topology" matrix elements in Eq. (18), estimating the dominant part of the decay amplitudes from experimental measurements.
Most importantly, in the commonly adopted convention, since the part of decay amplitude containing weak phase is suppressed by a very small CKM coefficient λ b , it is not feasible to extract this part from the experimentally observed decay rates, but rather calculate it directly, as it is done in this letter.

III. HADRONIC MATRIX ELEMENTS FROM LCSRS
Here we adapt for nonleptonic D-decays the approach to compute hadronic amplitudes for B → ππ decays suggested in [11] and based on the method of QCD LCSRs [17]. In particular, we will readily use the calculation of hadronic matrix elements of B → ππ decays with charm penguin topology performed in [12]. We start from the D 0 → π − π + transition aiming at estimating the s-quark "penguin" contribution P s ππ to this decay amplitude. Following [12], the starting object is the correlation function where j (π) α5 =dγ α γ 5 u and j (D) 5 = im cc γ 5 u are, respectively, the pion and D-meson interpolating currents, sandwiched together with the four-quark operator Q s 1 between the on-shell pion and vacuum state. The ellipsis denote the kinematical structures we do not use. Note that performing a Fierz transformation of the operator Q s 1 , the combination of current-current operators entering the effective Hamiltonian of Eq. (9) is transformed: so that the color-octet operator is the only one contributing to the penguin amplitude in the adopted approximation. Hence, we hereafter replace Q 1 → Q s 2 in the correlation function. The resulting hadronic matrix element obtained from LCSR below is related to the penguin matrix element: Following [11] we introduce an auxiliary 4-momentum k flowing from the vertex of the weak interaction in the correlation function (20) and assume k 2 = 0, and, for simplicity, p 2 = 0.
We adopt massless u and d quarks and a massless pion (q 2 = 0) approximation, so that the invariant amplitude determining the correlation function depends on three invariant variables As in [12], the essential OPE diagrams are the ones shown in Fig. 1. We remind the reader that these diagrams stem from the light-cone OPE of the correlation function. Each diagram contains a coefficient function (hard-scattering amplitude) calculated perturbatively and convoluted with the pion distribution amplitudes (DAs) of growing twist and multiplicity. Therefore, only the highly-virtual quark and gluon lines corresponding to near-light-cone separation (large spacelike momentum transfer) between the two external currents with momenta p − q and p − k are included explicitly. In particular, the s-quarks in the loops also have large virtualities. The small-virtuality quarks and gluons are by default included in the pion DAs. Hence, for example, the diagrams with gluons emitted from the s-quark loop and absorbed in the pion DA should not be included in OPE 3 .
Here, as in [12], we only keep the contributions of two-particle twist-2 and twist-3 DAs and neglect all contributions of multiparticle pion DAs, since they correspond to highertwist contributions to OPE and are suppressed by powers of characteristic large scales with 3 Similar diagrams with heavy c-quark loop in B decay LCSRs as discussed in Ref. [12] remain the part of OPE, albeit being the part of higher-twist power corrections. respect to the lowest-twist contributions. In particular, we neglect contributions of diagrams where low virtuality ("soft") gluons are emitted from the virtual quarks and form quarkantiquark-gluon DAs of twist 3 and 4. On the other hand, following [12] we take into account the diagram (c) containing the factorizable four-quark component of the pion DA, expressed via two-particle DA and vacuum quark condensate density. In the case of heavy charm-quark loop in B → ππ-decays considered in [12] there is a second diagram of this type with thē dd-condensate. In our correlation function such a diagram would only have small-virtuality s-quarks which cannot be resolved from the four-quark DAs, hence, it is absent in the context of OPE (see also discussion in [12]). Note that the s-quark pair in the correlation function, being generated from the V-A current, cannot itself form thess condensate. Hence, the small virtuality s-quark pairs can only form "genuine" nonfactorizable four-particle pion DAs originating from the nonlocal matrix element 0|d(y)s(0)s(0)u(x)|π + and described by the diagram (d). As mentioned above, we neglect these contributions.
The expressions for the diagrams in Figs. 1(a)-(c) are then taken from [12], replacing b-quark by c-quark and c-quark in the loop by s-quark. Derivation of LCSR follows then the same three-step procedure as the one developed and discussed in detail in [11] and used in [12]. The first step is to use the hadronic dispersion relation for the invariant amplitude F ((p − k) 2 , (p − q) 2 , P 2 ) in the variable (p − k) 2 = s of the pion channel, keeping the variables (p − q) 2 and P 2 fixed and spacelike, applying quark-hadron duality approximation in this channel with the effective threshold s π 0 and performing the Borel transformation (p − k) 2 → M 2 1 . The second step is the transition from spacelike to large timelike P 2 = m 2 D , assuming local quark-hadron duality. At this stage we obtain the relation: for the matrix element with two-pion final state, produced by the product of the operator Q s 1 and D-meson interpolating current with fixed (p − q) 2 < 0; f π in the above is the pion decay constant. The final third step is to apply to the above matrix element the hadronic dispersion relation in the variable (p − q) 2 = s , combined with the quark-hadron duality approximation in the D-meson channel with the effective threshold s D 0 and followed by the Borel transformation (p − q) 2 → M 2 2 . After that the auxiliary 4-momentum k vanishes for the D-meson pole term in the dispersion relation, leading to the sum rule for the D 0 → π − π + hadronic matrix element: where f D is the D-meson decay constant. The right-hand side of this expression is obtained by computing the double imaginary part of the sum of OPE diagrams, resulting in the final form of LCSR, where u D 0 = m 2 c /s D 0 . The above expression is valid at large positive P 2 s π 0 , and the power corrections O(s π 0 /P 2 ) are neglected. In Eq. (26) the quark-loop function is defined [12] as Note that the contributions from diagrams in Fig. 1a (Fig. 1b) are presented in the second and third (fourth) lines of Eq. (26), and the one from Fig. 1c occupies the fifth and sixth lines of Eq. (26). The pion DAs entering the LCSR are ϕ π (u) and ϕ π 3p (u), ϕ π 3p (u), of twist-2 and twist-3, respectively. They depend, as usual, on the fraction of the longitudinal momentum of the pion u (ū = 1 − u) carried by the quark (antiquark). We define φ σ 3π (u) = dφ σ 3π (u)/du and adopt the truncated expansion in Gegenbauer polynomials C (α) n (u −ū) for these DAs. We take for the twist-2 DA, and for the twist-3 DAs. Note that the Gegenbauer moments a π 2,4 entering the twist-2 DA, the normalization parameter µ π = m 2 π /(m u + m d ), and the parameters f 3π , ω 3π of the nonasymptotic parts in twist-3 DAs represent nonperturbative inputs determined from various sources. Their scale dependence is not shown for brevity and taken into account at leading order (see e.g., [18] for more details). In Eq. (26) , the operator Q s 1 → Q d 1 , and the final state |π + → |K + . Correspondingly, the diagrams in Fig. 1 will change their flavor content, in particular, the s-quark loop will be replaced by the d-quark loop, which is easily taken into account by putting the mass of the internal quark in Eq. (27) to zero. In the rest of the diagrams we intend to include O(m s ) and, correspondingly, O(m 2 K ) terms in order to assess the flavor SU (3) F violation in the hadronic matrix elements 4 .
Summarizing, the LCSR for the hadronic matrix element K + K − | O d 2 |D 0 is then obtained from Eq. (26) by the following replacements, as well as by replacing qq → ss , µ π → µ K and f π → f K . In the interest of succinctness we only quote the kaon DA of twist-2, where the Gegenbauer moment a K 1 reflects the SU (3) F -violating asymmetry of thes and u-quark average momentum fractions in the DA. The expressions for kaon twist-3 DAs can be found, e.g., in [18,19]. Apart from the parameters µ K , f 3K , ω 3K , analogous to the pion ones in Eq. (29), these DAs also contain certain O(m s ) corrections and an additional SU (3) F -asymmetry parameter λ 3K .
Note that here, similarly to the LCSR analysis of B → 2π decays [12], we neglect the penguin-annihilation contribution, which is expected to be α s -and power-suppressed. In principle, this contribution can be separately estimated using the same approach. That evaluation is however technically more involved, as it contains multiloop contributions.
In conclusion of this section we emphasize that the pion and kaon DAs of the lowest twist, as well as the interpolating currents of pion, kaon and D-meson all have a valence quark content of the corresponding hadrons, whereas the hadronic matrix elements P s(d) ππ (KK) are obtained from the diagrams where the quark pair in the relevant operator has a flavor different from the valence content. Hence, using the quark-hadron duality ansatz employed in the LCSRs, a hadronic matrix element with penguin topology is unambiguously identified, being "protected" at the level of correlation function from additional quark-antiquark insertions. Indeed, in the OPE such insertions in DAs or in interpolating currents would produce α s -suppressed and/or higher-twist (power suppressed) contributions.

IV. NUMERICAL RESULTS
In order to estimate the size of the computed hadronic matrix elements we need to provide numerical inputs for various parameters used in this calculation. We make conventional choices for the renormalization and factorization scale µ in LCSRs, adopting µ = 1.5 ± 0.5 GeV. The intervals of the (universal) Borel parameter M 2 1 = 1.0 ± 0.5 GeV 2 in the π, K -meson channels, and the corresponding thresholds s π 0 = 0.7 ± 0.1 GeV 2 , s K 0 = 1.2 ± 0.1 GeV 2 are inspired by the analysis of LCSRs for the pion and kaon electromagnetic form factors [20]. The corresponding parameters for D meson channel, M 2 2 = 4.5 ± 1.0 GeV 2 , s D 0 = 7.0 ± 0.5 GeV 2 are taken following the LCSR analysis of D → π, K form factors [18]. We display our choices for the remaining input parameters in Table I. They include the strong coupling, quark masses, quark-condensate densities and the parameters of pion and kaon DAs, all rescaled to the adopted scale. Finally, we use the value of f D = 201 ± 13 MeV for the D-meson decay constant obtained from the 2-point QCD sum rule analysis in [24], and the values f π = 130.5 MeV and f K = 155.6 MeV respectively [6] for the pion and kaon decay constants.
With the chosen input, the results for the hadronic matrix elements calculated from the sum rule in Eq. (26) and from its analogue for the kaon channel are π + π − | Q s 2 |D 0 = (9.50 ± 1.13) × 10 −3 exp[i(−97.5 o ± 11.6)] GeV 3 , where the imaginary parts generated by the quark loops should, in the quark-hadron duality approximation, reproduce the strong phases of these amplitudes 5 .
The values in Eqs. (33) and (34) will be used in the next section to estimate the direct CP-violating asymmetries and their difference in kaon and pion channels.

V. DIRECT CP-VIOLATING ASYMMETRY
Neglecting D 0 D 0 mixing, the partial rates for P = π, K can be written as where p * P is the decay 3-momentum in the D-meson rest frame. In terms of the amplitude parametrization in Eq. (16), the direct part of the CP-asymmetry for the decay of a D-meson into kaons is while the same asymmetry for the decay of a D-meson into pions is a dir CP (π − π + ) = 2r b r π sin δ π sin γ 1 + 2r b cos γ(1 + r π cos δ π ) + r 2 b (1 + 2r π cos δ π + r 2 π ) .
In Eqs. (37) and (38) it was convenient to represent the ratio of CKM parameters as where, in terms of Wolfenstein parameters, the modulus is while the phase γ = arg(ρ + iη) coincides with the angle γ of the unitarity triangle.
Due to the interplay of CKM coefficients in SM, the individual direct CP -violating asymmetries a dir CP (K − K + ) and a dir CP (π − π + ) have opposite signs, as expected. Taking the difference of Eqs. (37) and (38), expanding the result in r b , and keeping only the linear piece we obtain ∆a dir CP = a dir CP (K − K + ) − a dir CP (π − π + ) = −2r b sin γ(r K sin δ K + r π sin δ π ) + O(r 2 b ) . (41) In order to perform a numerical analysis of the above formula we take the central values of  (6). By substituting r b sin γ into Eq. (41) and neglecting the corrections of O(r 2 b ), we obtain the following interval for the penguin contribution parameters, 0.12 ≤ (r K sin δ K + r π sin δ π ) ≤ 1.46 , where we allowed for one standard deviation for the experimental data adding in quadrature the errors quoted in Eq. (6). We see that the combination of amplitudes extracted from experiment is quite uncertain. Still, within 1σ-deviation, quite large values of the relative magnitudes of penguin effects are allowed, especially if the phases are small.
To predict the ratios r K and r π from our calculations we extract the amplitudes |A ππ | and |A KK | relating them via Eq. (16) (in which we neglect the small O(λ b ) terms on the right-hand side) to the decay amplitudes. For determination of the latter, we use the experimentally measured branching fractions [6], B(D 0 → π − π + ) = (1.407 ± 0.025) × 10 −3 , B(D 0 → K − K + ) = (3.97 ± 0.07) × 10 −3 . (43) Inverting Eq. (36), and using the above values together with the lifetime τ D0 = 0.4101 ps, we obtain |A ππ | λ −1 s A(D → π − π + ) = (2.10 ± 0.02) × 10 −6 GeV , Finally, we use the estimated penguin hadronic matrix elements in Eqs. (33) and (34) to predict the ratios: Note that predicting the relative phase δ π or δ K of the total amplitude vs. penguin contribution is a difficult task which is beyond the calculation performed here. Although we obtained a certain prediction for the phase of the penguin contribution, the phase of the main amplitude still remains obscure. One substantial complication is a possible influence of nearby light-quark scalar resonances on the strong phases (see, e.g. [27]). Those resonances, however, are known to be rather broad, overlapping in the energy region of the D-meson mass. This gives us confidence that quark-hadron duality ansatz provides a reasonable approximation to the final result.
Substituting our estimates for r π and r K in Eqs. (37), (38) and (41) taken to O(r b ) and allowing for arbitrary strong phases δ π,K , we obtain the upper bounds, Alternatively, assuming that the phases δ π and δ K are given by the phases of P s ππ and P d KK calculated above and presented in Eq. (32), we obtain a dir CP (π − π + ) = −0.011 ± 0.001%, a dir CP (K − K + ) = 0.009 ± 0.002%.
This is equivalent to assuming that the dominant parts of the decay amplitudes parametrized according to Eq. (16) as λ s A ππ and λ s A KK both have a small phase relative, respectively, to P s ππ and P d KK . Such a situation is realized, for example, in a very simplified scenario when the decay amplitudes are dominated by the factorization ansatz. Yet, this might not be a very reliable approximation, as the decompositions (19) of A ππ and A KK contain hadronic matrix elements with different topologies, including the penguin-topology ones.
Our predictions (46) and (47) are consistent with the experimental results quoted in Eqs. (5), (6), and (7). Note however, that the predicted upper bound on the SM contribution to ∆a dir CP is about factor of five smaller than the magnitude of the central value of the currently available experimental interval (6).

VI. CONCLUSION
In this letter, we presented a new method to estimate the key hadronic matrix elements determining the direct CP asymmetries and their difference in D 0 → K − K + and D 0 → π − π + nonleptonic decays. The method is a variant of the QCD-based LCSR technique adopted from the computations of B → ππ decay amplitudes. A nontrivial strong rescattering phase emerges in the calculated hadronic matrix elements. Our results do not rely on any flavorsymmetry and/or model-inspired amplitude decomposition. They do, however, rely heavily on the assumption of quark-hadron duality, which introduces a yet unaccounted systematic uncertainty.
An interesting question, directly related to the duality violation is the influence of intermediate scalar-isoscalar f 0 resonances on the decay amplitudes. One possibility to address this question is to modify our method in the following way. Instead of directly continuing the calculated correlation function to the physical point P 2 = m 2 D , one can match the result to the hadronic dispersion relation in the variable P 2 at spacelike P 2 < 0, adopting a certain pattern of resonances in this dispersion relation. The hadronic matrix element is then obtained by setting P 2 = m 2 D in the fitted dispersion relation. This, more involved version of the LCSR method is postponed to a future study.
Our main results are the ratios (45) of the calculated "penguin topology" matrix elements to the total D 0 → P − P + decay amplitudes, extracting the absolute values of the latter amplitudes from experimental data on D 0 -decay rates to charged pions and kaons.
The upper bounds (46) and estimates (47) obtained here for the direct CP-asymmetry in both pion and kaon modes and their difference quantitatively assess the expected amount of direct CP violation in the charm sector of Standard Model. We believe that our results will become useful for the interpretation of current and future measurements of this elusive effect.