The New Physics Reach of Null Tests with $D \to \pi \ell \ell$ and $D_s \to K \ell \ell $ Decays

$|\Delta c|=|\Delta u|=1$ processes are unique probes of flavor physics in the up-sector within and beyond the Standard Model (SM). SM tests with rare semileptonic charm meson decays are based on an approximate CP--symmetry, a superior GIM--mechanism, angular distributions, lepton-universality and lepton flavor conservation. We analyze the complete set of null test observables in $D \to \pi \ell \ell^{(\prime)}$ and $D_s \to K \ell \ell^{(\prime)}$ decays, $\ell^{(\prime)}=e, \mu$, and find large room for new physics safely above the SM contribution. We identify signatures of supersymmetry, leptoquarks and anomaly-free $U(1)^\prime$-models with generation-dependent charges, for which we provide explicit examples. $Z^\prime$-effects in $c \to u \ell \ell^{(\prime)}$ transitions can be sizable if both left-handed and right-handed couplings to quarks are present.

The plan of the paper is as follows: In Sec. II we review the predictions in the SM, taking into account improved form factor computations. In Sec. III constraints on Wilson coefficients are obtained, and null test observables for D → π ( ) and D s → K ( ) decays, ( ) = e, µ, are analyzed.
Signatures of concrete BSM models in c → u ( ) transitions are worked out in Sec. IV. In Sec. V we conclude. In App. A we provide details on the D → π form factors. Constraints from D 0 − D 0 mixing are discussed in App. B. In App. C we give details on and explicit examples of anomaly-free Z -models with generation-dependent charges.

II. STANDARD MODEL PREDICTIONS
We discuss exclusive rare charm meson decays in the SM. In Sec. II A we review the SM contribution in an effective field theory framework at the charm mass scale. We discuss resonant contributions from intermediate mesons in Sec. II B. In Sec. II C we present the D → π + − and D s → K + − differential branching ratios and discuss SM uncertainties.
A. An effective field theory approach to charm physics Rare c → u + − processes can be described by the effective Hamiltonian, i=7,9,10,S,P where the dimension 6 operators which can receive BSM contributions are defined as follows: The operators O i are obtained from the O i by interchanging left-handed (L) and right-handed (R) chiral fields, L ↔ R. As in B-decays, in the SM, C S,P,T,T 5 = 0, and all C i are negligible. Unlike in B-physics, the GIM-mechanism kills penguin contributions to C 7,9,10 at the W -scale µ W in H eff .
Therefore, the SM contributions to O 7,9,10 are induced by the charged-current operators by renormalization group running to the charm scale µ c . After two-step matching at µ W and the b-mass scale, see [5,18] for details, the effective coefficients C eff 7/9 at µ c can be estimated as [5] C eff 7 (q 2 ≈ 0) −0.0011 − 0.0041 i , C eff 9 (q 2 ) −0.021[V * cd V ud L(q 2 , m d , µ c ) + V * cs V us L(q 2 , m s , µ c )] , where L(q 2 , m q , µ c ) = with x = 4m 2 q /q 2 and the dilepton invariant mass squared q 2 . Taking µ c = m c = 1.275 GeV one obtains |C eff 9 | 0.01. Im[C eff 7 ] increases from -0.004 at q 2 = 0 to -0.001 at high q 2 at NNLL order [18]. Importantly, C SM 10 = 0, which implies that at the charm scale effects from the V-A structure of the weak interaction are shut off. Numerically, the short-distance SM contributions are negligible in the D → P + − decay rates compared to the resonance-induced effects, discussed in Sec. II B.  The dominant process inducing the O 9,P operators is D → P γ * with γ * → M → + − , which can be parametrized by a phenomenological shape, 0.27, consistent with the isospin limit. The strong phases δ ρ, φ, η remain undetermined by this procedure and are varied within −π and π in the numerical analyses. Fig. 1 illustrates the impact of the uncertainties from the unknown phases on real and imaginary parts of C R 9 in the high q 2 -region ( q 2 ≥ 1.25 GeV) for D → πµ + µ − decays. The phases δ ρ, φ give rise to huge uncertainties, shown by the yellow wider bands. As a result, even the sign of C R 9 cannot be predicted. Once the phases are fixed, the residual uncertainties from a ρ , a φ , shown by the blue bands, are small. In addition, they could be reduced by improved measurements of D → P M and M → + − branching ratios. Darker shaded solid lines correspond to central values of input with strong phases fixed to δ ρ = 0, δ φ = π, consistent with the SU (3) F limit δ ρ − δ φ = π.
FIG. 1: Resonance contributions to the real part (plot to the left) and the imaginary part (plot to the right) of C R 9 in the high q 2 -region using the D + → π + parameters from Tab. I. The wider yellow bands correspond to the uncertainties from strong phases δ ρ , δ φ varied within [−π, +π], whereas the smaller blue bands correspond to fixed strong phases δ ρ = 0, δ φ = π. Additional uncertainties arise from a ρ , a φ , and are included in the plots. The dashed curves illustrate the uncertainties in the case of D 0 → π 0 parameters.

