Disentangling QCD and new physics in D → πℓ+ℓ−

In this paper we consider the decay D+ → π+ℓ+ℓ−, addressing in particular the resonance contributions as well as the relatively large contributions from the weak annihilation diagrams. For the weak annihilation diagrams we include known results from QCD factorisation at low q2 and at high q2, adapting the existing calculation for B decays in the Operator Product Expansion. The hadronic resonance contributions are obtained through a dispersion relation, modelling the spectral functions as towers of Regge-like resonances in each channel, as suggested by Shifman, imposing the partonic behaviour in the deep Euclidean. The parameters of the model are extracted using e+e− → (hadrons) and τ → (hadrons) + ντ data as well as the branching ratios for the resonant decays D+ → π+R(R → ℓ+ℓ−), with R = ρ, ω, and ϕ. We perform a thorough error analysis, and present our results for the Standard Model differential branching ratio as a function of q2. Focusing then on the observables FH and AFB, we consider the sensitivity of this channel to effects of physics beyond the Standard Model, both in a model independent way and for the case of leptoquarks.


Introduction
The anomalies in measurements of b → s transitions remain unresolved, and whether these can be interpreted via new particles, e.g. leptoquarks, or an underestimation of uncertainties, i.e. experimental errors or in the treatment of charm loops, remains uncertain. If these are indeed a sign of physics beyond the Standard Model (BSM), it is plausible that such physics could affect other flavour changing neutral current (FCNC) processes, amongst JHEP04(2021)158 which the c → u transition is the least constrained (for a recent study of K → π ( ) see ref. [1]). At the same time LHCb is currently producing unprecedented numbers of D mesons, the charm production cross-section, σ(cc) = 1200 µb [2], exceeds the bottom production cross section of σ(bb) = 75 µb [3] by far.
The lack of constraints on the c → u transition is due to the relative difficulty with respect to b → s decays in making the necessary predictions. This is in part due to the reduced hierarchy between the charm mass and Λ QCD compared to the corresponding hierarchy for B decays, which makes expansions less effective. In addition, it is due to the fact that for D decays the resonances affect a larger portion of the phase space. In the case of B meson observables, the contributions due to the light resonances (ρ 0 , ω 0 and φ) can be circumvented as the resulting effects in binned observable are negligible since the typical bin size is large compared to the width of the states. Nonetheless, polluting resonant effects due to the c-quark loop, are much larger and the kinematic regime where charmonium resonances are produced is ignored and often vetoed in the experimental analyses. However, in the case of D decays the resonances affect a larger fraction of the available phase space, such that predictions are required in regions not sufficiently far from the resonance tails and modelling the structure of the hadronic resonances becomes crucial.
In a first approach to this problem in D + → π + + − decays, the resonances were added "by hand", for example see refs. [4,5], by means of Breit-Wigner functions, on top of a non-resonant background described by the partonic result for the quark vacuum polarisation. In ref. [6], an alternative approach is advocated, where a subtracted dispersion relation is used to reconstruct the s-and d-quark vacuum polarisations from the respective imaginary parts, which must be modelled. In this approach, the partonic result is recovered asymptotically and the resonances are described following Regge trajectories with a number of simplifications, the main one being that the isospin 1 and isospin 0 channels are not treated separately, but in terms of a single tower of resonances with "effective" parameters.
We aim to improve upon the model of [6] by adapting the strategy applied to B → K in ref. [7], where the authors investigate the charm resonance contribution, extracting the charm vacuum polarisation from e + e − → (hadrons) data by means of a dispersion relation and the optical theorem, an approach pioneered in the context of b → s transitions in ref. [8]. Here we extract the u, d and s vacuum polarisations from e + e − → (hadrons) and τ → hadrons + ν τ data modelling the resonances present in the different channels with a Regge-based description similar to that of ref. [6]. The fact that in the low-energy region the three quark flavors are active makes the interplay between the different resonances more intricate and disentangling the contributions from the different quark flavours more challenging. The uncertainties on our final vacuum polarisations are sizeable; this stems from the uncertainties on the data sets and from the inherent limitations in modelling hadronic resonances. Given the challenges arising in obtaining a theoretical description of rare FCNC D decays, it is necessary to remain conservative and seek observables which evade these limitations while being sensitive to beyond the Standard Model (BSM) physics.

JHEP04(2021)158
Here we propose to go beyond the existing work using the above-described novel treatment of resonances in combination with the treatment of the perturbative contribution including QCD factorisation (QCDf) corrections derived in ref. [6] and the operator product expansion (OPE) weak annihilation corrections adapted from ref. [14]. Note that this is the first phenomenological study of D + → π + + − where these QCDf and OPE corrections are implemented which, due to the large contribution of weak annihilation in this channel, should not be neglected. Note however that weak annihilation effects for D → V γ and three-body hadronic D decays have been studied in the QCDf framework in refs. [15,16]. Our results are further ameliorated by the recent advances on the form factors from the Lattice [17,18] and Wilson coefficients from ref. [19] which we incorporate. Our aim is to, using the theoretical framework just described, provide accurate and realistic predictions of the observables sensitive to BSM physics with little dependence on the non-perturbative physics, along with conservative error estimates.
In this paper, we first lay out the theoretical framework necessary for the perturbative part of the calculation in section 2, most importantly introducing the operator basis and the calculation of the matrix element. This is followed by a detailed description of our formalism for the resonant contribution in section 3, providing details of the model and the fit to the experimental data. In section 4 we define a set of observables which are both sensitive to BSM physics and minimally affected by the hadronic uncertainties, and explore the related phenomenology both in the SM and for model-independent effects, keeping in mind the existing experimental constraints. Our conclusions can be found in section 5.

Naive amplitude and annihilation contributions
In this section we present the framework for the calculation of the amplitude for D + → π + + − , ignoring the effect of resonances, which can be divided into a naively factorisable part and a part containing non-factorisable corrections, where by non-factorisable we mean that the diagram cannot be factorised into a perturbative part and a form factor.

Theoretical framework
For c → u transitions, between the scale where the W boson is integrated out (µ W ∼ M W ) and the scale µ b ∼ m b , only operators O 1 and O 2 are present in the SM effective Hamiltonian. However, integrating out the b quark the penguin operators C 3−9 are generated. We adopt the basis of operators and effective Hamiltonian used in ref. [5] and write: where G F is the Fermi constant, the combination of CKM matrix elements, λ q , is given by λ q = V * cq V uq , and

JHEP04(2021)158
However, note that since λ b λ d , all contributions entering H (b) eff are heavily CKM suppressed. The complete set of operators used in this paper is then composed of the current-current operators where p = d, or s, and the operators O 3−10 , which read: (pγ µ T a p), where T a are the SU(3) c generators, q L/R = (1 ∓ γ 5 )q/2 denote the left/right-handed quark fields, m u is given in the MS scheme at the scale µ c ∼ m c , g e = √ 4πα e is the electromagnetic coupling, α e the fine structure constant, and F µν and G a µν are, respectively, the electromagnetic and chromomagnetic field strength tensor. In the semileptonic operators, represents the lepton field. To be exhaustive, the SM basis also contains the chirality-flipped operators O i identical to the O i up to the transformation q L/R → q R/L . As the Wilson coefficients of these (and also that of O 10 ) are negligible in the SM, these will only be included amongst the BSM contributions to the Hamiltonian. We calculate the Wilson coefficients in the SM at nextto-next-to-leading logarithmic order (NNLL) following ref. [19], details can be found in appendix B.
This SM Hamiltonian for D + → π + + − can be extended to receive BSM dimension-six contributions, which give rise to additional Lorentz structures where the following operators are defined:  In the SM, the leading contributions in an expansion in the strong coupling α s and 1/m c for the D + → π + + − decay arise from: • interactions via the semileptonic operator O 9 , • interactions where the charged lepton pair originates from a virtual photon γ * emitted via the EM dipole operator O 7 (as shown in figure 1(a)) or via the 4-quark operators (as in figure 1(b)), • weak annihilation as seen in figure 1(c).
In naive factorisation, only the first two are taken into account. Annihilation topologies fall under the set of so-called non-factorisable contributions. However, in contrast to the corresponding B decays, the annihilation topology is not CKM suppressed, and therefore, as observed first in [6], the decay amplitude turns out to be dominated by non-factorisable dynamics.
The closed fermion loop in figure 1(b) is calculable perturbatively as long as the qq pairs remain off-shell and for invariant dilepton masses away from hadronic resonances. In the region where the invariant dilepton mass squared s corresponds to that of hadronic resonances, non-perturbative methods are required.

