CP violation in quasi-two-body $D\to VP$ decays and three-body $D$ decays mediated by vector resonances

We re-examine direct CP violation in the quasi-two-body $D\to V\!P$ decays and study CP asymmetries in three-body $D$ decays proceeding through intermediate vector resonances within the framework of the topological amplitude approach for tree amplitudes and the QCD factorization approach for penguin amplitudes. Using the same mechanism of incorporating long-distance penguin-exchange amplitude induced from final-state rescattering that nicely accounts for the CP asymmetry difference between $D^0\to K^+K^-$ and $\pi^+\pi^-$ modes, we find that CP asymmetry can also occur at the per mille level in many of the $D\to V\!P$ channels or otherwise be negligibly small. We point out six golden modes which have sufficiently large branching fractions and direct CP violation of order $10^{-3}$. In particular, the direct CP asymmetry difference $\Delta a_{CP}^{V\!P}$ between $D^0\to K^+K^{*-}$ and $\pi^+\rho^-$ is predicted to be $(-1.61\pm0.33)\times 10^{-3}$, very similar to the counterpart in the $P\!P$ sector. Besides, we find a positive and large CP asymmetry of order $10^{-3}$ in $D^0\to K_S K^{*0}$, an example of the asymmetry induced at the tree level. We take into account the flavor-singlet QCD-penguin contributions to calculate their effects on the modes such as $D\to\pi\phi$ and $D\to \eta \phi$. We compare our results with the other approach in detail. For three-body $D$ decays, we perform the Dalitz plot analysis of those decays mediated by several vector resonances with the same mass but different electric charges. Local CP asymmetry varies in magnitude and sign from region to region and can reach the percent level in certain local regions of phase space due to interference.


I. INTRODUCTION
In 2019 the LHCb Collaboration has announced the measured difference between the timeintegrated CP asymmetries of the decays D 0 → K + K − and D 0 → π + π − based on pion and muon tagged analyses [1]. Combining these with previous LHCb results in 2014 [2] and 2016 [3] leads to a result of a nonzero value of ∆A CP [1] ∆A CP ≡ A CP (K + K − ) − A CP (π + π − ) = (−1.54 ± 0.29) × 10 −3 . (1.1) The time-integrated asymmetry can be further decomposed into a direct CP asymmetry a dir CP and a mixing-induced indirect CP asymmetry characterized by the parameter A Γ which measures the asymmetry between D 0 → f andD 0 → f effective decay widths where t is the average decay time in the sample. To a good approximation, A Γ is independent of the decay mode and hence, where ∆ t is the difference in the mean decay times of the two decays. Based on the LHCb average of A Γ , it follows that the direct CP asymmetry difference is given by [1] ∆a dir CP = (−1.57 ± 0.29) × 10 −3 . (

1.4)
This is the first observation of CP violation at 5.3σ in the charm sector. The most important question to be addressed is whether or not the observed CP violation in the charm sector is consistent with the standard model. This is the main issue discussed in most of recent papers appearing after the 2019 LHCb measurement [4][5][6][7][8][9][10][11][12][13][14][15][16][17]. 1 One may naïvely argue that it is difficult to accommodate the LHCb's measurement of CP asymmetry within the standard model. To see this, consider the tree T and penguin P contributions to D 0 → K + K − and D 0 → π + π − . Their interference gives rise to ∆a dir CP . A simplified expression of the CP asymmetry difference between them is given by 2 where the factor of 1.3 × 10 −3 comes from the imaginary part of a certain combination of CKM matrix elements, 2Im(λ d λ * s )/|λ d | 2 , with λ p ≡ V * cp V up for p = d, s, b, and θ KK is the strong phase of (P/T ) KK and likewise for θ ππ . Since |P/T | is naïvely expected to be of order (α s (µ c )/π) ∼ O(0.1), 1 Various models including many new physics ones (see, for example, references  listed in [9]) had been proposed in the past to explain the large ∆A CP = −(0.82 ± 0.21 ± 0.11)% reported by the LHCb in 2012 [18]. However, this large ∆A CP was not seen anymore in subsequent LHCb analyses [2,3,19]. 2 A complete expression of ∆a dir CP will be given in Eq. (1.9) below.
it appears that ∆a dir CP is most likely of order 10 −4 even if the strong phases are allowed to be close to 90 • . Indeed, using the results of |P/T | obtained from light-cone sum rules [20], Before the LHCb observation of ∆A CP in 2019, there existed two independent studies in 2012 in which direct CP violation in charmed meson decays was explored based on the topological diagram approach for tree amplitudes and the QCD-inspired approach for penguin amplitudes [21][22][23]. Interestingly, both works predicted a ∆A CP of order 0.1% with a negative sign seven years before the LHCb observation of CP asymmetry in the charm sector. It followed that were found in the QCD factorization (QCDF) calculation [9,21,22], while were obtained in the so-called factorization-assisted topological amplitude (FAT) approach [23]. Contrary to the sum-rule calculation, the magnitude of P/T turns out to be of order (0.22 ∼ 0. 30) in the QCD-inspired approaches. However, it does not suffice to explain the measured value of ∆a dir CP . The complete expression of ∆a dir CP is given by [22] ∆a dir CP = −1.30 × 10 −3 P d + PE d + PA d T + E s − ∆P KK sin δ KK + P s + PE s + PA s T + E d + ∆P ππ sin δ ππ , (1.9) where E is the W -exchange amplitude, PE is the QCD-penguin exchange amplitude and PA the QCD-penguin annihilation amplitude, the superscript d or s refers to the quark involved in the associated penguin loop, and δ KK is the strong phase of (P d + PE d + PA d ) KK relative to (T + E s − ∆P ) KK and likewise for the phase δ ππ . The parameter ∆P is defined by ∆P ≡ (P d + PE d + PA d ) − (P s + PE s + PA s ). The observation of the decay D 0 → K 0 K 0 indicates that SU(3) symmetry must be broken in the topological amplitude E. A fit to the data yields two possible solutions [9] I : E d = 1.10 e i15.  (1.12) in FAT [23]. Now we see a crucial difference between the two calculations. The calculated phases δ ππ and δ KK are close to 180 • in QCDF, while they deviate from 90 • by only around 40 • in FAT. As a consequence, Eq. (1.12) leads to ∆a dir CP ≈ −1.0 × 10 −3 , but not so for QCDF, where ∆a dir CP ≈ 6.8 × 10 −5 and −4.9 × 10 −5 for Solutions I and II of W -exchange amplitudes, respectively. One of the most salient features of the topological approach for charm decays is that all the topological amplitudes except the tree amplitude T are dominated by nonfactorizable long-distance effects. That is, final-state interactions play an essential role in describing the hadronic decays of charmed mesons. In [21] we have pointed out the importance of a resonant-like final-state rescattering which has the same topology as the QCD-penguin exchange topological graph PE. Hence, penguin annihilation receives sizable long-distance contributions through final-state interactions. Making the ansatz that (PE) LD is of the same order of magnitude as E and flavor independent, we obtain (1.13) Comparing this with Eq. (1.11), we see that, in the presence of the long-distance penguin exchange, not only the magnitudes of the ratios are enhanced but the strong phases are also closer to 90 • rather than to 180 • . We have predicted that ∆a dir CP is about (−0.139 ± 0.004)% and (−0.151 ± 0.004)% for the two solutions of W -exchange amplitudes [22]. Those were the main predictions among others made in 2012. Nowadays, we know that the LHCb's new measurement (1.1) almost coincides with our second solution. This implies that one does not need new physics at all to understand the first observation of ∆a dir CP by the LHCb. Although both QCDF and FAT approaches lead to ∆a dir CP at the per mille level [21][22][23], they can be discriminated in the D → VP sector. While our approach based on QCDF predicts CP asymmetries of D 0 → π + ρ − and D 0 → K + K * − similar to that of the corresponding D → P P decays [9], the FAT approach leads to very tiny values of CP violation for D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ [24]. In this work, we shall examine the underlying reason for the discrepancy between the two approaches for D → VP decays and present a more detailed analysis. We will also study the implications for CP asymmetry distributions in three-body D decays. For example, the experimental Dalitz plot analysis of D 0 → π + π − π 0 indicates that this channel is predominated by the ρ ± and ρ 0 resonances (see Sec. III B below). The interference between these resonances will yield rich information on local CP asymmetries.
The layout of the present work is organized as follows. In view of updates in several branching fractions, we will perform a re-fit to the Cabibbo-favored modes to extract various topological tree amplitudes in Sec. II. Then we elaborate the calculation of direct CP violation in D → VP decays in detail and compare our results with those of FAT. In Sec. III we discuss CP asymmetry distributions in the following 3-body D decays: D 0 → K + K − π 0 , D 0 → π + π − π 0 , D + → K + K S π 0 and D + s → K 0 π + π 0 . We focus on the phase-space region where these three-body decays are governed by the quasi-two-body D → VP ones. Our conclusions are presented in Sec. IV.