C. Phenomenological q 2 -distributions
The differential decay distribution of D → P + − can be written as [19] To ease notation, here and in the following, all Wilson coefficients C i , with the exception of the tensor ones, are understood as FIG. 2: The differential SM branching fractions of D + → π + µ + µ − (plot to the left) and D + s → K + µ + µ − (plot to the right) decays. Yellow (blue) bands show pure resonant (short-distance) contributions. The band widths represent theoretical uncertainties of hadronic form factors, resonance parameters and µ c . Darker shaded thin curves represent all parameters set to their central values including δ ρ = 0 and δ φ = π (solid) and δ ρ = δ φ = 0 (dashed). For the experimental limit by LHCb [20] see Eq. (9).
We neglect the up-quark mass. The corresponding D s → K + − decay distribution can be obtained by replacing m D → m Ds , using D s → K form factors f +,0,T and resonance parameters a ρ, φ, η, η from Tab. I. A detailed discussion of the form factors can be found in App. A.
In Fig. 2 we show the pure non-resonant q 2 -spectra determined by the SM short-distance C eff 7,9 coefficients (lower, blue curves) and the pure resonant ones determined by C R 9,P (yellow, wiggly curves). The non-resonant contribution is suppressed by several orders of magnitude with respect to the resonant one and negligible in the SM. The band width represents theoretical uncertainties of the form factors, the resonance parameters a M with corresponding phases δ M , see Eq. (6), and the scale µ c . The largest uncertainty stems from the unknown phases δ M , an effect which is pronounced through resonance interference.
The uncertainties from the strong phases δ M are highly dependent on the values of the a M . One can infer from Fig. 2 that the D-mode has huge uncertainties at low q 2 and around the ρ/ω peak, and it has much smaller uncertainties and more stable behavior above the φ mass, while for the D s -mode the situation is the opposite. The reason for this difference is the numerical values and hierarchy of a ρ and a φ (see Tab. I). Although the tiny uncertainty band in the low q 2 -region of the D s -mode looks promising for testing new physics, we stress that the uncertainty band width can change with future measurements of the D s → KM and M → + − branching ratios.

III. BEYOND THE STANDARD MODEL
BSM effects in the c → u ( ) transition are discussed in a model-independent way. In Sec. III A we update the constraints on the Wilson coefficients coming from D + → π + µ + µ − and D 0 → µ + µ − branching ratios and discuss the BSM sensitivity of B(D → P + − ) at high q 2 . In Sec. III B lepton universality tests are considered. Null tests of the SM based on the angular distributions are discussed in Sec. III C. CP-asymmetries and their sensitivity to BSM benchmark values is analyzed in Sec III D. In Sec. III E lepton flavor violating c → ue ± µ ∓ decays are discussed.