The naive amplitude
Based on the effective Hamiltonian described above, we can write the amplitude M(D + (p) → π + − + ) as in ref. [5]: The function β(s) and the Källen function λ(s) entering the expression of the tensor functions F T (s) and F T 5 (s) are given by In the SM, the only non-vanishing structure is the contribution of f + to F V (s), where the functions C (b) 9 (s) and C (d) 9 (s) can be expressed as [6] C (q) (2.8) in terms of the "hard kernels" T (q) . In naive factorization these hard kernels are given by where the coefficient C (0,q) contains only factorizable and leading contributions: (2.10) The 1-loop functions Y (q) combine the contribution from the four-quark operators

JHEP04(2021)158
whereh q (s) is the closed-quark loop function, discussed in detail the next section. Note that our closed fermion looph q (s) is related to that defined in [6], denoted h(s, m q ), byh q (s) = 9 4 h(s, m q ). For completeness, we give the explicit expression ofh q (s) in perturbation theory, h (pt) (s), at leading order where m q is the mass of the quark q and ζ = 4m 2 q /(s + i ). The quark-loop functions h(s, m q ) are needed in the resonance region. The perturbative description of these functions does not include any long-distance hadronic effects. In the next section, we will include these by reconstructing the functionsh(s, m q ) from their imaginary part by means of a dispersion relation, where the imaginary part is obtained by constraining the parameters of a hadronic model for the spectral functions [6]. Finally, we remark that there is a partial cancellation in the combination of Wilson coefficients that participate in eq. (2.11) which leads to an accidental suppression of the LO contribution.
Combining eqs. (2.8) and (2.9) and (2.10), the functions C (q) 9 in naive factorisation, denoted C (q) 9 Naive , can be summarised by As discussed previously, the naive result is not sufficient, primarily due to the weak annihilation diagrams. We will now consider corrections, calculated within QCD factorisation in the low q 2 region and in the OPE framework in the high q 2 region (above the φ pole).

Annihilation contributions
The expressions for the QCDf corrections were calculated in [6], where the authors adapted expressions for the B → K * + − from ref. [20] to the case of D → ρ(π) + − . Following [20], in QCDf the non-factorisable contributions can be classified under two categories: • A first category where the spectator quark participates in the hard scattering; annihilation topologies enter in this category. Here the spectator quark participates in the FCNC process via hard gluon exchange. The calculation of these processes leads to the so-called hard spectator scattering corrections.
• A second category where the spectator quark is connected to the hard process only through soft interactions and the hadronic transition can be described by the form factors. The calculation of which leads to the so-called form-factor corrections.
QCDf is applicable in the combined heavy-quark and large-energy (recoil) limit, where energy E refers to that of the final state meson, related to the dilepton invariant mass s via: (2.16)