A. Topological tree amplitudes
For D → VP decays, there exist two different types of topological diagrams since the spectator quark of the charmed meson may end up in the pseudoscalar or vector meson. For topological color-allowed tree amplitude T and color-suppressed amplitude C in D → VP decays, the subscript P (V ) implies that the pseudoscalar (vector) meson contains the spectator quark of the charmed meson. For the W -exchange amplitude E and W -annihilation A with the final state q 1q2 , the subscript P (V ) denotes that the pseudoscalar (vector) meson contains the antiquarkq 2 .
The partial decay width of the D meson into a vector and a pseudoscalar mesons is given by with p c being the center-of-mass (c.m.) momentum of either final-state particle. This expression can be simplified by the replacement ǫ ·p D → p c (m D /m V ) after the polarization states of the vector meson are summed over. Hence, with M =M(ǫ · p D ). By performing a χ 2 fit to the Cabibbo-favored (CF) D → VP decays, we have extracted the magnitudes and strong phases of the topological amplitudes T V , C V , E V , A V and T P , C P , E P , A P from the measured partial widths through Eq. (2.2) and found many possible solutions with local χ 2 minima [9,25]. Previously, we found six best χ 2 -fit solutions (S1)-(S6) as listed in Table IV of [9], when we restricted ourselves to χ 2 min < 10. Although Solutions (S1)-(S6) generally fit the CF modes well, they led to very different predictions for some of the singly Cabibbo-suppressed (SCS) decays. Especially, the D 0 → π 0 ω, D + → π + ρ 0 and D + → π + ω decays were very useful in discriminating among different solutions. For example, Solutions (S1), (S2), (S4) and (S5) implied that the branching fraction of D + → π + ρ 0 should be smaller than that of D + → π + ω, in contradiction with the experimental data. As shown in [9], Solution (S6) gave a better description of the SCS modes.
In light of the recent new measurements of D 0 → K * 0 η by Belle [26], [27,28], we have performed a new fit to the data given in Table I and now find five solutions (S1')-(S5'), listed in Table II. It is noted that there is a sign ambiguity if all the strong phases change their sign simultaneously. The topological amplitudes of all these solutions respect the hierarchy pattern: Among them, Solution (S3') is very similar to Solution (S6) and provides a best description of the SCS modes. The topological amplitude sizes and the strong phases for this solution are presented in Table II, where we have chosen the η − η ′ mixing angle φ = 43.5 • [29], which is defined in the flavor basis  [30], Belle [26] and BESIII [27,28]. Here λ sd ≡ V * cs V ud . The column of B theory (S3 ′ ) shows predictions based on Solutions (S3') presented in Table II. All branching fractions are quoted in units of %.
b The new measurement from BESIII [27] is taken into account in the world average. c The new measurement from BESIII [28] is taken into account in the world average.
with η q = 1 √ 2 (uū + dd) and η s = ss. Note that (S3') in Table II here and (S6) in Table IV of [9] are quite similar except for E V and the phase of T P . It is not surprising that the magnitude of E V (S3 ′ ) = 0.58 ± 0.06 is about twice bigger than that of E V (S6) = 0.26 ± 0.04 [9] because the new average of B(D 0 → K * 0 η) = (1.352 ± 0.115)% is increased by 35% from the old one (1.02 ± 0.30)% [30], recalling that its amplitude is proportional to 1 √ 2 (C P + E P ) cos φ − E V sin φ . As a consequence, the phase of E V is also changed significantly from (224 +22 −40 ) • [9] to (283 ± 5) • with a much improved accuracy. Since A(D 0 → K − ρ + ) ∝ (T P + E V ), T P is affected by a change of E V accordingly. 2) and φ = 43.5 • . We will take Solution (S3') in subsequent analyses. The amplitude sizes are quoted in units of 10 −6 (ǫ · p D ) and the strong phases in units of degrees.  In the PP sector, we need SU(3) breaking in the W -exchange diagrams [see Eq. (1.10)] in order to induce the observed D 0 → K S K S decay and to explain the large rate difference between D 0 → K + K − and D 0 → π + π − decays. Likewise, we also need SU(3) breaking in the W -exchange amplitudes in the VP sector because (a) the ratios, Γ(D 0 → K + K * − )/Γ(D 0 → π + ρ − ) = 0.32±0.03 and Γ(D 0 → K − K * + )/Γ(D 0 → π − ρ + ) = 0.45±0.03 [30], deviate substantially from unity expected in the SU(3) limit and (b) the predicted rates of D 0 → K 0 K * 0 and D 0 → K 0 K * 0 modes are too large [9]. Writing we are able to determine the eight unknown parameters e d V , e d P , e s V , e s P and δe d V , δe d P , δe s V , δe s P from the branching fractions of the following eight modes: D 0 → π + ρ − , π − ρ + , π 0 ρ 0 , π 0 ω and D 0 →  solutions for the parameters e d,s V,P and the associated phases (see Table VII of [9]). As we shall show in Sec. II C below, the four different solutions of SU(3) breaking in W -exchange amplitudes can be discriminated using the SCS mode D 0 → ηφ. It turns out that Solution (iv) is preferred.