A. Updated constraints on Wilson coefficients
Using the experimental limits on the branching fraction of D + → π + µ + µ − in high and full q 2 -regions at 90% CL [20], and neglecting the SM contributions, we obtain the following constraints on the BSM Wilson coefficients in the high q 2 -region, and in the full q 2 -region, We find good agreement between Eq. (10), in the limit when only one C i is present, and Ref. [6].
The limits (10), (11) are also consistent with the fact that in all D → P + − distributions the leptonic vector current is probed only through the combination with γ = 2mc m D +m P f T f + , which numerically is around one in the full q 2 -region, γ ≈ 1, as shown in App. A. Present bounds on dipole operators from D → ργ [21,22] are in agreement with (10), (11).
The numerical coefficients in Eqs. (10) and (11) are smaller than the corresponding ones in Eqs. (29) and (30) of Ref. [5]. This difference is caused by the D → π form factors, for which we use new lattice results [11,12] (for details, see App. A). This in particular affects f T and hence the available space for NP in C T,T 5 and C ( ) 7 . Note that in Eqs. (10) and (11) we neglected the contributions from C R 9, P . Taking them into account, we obtain one additional constant term plus several interference terms proportional to Re[C i ] and Im[C i ] on the left-hand side of Eqs. (10) and (11). When varying strong phases, the constant term turns out to be smaller than 0.1. The largest numerical coefficients of the interference terms can vary in the full q 2 -region within [-0.1, 0.1] and [-0.2, 0.2] for C 9 and C 7 , respectively, whereas at high q 2 the corresponding ranges turn out to be larger, about a factor ∼ 2 to 5. Despite being less tight we use in the following the constraint from the D + → π + µ + µ − branching fraction in the full q 2 -region as it is more stable than the high q 2 one with respect to unknown resonance parameters .
One may wonder whether Wilson coefficients as large as order one of |∆c| = |∆u| = 1 fourfermion operators are constrained from high-p T searches. While a dedicated analysis for upsector flavor changing neutral currents (FCNCs) is currently not available, constraints derived from 36.1 fb −1 ATLAS data on semileptonic operators with two charm singlets or two second generation quark doublets [23] are close to probing the range allowed by charm decays into muons, and are somewhat stronger for electrons.
Further constraints on C ( ) 10,S,P can be obtained from the D 0 → + − branching ratio The upper limit B(D 0 → µ + µ − ) < 6.2 × 10 −9 at 90% CL [24] gives In addition, tensor Wilson coefficients are constrained by the leptonic decays as they induce scalar ones from renormalization group running [25]. In the subsequent analysis we use |C S,P | 0.1 and |C i | 0.5 for all other NP Wilson coefficients as benchmark values, except in Sec. III D on CPviolation effects, which are subject to model-dependent, generically stronger D-mixing constraints.
In the following C i(j) stands for "C i or C j ".  Fig. 3. In the third to sixth column, upper entries correspond to NP-only branching ratios while for the lower entries the resonance contributions are taken into account.

B. Testing lepton universality
Lepton universality can be probed in semileptonic decays with the ratios [6,7] Here, q 2 min (q 2 max ) denotes the lower (upper) dilepton mass cut; to ensure maximal cancellation of hadronic uncertainties and hence a controlled SM prediction near unity it is crucial that the cuts for electron and muon modes are identical [26]. Assuming that NP contributes only to the muon mode, we present in Tabs. III and IV the predicted ranges of R D π and R Ds K , respectively, in the full, ×10 −9 FIG. 3: The differential branching fractions in the SM and in three BSM scenarios in the high q 2 -region for D + → π + µ + µ − (plots to the left) and D + s → K + µ + µ − (plots to the right). Orange bands correspond to the pure resonant contributions. The band width represents theoretical uncertainties, where we distinguish between pure BSM contributions (upper plots) and BSM benchmarks including interference with the SM resonance contributions, such that the variation of the strong phases also effects the BSM predictions (lower plots). Predictions for BSM coefficients C low and high q 2 -integrated intervals. As for the resonance parametrization, we use the same C 9,P for electron and muon modes. Due to poor knowledge of the resonances' behavior, the predictions for the low q 2 -region ( q 2 ∈ [250, 525] MeV) have sizable uncertainties and we only give the order of magnitude of the largest values found. The main uncertainty comes from the phases δ ρ, φ, η which are varied within −π and π, while the uncertainties due to hadronic form factors are of the order of few percent. Integration over the full q 2 -interval gives ratios R D π and R Ds K which are insensitive to NP. On the other hand, in the high q 2 -region NP can induce significant effects. BSM effects in the low q 2 -region can be huge in the D-mode, however, there are large uncertainties.  (15) in the SM and in NP scenarios for different q 2 -bins. Ranges correspond to uncertainties from form factors and resonance parameters. Due to large uncertainties at low q 2 in some cases only the order of magnitude of the largest values is given.

