$D \to \rho \,\ell^+\ell^-$ Decays in the QCD Factorization Approach

We consider rare semileptonic decays of a heavy $D$-meson into a light vector meson in the framework of QCD factorization. In contrast to the corresponding $B$-meson decays, the naive factorization hypothesis does not even serve as a first approximation. Rather, the decay amplitudes appear to be dominated by non-factorizable dynamics, e.g. through annihilation topologies, which are particularly sensitive to long-distance hadronic contributions. We therefore pay particular attention to intermediate vector-meson resonances appearing in quark-loop and annihilation topologies. Compared to the analogous $B$-meson decays, we identify a number of effects that result in very large theoretical uncertainties for differential decay rates. Some of these effects are found to cancel in the ratio of partially integrated decay rates for transversely and longitudinally polarized $\rho$ mesons. On the phenomenological side this implies a very limited potential to constrain physics beyond the Standard Model by means of these decays.


Introduction
In view of the persistent non-observation of any direct signals for new particles or interactions at the high-energy frontier, which is currently explored at the Large Hadron Collider (LHC), indirect probes of physics beyond the Standard Model (BSM) from low-energy observables gain ever more importance. In particular, depending on the specific model, precision measurements of rare flavour decays together with reliable theoretical predictions constrain the BSM parameter space for masses and coupling constants (for reviews and further references, see e.g. the according sections in [1][2][3]). This is particularly true for rare decays induced by flavour-changing neutral currents in the down-quark sector, i.e. b → s, d and s → d transitions, as well as B-B and K-K-mixing. On the other hand, rare c → u transitions are known to be plagued by serious theoretical uncertainties related to longdistance hadronic effects which are prominent because, due to the small Yukawa coupling y b y t , the GIM cancellation is more efficient for b-quarks in the loops than for top quarks. Nevertheless, there have been a number of phenomenological studies on the new-physics (NP) sensitivities of rare charm decays, including radiative and rare semileptonic D-meson decays induced by c → uγ and c → u + − transitions, see e.g. [4][5][6][7][8][9][10][11][12][13][14][15][16].
On the theoretical side, as a first approximation, one may employ the naive factorization hypothesis which expresses the decay amplitudes in terms of perturbatively calculable Wilson coefficients for c → uγ and c → u + − transitions, multiplied by hadronic form factors for D → π or D → ρ transitions. A simple model to estimate the non-factorizable long-distance effects is to assume vector-meson dominance, i.e. to describe the radiative decays via D → π(ρ) V (→ γ/ + − ) with suitable vector mesons V that couple to the corresponding hadronic current in the weak effective Hamiltonian and decay into a charged lepton pair via electromagnetic interactions. 1 In such an approach, however, the separation of short-and long-distance dynamics is no longer manifest.
To proceed, the systematic inclusion of strong-interaction effects in the relevant hadronic amplitudes requires additional approximations. In particular, one may consider an expansion in inverse powers of the charm-quark mass, which -however -is expected to work less effectively than in the corresponding B-meson decays because m c < m b . In such an approach, short-distance corrections in Quantum Chromodynamics (QCD) related to distances ∆x ≤ 1/m c will be calculated perturbatively. Radiative corrections between the electroweak scale and the charm-quark mass will be included in Wilson coefficients of the effective Hamiltonian for |∆C| = |∆U | = 1 transitions [17]. Recent next-to-leading order calculations for the relevant Wilson coefficients can be found in [18].
It is known from the analogous radiative and semileptonic b → s, d decays that corrections to the hadronic matrix elements from higher orders in a simultaneous expansion in the strong coupling and the inverse heavy-quark mass lead to sizeable effects, in particular in the region of large hadronic recoil (i.e. small invariant lepton mass q 2 ) [19][20][21][22][23][24][25][26][27]. Furthermore it seems a rather non-trivial task to construct realistic models for the effect of intermediate vector-meson resonances contributing to the q 2 spectrum in B → K ( * ) + − decays above and below the cc threshold, see e.g. the discussions in [28][29][30][31][32]. 2 Of course, we expect these issues to be even more pronounced in radiative and rare semileptonic D-meson decays. On the one hand, this severely limits the sensitivity to generic NP scenarios. On the other hand, the investigation of exclusive c → uγ * transitions may help to better understand the hadronic uncertainties in exclusive b → s(d)γ * decays.
The aim of this paper therefore is to critically assess the theoretical control on D → ρ + − (and also D → π + − ) decays within the framework of QCD factorization (QCDF) [34,35], closely following the analyses for the analogous B-meson decays in [19,24]. As a new ingredient we propose a simplified approach which allows to estimate the potential effect of light vector resonances on the level of the individual decay topologies that appear in QCDF (for simplicity, we restrict ourselves to the leading-order expressions in the strong coupling). To this end we model the tower of vector-meson resonances as in [36][37][38] and connect it to the QCDF expressions via dispersion relations such that the asymptotic behaviour (i.e. far away from the resonances) of the perturbative result is reproduced. Clearly, in this way we ignore additional (non-factorizable) hadronic rescattering effects that would potentially lead to a more complicated decay spectrum. As our main goal is not to get a fully realistic description of the differential decay width, our simple procedure should be sufficient to get a rough estimate on the associated hadronic uncertainties.
Our paper is organized as follows. In the next section, we give a brief overview over the theoretical framework, specifying the operator basis in the weak effective Hamiltonian and providing the definitions and factorization formulas for the generalized form factors and coefficient functions appearing in the QCDF approach. In the following Section 3 we provide detailed formulas for the contributions from the different decay topologies within QCDF. Here, we remind the reader that in the "naive factorization approximation", only contributions from the electromagnetic dipole operator O 7 and the semi-leptonic operators O 9,10 appear. Non-trivial contributions from the hadronic operators can be obtained by either closing two quark lines to a loop, or annihilating/pair-creating the valence quarks in the initial-and final-state mesons. Radiative corrections at first order of the strong coupling α s are included by adapting the results in [19,24] to the corresponding c → u transitions; in particular, this includes non-factorizable spectator scattering effects. Notice that the α s corrections to the annihilation topologies are presently unknown. In Section 4 we provide some numerical results for the individual contributions to the coefficient functions describing the decay amplitudes for transversely and longitudinally polarized ρ mesons for neutral and charged decay modes, respectively. On that basis we give numerical estimates for the central values and the dominating theoretical uncertainties for the partially integrated transverse and longitudinal decay rates, as well as for their ratio. We summarize and conclude in Section 5. In Appendix A we provide the explicit formulas that allow to deduce the decay amplitudes for D → π + − decays from the corresponding expressions with longitudinally polarized ρ mesons. Appendix B gives a detailed derivation of the hadronic model that we have used to estimate the effect of light vector resonances. Finally, Appendix C specifies a number of input parameters that have been used in the numerical analysis.