B. Penguin amplitudes in QCD factorization
Although the topological tree amplitudes T, C, E and A for two-body hadronic D decays can be extracted from the data, information on penguin amplitudes (QCD penguin, penguin annihilation, etc.) needed for inducing CP violation in the SCS decays is still absent. To consider the penguin contributions, we shall take the QCD factorization (QCDF) approach [31,32] to evaluate the hadronic matrix elements. Although the QCD-inspired approaches such as QCDF and pQCD [33] have been applied rather successfully to describing weak hadronic decays of B mesons, it should be kept in mind that these approaches provide only a crude estimate of the penguin amplitudes in charm decays because the charm quark mass is not heavy enough and 1/m c power corrections are so large that a sensible heavy quark expansion is not allowed.
For penguin effects in D → M 1 M 2 decays, we consider the QCD-penguin P , flavor-singlet QCDpenguin S, QCD-penguin exchange PE and QCD-penguin annihilation PA as depicted in Fig. 1 (see [34] for details). 3 In QCDF, the amplitudes of these QCD-penguin diagrams read [32] where p = d, s, p c is the c.m. momentum, and η = 1 (−1) for M 1 M 2 = P V (VP ). We have inserted an additional factor of m V /(m D p c ) in the PE and PA amplitudes so that they are expressed in units of ǫ · p D . In our convention, M 2 is an emitted meson and M 1 shares the same spectator quark with the charmed meson. For the PA diagram with qq produced from the exchanged gluon, M 1 refers to the meson containing the antiquarkq. The factorizable matrix element X reads where we have followed the conventional Bauer-Stech-Wirbel definition for form factors [35]. The chiral factors r M 2 χ in Eq. (2.6) are given by with N c denoting the number of colors, and contributions from vertex corrections V i , hard spectator interactions H i and penguin contractions P i have been taken into account [32]. Note that a i (M 2 M 1 ) = a i (M 1 M 2 ) in general. The annihilation operators b p 3,4 are given by where the annihilation amplitudes A i,f 1,2,3 are defined in [32]. Flavor amplitude decompositions of the SCS D → VP decays are updated in Table IV. Comparing with Table VI in [9], we have made the following improvements: (i) The notation of PE P and PE V in our previous publications [9,21] should be interchanged, i.e. PE P ↔ PE V . 4 (ii) Both QCD-penguin annihilation diagrams PA P and PA V contribute to D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ decays. (iii) Singlet QCD-penguin contributions denoted by S V and S P for SU(3)-flavor singlets η, η ′ , ω and φ are included. (iv) SU(3) breaking in W -exchange amplitudes is taken into account through Eq. (2.5). Values of SU(3)-breaking effects are explicitly shown in the following modes: Calculations based on QCDF give the hierarchical relation for QCD-penguin amplitudes (see e.g., Table VII below) (2.11) Roughly speaking, |PE/P | ∼ O(10 −1 ) and |PA/P | ∼ O(10 −2 ).

C. Branching fractions and SU(3) breaking
In the topological-diagram approach, we apply the topological amplitudes extracted from the CF modes to the SCS decays to predict their rates. It is well known that there exists significant symmetry breaking effects in some of the SCS modes from the flavor SU(3) symmetry limit. In principle, SU(3) breaking effects in the color-allowed T and color-suppressed C tree amplitudes can be estimated provided they are factorizable (in units of ǫ · p D ) (2.12) The topological amplitude sizes of T V , C P , T P and C V given in Table II are extracted from a fit to the data of Cabibbo-allowed D → VP decays. They are dominated by D 0 → K * − π + , K * 0 π 0 , K − ρ + and K 0 ρ 0 , respectively. SU(3) breaking effects in T V,P and C V,P in the SCS modes can be estimated based on factorization. However, as explained in [9], we found that the consideration of SU(3) breaking in T V,P and C V,P alone would render even larger deviations from the data. For example, we They will lead to too small branching fractions for D 0 → π + ρ − and D 0 → π − ρ + and too large B(D 0 → K + K * − ) and B(D 0 → K − K * + ) compared to experiment. Hence, the discrepancy becomes even worse. To circumvent this dilemma, we instead focused on SU(3) breaking in the 4 We see from Fig. 1 that PE P (PE V ) is always accompanied by P P (P V ). In our previous publications, we had, for example, Table II of [21]). Hence, in order to be consistent with the notation in Fig. 1, we need to interchange PE P and PE V . This is done in Table IV.  implies a singlet amplitude with M 1 = V and M 2 = η q (η s ). Values of SU(3) breaking effects in D 0 → π 0 φ and some of D + and D + s decays are shown.