C. Angular observables
The double differential distribution of D → P + − decays [19] where θ denotes the angle between the − -momentum and the P -momentum in the dilepton rest frame, offers the measurement of two angular observables: The lepton forward-backward asymmetry, and the "flat" term, both with normalization to the q 2 -integrated decay rate, with integration limits depending on the q 2 -bin. Since the scalar and pseudotensor operators are and both are even smaller at high q 2 , see Fig. 4. F H (D (s) → π(K)e + e − ) in the SM is even further suppressed. We learn that F H is

FIG. 5:
The forward-backward asymmetry A FB (upper plots) and F H (lower plots) in the high-q 2 region in different BSM scenarios for D + → π + µ + µ − (plots to the left) and D + s → K + µ + µ − (plots to the right).
yet another a very promising null test of the SM, sensitive to NP in (pseudo)-scalar and tensor operators.
In Fig. 5 we show A FB and F H in different NP scenarios for D + → π + µ + µ − (plots to the left) and D + s → K + µ + µ − (plots to the right) at high q 2 , that is, q 2 min = (1.25 GeV) 2 and q 2 max = (m D −m P ) 2 . One can see that A FB is mostly sensitive to the combinations of tensor and (pseudo)scalar operators, while it is significantly suppressed if only scalar or pseudotensor structure is present. The latter effect comes from interference terms, C R 9 C * S and C R P C * T 5 , which are suppressed by the light lepton mass and small η-resonance contribution to C R P , respectively. The flat term turns out to be sensitive to the (pseudo)scalar and especially to the (pseudo)tensor operators. The effect of simultaneous presence of scalar and tensor contributions (in cyan) is essentially indistinguishable from the one induced by pure (pseudo)tensor (in red) structure. However, large effects in both F H and A FB would be a signal of tensor and scalar nature of NP. The normalization of angular observables has significant impact on the shape of the distributions. Choosing dΓ/dq 2 , as in [27], instead of Γ (19), results in markedly different behavior for F H near the endpoint q 2 → (m D − m P ) 2 , where cancellations can occur.

D. CP-asymmetry
Another promising observable to probe NP is the CP-asymmetry, defined as [3,5] where Γ denotes the decay rate of the CP-conjugated mode, defined with q 2 -bin dependence as Γ in Eq. (19). The difference of the differential rates can be written as

(21)
A SM CP is determined by the first term in Eq. (21) and remains tiny due to small phases of the CKM factors in C 9 , see Eq. (4). Existing bounds from A CP (D 0 → ρ 0 γ) = 0.056±0.152±0.006 [21] do not provide further constraints on C ( ) 7 beyond the ones from branching ratio measurements [22]. Naïve T-odd CP-asymmetries from angular distributions in D → ππµ + µ − decays can probe CP-phases even for vanishing strong phases. First experimental studies are at the O(10 %) level [28] which is about where sensitivity to BSM physics starts [7].
For a sizable (T-even) CP-asymmetry, enhanced strong phases and therefore resonance effects are instrumental [3]. In Fig. 6 we show A CP in the φ-region (upper plots) and at high q 2 (lower plots) for C 9 = 0.1 exp(i π/4) and different values of δ φ . The band width corresponds to the 1σ uncertainty due to m c , form factors and resonance parameters (δ ρ,η , varied within −π and π).
NP effects in C T and C P,10 are suppressed by the light lepton mass and the completely negligible respectively. We learn from Fig. 6 that irrespective of the value of δ φ , sizable BSM effects occur making this observable promising for NP searches. Except for D s → K + − at high q 2 , A CP has rather small uncertainties and is useful to extract strong and weak parameters.
Note that A CP can change its sign around q 2 ∼ m 2 φ and hence, to avoid a vanishing integrated asymmetry, binning is required. Due to the way Eq. (12) with which leptonic vector contributions enter D → P + − decays, A CP has similar sensitivity to BSM effects from C ( ) 7 than from C ( ) 9 . Similar to the behavior observed in the angular observables F H and A FB , BSM effects in A CP in the D s → K + − mode are enhanced relative to the D → π + − one, due to the smaller decay rate of the former, caused predominantly by kinematics.