Theoretical Framework
In the following section, we briefly summarize the theoretical framework and fix the notation used for the computation of the decay amplitudes.

Effective Hamiltonian
The low-energy effective Hamiltonian for c → u transitions is written as follows, and 3) The operators are defined in the CMM basis [39]. The current-current operators read where q = d, s. The strong penguin operators are written as (2.5) The electro-and chromomagnetic penguin operators are given by and, finally, the semi-leptonic operators are chosen as The higher-order QCD corrections to the various short-distance Wilson coefficients C i have recently been calculated in [18]. For convenience, we have summarized in Table 1 the SM values for the Wilson coefficients at LL and NLL (NNLL for C 9 ) at a reference scale µ = µ c = 1.5 GeV. Notice that, contrary to B-meson decays, the Wilson coefficient C 10 vanishes for c → u transitions because of the perfect GIM cancellation at the electroweak matching scale (using m 2 b /M 2 W → 0).

Generalized form factors
As has been shown in [19], the hadronic matrix elements of the weak effective Hamiltonian simplify in the large-recoil limit, 3 where In the following, we will adopt the notation used in [24,40] (see also references therein). For radiative decays of a D-meson into light vector mesons, we express the hadronic transition matrix elements in terms of generalized form factors T Here for each set of operators (i = b, d) only two independent generalized form factors appear, referring to transversely or longitudinally polarized vector mesons (a =⊥, ), respectively. They can be further factorized in the form In the first term, the functions ξ a denote the universal ("soft") form factors for D → ρ transitions in the large-recoil limit [41], and the coefficient functions C (i) a contain the QCD corrections to the partonic c → u + − amplitude, dubbed "form factor corrections" in the following. The second term contains the spectator effects (including annihilation topologies) and is proportional to the meson decay constants (f D and f ρ,⊥( ) ). It contains a convolution integral of the relevant light-cone distribution amplitudes (LCDAs), φ D,± (ω) and φ ⊥, (u), and a hard-scattering kernel denoted as T (i) a,± (u, ω). It is to be stressed already at this point that the spectator interactions involve typical gluon virtualities of order m D Λ QCD ∼ 1 GeV, and therefore the corresponding value of the strong coupling constant α s and the associated perturbative uncertainty from scale variations will be large. Furthermore, we define the following coefficient functions, which are independent of the conventions chosen to renormalize the weak effective Hamiltonian, In terms of these, we obtain the double-differential decay rate as [24] with S = 1 for ρ − and S = 1/2 for ρ 0 . Here is the standard kinematic prefactor. Note, that we have used C 10 = 0, which also implies that the forward-backward asymmetry with respect to the angle θ vanishes. 4 For decays into light pseudoscalar mesons, analogous functions T (i) p can be defined [19,42]. The explicit formulas will be provided in Appendix A for completeness.

Naive factorization
In naive factorization, only the Wilson coefficients associated to the operators O 7,9,10 contribute, multiplied by the corresponding vector, axial-vector and tensor form factors for D → ρ transitions. In the large recoil limit, one could further reduce these form factors to the two "soft" form factors, appearing in (2.8). The normalization of the form factors at zero momentum transfer has been measured experimentally by the CLEO collaboration [43], as summarized in Table 2 below. For the q 2 dependence, we adopt the parametrization that has been proposed in [44], where , (3.2) with x = q 2 /M 2 D * and a = 0.55, b = 0.69. The q 2 -dependence of the soft form factors is shown in Fig. 1 and compared to the naive scaling behaviour which follows in the framework of SCET/QCDF at tree level [41]. As one can observe, the latter systematically leads to a somewhat too steep q 2 -dependence. We will therefore stick to the parametrization (3.2) in our numerical analysis below.

Corrections to the D → ρ form factors in the large recoil limit
Adapting the terminology of [19,24], factorizable form-factor corrections are accounted for by [24] C whereas C (f,d) and ∆M depends on the convention to define the charm-quark mass. In this work, we choose the MS scheme for simplicity, which amounts to setting ∆M = 0. With our convention for defining the "soft" D → ρ form factors, the factorizable spectator corrections to the form factors are reflected in the contributions , a,± = 0 in the large-recoil limit.