Mode
Representation  Table III). Experimental branching fractions taken from [30] in conjunction with the narrow width approximation are also quoted in units of 10 −3 .
Mode [38] is taken into account in the world average.
b Data from BESIII [38]. c Apparently, the published result of the BESIII measurement of D 0 → ηφ [39] was not taken into account in the world average in the 2020 version of PDG [30]. W -exchange amplitudes as we had done in [9]. By the same token, our predictions on the branching fractions of SCS VP decays of D + and D + s in [9] were made without taking into account SU(3) breaking amplitudes due to the absence of W -exchange contributions.
For the experimental branching fractions of the quasi-two-body D → VP decays listed in Table V, we have applied the narrow width approximation (NWA) to extract B(D → VP ). Notice that this relation is valid only in the narrow width limit, namely, Γ V → 0. Corrections to the NWA have been recently studied in [36,37] for B decays. We shall return to this issue in Sec. III E.
Very recently, BESIII [38] has reported the analysis of D + s → K S π + π 0 decay and obtained the branching fractions of various modes (see Table V). It turns out that our predictions of B(D + s → π + K * 0 ) = (3.65 ± 0.24) × 10 −3 and B(D + s → K 0 ρ + ) = (11.47 ± 0.48) × 10 −3 in [9] are too large compared to the BESIII new measurements listed in Table V. Especially, the latter is larger by a factor of 2. Also, it was already noticed in [9] that B(D + s → K + ω) = (2.12 ± 0.10) × 10 −3 was two to three times bigger than the BESIII result of (0.87 ± 0.25) × 10 −3 [41]. This calls for the necessity of incorporating the SU(3)-breaking effects in the tree amplitudes T V,P and C V,P of the SCS D + s decays.
There was a poorly measured branching fraction in the D + sector, namely, [30]. Thanks to BESIII, a new measurement of this mode with a significantly improved precision is just available [40]. It is based on the amplitude analysis of The new world average predominated by BESIII now becomes Our prediction of (9.80 ± 0.41) × 10 −3 in [9] is too small by about 4σ.
In light of the discrepancy between theory and several new data of VP decays of D + and D + s , it calls for SU(3) breaking in SCS decays of both D + and D + s . It appears that we have a rule of thumb in the absence of W -exchange: it is necessary to consider its SU(3)-breaking effects if only one of the T V,P and C V,P topological amplitudes appears in the decay amplitude. According to this simple rule, we need to account for SU(3) breaking in the color-allowed or color-suppressed tree topological amplitude only in the following SCS modes: 14) The SU(3) breaking effects in SCS D s → VP decays can be estimated from Eq. (2.12) and factorization and where we have assumed that a 1 and a 2 are process-independent. Likewise, for SCS D + → VP decays we have and . (2.18) Using the decay constants and form factors together with their q 2 dependence evaluated in the covariant confining quark model [42], we show the SU(3) breaking effects in the decay amplitudes in Table IV. Now we see from Table V that the agreement with the new BESIII measurements of D + s → π + K * 0 and K 0 ρ + is substantially improved. Likewise, the issue with B(D + s → K + ω) as noticed in [9] is also resolved. However, the predicted B(D + s → K + ρ 0 ) now becomes smaller than experiment. In the topological diagram approach, D + s → K + ρ 0 and D + s → K + ω should have similar rates as their amplitudes are mainly governed by λ d C P / √ 2. Their SU(3) breaking effects are also similar. It is thus desirable to have an improved measurement of B(D + s → K + ρ 0 ). For the D + s → K + φ mode, if we consider its SU(3) breaking in both T V and C P amplitudes, we will have T V (K + φ) = 1.45 T V and C P (K + φ) = 1.20 C P . This will lead to B(D + s → K + φ) = (0.46 ± 0.07) × 10 −3 which is too large compared to the current value of (0.182 ± 0.041) × 10 −3 [30]. This is why we mention in the rule of thumb that if two of the T V,P and C V,P topologies or more appear in the decay amplitude, we should not consider their SU(3)-breaking effects.
In the D + sector, the predicted B(D + → π + φ) = (3.22 ± 0.17) × 10 −3 in the absence of SU(3) breaking is smaller than the experimental value of (5.59 ± 0.10) × 10 −3 by a factor of 1.7. Thus a SU(3) breaking in C P is welcome, though it is not large enough in our calculation. The new prediction of B(D + → K 0 K * + ) = (16.3 ± 0.6) × 10 −3 is now in good agreement with the BESIII measurement. Unfortunately, the original nice agreement between theory and experiment for B(D + → K + K * 0 ) (see Table VIII of [9]) is spoiled when the SU(3) breaking in T V is taken into account. Our result of B(D + → K + K * 0 ) = (5.92 ± 0.18) × 10 −3 is too large compared to the measured value of (3.71 ± 0.16) × 10 −3 . Thus, we face a dilemma that not both D + → K 0 K * + , K + K * 0 modes can be accounted for simultaneously in the same SU(3) breaking scheme.
For this, we need to wait for further experimental justification. We have shown in Table III four sets of solutions for SU(3) breaking in the W -exchange amplitudes extracted from a fit to eight SCS channels. All the solutions are the same for E d V and E d P but differ in E s V and E s P . It turns out that D 0 → ηφ is the only channel which receives contributions from both E s V,P : Table IV). Hence, it can be used to discriminate these solutions. Explicitly, we find (in units of 10 −3 ) B(D 0 → ηφ) = 0.12 ± 0.01, 0.23 ± 0.01, 0.008 ± 0.005, 0.19 ± 0.02 (2.19) for solutions (i), (ii), (iii) and (iv), respectively. It is obvious that the calculated branching fraction based on solution (iv) agrees better with the measured value of 0.18 ± 0.03 than the other solutions. Therefore, we will stick to solution (iv) for SU(3) breaking in the W -exchange amplitudes hereafter.