E. Lepton flavor violation
To discuss LFV in c → u − + ( = ) decays we introduce the following effective Hamiltonian, where the K with other operators from Eq. (2) defined in similar way by changing flavor in lepton currents. Note that there is no O ( ) 7 contribution since the photon does not couple to different lepton flavors. The differential distribution for the LFV decays D → P e ± µ ∓ is given as where K i = K Similarly to Eq. (11), we obtain the following constraints on the LFV Wilson coefficients using the 90% C.L. upper limits B(D + → π + e + µ − ) < 2.9×10 −6 and B(D + → π + e − µ + ) < 3.6×10 −6 [29]: Tighter constraints on K ( ) 9,10,S,P can be obtained from data on leptonic decays, [30], we obtain K S − K S ± 0.04 K 9 − K 9 2 + K P − K P + 0.04 K 10 − K 10 2 0.01 . The differential decay rate of D + → π + e ± µ ∓ (plot to the left) and D + s → K + e ± µ ∓ (plot to the right) decays in different LFV-BSM scenarios.
In Fig. 7 we present the differential branching fractions for LFV decays D + → π + e ± µ ∓ and D + s → K + e ± µ ∓ for various NP Wilson coefficients allowed by Eqs. (25), (27). Since the resonant contributions are absent in the LFV processes, the uncertainties are due to the form factors and m c , making the band widths in Fig. 7 significantly smaller than in Fig. 2. The difference in the shapes for vector and (pseudo)tensor/scalar structure allow to experimentally distinguish different operators and BSM scenarios.

A. Leptoquark signatures
In order to bypass the constraints from the kaon sector in a most straightforward manner, we consider only the scalar leptoquarks S 1(2) with right(left)-handed couplings and the vector ones V 1,2 , e.g. [5,31]. The subscript "1" denotes SU (2) L -singlets, whereas "2" denotes doublets. This precludes any (pseudo-)scalar or tensor operators, which otherwise would be induced by S 1,2 , as well as (axial-)vector ones with left-handed quarks. We stress, however, that even small Wilson coefficients involving doublet quarks can signal NP in A CP , if CP-violating [5]. The interaction Lagrangian contributing to c → u − + processes then reads [32] L where L L denote lepton doublet and l R , u R lepton and up-type quark singlets; i, j are the generation indices. Hypercharge-assignments of the leptoquarks can be read-off from (29); all leptoquarks are color-triplets. Charm signatures of S 2 have been studied in [6,27] and V 1 in [6]. Only S operators can be written as, e.g. [5], where M X , X = S 1,2 , V 1,2 denotes the leptoquark mass. For the singlets (S 1 , V 1 ) K 9 = K 10 , while for the doublets (S 2 , V 2 ) K 9 = −K 10 . The corresponding lepton flavor conserving contributions to C 9,10 can be obtained from Eq. (30) for = . Note that scenarios with larger coupling λ cµ S 2 3.5 [6,27], that could induce significant C S,P,T,T 5 together with even a suppressed coupling to doublet quarks, are excluded by high-p T studies, λ cµ S 2 0.4M S 2 /TeV [23]. Taking this constraint into account leptoquark effects in D → πµ + µ − decays in A FB (q 2 ) do not exceed permille-level. On the other hand, F H (q 2 ) 0.1 in the high q 2 -region which is small, however, above the SM value shown in Fig. 4.
Using Eqs. (13), (26) and neglecting the SM contribution, we obtain the following constraints from the upper limits on B(D 0 → µ + µ − ) [24] and B(D 0 → e ± µ ∓ ) [30], Constraints in a similar ballpark are obtained from B(D → πµ + µ − ) using Eq. (11), Leptoquark contributions to D 0 − D 0 mixing are induced by one-loop box diagrams. Yet, data on ∆m D 0 results in a somewhat more stringent constraint than the rare decays, [5], which can be eased for the individual, lepton-specific terms due to cancellations. Note that the corresponding mixing bound on the imaginary part is about a factor 0.2 stronger, see App. B. Also dipole operators are induced by leptoquark-lepton loops [22]. As we do not consider leptoquarks with couplings to left-handed quarks, no chirality-flipping τ -loops are available, and the resulting C ( ) 7 is below O(10 −3 ) for S 1,2 , but could reach O(10 −2 ) for V 1,2 , see [22] for details. Leptoquark signals can hence be observed in D → P ( ) decays in lepton-universality tests, LFV searches and, if CP-violating, in A CP . LFV rates from S 1,2 (cyan) together with SUSY predictions (brown, magenta) are shown in Fig. 8.