JHEP04(2021)158
In QCDf the decay amplitude is schematically expressed as The non-factorisable corrections enter the coefficients C (q) as well as the second term made of a convolution product between the so-called hard kernel T (q) where φ ± D is the LCDA of the D meson. More explicitly, the naive factorisation results in eq. (2.9) can be extended by where the ± subscript refers to the projection of the amplitude on the D meson LCDA.
The perturbative quantities C (q) and T (q) ± are given by where we remind the reader that a s ≡ α s /(4π) and C (0,q) is defined in eq. (2.10). The form-factor corrections are contained in C (1,q) , the annihilation corrections in T (0,q) and the hard spectator scattering corrections in T (1,q) , expressions for all these quantities can be found in ref. [6]. Finally, in the QCDf framework the function C 9 FF are the annihilation, the spectator scattering, and the form-factor corrections to the C (q) 9 functions, respectively. Note that unlike ref. [6] we have chosen to use the pole mass for the charm quark, and have therefore set the quantity ∆M accordingly the definition of which can be found in ref. [21]. We have studied these sets of QCDf corrections numerically in section 4.2 (see table 5, and find that, in agreement with [6], the only one numerically relevant to our analysis, given the size of the theoretical errors, is C (q) 9 Ann . We will therefore discuss these annihilation corrections in more detail.
The four annihilation diagrams shown in figure 1(c) contribute at different powers in the 1/m c expansion. With the convention that the π + meson momentum is nearly lightlike in the minus light-cone direction, the amplitude for the surviving contributions depend only on the minus component and T being u-independent, corrections to C (q) 9 are then given by with the s-dependent moment λ − D (s) given by We note that C (q) 9 Ann leads to a sizeable contribution as C 2 appears without any cancellation from other Wilson coefficients and, as a result, it turns out that annihilation gives very large contributions to the decay rate. It should also be mentioned that eq. (2.25) is purely partonic and does not include any effect due to the resonance spectrum. In our final results we will modify λ − D (s) to include those effects employing the ansatz of ref. [6]. We remark that the calculation of the annihilation contribution could also be performed in light-cone sum rules as applied to the B-decay case in refs. [22,23]. We have verified that the two results are approximately compatible. Given the dominance of and the large uncertainties on the contribution of the hadronic resonances, we don't expect significant changes in the results should the LCSR framework be adopted for this contribution.
The QCD factorisation framework is valid at small q 2 , as here the pion is energetic, and in the heavy-quark limit. To be precise, the energy of the pion E π should be large compared to Λ QCD . We would also like to be able to make predictions at large q 2 , away from the dominant hadronic resonances, where the theoretical uncertainties are under better control. Since the QCDf weak annihilation provides such a large contribution in the low-q 2 region one wonders whether this is also true at high q 2 and how to rigorously estimate these contributions in this regime. Several papers exist which tackle this issue for B → K ( * ) transitions [14,24]. Here we follow ref. [14] where one can perform an operator product expansion (OPE) exploiting the large size of q 2 compared to the pion energy ( q 2 E π , Λ QCD ) to expand the amplitude in powers of E π / q 2 . The amplitude can then be easily factorised into the form factor and a coefficient which can be calculated perturbatively to order α s . The results obtained in ref. [14] for the weak annihilation contribution (which arises at leading order in α s ) to C 9 can easily be adapted to the case of D + → π + + − , and we find Being proportional to C 2 , as opposed to C 4 + C 3 /3 for the case of B → K , annihilation also dominates in the high-q 2 regime (we neglect the other contributions which are Cabibbo suppressed). In our implementation, the QCDf and OPE corrections must be adopted only in the appropriate range in q 2 , i.e. for QCDf E π Λ QCD and for the OPE q 2 E π , Λ QCD . These conditions can be satisfied if we use the QCDf corrections only up to the φ pole,

JHEP04(2021)158
and for higher energies the corrections in the OPE framework. Note that the latter are not valid up to the lowest recoil point, and therefore in our phenomenological analysis we only consider q 2 up to 2.3 GeV 2 .

Hadronic resonances
As discussed in the previous section, in D + → π + + − decays, contributions with a closed quark-loop arise at leading order, where the d-and s-quark vacuum polarisations intervene. These contributions are encoded in theh(s, m q ) functions of eqs. (2.11) and (2.12). Here we discuss in detail how we obtain a description for these functions beyond perturbation theory, including effects of the hadronic resonances.
Due to phase space restrictions, any realistic phenomenological study of D + → π + + − decays will require a description of resonance effects. The resonances must be included, ideally, without losing contact with the first-principles calculation. In ref. [4], the resonances were added by means of Breit-Wigner functions, on top of a non-resonant background described by the partonic result for the quark vacuum polarisation (at leading order). A somewhat more sophisticated approach was advanced in ref. [6] in which the vacuum polarisations are modelled following ref. [25]. In this approach, the partonic result is recovered asymptotically, away from the resonant region, and the resonances are modelled following Regge trajectories with a number of simplifications, the main one being that the isospin 0 and isospin 1 channels are not treated separately, but in terms a single tower of resonances with "effective" parameters [6].
Here, we improve on the description of ref. [6] by implementing a strategy similar to that of ref. [7], where the charm vacuum polarisation was extracted from e + e − → (hadrons) data by means of a dispersion relation. The general framework is similar to ref. [7], but the fact that we are dealing with light quarks leads to a number of complications. Below the charm threshold, quarks u, d and s contribute to the vacuum polarisation probed in e + e − → (hadrons), with several resonances that interfere. In the case of charm, however, one can safely subtract the perturbative background arising from the light-quark contributions and work with data that are "pure" cc. Here, in order to be able to disentangle the different channels better we will make use of τ → (hadrons) + ν τ data as well, which are pure isospin 1, as opposed to the e + e − → γ * → (hadrons) data which are a mixture of isospin 0 and 1.
However, no treatment of the resonances in the vacuum polarisations needed for the description of D → π + − is completely sound from a theoretical perspective. The reconstruction of the vacuum polarisation functions using a dispersion relation requires knowledge about their imaginary parts from threshold to infinity. This description cannot be done within a unique framework from the low-q 2 region, where chiral symmetry plays a role, up to the high-q 2 regime, where perturbative QCD is the main contribution, passing by resonance peaks and residual resonance oscillations. In the end, non-factorisable corrections will also play a role and the quark-loop functions that appear in the D + → π + + − amplitude are expected to be modified with respect to the description of the quark vacuum polarisations. The different treatments of the resonant contribution should serve, in the end, to obtain a reasonable estimate of their effect, with a conservative error band, in JHEP04(2021)158 order to identify regions of the spectrum where room for BSM physics still exist and where theory is under control. With these caveats in mind, we discuss now how the framework of ref. [7], together with the model from refs. [6,25], can be adapted to our case.
In the calculation of D + → π + + − one needs the d and s quark vacuum polarisation functions, the imaginary parts of which intervene in the observable where the second (approximate) equality is valid at LO and for values of s for which the muon mass can safely be neglected. Below the charm threshold, experimental data for R(q 2 ) contains the entangled contributions from u, d and s quarks.
On the theory side, R(q 2 ) can be written in terms of the imaginary part of theqq contributions to the photon vacuum polarisation as where the vacuum polarisation is defined as with s = q 2 and the electromagnetic current J EM µ given by (where we omitted the sum over colour indices). Henceforth we will disregard the mixed quark flavour contribution which can occur through disconnected diagrams that are 1/N c suppressed [26]. Below charm threshold, we then write the R ratio as the sum of the three quark-flavour contributions We then define the correlators Π (q) (q 2 ) in terms of the flavour-singlet currents j (q) µ = qγ µ q, in complete analogy with eq. (3.3). To make contact with refs. [6,7] it is convenient to define the functionh(q 2 ) ash At sufficiently high momenta in the deep Euclidean region, q 2 → −∞, we have, after normalisation,h Upon analytic continuation one obtains, for q 2 > 0, in the chiral limit, Imh q (q 2 ) = π.
For the purpose of modelling the resonances it is appropriate to consider the different isospin channels separately. The electromagnetic current of eq. (3.4) contains an I = 1 and I = 0 light-quark part, as well as the strange-quark contribution, the respective dominant vector resonances being the ρ(770), the ω(782), and the φ(1020). Taking into account the quark-charge factors, the electromagnetic current can be written as

JHEP04(2021)158
where we made the I = 1 and I = 0 light-quark contributions explicit in the first and second terms in between parenthesis on the right-hand side, respectively, where J µ For the R ratio one has then 1 whereh 1/0 (q 2 ) represent the light-quark I = 1 and I = 0 contributions andh s (q 2 ) is the contribution from the strange quark. The imaginary part of these functions are proportional to the respective spectral functions. In the case of the I = 1 current, additional experimental data for the spectral function exists from τ → (hadrons) + ν τ decays [28][29][30].
We will use these data sets as a way to help disentangling the different contributions in R(q 2 ). Ultimately, in the application to D + → π + + − , within our assumptions, we need the functionsh d (q 2 ) andh s (q 2 ). The former can be obtained fromh 1/0 (q 2 ), which contribute equally to the d-quark vacuum polarisation. We then need a concrete model for the imaginary part of the differenth I (q 2 ) functions of eq. (3.8), with I = 1, 0 or s. Our description is based on the model suggested in appendix B of ref. [6] which, in turn, is based on a proposal by Shifman [25]. The model can be summarised as follows. The imaginary part of each channel is modelled with a dominant vector resonance plus the sum of an infinite tower of resonances, starting from the first excitation, with masses following Regge trajectories. The dominant resonance in each channel is modelled by a Breit-Wigner function, f (R) BW , while the sum over the infinite tower of resonances is performed analytically within Shifman's model [25]. Concretely, after summing over all the excited states, we have The first term on the r.h.s. corresponds to the dominant vector resonance and is discussed below, the second term represents the sum over the infinite tower of equally spaced resonances with masses M 2 I (n) = (n + a I )σ 2 I (3.10) and widths In eq. (3.9), Ψ(z) is the digamma function and for I = 1, 0 while, for the strange quark,

JHEP04(2021)158
(We are treating the light quarks and the pions as massless.) The term with the Ψ function corresponds to the description of ref. [25], reviewed in [6], to which we refer for more details about the model. The description of the tower of resonances adds three parameters per channel to the model (σ i , a i , and b i ), therefore nine in total.
The dominant resonances are the ρ(770) in the I = 1 channel, the ω(782) in I = 0, and the φ(1020) in the ss channel. For the description of the leading resonances we use the following Breit-Wigner function where M R and Γ t are the Breit-Wigner parameters related to the mass and the total width. The phase α ρ(770) is taken to be zero and is used as a reference for the phases of the ω(782) and φ(1020). The imaginary part in the denominator of the f to account for the non-zero kaon mass. We fix the parameters of our model from a comparison to e + e − → (hadrons) and τ → (hadrons) + ν τ data -a strategy similar to the one of ref. [7]. It is known that models related to the one we are employing here agree well with the data in an asymptotic region [26,30,31], where q 2 is large enough (in practice this means q 2 1.5 GeV 2 , but this value is channel dependent). In the case of the τ data, which is pure I = 1, this type of description has also been extended to lower energies with the inclusion of a Breit-Wigner for the ρ(770) [25].
Here we use the publicly available Particle Data Group compilation of R(q 2 ) data [32] supplemented with R(q 2 ) measurements from the BES and KEDR collaborations published recently in refs. [33][34][35]. No correlations are publicly available for these data, and the data are therefore treated as uncorrelated. These data are a mixture of the three channels I = 1, I = 0, and ss as described in eq. (3.8). To better disentangle the three channels, we also use the ALEPH [28] and OPAL [29] data for the vector-isovector spectral function from τ → (hadrons) + ν τ . In the case of OPAL data we use the updated version of ref. [30]. These data sets are also treated as uncorrelated, to be fully coherent with our treatment of the data for R(q 2 ).
Since the number of free parameters in the model is quite large 2 and the interplay between the different contributions to R(q 2 ) is intricate, we need to make a few assumptions in order to fix all the parameters in our description. The isospin 1 parameters a 1 and σ 1 from eq. (3.10) are fixed, from a fit to the τ -decay data, to be σ 2 1 = 2.476 GeV 2 and a 1 = 0.974, which is in the ballpark of values expected from Regge behaviour (2 GeV 2 and 1.0, respectively [6]). We also fix the parameter b 0 , related to the widths of the resonances in the . Additionally, we assume that σ 2 1 = σ 2 0 -which amounts to the assumption that the masses of higher resonances with I = 1 and I = 0 follow a similar pattern. We then build an uncorrelated χ 2 function from the R(q 2 ) and τ data combined. We use all the τ -decay data available, which gives 174 data points. From the R(q 2 ) data we exclude the low-energy data below 0.3 GeV 2 since the use of our Breit-Wigner functions is not fully reliable at such low energies. We also exclude 8 data points around 2.1 GeV 2 where the data shows a fluctuation downwards, clearly visible in figure 2, which cannot be described within our model. The inclusion of these data points do not change the values of the parameters significantly, but leads to a worse χ 2 value. In total we use 637 points from the R(q 2 ) data set, with the highest energy bin being at 4.0 GeV 2 . Little information is added if one includes data beyond this point (mainly from BES and KEDR), since data points become scarce. The fit gives then the value χ 2 min = 828.1 for 795 degrees of freedom. The resulting parameters from this minimisation are shown in table 1. In figures 2 and 3 we compare the results of the fit with the R(q 2 ) and τ -decay data. The errors given in our table 1 are quite small for some of the parameters. However, to accommodate further non-factorisable effects, in our practical use of these results we will allow for a variation of the resonance phases, α R of eq. (3.14). In the end, this variation will be one of the main sources of error in the description of resonances in D + → π + + − decays.
Our model is sufficient to achieve a reasonably good representation of the data. One should note that, at higher energies, the model is systematically below the data due to the lack of perturbative corrections, which are of the order of 15%. With this description, we obtain the imaginary parts of the three non-perturbative functionsh I (q 2 ).
From the imaginary part ofh I (q 2 ) we can reconstruct the full function using a dispersion relation. We follow the suggestion of ref. [6] and use a once-subtracted dispersion relation with the subtraction constant fixed from the perturbative description. Accordingly, the subtraction point is chosen in the deep Euclidean at q 2 = −s 0 and the functions h I (q 2 ) are given byh OPAL data ALEPH data our model Figure 3. Vector isovector spectral function from hadronic tau decay data [28][29][30] compared with our model with parameters given in table 1.

Parameter
Central value Relative error The subtraction constant is calculated from the perturbative description (without α s corrections) which is given by eq. (2.13). (For the light-quark contributions we take eq. (2.13) in the chiral limit.) When reconstructing the functionsh q (q 2 ) we have checked that using a dispersion relation with more subtractions leads to very similar results. We have also JHEP04(2021)158 checked that the results are stable upon variation of the subtraction point in eq. (3.15). For our final results we use s 0 = 10 GeV 2 and µ 2 = (1.5 GeV) 2 .

Branching ratios and the φ normalisation
The direct application of the results from the above fit to e + e − and τ data assumes that the quark loop can in fact be factorised from the mesonic D → π transition. While in ref. [7] it was argued that this is a good approximation for the analogous charm loop in B → K , it is prudent to verify the extent to which it can be trusted here. We cannot follow the strategy of ref. [7] for want of direct measurements of the resonance spectrum in the D + → π + + − decays. We do this by comparing our predictions for the branching ratio D + → Rπ + → π + µ + µ − , with R = φ, ρ or ω to the limited experimental results available. 3 In the case of D + → φπ + → π + µ + µ − , there is only one measurement [36] which gives B(D + → φπ + → π + µ + µ − ) = (1.8 ± 0.8) × 10 −6 , and is consistent with the estimate obtained from the product of the individual branching ratios of D + → φπ + and φ → µ + µ − . The calculation of this branching ratio within our description (in the SM), and using the parameters of table 1 for the functionsh I (s), gives 1.3 × 10 −8 , which is still compatible within a little more than 2σ with the experimental result, given the large uncertainties. 4 We decide, nevertheless, to re-scale the parameter n φ by a factor of 11.8, in order to account for this mismatch, presumably a result of the combination of an experimental fluctuation and non-factorisable effects. After this re-scaling, the central value for B(D + → π + φ → π + µ + µ − ) calculated with our description reproduces the central value of the experimental measurement. Integrating our decay distribution around the ρ − ω region we find a branching ratio of ∼ 10 −7 . The experimental counterpart to this value can be estimated from the sum of B(D + → Rπ + )×B(R → µ + µ − ) with R = ρ, ω [37], which gives (6 ± 1) × 10 −8 , not too far from the value we obtain. Therefore, we do not perform any further rescaling of the normalisation of the resonances. Upper bounds on the non-resonant branching fractions exist from LHCb [38]. The limits are shown in the table 2 in two regions, Region I (low q 2 ) and Region II (high q 2 ). After the re-scaling of the φ contribution (which affects the branching ratio in the non-resonant region through the tale of the resonance) we find, in the SM, 5 which fulfil the experimental bounds. (Our errors are 68% C.L. We postpone to section 4.3 a detailed discussion about the estimate of theoretical uncertainties.) These results appear 3 Note that we ignore the contributions from the narrow pseudoscalar mesons, e.g. η, η in our D + → π + + − description. 4 The values for all the parameters that enter this calculation can be found in table 3. 5 The calculation of the branching ratio in the high-q 2 region requires integration up to the kinematic end point of the spectrum. This is beyond the expected validity range of the OPE approach to the weak annihilation (q 2 ≤ 2.3 GeV 2 ) contribution. This uncontrolled systematic error should be relatively small since the main contributions to the integral arise from the region where we expect the OPE to be valid.  Table 2. Current best-world limits on B(D + → π + µ + µ − ) in two bins of the dilepton invariant mass s [38]. In the experimental search, the branching fraction excluding the resonant contributions is extrapolated assuming a phase space model.

JHEP04(2021)158
to be in line with the findings of refs. [4,5] although there they were not given in the same form. Note that the purely short-distance contribution to the functionsh(s) would lead to much smaller branching fractions, of the order of 10 −12 , again in agreement with the results of refs. [4,5]. Branching fractions of the order of 10 −9 can be at LHCb reach once the full Run2 data is analysed. We conclude that the re-scaling the φ normalisation is compatible with all known experimental bounds and, from now on, we will use the value of n φ given in table 1 multiplied by 11.8.

Numerical and phenomenological analysis
We will now perform a thorough phenomenological analysis of the decay D + → π + + − : we first define the observables of interest; we present our results in the SM, where the set of observables reduces to the differential Branching Ratio; allowing BSM contributions to the Wilson coefficients, we then assess the existing experimental bounds on such contributions and determine how the aforementioned observables are affected.

Definition of the observables
Following [5,39], the double differential D + → π + + − decay width with respect to s the dilepton invariant mass and θ, the angle between the three-momenta of D + and − in the rest frame of the lepton pair, is given by where the kinematic function β and λ are given in eq. (2.6), and the dependence on s is assumed to be understood. The normalisation factor N is given by and the three angular coefficients a , b and c by

JHEP04(2021)158
Note that these three angular coefficients form the basis of the observables we will construct, and given that there are three independent coefficients (in the absence of CP violation) there will clearly only be three independent observables: we choose the differential decay width, the flat term, and the forward-backward asymmetry. On the other hand, in the SM the three angular coefficients are not independent, one finds the expressions: where F SM V is the SM contribution to F V defined in eq. (2.5), depending on C q 9 (s), f + (s) and CKM factors. There is therefore only one independent observable in the SM, the decay rate, and all others will vanish.
Coming back to the general case, in order to obtain the first of these observables, the differential decay width, we integrate eq. (4.1) with respect to cos θ, Here it is convenient to introduce dΓ (s) as it arises in the definition of several observables: and depends on all the considered Wilson coefficients and on the D to π transition form factors, f + (s), f 0 (s) and f T (s).
Following [5], the next observable, the flat term, F H , can be defined as where the second line in eq. (4.7) is valid only in the limit m → 0 or in the high-s region where β ∼ 1. In the SM, up to O(m ) corrections the flat term is only a function of β: as the Wilson coefficients and form factors cancel, clearly resulting in a small uncertainty on the helicity suppressed SM distribution. We note that F H is numerically very small for q 2 m 2 . Effects of physics beyond the Standard Model would therefore stand out for this observable, which is primarily sensitive to C P , C S , C T and C T 5 [4,5].

JHEP04(2021)158
Noticing that b = 0 might no longer hold in BSM scenarios, it is interesting to build observables sensitive to the angular coefficient b , one example of which is the forwardbackward asymmetry A FB . Following [5], A FB can be defined as Finally, in order to explore the sensitivity to possible BSM CP phases, we consider the CP -asymmetry (A CP ) defined by [4]: where dΓ/ds is the differential decay rate of the CP-conjugate mode,  where the range if integration in the dilepton mass squared is between q 2 min and q 2 max .

Results in the Standard Model
In the SM, making use of eq. (4.4), the double differential decay rate can be written such that in the SM, beyond the kinematic functions (λ and β), the differential decay width only depends on the function F SM V (s). Integrating over cos θ, and again using eq. (4.4), we find dΓ(D + → π + + − ) It is therefore of interest to examine the different contributions to F SM V (s) in more detail. Before doing this we first need to discuss the numerical values of the parameters used in our work, summarised in table 3 along with the appropriate references. The first set of JHEP04(2021)158 parameters are the masses of the quarks, where we adopt the MS values of the strange, bottom and top masses and the pole mass of the charm quark as explained in section 2.3, and the on-shell mass of the W boson, taken from ref. [32]. For the scales, µ c is set to the pole mass of the c quark, for the scale µ b we use the MS mass m b (m b ), and for the scale µ W the on-shell mass of the W boson. These scales µ c , µ b , µ W are varied between µ/ √ 2 and √ 2 µ, to account for residual uncertainties related to renormalisation-scale variation. Next we come to the hadronic parameters characterising the decay constants of the charged pion and D meson, which we take from ref. [32], and the pion and D meson lightcone distribution amplitudes (LCDAs); definitions of ω 0 for the D meson LCDA and a 2 and a 4 for the pion LCDA can be found in appendix A.2. Unlike the LCDA for of the B meson, very little is known about that of the D meson. Therefore, following ref. [6], we adopt an ad-hoc range for ω 0 , keeping in mind the naive expectation from heavy-quark symmetry and choosing a sufficiently large uncertainty to remain conservative. As can be seen in eq. (A.4), the only numerical input required for the twist-2 π LCDA are the Gegenbauer moments. These can only be calculated via non-perturbative methods. While there has been a great deal of progress from the Lattice in calculations of a 2 [40], the most recent calculation being that of the RBC collaboration with N f = 2 + 1 flavours of dynamical Wilson-clover fermions [41], a 4 has not so far been calculated. We therefore adopt a 2,4 (1 GeV) from ref. [42], where the light-cone sum rules (LCSR) result for the pion electromagnetic form factor [43] is fitted to experimental data [44]. The extracted values, a 2 (1 GeV) = 0.17 ± 0.08 and a 4 (1 GeV) = 0.06 ± 0.10, where the errors reflect both experimental and theoretical uncertainties, are consistent with previous results from sum rules and Lattice QCD. 6 For the form factors, we adopt the recent calculation on the Lattice by the ETM collaboration [17,18] (again the details of the parametrisation can be found in appendix A.1). The parameters given in table 3 are taken from refs. [17,18]. Note that in the low-q 2 region, these form factors were also calculated in LCSR in ref. [46]. However, table 3 does not provide the values of the Wilson coefficients C 1−9 defined in section B. The numerical results for the Wilson coefficients, calculated as described in appendix B for the central value of these scale (as given in table 3), are summarised in table 4. We note that only C 1 , C 2 and C 9 have sizeable values, where however the signs of C 1 and C 2 are opposite, the Wilson coefficients related to the strong penguin and the dipole operators are numerically small. Further, we have checked that for the same set of input parameters, our results agree with those given in ref. [49].
Having provided the numerical values of all parameters that enter our analysis, we are now ready to examine the contributions to F SM V (s), where C (q) 9 = C (q) 9 Naive + C (q)  Table 3. Summary of the numerical input used in the study of D + → π + + − . In naive factorization, 0.028 − 0.033i 0.005 + 0.002i Table 5. Breakdown of individual contribution at NLO for the D + → π + + − decay at s = 0.5 GeV 2 .
where the QCDf corrections C (q) 9 Ann , C 9 SS and C (q) 9 FF are provided in section 2. The numerical results for these corrections at s = 0.5 GeV 2 are presented in table 5. We note that Y (q) (containing the resonance contribution) and C (q) 9 Ann are the only sizeable contributions. We further remind the reader that the (q = b) contributions are CKM suppressed since λ b λ d . Hence, F SM V (s) can be approximated by Ann f + (s). (4.17) As was first noticed in ref. [4], Y (d) is further suppressed due to the cancellation between C 1 and C 2 occurring at the scale m c in the factor (2/3C 1 + 1/2C 2 ) ∼ 0.34. In our final results, we further modify the q 2 -dependent moment, λ − D (s), of eq. (2.25) to take into account effects due to the hadronic resonances as where j d (s) = 1 π Imh d (s). The constant n d is introduced to ensure that the perturbative result is recovered for q 2 −m D Λ QCD , as discussed in appendix B of ref. [6]. The branching ratio distribution is related to the decay width distribution defined in eq. (4.5) via where Γ D is the total decay width of the D 0 meson. The final result for the differential branching fraction in the SM, taking into account all uncertainties, is given in figure 4. (Below, in section 4.3, we will discuss in detail how the uncertainty bands are obtained.) In the same figure, we also show the two separate contributions in eq. (4.17): in red the component due to Y (d) and in green the weak annihilation contribution, C (d) 9 Ann . The decay width is largely dominated by weak annihilation for q 2 0.7GeV 2 as well as by the contribution from the narrow η resonance, which is not included in our description [4,5]. For q 2 above the φ peak the contribution from Y (d) becomes dominant. (We remind the reader that at the φ peak we start using the OPE expression for the weak annihilation The effect of weak annihilation for the higher-q 2 part of the spectrum is not completely negligible, though. To quantify it we plot, in figure 5, the ratio of the branching ratio obtained with the weak annihilation contribution to the branching ratio without this contribution For our central values, including the weak annihilation results in a decrease in the branching fraction of about 10%. But the magnitude and sign of the effect depends on interferences that are sensitive to the resonance phases. Within 68% CL, this means that, for the differential branching ratio, either a decrease of up to 30% or an increase by up to 15% is possible, locally. Since the contribution of weak annihilation can be sizeable, we include it in our final results. Finally, as we mentioned in section 3.1, our integrated branching ratios are compatible with other results found in the literature, although our values tend to be a bit larger. The results of figure 4 can also be compared locally. For example, comparing our differential branching ratio at q 2 = 0.5 GeV 2 with the results of ref. [4] we find, in our description, values ranging from 7.1 × 10 −8 to 2.8 × 10 −7 , while in ref. [4] one has values ranging from 3.6 × 10 −9 to 1.1 × 10 −7 , which is fully compatible, although our upper boundary is larger by a factor of almost 3. At q 2 = 2.5 GeV 2 , we find a total branching ratio between 1.4 × 10 −10 and 1.5 × 10 −9 whereas in ref. [4] the corresponding result is between 3 × 10 −11 and 6.2 × 10 −10 . The two ranges are again compatible, but the upper bound of our result is a factor of 2.4 times larger.

Calculation of the uncertainty band
We have identified a set of seven parameters playing a major role in the uncertainties entering our predictions. These consist of the phases of the resonances (φ ω and φ Φ ), three parameters entering the φ resonance structure and the ss resonance excitations (n Φ , σ 2 Φ and a Φ ) as well as two of the three renormalisation scales which enter the Wilson coefficients (µ c and µ W ). The dependence on these parameters is highly non-linear and there is no reason to assume a Gaussian distribution for the scale variation. Therefore, we use a Monte Carlo method in order to calculate the uncertainty band. In the case of the SM decay width, for a list of carefully chosen q 2 values within the kinematic range, we calculate the value of the observable N = 1000 times varying the parameters listed above following a specific distribution (uniform for the phases and the scales and Gaussian for the other parameters). We then use the median value as the central value, the uncertainty band being obtained looking for the interval that contains 68% of the values around the median. (The upper and lower uncertainty bands shown in figure 4 are the result of an interpolation.) Note that on calculating the observables in BSM Physics scenarios, it is again impractical to apply this MC technique to a large number of points in the plane of BSM Physics contributions to the Wilson coefficients, due to the extensive computing time that it requires. We therefore perform the full Monte-Carlo error propagation (N = 1000) on the desired observable for a small subset of 9 points spread equally across the plane and extrapolate between these points.

Constraints on the BSM Wilson coefficients
Before investigating the possible size of the effect of BSM physics on our observables, we first need to verify the existing constraints on the Wilson coefficients from the D 0 → + − and D + → π + + − branching ratios. Note that here we do not consider the interplay with JHEP04(2021)158 direct searches at the LHC, studied in ref. [50]. We further compare these constraints with the possible sizes of these Wilson coefficients in certain BSM models.

Constraints from upper limits on D 0 → + −
The measurement of the D 0 → + − branching ratio provides constraints on the BSM Wilson coefficients C 10 , C S and C P . The expression for this branching ratio is given by [4]: where in β we replace s by m 2 D , and Note that by simply adapting the result for B s → + − for D → + − , one neglects the possibility of additional long-distance contributions which might arise at the scale m c , thereby increasing the uncertainty on this constraint, an estimate of these effects can be found in ref. [51]. The constraints obtained from the D 0 → e + e − decay are much weaker than those obtained from D 0 → µ + µ − . For D 0 → µ + µ − , using the current best upper limit at 90%C.L., 6.2 × 10 −9 [52], we find: These constraints are similar to that obtained in ref. [4]. We note that the constraint on C 10 is ten times weaker than those on C P and C S . This is due to the fact that the axial-vector contribution to the branching ratio is helicity suppressed.

Constraints from upper limits on D + → π + µ + µ −
As for D 0 → + − , the constraints obtained from the c → ue + e − transition are much weaker than those obtained from the c → uµ + µ − transition. Compared to the leptonic D 0 → µ + µ − decay, which is only sensitive to effects in C ( ) 10 , C ( ) S and C ( ) P , the D + → π + µ + µ − decay is sensitive to a more diverse range of BSM Physics effects, namely C ( ) 10 , C ( ) S , C ( ) P , C T , C T 5 but also to possible BSM contributions to C 7 and C 9 . In the SM, these Wilson coefficients are multiplied by λ b as visible in eq. (4.16) and (2.7). We follow ref. [5], defining possible BSM contributions to C 7 and C 9 as  Table 6. Constraints on the maximal values of the Wilson coefficients from the 90% C.L. limit on the low-q 2 (Reg. I) and high-q 2 (Reg. II) bins of D + → π + µ + µ − branching ratio [38]. Table 7. Summary of the maximal value allowed by the current experimental limits on D 0 → + − and D + → π + + − at 90% C.L. if we assume the BSM Wilson coefficients are real and that the chirality-flipped Wilson coefficients are all equal to zero.

BSM Wilson coeff. Maximal values
where we assume that C BSM 7 and C BSM 9 are not necessarily subjected to the same CKM suppression as the SM contributions are.
The constraints on the BSM Wilson coefficients from the upper limits on the binned branching ratios of D + → π + µ + µ − are determined using the results at 90% C.L. given in table 6. We use the central value of the parameters given in table 3 and our expression for the branching ratio distribution given in eq. (4.20), integrated over two bins at low and high q 2 to compute the constraints on the Wilson coefficients C BSM 7 , C BSM 9 , C 10 , C S/P/T /T 5 in the two regions Reg. I and Reg. II respectively. These constraints can be found in table 6.
It is relevant to note that while the upper limits lie very close to the SM prediction in Reg. I compared to Reg. II; bounds on BSM coefficients are not stronger in this region. This can be understood by the fact that this region is dominated by weak annihilation, depending on C 1 and C 2 for which we do not consider BSM contributions.

Implications of the constraints for BSM models
If we assume that the BSM Wilson coefficients are real and that the chirality-flipped Wilson coefficients are all equal to zero, the maximal values allowed by the current experimental limits on B(D 0 → + − ) and B(D + → π + + − ) are summarised in table 7. Several models generating c → u + − transition have been studied in the past, e.g. Minimal Supersymmetric Standard Models [9,11,51,53], two Higgs doublet models [9], Little Higgs model [10,54], vector-like quark singlet [11] and leptoquarks [4,5]. It is of interest to compare the size of the Wilson coefficients allowed by the current experimental constraints on B(D 0 → µ + µ − ) and B(D + → π + µ + µ − ) summarised in table 7 to the values obtained in concrete BSM models.

JHEP04(2021)158
We will focus on scalar and vector leptoquarks (LQs), due both to the significance of these models in light of the B anomalies [55][56][57][58] and the interesting effects that they can generate in c → u transitions [4,5]. Following ref. [4], the contributions of the LQs to the Wilson coefficients C  (4.29) where i, j = L, R, and M is a generic scale of the leptoquarks. We observe that, beyond (semi)leptonic D decays, the only observables constraining the couplings λ u /c L/R are the branching ratios of K + → π + νν for S 1L and K 0 L → µµ for S 2R , V 2 and V 3 , as explained in ref. [4]. As these constraints are very strong, the Wilson coefficients C 9 , C 10 , C S and C P affected can be neglected. In table 8 we provide details of the remaining scalar and vector leptoquarks after taking into account these constraints.
The Scalar LQs S 1 with quantum numbers (3, 1, −1/3) and S 2 with (3, 2, −7/6) are interesting as they contribute to all the WCs we consider. For S 1 , we assume that the lefthanded up-quark-muon coupling can be neglected, i.e. λ uµ L ∼ 0, such that the combination λ uµ R λ cµ R controls C 9 = C 10 , and λ uµ R λ cµ L controls C S = C P = C T /2 = C T 5 /2. Analogously for S 2 , we assume that the right-handed up-quark-muon coupling can be neglected, i.e. λ uµ R ∼ 0, such that the combination λ uµ L λ cµ L controls C 9 = −C 10 , and λ uµ L λ cµ R controls C S = −C P = C T /2 = −C T 5 /2. As for the Vector LQs, these only contribute to C ( ) 9 and C ( ) 10 . After taking into account the constraints from kaon decays, the only vector LQs which can give rise to non-negligible Wilson coefficients areṼ 1 with quantum numbers (3, 1, −5/3) with C 9 = C 10 andṼ 2 with (3, 2, 1/6) with C 9 = −C 10 . Large values of un-primed C 9 and C 10 JHEP04(2021)158 cannot be generated. In addition to our model independent analysis, we will also study these models in the following subsection.

Results for physics beyond the SM
In figure 4 we see that the resonances have a larger effect on the differential Branching Ratio in the lower half of the q 2 region. In this section, we therefore choose to study the sensitivity of the observables defined in eqs. (4.7), (4.9) and (4.10) in the high q 2 region to possible BSM effects, integrating these observables over a large bin in q 2 in order to reduce the impact of the OPE assumptions adopted here for the weak annihilation corrections. We first do this in a model independent way, considering generic BSM contributions to pairs of Wilson coefficients, and then study the specific case of leptoquark models.

The model independent case
In order to study the potential size of these observables in BSM physics models we start by considering F H , as defined in eq. (4.11), integrated from q 2 min =1.8 GeV 2 to q 2 max =2.3 GeV 2 . The range is chosen for the following reasons: the OPE is valid for the case q 2 E π , Λ QCD ; on approaching the φ resonance the uncertainties increase. Note that for those cases where the BSM contributions are real and the tensor operators vanish, this is the only observable apart from the Branching Ratio. We start by looking at the effect of BSM contributions to C 9 and C 10 in figure 6. We include the constraints from B(D → µ + µ − ) and B(D + → π + µ + µ − ). We find that for the allowed ranges of C 9 and C 10 , values of F H almost reaching 0.04. This small value can be attributed to the fact that C 10 is restricted from B(D 0 → µ + µ − ). Further it turns out that the errors on F H in this plane are too large to be able to differentiate between different values of the Wilson coefficients. To be more precise, given the small variation in F H in the C 9 -C 10 plane, the uncertainty, of order 20% when not very close to the axis C 10 = 0, is too large to distinguish between different points in this plane. Note that we have checked that our results for F H and A FB are compatible with those in ref. [12] (taking into account the different definition of these observables), in spite of differences in the treatment of the decay amplitude, including the inclusion of weak annihilation and the treatment of the resonances.
Let us now consider a case where the uncertainties are small enough such that different values of the Wilson coefficients can be distinguished between. Since the largest contribution to the uncertainties come from the resonance structure affecting C 9 , it turns out that observables can be predicted more accurately as a function of BSM contributions to other Wilson coefficients. As a first example in figure 7 we vary C 10 and C P , again including constraints from B(D → µ + µ − ) and B(D + → π + µ + µ − ). The uncertainty bands are shown around the contours of constant values of F H . We see that, for this case, the uncertainty bands are narrow, particularly in the upper-left and lower-right quadrants of the plane. Note that while the bounds from the di-muon D decay are highly constraining, on introducing a right-handed Wilson coefficient, i.e. making the replacements C 10 → (C 10 +C 10 )/2 and C P → (C P + C P )/2, the results for D + → π + µ + µ − will be unchanged and the NP contribution to the branching ratio B(D → µ + µ − ) will vanish. For this reason in the plots JHEP04(2021)158 Now, allowing non-vanishing tensor Wilson coefficients we can consider the observable A FB , also defined in eq. (4.11), which vanishes in the SM and is helicity suppressed unless the combination of Wilson coefficients C S and C T or C P and C T 5 are non-zero. In figure 8, we study this observable in the C P and C T 5 plane, again including constraints from B(D → µ + µ − ) and B(D + → π + µ + µ − ). We show the uncertainty band around the contours of constant values of A FB . We find that the uncertainty on the results is small, particularly for small C P or C T 5 . This seems to be a promising observable for studies of the Wilson coefficients considered, with allowed values ranging from −0.35 to 0.35.
We have not yet contemplated the possibility of complex contributions to the Wilson coefficients. The ideal observable to look for CP violation is the CP asymmetry, A CP defined in eq. (4.10). However, on studying this CP asymmetry we find that the theoretical uncertainties are much larger that the size of the asymmetry. Our findings slightly differ from those of ref. [12], this can be understood for the following reasons: unlike in ref. [12] we choose to study the asymmetry integrated over a region in q 2 above the resonances, to reduce the dependence on the non-perturbative contribution, the size of the asymmetry is very much dependent on the strong phases arising due to the treatment of the resonances, which is very different in the two approaches. Note that the size of the asymmetry is larger at lower q 2 . We therefore choose not to show plots of the CP asymmetry as it would be difficult to disentangle any BSM Physics from the hadronic physics.  . Contours (blue, red, green) of the observable A FB , as defined in the text, in the C P -C T5 plane. The bands around the contours signify the theory uncertainty. The grey lines and shaded areas indicate the exclusion from the LHCb bounds on the branching ratios B(D 0 → µ + µ − ) [52] and B(D + → π + µ + µ − ) [38] as labelled.

JHEP04(2021)158
To summarise, in figure 6 we see the following: not only is the size of F H very small when BSM contributions only affect C 9 and C 10 , but the errors on the contours are so large that it is difficult to differentiate between them. On the other hand in figure 7 we see that larger values of the observable F H are achievable in the C 10 -C P plane, where the uncertainties are smaller and an experimental sensitivity below the O(10%) level could already provide very interesting information. As for figure 8, we see that probing A FB down to the 10% level would also result in important constraints on the combination of Wilson coefficients C P and C T 5 , where a non-zero value of A FB would be a decisive sign that either C P and C T 5 or C S and C T are both non-zero. On the other hand, in order to perform a precise fit to the Wilson coefficients an experimental sensitivity of O(1%) would be preferable.

The model dependent case
As mentioned earlier, we are interested in the possible affects of leptoquarks for the observables F H and A FB . In section 4.4 we discuss two scenarios for scalar and vector LQs which we would now like to investigate further. It turns out that the vector LQ scenarios,Ṽ 1 with quantum numbers (3, 1, −5/3) with C 9 = C 10 andṼ 2 with (3, 2, 1/6) with C 9 = −C 10 , were already investigated in figure 6, where the diagonal dashed lines represent these two scenarios as indicated. However, the uncertainties on F H are too large to allow information about the Wilson Coefficients C 9 and C 10 to be extracted. We therefore turn to the scalar LQs S 1 and S 2 defined in section 4.4.
We begin with the case of the leptoquark S 1 , in which there are two free parameters, C 9 = C 10 ≡ C 910 , and C S = C P = C T /2 = C T 5 /2 ≡ C ST . In figure 9 we study the observable A FB in the C 910 -C ST plane. We see that the uncertainty bands on the three contours shown, 0.025, 0.1 and 0.25, do not overlap, and the bands on the first two contours are particularly narrow. Note that in this scenario the primed and unprimed Wilson coefficients are fixed, such that the constraint coming from the branching ratio of D → µµ cannot be avoided. This imposes strict constraints on C ST , and only values lying approximately between −0.1 and 0.1 are viable, restricting A FB to a maximum value of 0.1. Fortunately the uncertainties in this region are small, such that different scenarios within this region would be theoretically distinguishable. Now coming to leptoquark scenario S 2 , in figure 10 we plot contours of constant F H for the process D + → π + e + e − in the in the C 910 -C ST plane, where here C 9 = −C 10 = C 910 , and C S = −C P = C T /2 = −C T 5 /2 = C ST . The decision to explore the case of electrons in the final state here can be explained in terms of eqs. (4.7) and (4.3), where we see that the dependence of the observable on the Wilson coefficients is much cleaner in the case of vanishing lepton masses. Replacing the muon mass by the electron mass not only makes it easier to interpret the results as the relationship between the observable and the Wilson coefficients is simpler, but also means that certain terms containing F V and hence the resonances are highly suppressed. This reduces the theoretical uncertainty on our predictions. From eqs. (4.3) and (4.9) we see that this reduction in uncertainty in D + → π + e + e − is the case for the observable F H but only true for A FB when C T is present. When C T is zero (and therefore C T 5 is non-zero in order for A FB to be non-JHEP04(2021)158   We conclude that in order to probe these two leptoquark scenarios, uncertainties on A FB (for D + → π + µ + µ − ) and F H (for D + → π + e + e − ) at the percent level would be required. We will comment further on the experimental feasibility of such measurements below.

Experimental prospects for D + → π + + −
The future prospects for D + → π + + − from LHCb, Belle-II and BES-III are very promising. Rough estimates find that the sensitivity of LHCb to the angular asymmetries and branching ratio for D + → π + can be found in table 9. Note that 50 fb −1 will only be collected with the Upgrade I at end of Run 4 in ∼2030. Then with the Upgrade II, by the end of Run 4 in ∼2038, 300 fb −1 should be collected. Belle-II is particularly suited to the electron channel, interesting as the results can be interpreted more easily in terms of the Wilson coefficients. By scaling BaBar results [59], assuming a similar efficiency, one finds the projected upper limit of the order 10 −8 on the branching ratio. 7 We note that this sensitivity to the branching ratio is of the same order as that of LHCb with 50 fb −1 , such that extrapolating we can say that a full angular analysis with percent-level sensitivity should also be possible at Belle-II. Finally, interesting results from BESIII in ref. [60] show that it will also have an important role to play for D + → π + .
Looking back at our results in the context of these estimated future sensitivities, it seems that by the time that LHCb has 50 fb −1 of data, the experimental errors will have far overtaken theoretical uncertainties, such that the various contours shown in our plots can easily be distinguished between. However, even before 50 fb −1 is collected, for F H and A FB which are highly suppressed in the SM a sensitivity at the 10% level could be sufficient to provide evidence for BSM physics. An experimental sensitivity at the 1% level would be enough to perform a precise fit to Wilson coefficients, the primary hindrance JHEP04(2021)158 being theoretical uncertainties. Note that for leptoquark scenarios, angular observables are constrained to be at most O(10%), such that probing these scenarios would require sensitivities at the O(1%). If Belle-II carries out the angular analysis for the electron case this will further provide important complementary information.

Consequences of updated theoretical framework on phenomenology
In this paper we advocate a dispersive approach to describe the resonance structure, employing the Shifman model and fitting to e + e − and τ data, as well as including weak annihilation corrections in both the low and high-q 2 regimes. Here we would like to address the question of how this formalism has an impact on the phenomenology. In considering this question, the following two factors are relevant: in our formalism, the residual effect of higher-order resonances are also taken into account, leading to larger uncertainties in the high-q 2 region; further in the high-q 2 region we include the weak annihilation corrections calculated in the OPE, and this constrains us to q 2 E π , Λ QCD . The impact on the phenomenology is as follows: • As a consequence of the range of validity of the OPE and the fact that uncertainties are still large above the φ, we propose to integrate observables in the range q 2 = [1.8, 2.3] GeV 2 .
• We find on allowing BSM contributions to certain Wilson coefficients (mainly C 9 ), the observables are subject to large theoretical uncertainties as a result of the resonance structure, such that, given the existing experimental bounds, determining the Wilson coefficient from the observable becomes very difficult. However the relation between the BSM contributions to certain pairs of Wilson coefficients and the integrated observables is clean, e.g. (C 10 , C P ) and F H or (C P , C T 5 ) and A FB .
• Within our framework we find that the uncertainty on A CP is very large (an order of magnitude large than the central value), due to the large phase uncertainty. While the strong phase of the SM contribution is unknown, differentiating between possible BSM phases would be difficult. However, a measurable A CP would probably be a sign of complex BSM couplings.
Therefore a number of differences emerge in the phenomenology as a result of the theoretical framework we adopt.

Conclusions
Charm physics is gaining increased interest, due to the large numbers of charm mesons produced at LHCb and Belle-II. The decay D + → π + + − , being a neutral current process is a particularly interesting probe of BSM physics, especially in light of potential links with the Flavour Anomalies observed in b → s transitions. We have conducted a comprehensive analysis of this decay, tackling the increased complexity compared to its b → s counterpart resulting from the fact that the expansion in Λ QCD /m c is less effective than Λ QCD /m b , and that resonances affect larger regions in JHEP04(2021)158 phase space. Our calculation of D + → π + + − includes weak annihilation within QCDf at low q 2 and also in the OPE at high q 2 . To the best of our knowledge, this is the first time both QCDf corrections have been analysed in detail in a phenomenological analysis of D + → π + + − , and that the OPE-approach has been applied to the weak annihilation contribution for these decays. We were motivated to include these contributions due to the fact the strong hierarchy λ b λ d means that they are not suppressed as in the b → s case, see table 5. Our work benefits from recent Lattice QCD results for the form factors, as well as recent calculations of the Wilson coefficients at next-to leading order. Since the effect of resonances on the SM prediction extends over the entire kinematic regime, it is important that these are modelled carefully. We have employed a novel method, fitting the Shifman model for the resonance structure to e + e − → (hadrons) and τ → (hadrons) + ν τ data, and further using the experimental value for the D + → π + R (R → + − ), with R = ρ, ω, and φ, branching ratios. We have performed a thorough and conservative error analysis involving Monte Carlo error propagation, taking into account the uncertainties from the resonance model, which dominate, as well as from the renormalisation scales.
Our prediction for the differential branching ratio as a function of q 2 , in the entire kinematic region, along with the uncertainty band is provided in figure 4. This can serve as a conservative prediction of the branching ratio and uncertainties throughout the phase space, which would be an important input for backgrounds in experimental searches. An important result is that, due to the large weak-annihilation contribution as well as the residual resonance contributions away from the resonance peaks, the integrated non-resonant branching ratios could be of the order of 10 −9 . Note that the sensitivity of LHCb to the branching ratio will be ∼ 10 −8 with 50 fb −1 and ∼ 10 −9 with 200 fb −1 [61].
We studied CP conserving observables, focusing on the observables F H and A FB . These are of interest as they are close to or equal to zero in the SM but can be strongly enhanced by BSM physics; to be more precise any contributions arising from the Wilson coefficients which are non-zero in the SM are helicity suppressed. We focus on the regions at large q 2 where the uncertainty from the resonances is best under control, and the experimental constraints are weakest. This calls for the OPE estimate for the weak annihilation contribution, such that we integrate over a certain range in q 2 in order to obtain more accurate predictions.
We considered effects of both model-independent BSM and a specific model, i.e. leptoquarks, on F H and A FB . We observe that BSM contributions to the vector Wilson coefficient C 9 are subject to large theoretical uncertainties, and do not contribute significantly to the observables F H and A FB . We therefore focused on combinations C 10 -C P for F H in figure 7 and C P -C T 5 for A FB in figure 8, where the uncertainties are particularly small and different values of the Wilson coefficients could be distinguished between (theoretically) in the allowed region of parameter space. Of the leptoquark scenarios considered, we find that for vector leptoquarks A FB vanishes and F H suffers large theoretical uncertainties. On the other hand, in scalar leptoquark scenarios, while being strongly constrained by D meson decays to di-muons, these observables, in particular A FB , can be precisely predicted in the remaining parameter space.

JHEP04(2021)158
We therefore look forward to the upcoming results for D + → π + + − from LHCb, Belle-II and BES-III (see table 9), from which we will obtain much-improved bounds on the Wilson coefficients and the models discussed. We urge the experimental collaborations to measure the observables advocated in this paper F H and A FB in the range q 2 from ∼ 1.8 to 2.3 GeV 2 , both for D + → π + µ + µ − and D + → π + e + e − , stressing that an experimental sensitivity of O(10%) to these observables would already provide evidence for BSM physics in certain scenarios, and furthermore a sensitivity of O(1%) would make a precise fit to the Wilson coefficients achievable.

JHEP04(2021)158
The moments are required at the scale µ c , obtained via [66] c i (µ) = L γc i /β 0 c i (1 GeV), (A.7) where L = α s (µ)/α s (1 GeV) and γ c i is given by The LCDA of the D meson is defined by (see e.g. ref. [21]) where v is the velocity of the D meson. We choose to parametrise the D meson LCDA using a simple exponential model [67]: where ω 0 is the sole input parameter.

B Calculation of the Wilson coefficients
Here we recall the two-step running of the Wilson coefficients for the c → u + − transition [19]. These coefficients are computed at the electroweak scale µ W ∼ M W and then run down to the typical mass scale of the decay under consideration, here µ c ∼ m c . This running involves an intermediate scale µ b ∼ m b , where the bottom quark is "integrated out". The matching coefficients and anomalous dimension matrices are taken to the required order by generalising and extending results from b → d/s transitions [20]. The Wilson coefficients for the c → u + − transition at NNLL accuracy were calculated for the first time in ref. [19]. The full procedure to calculate the Wilson coefficients can be broken down in the following steps: • We determine the Wilson coefficients via matching at µ ∼ µ W to second order in a s (µ W ), of which only C 1/2 receive non-zero contributions [19].
• Secondly, C 1 and C 2 are run down to the scale µ b following ref. [20]. The full 8 × 8 anomalous dimension matrices, γ i where i = 0-2, required for the RGE running can be broken down in the following way where Q (i) 1 are the 6 × 6 three-loop anomalous dimension matrices describing the mixing of the four quark operators C 1−6 , they are taken from [68].
• We then perform the matching from the five-quark (n f = 5) to the four-quark (n f = 4) effective field theory at the scale µ b following ref. [19]. The matching matrix R is different from the unit matrix because the operators O b 1/2 are absent below the b-quark threshold. C 3/9 receive non-zero contributions only from the matching of the five-flavor effective theory above the scale µ b to the four-flavor EFT below that scale and from the mixing of O 1/2 into O 3/9 . The R matrix is given by: where the non-zero elements of R (1) ij can be found in ref. [19].
• For C 9 , the running procedure follows ref. [19] and in turn ref. [20], where the sixdimensional list of 4-quark operators at the weak scale, C (µ W ), is evolved down to the scale µ b using the six-dimensional anomalous dimension matrices Q (i) 1 , and matched using the six-dimensional R 6 .
• Below µ b the contribution of these operators to O 9 are calculated via the 1 × 6 matrix Here U n f =4 6 (µ c , µ b ) and R 6 are 6 × 6 sub-matrices of U n f =4 (µ c , µ b ) and R are 6 × 6 mentioned earlier for the C 1 -C 8 case. The solution of eq. (B.3) can be found in appendix C of [20]. The leading order contribution to the initial condition for C 9 at the scale m b is given by The full procedure for C 1 -C 8 can be summarised by the following equation: and for C 9 by Note that C 10 does not mix under renormalisation and thus is zero at all scales to leading order in the 1/M W expansion. Moreover, the framework introduced above does JHEP04(2021)158 not compute the coefficients C 7/8 but the renormalisation-scheme independent effective ones, defined as: with y (7) = Q d 0, 0, 1, 4 3 , 20, 80 3 , Q d = −1/3 and y (8) = 0, 0, 1, − 1 6 , 20, − 10 3 .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.