Quark-loop topologies (w/o spectator effects)
Closing two quark lines from the 4-quark operators and radiating a (virtual) photon from the quark loop results in form-factor corrections to naive factorization that contribute to the effective Wilson coefficients in (2.10) as follows, Here, the 1-loop functions Y (i) (s) can be decomposed as follows [18], and where the function h(s, m) can be found, for instance, in [45], and -with our normalization convention 5 -is given in (B.15) in the appendix. Notice that the contribution of the function Y (b) (s) to the decay amplitudes is strongly CKM suppressed. On the other hand, the function Y (d) (s) vanishes in the limit m s → m u,d 0. Furthermore, at scales of the order of the charm-quark mass one encounters numerical cancellations in the combination of Wilson coefficients (4/3 C 1 + C 2 ); the effect is illustrated in Fig. 2. In particular, one observes a large shift when going from the LL to the NLL result, indicating that this combination of Wilson coefficients is particularly affected by higher-order radiative corrections, see also the discussion around Fig. 4.
At this point it is to be stressed that, strictly speaking, the perturbative calculation of the loop-function h(s, m) is only justified in the deep Euclidean region, where −s = Q 2 Λ 2 QCD . In contrast, for time-like momentum transfer, s = q 2 > 0, the q 2 -dependence would rather be given by a hadronic spectral function that describes the effects of a tower of light vector mesons and multi-hadron final states with the appropriate quantum numbers. While in the corresponding B-meson decays, the annihilation topologies only provide a 5 Our definition of h(s, m) differs from the one in [18] by a relative factor (−1/2). We also remark that the definition of the constant piece in the function Y (b) (s) depends on the operator basis, and only the sum C9 + Y (b) (s) is basis-independent. The quoted results for C9 in Table 1 and Y (b) (s) in (3.8) refer to the CMM basis [39]. correction to the total decay amplitude, this is no longer the case for D → ρ + − decays, see also earlier estimates in Refs. [11,46]. In addition, for B-meson decays the region of small momentum transfer, 4m 2 ≤ q 2 1-2 GeV 2 , could simply be excluded from the phenomenological analysis. Due to the restricted phase space this is no longer possible in D-meson decays. For this reason one should carefully discuss the associated parametric and systematic hadronic uncertainties. Inevitably, this requires some hadronic modelling. In particular, the amount of GIM cancellation in the difference h(s, m s ) − h(s, m d ) in the perturbative calculation will be quite different from estimates based on hadronic models.
A straightforward approach, which has already been extensively used in the past [7,9,12,13,15], assumes that the long-distance hadronic effects are completely dominated by the lowest-lying narrow vector states (ρ, ω, φ etc.) which are then modelled by Breit-Wigner resonances. In this work, we propose a more involved (but still oversimplified) picture, following [37,38], 6 where (i) a model for the infinite tower of higher resonances is included, and (ii) the shape of the Breit-Wigner resonances for the lowest states is modified to be in accordance with simple (leading-order) analyticity arguments. As explained in more detail in Appendix B, our idea amounts to replacing where the precise form of the spectral function j q (s), which models the effect of hadronic vector resonances with the corresponding quantum numbers, 7 and the meaning of the parameter σ 2 can be found in Appendix B. In Fig 3 we illustrate the effect of the resonance in the low-q 2 region, relevant for D → ρ decays. One observes that neither the sign nor the order of magnitude nor the shape of the hadronic model can be reproduced by the perturbative result (despite the fact that - 6 This model has also been used to estimate the effects of higher charmonium resonances in B → K + − decays at large values of q 2 [30]. 7 For simplicity, we do not distinguish between ρ and ω mesons and set mu = m d = 0. Moreover, as by construction -for q 2 < 0 one has almost perfect numerical agreement, see Figs. 9,10 in the appendix).
It is important to note that the numerical effect of the vector resonances in rare semileptonic D-decays will be quite different from B-decays or e + e − annihilation. While in the latter case, the cross section for e + e − → hadrons will be given by the imaginary part of the loop function, in B-decays the dominating correction to the decay width will arise from the interference with the short-distance terms (which do not yield strong phase differences), and therefore one will probe the real part of the loop function (see e.g. [30]). In D-meson decays, however, the short-distance part is heavily suppressed, and therefore we expect the dominating effect to be given by the absolute value of the loop function, see

Non-factorizable c → u form-factor corrections
The remaining non-factorizable contributions can be expressed in terms of q 2 -dependent short-distance functions F (j) i , which enter as follows (see [19,24]) At this point, a few comments are in order on how to obtain the functions F   can be obtained by changing the corresponding charge factors (which amounts to multiplying by Q u /Q d , where Q u,d are the charge factors for up-and down-quarks, respectively).
• In order to get the functions F 2,d and F (9) 2,d , one has to go back to the bare (i.e. unrenormalized) functions from b-decays [40], replace the charge factors and renormalize the functions. Here, the renormalization constant with general charge factors can be derived from the results given in [47] and reads • The results for the functions F 1 , F 1 , F 2 and F 2 with general charge factors have been reconstructed from Mathematica notebooks that have been kindly provided by Christoph Greub (related to the work in [48]).
• For the limit q 2 → 0, the relevant functions F can be directly extracted from [49].
As explained in Section 3.2, the leading contributions from quark-loop topologies in D → ρ transitions suffer from a renormalization-scale uncertainty due to partial numerical cancellations in the combination of Wilson coefficients ( 4 3 C 1 + C 2 ). The higher-order corrections encoded in (3.11) are expected to reduce this uncertainty. To numerically investigate the effect, we plot in Fig. 4 the contribution to the partially integrated branching ratio, 0.5 GeV 2 ≤ q 2 ≤ 0.7 GeV 2 , of the quark-loop topologies as a function of the renormalization scale µ (where we have restricted ourselves to the CKM-favoured contributions proportional to λ d ). We observe that the NLO corrections actually dominate over the LO contribution. On the one hand, this removes the issue with the accidental numerical cancellations at LO. On the other hand, it implies a very slow convergence of the perturbative series, i.e. large renormalization-scale uncertainties even at NLO. (We remind the reader that annihilation and spectator-scattering topologies will induce additional and formally independent scale uncertainties.)