B. Charm reach of SUSY models
Supersymmetric SM extensions offer two ways for BSM signals in rare charm decays. One is through enhanced dipole couplings C ( ) 7 from scalar quark mixing; corresponding contributions to c → uγ decays have been studied recently in [22], to which we refer for details. The sensitivity of rare semileptonic D → P + − decays to C ( ) 7 can be inferred from Fig. 3 for the branching ratio and in Fig. 6 for the CP-asymmetry using Eq. (12). Note that the D → P + − distributions are sensitive to the sum of dipole coefficients, C 7 + C 7 , only. This ambiguity can be resolved with polarization studies in D → V γ, V = ρ, K * , φ, K 1 [33,34], which probe the fraction of righthanded to left-handed photons. While the sensitivity window to NP in the branching ratio is rather
small, there is large room in A CP for CP-violating dipole contributions from supersymmetric flavor violation.
Another possibility to probe SUSY in c → u FCNCs is with R-parity violating terms λ LQD, which induce (axial)-vector couplings with C 9 = −C 10 through tree level exchange of a scalar partner of a singlet down quark, see [1] for details; also [35]. However, this unavoidably involves doublet up-type quarks, hence is subject to constraints from rare kaon decays, specifically K → πνν.
This situation resembles the one of the coupling to left-handed leptons of the scalar leptoquark S 1 , which has the same quantum numbers as a singlet down quark. In SUSY additionally even stronger constraints from K L → e ± µ ∓ apply if the squark doublet masses are smaller or close to the downsinglet ones. We illustrate maximal LFV rates in Fig. 8 for K 9 = −K 10 0.009 (benchmark I) from K → πνν [36] for decoupled doublet squarks and K 9 = −K 10 0.001 (benchmark II) from K L → e ± µ ∓ assuming degenerate squarks [14,37].
To summarize, the strongest opportunity for SUSY to signal NP in D → P + − decays is with CP-violating contributions in A CP , and in LFV decay rates in R-parity violating models.