D. Direct CP violation
It has been noticed that the QCD-penguin exchange diagrams receive sizable long-distance contributions from final-state rescattering [21]. We shall assume that the long-distance PE V and denotes CP asymmetry arising from purely tree amplitudes. The superscript (t+p) denotes tree plus QCD-penguin amplitudes, (t+pe+pa+s) for tree plus PE, PA and S amplitudes, (t+pe LD ) for tree plus long-distance PE amplitude induced from final-state rescattering and "tot" for the total amplitude. As explained in the text, we use solution (iv) for the SU(3) breaking effect in the W -exchange amplitudes (see Table III). The predictions from [24] in the FAT approach with the ρ − ω mixing are listed in the last column for comparison.
Mode a  PE P are of the same order of magnitude as E P and E V , respectively. 5 For concreteness, we follow [9] to assign by choice the same magnitude and phase as the W -exchange amplitudes with 20% and 30 • uncertainties, respectively, so that (in units of 10 −6 (ǫ · p D )) For simplicity, we shall assume its flavor independence, that is, (PE) LD d = (PE) LD s . CP asymmetries of the SCS D → VP decays are updated in Table VI. The improvements over our previous work (see Table IX of [9]) are as follows: (i) In our previous work, the factor of ǫ·p D term in the QCD-penguin amplitudes has been replaced by p c (m D /m V ), while the tree amplitudes are expressed in terms of ǫ · p D . This inconsistency is corrected in this work. (ii) Typos in the code for D + → π 0 ρ + and D + s → K + φ are corrected. The resultant CP asymmetries a dir CP (D + → π 0 ρ + ) = (−0.44±0.52)×10 −3 and a dir CP (D + s → K + φ) = (−1.33 ± 1.59) × 10 −3 are quite different from the previous values of (0.08 ± 0.11) × 10 −3 and 0, respectively. (iii) Singlet QCD-penguin contributions as well as weak penguin annihilation (PE and PA) effects are included under the column denoted by a (t+pe+pa+s) dir . The D 0 → π 0 φ and D + → π + φ decays proceed only through the tree diagram C P . Nevertheless, they receive a small flavor-singlet QCD-pengion contribution S P . Owing to the interference between the tree and singlet QCD-penguin, CP asymmetries no longer vanish, though they are very small of order 10 −6 .
From Table VI we identify several golden modes which have large branching fractions and sizeable CP asymmetries at the order of 10 −3 : It is interesting to notice that the CP asymmetry difference defined by ∆a VP CP ≡ a CP (K + K * − ) − a CP (π + ρ − ), (2.22) in analogy to ∆A CP defined in Eq. (1.1) for the corresponding PP final states, is predicted to be (−1.61 ± 0.33) × 10 −3 , which is very similar to the recently observed CP asymmetry difference between D 0 → K + K − and D 0 → π + π − . This is an attractive and measurable observable in the near future. It is thus desirable to first search for CP violation in the aforementioned golden modes.
Since the quasi-two-body decays D 0 → π + ρ − (D 0 → π − ρ + ) and D 0 → K + K * − (D 0 → K − K * + ) are connected by an interchange of all d and s quarks, there exists a general U -spin relation between corresponding CP -rate differences [43] |A where U denotes a U -spin transformation (d ↔ s) and the overall minus sign on the right-hand side comes from a change in the CKM factors. Hence, Since the CP -averaged rates of D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ are identical in the U -spin limit, we are led to the U -spin relations (2.25) CP asymmetries for the above-mentioned two modes have equal magnitudes and opposite signs. It is evident from  Table VI. It is obvious that CP asymmetries of the following modes are dominated by the long-distance QCD-penguin exchange: The six golden modes listed in Eq. (2.21) all belong to this category. To see this, we consider the four decay modes D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ . From Table IV, their amplitudes read , s, b). Direct CP asymmetries are expressed by where θ πρ is the phase of (P s V + · · · + PA d P ) relative to (T V + E d P + ∆P ) and likewise for the other phases.
In Table VII we show the numerical values of the topological tree and penguin amplitudes of D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ . The magnitudes of T V,P and E d,s V,P are taken from Tables II and III, while QCD-penguins P, PE and PA are calculated in QCDF through Eq. (2.6). It follows that which are consistent with the results shown in Table VI. 6 Obviously, the phase of (P d V + PE d V + PA d V + PA s P ) is almost compensated by that of (T V + E s P − ∆P ), that is, θ K + K * − = 1.2 • . so that the resulting CP asymmetry in D 0 → K + K * − is vanishing small of order 10 −6 . Likewise, θ π + ρ = −8.9 • also implies a very small direct CP violation in D 0 → π + ρ − , of order 10 −5 . For completeness, we give the results for π − ρ + and K − K * + modes a (t+p+pe+pa) dir (π − ρ + ) = 2.1 × 10 −5 , a (t+p+pe+pa) dir After including the long-distance contributions to PE via final-state rescattering (see Eq. (2.20)) we have (2.33) It follows that a dir CP (π + ρ − ) = 0.89 × 10 −3 , a dir CP (K + K * − ) = −1.00 × 10 −3 , (2.34) 6 In the absence of tree contributions to CP violation (i.e. a t dir = 0), a (t+p+pe+pa) dir will be the sum of a  Table VI. Therefore, a  Table VI. 7 Two remarks are in order: (i) In the topological amplitude approach, the magnitude and the relative strong phase of each individual topological tree amplitude in charm decays can be extracted from the data. Consequently, direct CP asymmetries in charmed meson decays induced at the tree level can be reliably estimated. Since the contribution from PA to D 0 → K 0 K * 0 and D 0 → K 0 K * 0 are very small compared to the W -exchange amplitudes, CP asymmetries in these two modes are induced at the tree level, in analogue to D 0 → K S K S . We find a large CP asymmetry at the per mille level for D 0 → K S K * 0 but not for D 0 → K S K * 0 : This is consistent with the upper limit of 0.3% found in [44]. In view of the current efforts in search of CP asymmetry in D 0 → K S K S , it is worthwhile to pursue that in the D 0 → K S K * 0 decay as well. (ii) Although the magnitude of PE and PA is smaller than P (see e.g., in some channels. It depends on the relative phases between T and P as well as between T and PE.

E. Comparison with the FAT approach
From Table VI it is evident that the predicted CP asymmetries given in [24] based on the FAT approach are generally smaller than ours by one to two orders of magnitude. To see the underlying reason, we consider QCD-penguin exchange PE and annihilation PA amplitudes. It follows from Eqs. (2.6) and (2.10) that where the superscripts 'i' and 'f ' refer to gluon emission from the initial and final-state quarks, respectively. The subscript 'k' on A i,f k refer to one of the three possible Dirac structures: k = 1 for [32]. The amplitudes A k have the expressions As stated in [9], predictions in Tables V and VI are Table VI. with (q 1 q 2 ) S±P ≡q 1 (1 ± γ 5 )q 2 and (q 1 q 2 ) V ±A ≡q 1 γ µ (1 ± γ 5 )q 2 . Now A f k corresponds to the factorizable contribution, while A i k to the nonfactorizable contribution of A k . It turns out A f 1 = A f 2 = 0 because of helicity suppression, but not so for A f 3 owing to its (S − P )(S + P ) structure for the four-quark operator; that is, In [24] the factorizable PE f V (P ) amplitude has the expression It was evaluated in the pole model by assuming its dominance by resonant pseudoscalars. Explicitly, it reads 8 where P * represents the pole resonant pseudoscalar meson and H s is the corresponding strong Hamiltonian. Since P V |H s |P * = − VP |H s |P * [24], this leads to the relation A f 3 (P V ) + A f 3 (VP ) = 0. Take D 0 → π + ρ − as an example where P V and PE f V are given by [24] It was claimed in [24] that there was a numerical coincidence that P V and PE f V canceled each other. As a consequence, CP asymmetries in D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ decays are very small of order 10 −5 .
In QCDF, A f 3 is expressed in terms of the twist-2 light-cone distribution amplitude (LCDA) Φ M and the twist-3 one Φ m . A direct evaluation of the weak annihilation diagram with the four-quark operator (ūq) S+P ⊗ (qc) S−P yields [32] in the QCDF approach, contrary to the aforementioned pole model assumption. It is known that the integrals in Eq. (2.42) involve endpoint divergences. We shall follow [31] to model the endpoint divergence X ≡ 1 0 dx/x in the penguin annihilation diagram as (2.43) with Λ h being a typical hadronic scale of 0.5 GeV. The expressions of A f 3 can be further simplified by using the asymptotic distribution amplitudes Φ P,V (x) = 6xx, Φ p (x) = 1 and Φ v (x) = 3(x −x) [32]: (2.44) In principle, one can also add the superscripts 'VP ' and 'P V ' to distinguish penguin annihilation effects in D → VP and D → P V decays. Unfortunately, unlike B decays we do not have any knowledge about the parameters ρ A and φ A for charm decays. Therefore, we shall consider the "default" value X A = ln(m D /Λ h ) in this work. The numerical values of PE and PA shown in Table VII are obtained using the default value of X A . We see from this table that −λ b PE V is smaller than λ d P d + λ s P s and they are of the same sign. Therefore, there is no cancellation between P V and PE V in QCDF! In short, before the consideration of long-distance contribution to PE through final-state rescattering, the predicted CP asymmetries in D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ in QCDF are also very small, of order 10 −5 (K + K * − is further suppressed; see Eqs. (2.31) and (2.32)). However, they are very small in QCDF for a reason quite different from FAT: the phase angles θ πρ , θ ρπ , θ KK * and θ K * K in Eq. (2.28) become smaller or even close to zero after including the contributions from PE, PA and W -exchange to the ratio of P/T . It is the long-distance QCD-penguin exchange that explains why our predictions of CP asymmetries in D 0 → π ± ρ ∓ and D 0 → K ± K * ∓ are much bigger than those in the FAT approach.

III. CP VIOLATION IN THREE-BODY D DECAYS THROUGH VECTOR RES-ONANCES
In this section we would like to study the Dalitz plots of CP asymmetry distributions in some of the SCS D → P 1 P 2 P 3 decays. In particular, we shall focus on the three-body D 0 → K + K − π 0 , D 0 → π + π − π 0 , D + → K + K S π 0 and D + s → K 0 π + π 0 decays as they receive contributions from quasi-two-body D 0 → K ± K * ∓ , D 0 → π ± ρ ∓ , π 0 ρ 0 , D + → K S K * + , K + K * 0 and D + s → π +,0 K * 0,+ decays, respectively. Due to the interference between various vector resonances carrying the same mass but different electric charges, the magnitude and sign of local CP asymmetry vary from region to region. Large CP asymmetries can occur in some localized regions of phase space. For previous studies of CP violation in three-body D decays, see [17,45]. For a study of three-body hadronic D decay amplitudes within the framework of QCDF, the reader is referred to [46]. For the proper three-body approaches using more complex theoretical background in Fadeev techniques, Khuri-Treiman and Triangle singularities, see [47][48][49][50].
In principle, CP asymmetry may also arise from the interference between P -and S-wave contributions. For example, a significant CP asymmetry in the ρ(770) region coming from the interference between the ρ(770) and S-wave resonances has been observed by the LHCb [51,52]. See also the discussions in [46]. Owing to the lack of knowledge on scalar resonances, we will not pursue this direction in this work.
Consider the Dalitz plot analysis of D 0 → K + K − π 0 in the overlapped region of the vector resonances K * + (892) and K * − (892) so that The explicit expression of A K * + reads where M(D 0 → K − K * + ) has been given in Eq. (2.27) and we have considered the relativistic Breit-Wigner line shape for K * (892), and a mass-dependent width with q = | p 1 | = | p 3 | being the c.m. momentum in the rest frame of K * + , q 0 the value of q when s 13 is equal to m 2 K * , and X 1 a Blatt-Weisskopf barrier factor given by with r BW ≈ 4.0 GeV −1 . In Eq. (3.5), Γ 0 K * is the nominal total width of K * with Γ 0 K * = Γ K * (m 2 K * ). Likewise, the amplitude A K * − is given by To derive Eq. (3.3), we have taken the Lagrangian to obtain K + π 0 |K * +µ K − ↔ ∂ µ π 0 |K * + = iǫ * · (p 3 − p 1 ). The form factor F (s, m K * ) in the amplitudes A K * ± is introduced for the following reason. The coupling |g K * ± →K ± π 0 | = 3.15 is extracted from the measured K * (892) width through the relation When K * (892) is off the mass shell, especially when s 13 is approaching the upper bound of (m D − m K ) 2 , it is necessary to account for the off-shell effect. For this purpose, we shall follow [37] to introduce a form factor F (s, m R ) parameterized as with the cutoff Λ not far from the resonance, where the parameter β is expected to be of order unity. We shall use n = 1, Λ QCD = 250 MeV and β = 1.0 ± 0.2 in subsequent calculations. The sum over the polarizations of K * yields the familiar relation However, when the transversality condition ǫ µ k µ = 0 is imposed, the mass in the denominator of the second term should be replaced by the invariant mass m 13 = √ s 13 or m 23 = √ s 23 [53], where , it is straightforward to show that where M =M(ǫ · p D ). Notice that the term in the big parentheses can be recast to where | p 1 | and | p 2 | are the momenta of K + and K − , respectively, in the rest frame of K + and π 0 . Hence, the angular distribution is proportional to the Legendre polynomial P 1 (cos θ). From Eqs. (3.32) and (2.27), we obtain and Note that the branching fractions given above are consistent with that in Table V in conjunction with the narrow width approximation. Likewise, the calculated CP asymmetries are consistent with Table VI. We show in Fig. 2 the Dalitz plot of the CP asymmetry distribution in the overlapped regions of K * (892) + and K * (892) − . Owing to the interference between D 0 → K + K * − → K + K − π 0 and D 0 → K − K * + → K + K − π 0 , local CP violation varies from region to region. For example, a CP is of the percent level in the region around (s K + π 0 , s K − π 0 ) = (0.75, 0.80) GeV 2 (the pink area in Fig. 2) and it becomes negative of order −2.0 × 10 −3 in the vicinity of (s K + π 0 , s K − π 0 ) = (0.70, 0.77) GeV 2 .
Consider the Dalitz plot of D 0 → π + π − π 0 in the overlapped region of the vector resonances ρ(770) ± and ρ(770) 0 so that where we have used the ρππ interaction Lagrangian Numerically, the on-shell coupling is |g ρ→ππ | = 6.0. Hence, Since the ρ(770) resonance is broad, a popular choice for describing its line shape is the Gounaris-Sakurai (GS) model [56] given by . (3.24) In this model, the real part of the pion-pion scattering amplitude with an intermediate ρ exchange calculated from the dispersion relation is taken into account by the f (s) term in the propagator of T GS ρ (s). Unitarity far from the pole mass is thus ensured. Explicitly, and (3.26) The constant parameter D is given by We find and a CP (D 0 → π 0 ρ 0 → π + π − π 0 ) = 0.295 × 10 −3 .
The D + → K + K − π + mode is dominated by the quasi-two-body decays π + φ, K + K * 0 (892) and K + K * 0 (1430) with the fit fractions of order 28%, 26% and 19%, respectively [58]. In this work we find CP asymmetries with the values of −3.7 × 10 −6 and (−0.80 ± 0.24) × 10 −3 , respectively, for π + φ and K + K * 0 (892). In [6] it was claimed that a sizable CP violation could occur in D + → K + K * 0 (1430) within the FAT approach. A search of CP violation in D + → K + K − π + has been carried out by LHCb without any evidence [59]. In short, among the aforementioned four D + → VP modes, analysis of π 0 ρ + and ηρ + is not available yet and the branching fraction of π + ρ 0 is too small. Finally, we turn to the decay D + → K + K 0 π 0 whose amplitude analysis was recently performed by BESIII [40]. It is dominated by the quasi-two-body decays K S K * + and K + K * 0 with the fit fractions of order 57% and 10%, respectively, Following the similar analysis for D 0 → K + K − π 0 , we obtain with M =M(ǫ · p D ). Numerically, we obtain and The CP asymmetry distribution of D + → K + K 0 π 0 in phase space is exhibited in Fig. 4. Regions with positive and negative a CP at the percent level are clearly seen in the plot.
Considering the decay amplitude of D + s → K 0 π + π 0 in the overlapped region of the K * + and K * 0 resonances, we have with where use of the K * Kπ interaction Lagrangian (3.8) has been made and FIG. 5: CP asymmetry distribution of D + s → K 0 π + π 0 in the overlapped regions of K * (892) + and K * (892) 0 .
with the SU(3) breaking effects in T V and C V included (cf . Table IV). Hence,
(3.41) Fig. 5 shows the Dalitz plot of CP asymmetry distribution of D + s → K 0 π + π 0 . Since CP violation in D + s → π 0 K * + and D + s → π + K * 0 is large and positive, CP asymmetry in most of the phase space is positive, though it can become negative in some regions. For example, it is of order of order −0.35 × 10 −3 at (s K 0 π 0 , s K 0 π + ) = (1.97, 1.10) GeV 2 .

E. Finite-width effects
Since all the framework has been set up, we are ready to discuss the finite-width effects. First of all, we need to show that the relation Eq. (2.13) is valid in the narrow width limit. For illustration we take D 0 → π + ρ − → π + (p 1 )π − (p 2 )π 0 (p 3 ) as an example. Since s 12 − s 13 = 4 p 1 · p 2 = 4| p 1 || p 2 | cos θ, (3.42) where | p 1 | and | p 2 | are the momenta of π + and π − , respectively, in the rest frame of π − and π 0 , and withp c being the c.m. momentum of π + (p 1 ) in the D 0 rest frame, then the amplitude A ρ − in Eq. (3.23) can be recast to The decay rate is given by The angular distribution part is given by (see Eq. In the narrow width limit, we have [37] m ρ Γ ρ (s)(1 As a result of the δ-function, s 23 → m 2 ρ . We then obtain the desired NWA: where use of the relations has been made with q 0 and p c being q andp c , respectively, except for the replacement of s 23 by m 2 ρ . Numerically, we have also verified the NWA (3.48). The finite-width effect is accounted for by the quantity η R defined by [36,37] η R ≡ Γ(D → RP 3 → P 1 P 2 P 3 ) Γ R →0 Γ(D → RP 3 → P 1 P 2 P 3 ) = Γ(D → RP 3 )B(R → P 1 P 2 ) Γ(D → RP 3 → P 1 P 2 P 3 ) = 1 + δ , (3.50) so that the deviation of η R from unity measures the degree of departure from the NWA when the resonance width is finite. It is naïvely expected that the correction δ will be of order Γ R /m R . The calculated η V parameters for vector resonances ρ and K * produced in the three-body D decays VIII: A summary of the η R parameter for vector resonances produced in the three-body D decays. The parameter η ρ + extracted from D 0 → π − ρ + → π + π − π 0 is the same as η ρ − and likewise for η ρ 0 . K * (892) D + s → K * + π 0 → K 0 π + π 0 47.3 ± 0.5 0.053 1.099 (BW) are summarized in Table VIII. Since ρ is three times broader than that of K * (892), it is naïvely expected that the deviation of η ρ from unity should be larger than that of η K * . However, our calculation shows that they are close to each other. Our results are to be compared with η ρ = 0.93 in the GS line shape model and 1.11 in the BW line shape for B + → ρ 0 π + → π + π − π + , and η K * = 1.067 ± 0.002 for B + → K * 0 π + → K + π + π − [36,37]. When the resonance is sufficiently broad, it is necessary to take into account the finite-width effects characterized by the parameter η R . Explicitly, where B(D → RP ) NWA denotes the branching fraction obtained from Eq. (2.13) valid in the NWA. Take D 0 → πρ as an example. Their branching fractions have been extracted from D 0 → π + π − π 0 by BaBar using the Breit-Wigner line shape [55]. According to our calculation, η BW ρ = 1.116. Hence, the PDG values of (5.15 ± 0.25), (10.1 ± 0.4) and (3.86 ± 0.23) (in units of 10 −3 ) [30] for the branching fractions of D 0 → π + ρ − , π − ρ + and π 0 ρ 0 , respectively, should be corrected to (5.75 ± 0.28), (11.3 ± 0.5) and (4.31 ± 0.26).

IV. CONCLUSIONS
In this work we have re-examined direct CP violation in the quasi-two-body D → VP decays and studied CP asymmetries in three-body D decays proceeding through intermediate vector resonances within the framework of topological amplitude approach for tree amplitudes and the QCD factorization approach for penguin amplitudes. As we have pointed out in 2012, the longdistance penguin-exchange through final-state rescattering gives the major direct CP violation to both D 0 → K + K − and D 0 → π + π − . It accounts for nicely the CP asymmetry difference between the aforementioned two modes observed by the LHCb in 2019. The same mechanism implies that CP asymmetry can occur at the per mille level in many of the D → VP channels.
Our main results are: • In light of several new measurements of Cabibbo-favored modes, we have performed a re-fit to the data to extract topological tree amplitudes. The topological W -exchange E V has a size twice larger than the old value and its phase is significantly different from the old one.
• In the U -spin limit, there are the relations a dir CP (D 0 → K ± K * ∓ ) = −a dir CP (D 0 → π ± ρ ∓ ). We have found that the U -spin relation is approximately respected for D 0 → K + K * − but not for D 0 → K − K * + .
• Analogous to the PP sector, SU(3) breaking in E V,P is obtained from a fit to eight singly Cabibbo-suppressed D 0 decays. There are four different solutions, which can be discriminated by the decay mode D 0 → ηφ to favor solution (iv). An important consequence is that D 0 → K S K * 0 is predicted to have a positive CP asymmetry at the per mille level. It is an example that the asymmetry is induced at the tree level.
• The consideration of SU(3)-breaking effects in the VP sector is more subtle than the PP case. It appears that SU(3) breaking can only occur in one of the topological tree amplitudes in each decay channel. We are able to account for the recent new measurements from BESIII except for D + → K + K * 0 . For this, we need to wait for further experimental justification.
• We have included the flavor-singlet QCD-penguin contributions, which were missing in our previous studies, to calculate their effects on modes such as D → πφ and D → ηφ. As shown in Table VI, the inclusion of the singlet penguins in these decay modes results in nonzero a dir , albeit not to a level that can be easily measured experimentally.
• We compare our approach with the FAT approach in detail. The main difference lies in the treatment of the factorizable part of the penguin-exchange amplitude PE characterized by A f 3 and the consideration of long-distance contribution to PE through final-state rescattering. In QCDF, A f 3 is evaluated in terms of twist-2 and -3 light-cone distribution amplitudes and it does not have near cancellation with PE as claimed by the FAT analysis.
• For three-body D decays, we show the Dalitz plot of CP asymmetry distribution of D 0 → K + K − π 0 in the overlapped regions of K * + and K * − resonances and likewise for D 0 → π + π − π 0 , D + → K + K S π 0 and D + s → K 0 π + π 0 . Regional asymmetry varies in magnitude and sign from region to region and can reach the percent level in certain invariant mass regions.