Annihilation
As already mentioned, the so-called annihilation topologies will turn out to give large contributions to the decay rate, and therefore the associated hadronic uncertainties will be essential for phenomenological studies. The fact that the quark propagators in the annihilation diagrams involve time-like virtualities implies a particular sensitivity to the modelling of hadronic resonance effects. Furthermore, annihilation diagrams with the photon radiated from the quarks in the final-state meson are formally power-suppressed in the 1/m c expansion, but may still be phenomenologically important. In particular, radiative corrections to these topologies have not been systematically computed in QCDF so far. As a consequence, the ambiguities related to the renormalization-scale setting in the relevant Wilson coefficients will remain a major source of theoretical uncertainties. In the heavy-quark limit, the leading contributions from quark annihilation topologies originate from photon radiation off the light-quark in the D-meson. The results thus depends on the charge factor e q of the spectator quark. Translating the results from b → s, d to c → u transitions, we obtain Again, the contributions from T (0,b) ,− will be strongly CKM suppressed. Therefore, for charged D-meson decays (q = d) the annihilation contribution is triggered by the Wilson coefficient C 2 , while for neutral D-meson decays (q = u) it is proportional to the same combination 4C 1 3 + C 2 as appearing in the quark-loop topologies. As a consequence, we expect that charged D-meson decays will actually be dominated by annihilation topologies, since the partial numerical cancellation in Wilson coefficients does not occur here. In contrast, in neutral D-meson decays annihilation and quark-loop topologies will enter with similar magnitudes.
The convolution of the expressions in (3.13) with the D-meson LCDA leads to the q 2 -dependent "moment" 14) whose analytic properties are further discussed in appendix B.2. Notice that the limit q 2 → 0 does not exist in (3.14) which limits the applicability of QCDF for that part of the amplitude to hard-collinear values of the momentum transfer that formally scale as Moreover, as we have already discussed for the quark-loop contributions to the form-factor-like terms in the previous subsection, the physical q 2 spectrum in that region will be significantly influenced by light vector-meson resonances and looks quite different from the partonic result following from (3.14). Applying the same kind of model for the hadronic spectrum, as explained in appendix B.2, we end up with an estimate for the hadronic effects in λ − D (q 2 ). Here, the q 2 -spectrum (given by the imaginary part of (λ − D ) −1 ) is assumed to factorize into the D-meson LCDA and a hadronic model for the spectrum associated to the light vector current, The parameters in the function j q (s) are adjusted to reproduce the perturbative result in the limit q 2 −m D Λ. To illustrate the numerical effect, we plot in Fig. 5 the absolute value of λ − D (q 2 ) −1 (the real and imaginary part are plotted in Fig. 11 in the appendix).
Here we have taken a simple exponential model [50], for the LCDAs φ ± D (ω). The resulting picture looks quite similar as for the quark-loop functions h(s, m). In particular the asymptotic behaviour for large values of |q 2 | is unchanged (by construction), while the typical modifications from the resonances occur in the region below 5 GeV 2 .

Power corrections of order 1/m D
The contribution of annihilation topologies to the amplitudes for transversely polarized vector mesons only start at relative order 1/m D . It is known from the analysis of the analogous B-meson decays that these terms should not be neglected, in particular for observables that are sensitive to the transverse decay amplitudes or isospin asymmetries [20][21][22][23][24]. Translating again the known results from the B-meson sector, we end up with the following expressions, (3.17) Here, the first term arises from photon radiation off the constituents of the ρ-meson 8 and therefore contains a non-trivial convolution with respect to the corresponding momentum fractions u andū = 1 − u. The other terms stem from sub-leading contributions from photon radiation off the D-meson constituents, which now give rise to the q 2 -dependent "moment" As before, our hadronic model for the vector resonances can be implemented by the replacement The numerical effect is illustrated in Fig. 6 (see also Fig. 12 in the appendix). Notice that at q 2 = 0 the function λ + D (0) describes the "soft" contribution to the D → γ transition form factor. The hadronic effects in our simple-minded model lead to a reduction of the form factor by 20 − 30% which is in qualitative agreement with the findings in [51] for the B → γ form factor.

Non-factorizable spectator scattering
Non-factorizable spectator scattering effects arise from matrix elements of the hadronic operators O 1−6,8 where the (would-be) spectators in the D → ρ transitions take part in the short-distance scattering process. The various contributions that appear to order α s , again, can be adapted from [19,24] with the appropriate modifications for c → u transitions. 9 We find and together with and Here the quark-loop functions describing the relevant sub-diagrams are given by while the functions t ⊥, (u, m q ) and the definition of C eff 8 can be found in [19]. Notice that the contribution to the CKM-favoured amplitudes for transversely polarized vector mesons, T (nf ),d ⊥,+ in (3.23), are again GIM-suppressed. We therefore also include power corrections of relative order 1/m D to the transverse amplitudes, which again can be adapted from the corresponding expressions for B-meson decays. We find (3.27)