C. Flavorful Z -models
The effective Z -interaction Hamiltonian part for c → u − + processes can be written as The following Wilson coefficients are induced at tree-level Using Eqs. (13), (26) and (11), respectively, we obtain the following constraints on the Z -model where M Z denotes the mass of the Z . Stronger bounds, however, arise from the Z tree level contribution to D 0 − D 0 mixing, see App. B for details (Here and in the following the superscript 'uc' has been dropped.), In the following we show that these conditions can be met in flavorful Z -models by flavor rotations without introducing unnatural hierarchies. Here, FCNC c → u-transitions arise from non-universal charges, denoted by F . Specifically, inducing g L (g R ) requires the charges of the doublet (singlet) charm and up-quarks to be different from each other i.e., ∆F L(R) = 0. Another requisite ingredient is flavor mixing. Four such rotations between flavor and mass bases exist, those for up-singlets, U u , for down-singlets U d , and those for up-and downdoublets, V u and V d , respectively. U d , U u , V u and V d are unitary matrices; in absence of a theory of flavor they are unconstrained, except for V † u V d = V CKM . Here we consider rotations residing in the up-sector, U d = 1, V d = 1, hence V u = V † CKM , as flavor rotations in the down sector are subject to strong constraints from the rare kaon decays, e.g. [38]. Small mixings, (V d ) 12 /(V u ) 12 10 −3 , however, would still be consistent with sizable effects in charm. To understand the interplay between rotations and charges in order to satisfy Eq. (39), we simplify the analysis by assuming the third generations to be sufficiently decoupled. This allows us to work with 2 by 2 orthogonal matrices, parametrized by a single angle each, θ u and Φ u for the up-singlets and -doublets, respectively. Then, with g 4 the U (1) gauge coupling strength and, since V u = V † CKM , we obtain where λ ∼ 0.2 denotes the Wolfenstein parameter.
We obtain g R /g L = X for and g R /g L = 1/X for The former case, in the following termed right-handed (RH)-dominated, requires some mild hierarchy between the left versus right charges, whereas the latter case, termed left-handed (LH)dominated, can be accommodated with mixing alone θ u ∼ 10 −2 . In either case, the ratio of leftand right-handed couplings is fixed, C 9/10 C 9/10 , K ( ) 9/10 K ( ) ( ) 9/10 ∼ X (LH) or ∼ 1/X (RH) .
In Tab. V we present sample scenarios of Z -models with sizable BSM coefficients, which obey We work out constraints on the parameters of the concrete Z -models. The constraints from D → πµ + µ − branching ratios in Eq. (11) can be written as In Fig. 9 the BSM coupling g 4 is shown as a function of the Z mass. Each line corresponds to an upper limit in the scenarios from Tab. V with one specific choice for the charges F L 1 (F L 2 ) and F e 1 (F e 2 ) for electrons (muons). Constraints for RH cases are stronger than for corresponding LH ones in Eq. (47). The black region shows the exclusion region due to resonance searches in dilepton searches of M Z > 4.5 TeV [14]. The true lower bound for the Z mass is different for every solution due to the different quark charge assignments and overall coupling strength and in general different for electrons and muons. Therefore, part of the region g 4 < 0.5 and M Z < 4.5 TeV might still be viable or constrained by other searches [39].
Z -contributions to charged LFV arise by a similar misalignment between flavor and mass bases as for the quarks. Likewise, there is no charged LFV if the PMNS-matrix is only due to rotations in the neutrino sector. In the presence of charged lepton rotations (i) τ → (µ, e) with = e, µ , as well as (ii) µ → eee and µ → eγ are the most stringent ones [14]. Charm constraints continue to be more restrictive for the models in Tab Concrete models with potentially large (pseudo)scalar and tensor contributions, which could be signaled with angular observables, are scalar leptoquarks with representations S 1 and S 2 . However, it is exactly those representations that allow for chirality flipping operators that are unavoidably subject to tight constraints from the kaon sector. The largest effects in charm decays within leptoquark models are therefore from couplings to right-handed up-quarks, inducing primed (axial-)vector operators. Both right-handed and left-handed (axial-)vector couplings C   [11,12] and f + from HFLAV [40] (gray band) with 1σ theoretical uncertainties.