Numerical Results
In this section we present some numerical estimates following from our theoretical analysis in the previous section. The theoretical predictions for the differential decay rates depend on a number of parameters. The most important hadronic input parameters are listed in Table 2. The remaining input parameters that have been used for the numerical analysis are listed in Table 11 in the appendix for completeness. Here a comment is in order about our choice for the parameter ω 0 which determines the average value for the light-cone momentum of the spectator quark in the D-meson. While the analogous parameter for B-meson decays has been studied in some detail in the past (see e.g. [51,52] and references therein), there is practically no theoretical or phenomenological information on the Dmeson LCDAs. We have therefore used an ad-hoc range for ω 0 which reflects the naive expectation from heavy-quark symmetry with a sufficiently conservative uncertainty. a 2 (ρ) ⊥, 0.15 ± 0.07 [54] f D (209 ± 3) MeV [55] ω 0 (450 ± 150 MeV (ad-hoc) 4.1 Detailed breakdown of contributions to C 9,⊥ and C 9, We first summarize the individual contributions to the coefficient functions C 9,⊥ (q 2 ) and C 9, (q 2 ) as defined in Eq. (2.10) at NLO, 10 at a benchmark value q 2 = 0.5 GeV 2 for the momentum transfer, see Tables 3 and 4. Compared to the situation in the analogous Bmeson decays (see the discussion in [24]), one observes a number of important differences: • As is well known, the purely short-distance contribution from the Wilson coefficient C 9 is heavily suppressed by two effects: (i) the GIM cancellation between downtype quarks running in the loops, leading to a factor-10 smaller value for C 9 in D-decays compared to B-decays; (ii) the CKM suppression for c → u transitions, reflected by λ b λ d . The decay amplitudes are therefore dominated by long-distance contributions proportional to the CKM structure λ d .
• Among these, it turns out that for the neutral decay mode -at least at the considered value of q 2 -non-factorizable form-factor corrections (FFnf) and annihilation topologies (Ann) enter with the same order of magnitude, and also non-factorizable spectator effects (Specnf) give a non-negligible contribution. In the case of transverse vector mesons this requires to include (numerically unsuppressed) power corrections to annihilation spectator topologies (1/M Ann) as well.
• The charged decay modes are completely dominated by contributions from annihilation topologies, where again the case of transverse vector mesons requires to include terms that are formally suppressed by 1/M D .
• We stress again that in each case, the estimates for the relevant long-distance contributions within the QCDF approach are very sensitive to hadronic resonance effects which adds to the systematic theoretical uncertainties.

Decay
Contr.  Table 3. Breakdown of individual contributions to coefficient functions C 9,⊥ and C 9, at NLO for the neutral decay mode, D 0 → ρ 0 + − at q 2 = 0.5 GeV 2 . The renormalization scale is set to µ = 1.5 GeV. All other input parameters are set to their default values. Here C 9 denotes the purely short-distance contribution. We further list factorizable (FFf) and non-factorizable (FFnf) formfactor corrections, annihilation (Ann) at leading power, factorizable (Specf) and non-factorizable (Specnf) spectator interactions at leading power, as well as the included power corrections (1/M ) to annihilation and spectator topologies.

Differential decay rates
We next turn to the differential decay rates, where we distinguish between the contributions of transverse and longitudinal ρ mesons, dΓ T,L , which are obtained by projection onto the terms with (1 ± cos 2 θ) in Eq. (2.11). In Figs. 7 and 8 we show our result for the differential branching fractions in the case of neutral and charged mesons, respectively. Here, we compare the LO and NLO results, displaying only the uncertainties from scale variation for simplicity. The following comments can be made: • The difference between the central values for the LO and NLO predictions (for our default choice of µ) is not very pronounced.
• Still, at least for the neutral decay mode, we observe a large renormalization-scale dependence. This can be traced back to the issues with the combination of Wilson coefficients (4/3C 1 + C 2 ) (see the discussion around Fig. 2) appearing in the relevant annihilation contributions in (3.13) and (3.17). Such cancellations do not occur in the charged decay mode.
• One should be aware that (presently unknown) NLO corrections to annihilation are not included. This particularly concerns the charged decay modes which are completely dominated by annihilation topologies.
• The hadronic-resonance model mainly leads to an enhancement of the differential rates. This is related to the fact that the rates are sensitive to the absolute values of the modelled complex functions, which has already been pointed out above. As a consequence, in contrast to the well-known R-ratio, the pseudo-realistic q 2 -spectrum does not show the naively expected oscillations around the perturbative result.
• There is, however, a small region around q 2 0.75 GeV 2 , where the perturbative result does not seem to be very much affected by resonance effects, at least within our simplified hadronic model. However, in that region the differential branching fractions are small, and the relative uncertainties are large.
Apart from the renormalization-scale dependence, the main parametric uncertainties stem from the LCDA of the D-meson, while the modelling of the hadronic resonances yields an estimate for the systematic hadronic uncertainties. These are summarized in Tables 5-8, where we display our predictions for partially integrated branching ratios in specific q 2 bins. As one can see, the uncertainties in the chosen q 2 bins can easily exceed 100%. In addition, we also expect substantial contributions from higher orders in the 1/m D expansion which, however, are difficult to quantify.   Table 6. Same as Table 5 for longitudinally polarized vector mesons.  Given the large parametric and systematic uncertainties, we do not expect to obtain reasonably reliable predictions for the D → ρ + − differential decay rates within the QCDF framework. We will therefore briefly investigate to what extent at least some of the theoretical uncertainties might be reduced in ratios of decay widths. To this end we consider the ratio of partially integrated rates for transversely and longitudinally polarized vector mesons, In a generic NP scenario with sizeable short-distance contributions this ratio would be sensitive to the following combinations of Wilson coefficients in the large-recoil limit (see e.g. the discussion in [56]), where the primed coefficients refer to the operators with flipped quark chiralities,m c = m c /M D ,ŝ = q 2 /M 2 D , and possible contributions from scalar or tensor operators have been neglected. We present estimates for the quantity R T L in specific regions of momentum transfer ∆q 2 in Tables 9 and 10. We observe that, indeed, the parametric and systematic uncertainties are somewhat reduced, at least for the lower two q 2 bins: • For the neutral decay modes the scale dependence at NLO still can reach several 10 percent.
• The uncertainties induced by the parameter ω 0 in the D-meson LCDA typically reaches up to 30%. Notice that this estimate has been obtained using a particular (simple) model function (3.16) for the two 2-particle D-meson LCDAs.
• As our model correlates the hadronic-resonance effects in the different decay topologies, we also found a slight reduction of the associated uncertainty in R T L . We emphasize again that the hadronic model is only meant for illustration, and therefore one should not draw specific conclusions about the resonance effects from this.
We have also investigated to what extent hadronic uncertainties may cancel in the leptonic forward-backward asymmetry. Normalizing to the transverse rate, depending on the relative phase of a possible NP contribution to the Wilson coefficient C 10 , this observable is determined by the ratio Re C 9,⊥ (q 2 ) C * 10 /|C 2 9,⊥ (q 2 )|. As the real or imaginary part of the resonant contributions to C 9,⊥ look rather different from the square of the absolute value, we did not observe a significant reduction of the associated theoretical uncertainties in that case.

Conclusions
In this paper we have investigated the rare semileptonic decays D → ρ + − in the framework of QCD factorization (QCDF). Here, our primary goal was not to obtain very precise predictions for the SM decay rates, but rather to achieve a sufficiently realistic (i.e. conservative) estimate of hadronic uncertainties related to long-distance QCD effects.
The QCDF framework allows one to separate contributions from various decay topologies via a simultaneous expansion in the strong coupling constant and inverse powers of the heavy charm-quark mass. This includes factorizable and non-factorizable effects from quark loops, annihilation topologies, as well as spectator-scattering effects. Within the perturbative analysis, we have found that in general the contributions from non-factorizable effects dominate. In the case of neutral mesons form-factor corrections, annihilation and spectator effects enter with similar magnitude, whereas the decay modes with charged mesons are dominated by annihilation topologies alone.
We have also seen that -not surprisingly -the convergence of the QCDF predictions for the considered decay rates is relatively poor, which is clearly to be attributed to the relatively small charm-quark mass, but also -in the case of neutral mesons -due to accidental cancellations between the Wilson coefficients multiplying the dominant expressions. Even more importantly, the restricted phase space for D → ρ transitions together with the suppression of the purely short-distance effects (i.e. the ones encoded in the semileptonic and electromagnetic operators) in the Standard Model imply a much stronger sensitivity to hadronic-resonance effects as compared to the analogous B-meson decays. In order to get a quantitative estimate of these effects, we have extended a model, that has been originally designed to describe the effect of vector resonances in e + e − annihilation or hadronic τ decays. This allows us to study aspects of quark-hadron duality within the QCDF approach via dispersion relations and to estimate the related systematic uncertainties in the q 2 spectrum. (We repeat that we do not claim a realistic description of the q 2 -spectrum itself.) Together with the remaining parametric uncertainties coming from non-factorizable contributions and (yet unknown) higher-order effects in the QCDF approach (notably NLO corrections to the annihilation topologies), we have found that reliable theoretical predictions for the differential decay widths in D → ρ + − within the QCDF approach are almost impossible. On the other hand, part of the systematic and parametric uncertainties tend to partially cancel in ratios of decay widths, at least in certain regions of phase space. As an example, we have studied the ratio of transverse and longitudinal decay rates. Still, our conclusion about possible new-physics sensitivity of rare semileptonic D-meson decays tends to be somewhat pessimistic, unless one concentrates on bins in the lepton invariant mass far below (or, in case of radiative D → π + − decays, also far above) the light resonances, or one considers observables that (practically) vanish in the Standard Model. In this way one can still look for new-physics scenarios with large Wilson coefficients which are not excluded by current exprerimental data, see e.g. the discussion in the recent literature [13,14]. (To a certain extent, this also applies to the purely radiative decays D → ργ which have recently been seen by the Belle experiment [57].) On the other hand, one may exploit the fact that non-factorizable hadronic effects are more pronounced in D-meson decays, and use D → ρ + − as a playground for QCD studies that are also relevant to estimate sub-leading effects in the corresponding B-meson decays.
A Explicit Formulas for D → π + − Decays Defining generalized form factors for D → πγ * transitions as γ * (q, µ) π + (p )|H following [19], we obtain the factorization formula where the same functions T (i) ,± appear as in the case of decays into longitudinally polarized vector mesons, and ξ P (q 2 ) is the soft form factor for D → π transitions as defined in [19]. If we define again the coefficient function 11 the twofold differential decay rate can be written as (A.5)

B Semi-naive Model for Duality violation
In the perturbative analysis of non-factorizing effects from hadronic operators, the virtual photon is treated as a point-like particle. However, at time-like virtualities, q 2 > 0, the coupling of the photon to a (perturbative) sum of intermediate partonic states has to be replaced by the coupling to an (infinite) sum over physical hadronic states. Partonhadron duality (PHD) is expected to hold if one averages over a "sufficiently" large phasespace region of q 2 . In order to estimate the numerical effect of violation of PHD, we will construct a simple model which reflects the main theoretical features and phenomenological constraints, following the ideas discussed in [36][37][38] (see also [30,[58][59][60][61]).

Preliminaries: (a) Modelling of vector resonances
The standard procedure to model a vector-meson resonance is to consider a Breit-Wigner (BW) ansatz, In the narrow-width approximation, one would further neglect Γ V /M V 1. The above form gives a successful description of the resonance shape in the vicinity of q 2 M 2 V ; however, for our purposes we would also need to control the effect away from the resonance peak. A more sophisticated ansatz for a modified BW-like shape has been suggested by Shifman [37,38], where Among others, this form has the correct analytic behaviour, i.e. a branch cut at q 2 = 0 (we neglect the pion mass in the following). On the other hand, it reproduces the BW-ansatz when q 2 M 2 V ≈ σ 2 V and b V 1. The imaginary part of the modified BW ansatz can then be approximated as , where in the last line we have takenσ 2 V = σ 2 V /M 2 V = 1 for simplicity and used b V 1. For vector resonances in the ss channel, the definition of z V has to be modified to include the threshold for KK production. Neglecting mixing effects, this amounts to setting The imaginary part of the modified BW ansatz now reads (forσ 2 Notice that in this form, the original BW ansatz is recovered in the limit m 2 K → −∞.

Preliminaries: (b) Resumming an infinite tower of vector resonances
Shifman [37,38] also gives a simple model to describe the effect of an infinite tower of equidistant vector resonances with masses M 2 n = (n + a 0 ) σ 2 and widths Γ n = bM n , for n = {0, 1, 2, . . .} and a 0 = const. To this end, one considers the function Each individual resonance contributes with a modified BW term as discussed above. However, the infinite sum over all resonances reproduces the asymptotic result In [37,38], a simplified model has been constructed, taking a 0 = 1. In this case, one has ψ(z + 1) = ψ(z) + 1 z = ψ(−z) − π cot(πz) , (B.9) and for time-like momentum transfer, q 2 > 0, the imaginary part of π(q 2 ) receives an oscillatory contribution from the second term, with 1 π Im π(q 2 ) osc. Im [cot(πz)] where the last approximation is valid for q 2 > s 0 ≡ σ 2 2π b (and b π ln q 2 σ 2 1). The general idea is then to start from a perturbative result in the OPE/factorization approach, where the leading q 2 dependence is logarithmic, such that the spectrum is given by j OPE (q 2 ) = θ(q 2 ) + . . . and to replace it by the hadronic model (B.7) -which reproduces the oscillatory behaviour in (B.10) -and to add a finite number of vector resonances in the region 0 < q 2 < s 0 . In [37,38] a good fit to τ decay data in the vector channel has been obtained for σ 2 = 2 GeV 2 and b = 1/6. (More sophisticated analyses with qualitatively similar parameter values can be found in [58][59][60][61]).
Our model ansatz for the hadronic spectrum with up/down or strange quarks in the loop therefore looks as follows. 12 and The real part of the function under consideration can then be recovered from an appropriate dispersion relation. For the numerical illustration, we fix the parameters in the above ansatz as follows (Notice that for σ 2 and a we have taken the same values as for the simplified model that has been fitted to data; the assumed flavour-invariance to fix the parameter σ 2 s is consistent with the mass splitting m 2 φ(2175) − m 2 φ(1680) 1.9 GeV 2 .) With this input, the normalization factors n V and n φ will be fitted by requiring that the relevant integrals over j (i) (s) are identical for the OPE and the hadronic result. We emphasize again that our aim is not to provide a sophisticated and fully realistic model for the hadronic resonances, but rather to illustrate the systematic uncertainties associated to our ignorance about long-distance QCD effects.

B.1 Application to quark-loop topology
As explained above, the leading effect of hadronic operators with closed quark loops can be described in terms of the loop functions We thus replace the perturbative function through a once-subtracted dispersion relation that involves the above model for the hadronic spectral function, where the hadronic parameters in j q (s) depend on the light quark flavour, as indicated.
The default renormalization scale in the functions h(s, m q ) is chosen as µ = 1.5 GeV. For the ansatz (B.11) and (B.12), the parameters n V and n φ are tuned by demanding that the OPE result for h(q 2 → −∞, m q ) is reproduced, which implies This yields n V 1.94 and n φ 2.37, if the spectral parameters are fixed as explained in the previous subsection. We stress that at this point our goal is not to get a precise phenomenological determination of n V and n φ itself, which would better be achieved by standard sum-rule techniques applied to 2-point correlators (see also the discussion below).
Rather, the so-obtained numerical values allow for a consistent comparison with the partonic prediction in the context of QCDF at LO. The numerical comparison between the original (partonic) functions and the result of the simple hadronic model is illustrated in Figs. 9,10. Here we plot the real and imaginary part of the functions h(q 2 , m q ) and the above hadronic modification at negative and positive values of q 2 . Zooming into the region of small positive values of q 2 -which is relevant for D → ρ decays -we also show the contribution of a simple Breit-Wigner model for the low-lying resonances for comparison. We observe that • For negative values of q 2 and for q 2 5 GeV 2 , the perturbative result approximates the hadronic model very well.
• For small positive values of q 2 the hadronic model exhibits the expected oscillations around the perturbative estimate.
• In case of the real part (and also for the absolute value) of h(q 2 , m u,d ) our model leads to somewhat larger values compared to a simple Breit-Wigner ansatz, which can be traced back to the additional contributions from the higher resonances in the dispersion integral.
• By construction, the imaginary part of h(q 2 , m q ) in the vicinity of q 2 = M 2 V looks similar in our model and for a Breit-Wigner resonance. However, for q 2 → 0, our model requires the imaginary part of h(0 + , 0) to vanish, while the Breit-Wigner resonance formula yields a finite result proportional to Γ ρ . 13 As discussed above the leading contribution from the quark-loop topologies enters through the function Y (d) (s) and involves the incomplete GIM cancellation between strange and up/down quarks in the difference h(s, m s )−h(s, m d ), for which we show the comparison between the perturbative result and our hadronic model in Fig. 3 in the main body of the manuscript.
As an aside, it is also instructive to perform a QCD-sum-rule analysis of our ansatz. In the spirit of the standard duality argument used in QCD sum rules, one replaces the "true" hadronic spectrum (i.e. in our case the model for j q (s)) by a single resonance plus continuum, leading to Applying the standard Borel transformation, introducing the Borel mass parameter M , this yields In [13] the rho-meson width is multiplied by hand with a factor q 2 /mρ in the BW formula such that the imaginary part again vanishes for q 2 → 0. However, we remark that such an ansatz has wrong analytic properties for q 2 < 0.
where we have identified the last term as the leading-order sum rule for the vector-meson decay constant, see e.g. [62]. Using f ρ = 205 MeV, this yields n V 2.75 which is of the same order of magnitude as the values found for n V and n φ by the procedure described in the previous paragraph. [The numerical difference may be taken as a rough estimate for the intrinsic uncertainties of our simplified model.] Figure 9. Comparison of the perturbative result for the real and imaginary part of the function h(q 2 , 0) (gray solid line) and the model (3.10) using (B.11) for its hadronic modification (black dashed line) as a function of q 2 (in units of GeV 2 ). On the right-hand side, we also compare to the simple Breit-Wigner ansatz (blue dotted line). The parameter values are chosen as follows: σ 2 = 2 GeV 2 , a = 1, b = 1/6, µ 2 = 1.5 GeV 2 . The value of n V = 1.94 is tuned to reproduce the perturbative result in the limit q 2 → −∞.

B.2 Application to annihilation topology
In the QCDF approach, the leading contributions to the annihilation topology are determined by the functions To the level of accuracy that we are working with, it is sufficient to neglect 3-particle LCDAs in the D-meson. In that case the LCDAs φ ± D (ω) are not independent, but fulfill a so-called Wandzura-Wilczek relation [41], which for the above moment implies Moreover, an alternative representation of the moments can be achieved in terms of the "Wandzura-Wilczek" wave function ψ D (x) as defined in [63], where and x is twice the energy of the spectator quark in the D-meson. With this, one has Looking at the analytic properties of these expressions, one easily verifies that and also (B.22) holds by simple differentiation. As a naive ansatz we may simply replace θ(q 2 ) on the right-hand side of the above equations by the same function j q (q 2 ) that has been used to model the hadronic spectrum for the function h(q 2 , m q ) in the previous subsection. This corresponds to the naive factorization into the subprocesses D → ρqq (described in QCDF in terms of φ ± D ) followed by qq → γ * (described by the resonance model). 14 With this, we take 14 Obviously, this simple recipe would not work at higher orders in the strong coupling. For a related discussion of the analytic properties of the spectator contributions from the operator O g 8 in the context of light-cone sum rules can be found in [64].
As before, we tune the parameter n V (independently for both cases) to reproduce the asymptotic limit q 2 −m D Λ, i.e. we require Using an exponential model (3.16) for the D-meson LCDAs with ω 0 = 0.45 GeV, we obtain n V = 2.40 for λ − D (q 2 ) and n V = 1.75 for λ + D (q 2 ), respectively. The comparison between the perturbative result and the hadronic model is shown in Figs. 11,12, see also the discussion around Figs. 5,6 in the main text.
We may again analyze our model for λ ± D (q 2 ) in the framework of QCD sum rules. With the analogous steps as in (B.18) this leads to and after Borel transformation one has Again, within the intrinsic uncertainties, the numerical values for n V fitted to the asymptotic behaviour of λ ± D (q 2 ) are consistent with the above findings (without going into details of the sum-rule analysis). In particular, by comparing Eqs. (B.29) and (B.19), the value of n V -to first approximation -is expected to be universal 15 for the modelling of h(s, m q ) and λ − D (q 2 ). 15 A similar argument has been used in [65,66] to show the equivalence of the QCD factorization approach and the light-cone sum rule approach for the description of radiative QCD corrections to heavy-to-light form-factor ratios.

C More Input Parameters
In Table 11 we summarize a number of input parameters that have been used in the numerical analyses. The parameters describing the q 2 -dependence of the the D → ρ form factors are quoted after Eq. (3.2) in the main body of the text.  Table 11. Summary of input parameters not quoted in Table 2. All values taken from [53]. (Here λ, A,ρ,η are the Wolfenstein parameters for the CKM matrix.)