Appendix A: Form Factors
The hadronic matrix elements for D → P transitions are parametrized in terms of three form factors, f + ,0 ,T (q 2 ), defined as 3 In this work we use the D → π form factors computed recently by Lubicz et al. [11,12] using lattice QCD (LQCD), parametrized in terms of the z-expansion as, i = +, 0, T , where The numerical values of f i (0), c i and P i parameters together with their uncertainties and covariance matrices are given in [11,12] 4 . In particular, P 0, T and c T have large uncertainties, about 40% 3 The normalization of the tensor matrix element in Ref. [5] is not mD + mP but mD. This in turn introduces different denominators in Eqs. (7), (17) and (18) in the terms involving fT with respect to Eqs. (D1)-(D3) of Ref. [5]. 4 By comparing the diagonal covariance matrix elements of Tab. 7 in Ref. [11] with σ 2 i from Tab. 6, we notice that the ordering of parameters in rows/columns of Tab. 7 must be {f (0), c0, c+, PS, PV } rather than {f (0), c+, PV , c0, PS}. We thank Vittorio Lubicz for confirmation. Here, P V /S corresponds to P +/0 in Eq. (A3). and 90%, respectively, such that naïve variation of parameters independently produces larger(huge) uncertainties for f +,0 (f T ) at high q 2 . Taking into account correlations, we present in Fig. 10 the D → π form factors within 1σ uncertainties. Also shown (gray band) is f + obtained from the HFLAV collaboration [40], parametrized as where the parameters r k = a k /a 0 are extracted from data on semileptonic D → π ν decays from several experiments. Both results for f + are in good agreement with each other up to q 2 2.3 GeV 2 .
Towards the low recoil end point, f HFLAV + /f LQCD + 1.4 at q 2 max = (m D −m π ) 2 . Symmetry breaking effects in the form factor relation f T = f + (1 + O(α s , Λ QCD /m c )) amount to about O(30%) in the LQCD results. We assume the D s → K and D → π form factors to be the same, supported by the lattice study [41].
In Fig. 11 we show the factor γ = 2mc (12) for the D → π (D s → K) transition in green (red). As noted previously, γ ≈ 1 is a good approximation in the full q 2 -range. where and employx |g R | 2 and 4/X 2 1.
Experimental constraints on CP-violation in D 0 −D 0 mixing, x 12 sin φ 12 2×10 −4 , are stronger than (B1) by about ∼ 0.04 [40,45]. Both SUSY and leptoquark models contribute to ∆C = 2 Wilson coefficients with the square of the couplings relevant for ∆C = 1 coefficients. For instance, C ( ) 7 in SUSY is proportional to one power of flavor-violating scalar quark mass terms whereas mixing is induced at second order. Constraints on the ∆C = 1 couplings from mixing allowing for order one phases are hence about ∼ 0.2 stronger than without CP-phases.
where we allow for the possibility of having three RH neutrinos ν. The charges F ψ i are subject to constraints from gauge anomaly cancellation conditions (ACCs). An excellent introduction to the subject is given in [46], for recent phenomenological applications, see, for instance, [47][48][49][50]. The ACCs read [48]: gauge-gravity: Since they are SM singlets, the RH neutrinos only appear in Eqs. (C5),(C7). In the SM+3ν R (the SM), the 18 (15) charges are constrained by 6 ACCs. Important features of the ACCs and their solutions are (see [48] and references therein for details): • Rational solutions: We assume that all U (1) -charges are rational numbers F ψ ∈ Q .
• Rescaling invariance: Any solution of the ACCs can be rescaled by any rational number k ∈ Q , F ψ → kF ψ , ∀ψ ∈ {Q i , u i , d i , L i , e i , ν i }, which constitutes another solution. As rescaling all charges is equivalent to a rescaling of the U (1) gauge coupling, these solutions are in the same equivalence class and, hence, are not independent from each other. We therefore assume integer solutions F ψ ∈ Z without loss of generality.
• Permutation invariance of fermions: The ACCs are invariant under the permutation of generation indices within each specific species ψ.
To extract concrete solutions of the non-linear Eqs. (C2)-(C7), we use computational algebraic geometry and perform a Gröbner basis computation [49] with the aid of Mathematica and obtain analytical expressions for the charges F ψ . We impose phenomenological constraints, for instance,   Sizable renormalization group coefficients can be induced in models with large U (1) -charges.
The study of these effects is beyond the scope of this work.