Theoretical and experimental status of rare charm decays

Rare charm decays offer the unique possibility to explore flavour-changing neutral-currents in the up-sector within the Standard Model and beyond. Due to the lack of effective methods to reliably describe its low energy dynamics, rare charm decays have been considered as less promising for long. However, this lack does not exclude the possibility to perform promising searches for New Physics per se, but a different philosophy of work is required. Exact or approximate symmetries of the Standard Model allow to construct clean null-test observables, yielding an excellent road to the discovery of New Physics, complementing the existing studies in the down-sector. In this review, we summarise the theoretical and experimental status of rare charm $|\Delta c|=|\Delta u|=1$ transitions, as well as opportunities for current and future experiments such as LHCb, Belle II, BES III and the FCC-ee. We also use the most recent experimental results to report updated limits on lepton-flavour conserving and lepton-flavour violating Wilson coefficients.


Introduction
Despite being consistent with an enormous amount of experimental results, there are undoubtedly phenomena that the Standard Model (SM) fails to explain and a more fundamental theory has to exist. Nowadays, the combined effort of theoreticians and experimentalists is to formulate extensions to account for the apparent shortcomings of the SM, collectively referred to as New Physics (NP) models, and to find hints for its breakdown.
Rare decays of flavoured mesons containing an s, c or b quark receive contributions from flavour-changing neutral-current (FCNC) processes and are sensitive probes to heavy degrees of freedom at mass scales much higher than the available center-of-mass energies in the most powerful particle colliders. New and yet unknown particles and interactions can modify the rate of such a process, change the angular distributions of the decay products, or introduce additional sources of CP violation. Precision measurements of FCNC processes have been proven to have the potential to indirectly point towards the existence of new particles long before their direct detection. Famous examples are the first hints for the necessity of the charm quark, 1 or the extremely heavy mass of the top quark, inferred from the suppression of neutral kaon decays, 2, 3 and the observation of neutral B 0 meson mixing, 4 respectively.
While in the past most of the experimental effort has focused on studying rare processes in the kaon and beauty sectors, investigations of rare charm decays, which are sensitive to |∆c| = |∆u| = 1 transitions, have only started. Due to the strong suppression of rare charm decays, a major part of experimental analyses has been restricted to setting upper limits on rare and forbidden decay processes. Since most decays are dominated by resonant contributions, which cannot be described in a consistent theoretical framework, rare charm decays have been considered as less promising for a long time. The presence of large uncertainties coming from effects of the strong interaction at low energies, such as hadronisation and the formation of light intermediate resonances, does not prevent clean searches for NP in rare charm decays. However, a different philosophy of work is required compared to beauty and kaon physics, aiming for optimised observables.
In this review, we highlight the unique phenomenology of rare charm decays and how this can be used to search for NP in the up-sector. The SM symmetries in the charm system offer the possibility to define null-test observables with small or negligible theory uncertainties in resonance-dominated semi-leptonic and radiative decays. The possibility to investigate angular distributions, CP asymmetries and tests for lepton universality experimentally in these decays with typical branching fractions of 10 −5 − 10 −7 has only opened recently, and precision measurements are expected to be possible in the near future at the current flavour experiments LHCb, Belle II and BES III.
In addition, in light of the persistent anomalies in rare B decays (e.g. see Ref. 5 and references therein), the charm systems offers a complementary opportunity whose potential has hardly been exploited so far.
The article is organised as follows: We start with an overview of the theoretical framework in Sect. 2, discussing the short-and long-distance description of rare charm decays. In Sect. 3, we briefly summarise the NP models that have been explored in the context of rare charm decays. We then present a summary of experimental searches for rare and forbidden decays in Sect. 4, and how experimental limits can be translated in model-independent bounds on Wilson coefficients in Sect. 5. In Sect. 6, we discuss strategies to test the SM with clean null tests in semi-leptonic and radiative decays. We conclude the review with a brief outlook to future prospects in Sect. 7 and closing remarks in Sect. 8.

Theoretical framework
In the SM, the leading contribution to |∆c| = |∆u| = 1 transitions appears at 1loop level, consequently its study provides an excellent window to test its quantum structure. Fig. 1 displays a possible 1-loop contribution via a W boson, with internal down-type quarks (d,s,b) flowing in the loop. The amplitude of this diagram A(c → u) can be written as where λ i ≡ V * ci V ui encodes the dependence on the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements and thus the CP violating phenomena in the SM. The loop  Using the unitarity condition of the CKM matrix, i=d,s,b λ i = 0, one can remove the dependence on λ d , and Eq. (1) becomes where ξ b ≡ λ b /λ s . Eq. (2) exhibits two important features that characterise rare charm decays in the SM: • Branching fractions for rare c → u transitions are suppressed via the Glashow-Iliopoulos-Maiani (GIM) mechanism, as shown in the first term of Eq. (2), resulting in the naïve estimation of A(c → u) ∼ O(10 −8 ).
Furthermore, the very effective GIM mechanism leads to specific angular distributions of the final states, as will become clearer later. In summary, short-distance contributions in the SM for |∆c| = |∆u| = 1 transitions are well below the current experimental precision due to both GIM-and CKM-suppressions. Hence, the study of observables like branching fractions, CP and angular asymmetries provides a fantastic road towards the discovery of physics beyond the SM, since large signals indicate clear signs of NP. The above statements are based on a naïve estimation of the SM contribution. For the benefit of more precise calculations, a robust theoretical framework is needed. In the following sections, we give a brief summary of the current tools used to describe both the short-distance and long-distance effects involved in rare charm decays.

Short-distance description
In the most general form, rare |∆c| = |∆u| = 1 transitions are described by the following effective Hamiltonian, [6][7][8] with the following local dimension-six operators, [9][10][11] where q L,R = 1 2 (1∓γ 5 ) q are chiral quark fields, T a are the generators of SU (3) C and g s is the strong coupling. Furthermore, G F is the Fermi constant and α e = e 2 /(4π) stands for the fine-structure constant with the electromagnetic coupling e. Finally, σ µν = i 2 [γ µ , γ ν ] and F µν , G µν a with a = 1, ..., 8 denote the electromagnetic and gluonic field strength tensor, respectively. Primed operators are obtained replacing L(R) → R(L) in Eq. (4). Note that Eq. (3) is a direct consequence of the operator product expansion (OPE), which allows the factorisation between the matrix elements of local operators O i and the Wilson coefficients C i , which parameterise the strength of each operator. The operators given by Eq. (4) are constructed with the light fields with masses below µ < m b , with µ being the renormalisation scale. The heavy fields with masses greater than µ m b are removed as dynamical degrees of freedom. Effects of the heavy fields are implicitly encoded in the Wilson coefficients. Hence, experimental deviations of these coefficients from their SM predictions indicate a signal of physics beyond the SM (BSM), since such a discrepancy requires additional dynamical degrees of freedom at high-energy scales.
The general OPE setup (i.e. definition of the operator basis, matching of the SM contributions onto the effective theory, renormalisation group (RG) evolution of Wilson coefficients from the high-energy to the low-energy scale) needed for the computation of the Wilson coefficients relevant in |∆c| = |∆u| = 1 transitions is almost analogous to |∆b| = |∆s| = 1 transitions. [9][10][11][12][13][14][15][16] However, due to the specific CKM and mass structure of charm FCNCs, the predictions for this sector can differ from each other by several orders of magnitude depending on which corrections are taken into account. Ref. 8 provides a consistent expansion of the complete SM computation for rare charm transitions to O(α s ) with α s = g 2 s /(4π). In the following, we give an overview of the steps required to achieve a perturbative nextto-(next-to-) leading order precision of the short-distance contributions. 8 These are: (i) matching of SM contributions onto Weak Effective Theory at µ = M W , (ii) RG-evolution of Wilson coefficients from M W to m b , (iii) integrating out the b quark and second matching at µ = m b , (iv) RG-evolution of Wilson coefficients from m b to the charm scale µ c .
In the SM, |∆c| = |∆u| = 1 transitions are driven via the exchange of a W boson between two weak charged quark currents as shown in Fig. 2 (a). At high energies, where quantum chromodynamics (QCD) corrections are small, however, still below the scale of electroweak symmetry breaking, the W boson can be integrated out, and the interaction can be well described in terms of a single operator, as shown Fig. 2 Fig. 2. Exemplary contributions to |∆c| = |∆u| = 1 transitions at different energy scales. Current-current topology shown in (a). Current-current topologies when the W boson is integrated out, with and without gluonic (denoted by g) corrections (c) and (b), respectively. Photon γ penguin contribution generated at the b quark mass threshold shown in (d).
With decreasing energies, µ < M W , the gluonic corrections become sizeable and have to be summed up using the OPE and the renormalisation group equations (RGEs). [15][16][17] This procedure requires an additional operator, see Fig. 2 The operators (5) and (6) are called current-current operators and are also responsible for non-leptonic D-meson decays. At µ > m b , the treatment of light quarks as massless leads to a fully effective GIM mechanism causing the cancellation of penguin contributions, 6 as already seen in Eq. (2). Penguin contributions appear when the b quark is integrated out as an effective degree of freedom at its threshold, see Fig. 2 (d). After the matching, the Wilson coefficients are evolved from m b to m c using the RGEs, which induces mixing between different |∆c| = |∆u| = 1 operators, given by Eq. (4). 18 The last step is the evaluation of the Wilson coefficients down to m c . The computation of the matrix elements associated to certain operators result in "effective" coefficients which account for effects of non-zero light quark masses. Strictly speaking, the effective coefficients are not Wilson coefficients, since they contain low energy effects. Wilson coefficients are process universal, independent of the operator basis, regularisation and renormalisation schemes. 19 The effective coefficients C eff 7,9 , which are the only non-vanishing ones, at the charm scale µ c read 18, 20 where q 2 is the dilepton invariant mass squared and L(q 2 , m q , µ c ) is a function that accounts the low dynamical effects due to m q = 0. The explicit form of L(q 2 , m q , µ c ) can be found in Ref. 18. The coefficients C eff 7, 9 are dominated by the matrix elements Especially C SM 10 = 0 distincts charm FCNCs from K or B physics. Since C 10 corresponds to an axial vector coupling for the leptonic part, effects on the V − A structure of the SM are shut off at the charm scale and charm physics within the SM is dominated by QCD and quantum electrodynamics (QED) effects.

Long-distance description
The OPE allows for the separation of short-distance and long-distance effects. In order to fully assess the decay amplitude of rare c → u + − transitions, besides the Wilson coefficients C i (µ), we also need to determine the hadronic matrix elements O i (µ) which encode the non-perturbative dynamics at low energies. Since measurable observables cannot depend on the renormalisation scale, the dependence on µ has to cancel in the product of C i (µ) and O i (µ) . A proper understanding of the long-distance dynamics with a solid effective field theory framework is still missing in the literature, mainly due to Λ QCD ∼ m c , which makes a perturbative expansion in powers of 1/m c slowly converging at best. In the following, we summarise the available techniques to deal with this challenging task.

Form factors
Form factors (FFs) parameterise our lack of knowledge about the hadronic effects in a hadronic transition. The hadronisation of the |∆c| = |∆u| = 1 operators given by Eq. (4) leads to the factorisation between the lepton and quark currents, where h c represents a charmed hadron b and F is the final state. H i α1,...,αn and L α1,...,αn i are the quark and lepton currents of O i , respectively, α 1 , ..., α n being shared Lorentz indices between both currents. The quantity 0|L α1,...,αn i | + − can be computed applying perturbation theory in QED, while effects contained in h c |H i α1,...,αn |F require non-pertubative techniques to deal with the dynamics of QCD at low energies. Even though there is no rigorous effective field theory to face these effects, we can reduce the problem using symmetries at low energies. In particular, requiring Lorentz structure and parity invariance of QCD, the vectorial V (n = 1) current in which a D meson decays into a pseudoscalar P can be written as where the functions F V ± are the FFs associated to the V current, and the fourmomenta q, p are given by q α = (p D − p P ) α and p α = (p D + p P ) α , respectively. Due to four-momentum conservation, FFs can only depend on the four-momentum transfer q. In contrast, for a tensorial T (n = 2) current, due to its antisymmetric Lorentz structure, only one FF appears Currently, FFs can only be computed using methods such as lattice gauge theory, QCD light-cone sum rules (LCSR) c or fitting theory/models d to experimental data. Mode 0.20 ± 0.04 0.111 ± 0.008 -- where a M is the resonance parameter with M = ρ 0 , φ, η, η , f and isospin has been employed in Eq. (12) to relate the ρ 0 to the ω. Here, m M and Γ M denote the mass and the total decay rate of M , respectively. The a M parameters can be extracted from measurements of branching fractions B(h c → F M ) and B(M → + − ), and are given in Table 2 for some examples studied recently. The strong phases δ ρ 0 , φ, η are the largest source of uncertainty. They can only be constrained if experimental input on branching fractions in resonance-dominated regions is provided, which is often not the case. In most numerical analyses, they are varied between −π and π. Resonance contributions in radiative decays D → V γ, with V a vector meson, can be included through phenomenological models. For example, the model considered in Refs. 43, 44 is a mixture of factorisation, heavy quark effective theory and chiral theory. The parameters of this model are unknown and have to be extracted from experimental data, allowing implicitly for a breaking of the SU (3) F flavour symmetry.
Since sensitivities of current experimental searches (grey shaded areas) are close to the orange/red resonant curves, searching for NP in branching fractions is a challenging endeavour, as interference effects between NP and long-distance contributions have to be taken into account and increase the theoretical uncertainties in the interpretation of NP contributions. Despite the challenges of discovering NP in branching fractions, upper limits can still serve to provide constraints on NP, as we will see in Sect. 5. Furthermore, branching fraction measurements of resonancedominated regions give insight into QCD at low energies and help to constrain strong phases, for instance.

New Physics models generating rare charm decays
New particles and interactions from various NP models are suitable to generate |∆c| = |∆u| = 1 transitions at tree-or loop level. Before discussing the experimental status of rare charm decays, in this section we briefly summarise the literature on model-dependent BSM studies. Examples are models with non-minimal Higgs sector, such as two Higgs doublet models, 47 28, 41, 49, 50, 61, 75, 76. In this review, we mainly focus on a model-independent description of rare charm decays.

Experimental searches for rare and forbidden decay modes
Searches for rare and forbidden decays have been carried out investigating purely leptonic final states, semi-leptonic multi-body decays of charged (D + , D + s ) and and neutral (D 0 ) mesons, as well as decays of charmed baryons. Searches for semileptonic decays are often restricted to regions in q 2 away from the resonances to enhance sensitivity to NP as shown in Fig. 3. Forbidden decays refer to those which violate conservation of charged lepton flavour and lepton number. As resonance pollution to these decays is absent, no q 2 binning is needed for forbidden decay modes, and their studies represent a clear null test, complementary to those presented in Sect. 6. Searches for decays of D 0 and D + mesons violating conservation of baryon number also exist, 78, 79 but will not be discussed further in this review.
The most recent experimental results have been obtained by the LHCb, BaBar, Belle and BES III collaborations. The LHCb detector 80, 81 is a single arm forward spectrometer designed to study decays of mesons containing a c or b quark, sited at the LHC (CERN, Switzerland). LHCb has been designed to study proton-proton (pp) collisions in its main operation mode. BaBar 82 and Belle 83 are cylindrical largesolid-angle detectors that operated at the PEP-II (SLAC National Accelerator Laboratory, USA) and KEKB (KEK, Japan) asymmetric-energy e + e − colliders, known as b-factories. BES III 84 is a general-purpose detector recording e + e − collisions in the double-ring collider BEPCII (IHEP, China). Thanks to the large production cross section at hadron colliders in the forward region, 85,86 LHCb can profit from the world's largest recorded data set of charm hadron decays to date. Conversely, the other experiments benefit from detectors with excellent capabilities to reconstruct final states including neutral particles and electrons.
Older results published by CLEO II and the Fermilab E653 and E791 collaborations still hold the most stringent limits for some decay channels, see Refs. 87, 88, 89, 90 for details. These measurements are not discussed in detail in this review, however, their results will be added to the summary tables for completeness. All experimental limits in this section are quoted at a 90 % confidence level (CL).

Searches for purely leptonic decays
Hadronic uncertainties on theoretical predictions are minimal in purely leptonic decays. However, their decay rates are subject to an additional helicity suppression, making them extremely rare. In SM, these decays are dominated by long-distance contributions from D 0 → γ * γ * → + − 58, 91, 92 with an estimated branching fraction of order 10 −11 for muons, see Sect. 5. Experimentally, the detection of final states consisting of two leptons is rather simple, generally allowing to set the more stringent limits on the branching fractions compared to decays to final states involving hadrons. The possibility of NP searches in rare charm decays triggered the attention of experimental particle physicists already in 1988 and first searches for the rare decays D 0 → e + e − , D 0 → µ + µ − and the lepton flavour violating (LFV) decay D 0 → µ ± e ∓ started by the CLEO and ARGUS collaborations. 93,94 The world's most stringent limit nowadays on D 0 → e + e − decays has been set by Belle in 2010, analysing a data set corresponding to an integrated luminosity of 660 fb −1 collected at a center-of-mass energy at or close to the Υ (4S) resonance. 95 The best limits on the decays D 0 → µ + µ − and D 0 → µ ± e ∓ have been published by LHCb in the years 2013 96 and 2016, 97 respectively, using data sets corresponding to integrated luminosities of 0.9 fb −1 and 3.0 fb −1 . The best upper limits on the branching fractions of purely leptonic rare charm decays are [95][96][97]

Searches for semi-leptonic decays
Searches for semi-leptonic decays of neutral D 0 and charged D + , D + s charm mesons into two leptons and additional hadrons cover a large variety of final states, and we will briefly discuss the most recent publications. Singly Cabibbo-suppressed decays with two oppositely charged leptons + − ( = e, µ) in the final state of the form D → F + − are sensitive to FCNC processes. Here, F can be one or several neutral (π 0 , K 0 S , η, ρ 0 , ω, K * 0 , φ) and/or charged (π + , K + , ρ + , K * + ) mesons. Final states comprising two leptons of different flavour correspond to LFV modes, while modes with two leptons carrying the same electrical charge are lepton number violating (LNV). To date, no indications of non-resonant short-distance contributions to rare decay modes or hints for forbidden modes exist.
In autumn 2020, the LHCb collaboration has published a search for 25 rare and forbidden decays of D + and D + s mesons into two leptons and a charged kaon or pion, analysing a data set of pp collisions corresponding to an integrated luminosity of 1.6 fb −1 . 45 The achieved limits in the range (1.4 − 640) × 10 −8 improve the previous ones in most cases by at least one order of magnitude. Resonant contributions are minimised by vetoing the region [525, 1250] MeV/c 2 in dilepton mass and extrapolating the signal yields to the vetoed regions assuming a uniform distribu- Table 3. Most stringent upper limits (UL) at a 90 % CL on the branching fractions of (left) rare and (right) forbidden decays of D + and D + s mesons.
tion of the particles across the phase space. The resonant decay modes proceeding via an intermediate φ meson D + (s) → π + [ + − ] φ are used as calibration and normalisation. The most stringent results are obtained analysing final states involving muons. Decay channels of D + mesons including negatively charged kaons such as D + → K − e + e + , D + → K − µ + µ + , D + → K − µ + e + have not been investigated due to a large amount of background coming from misidentified D + → K − π + π + decays. A listing of the observed limits can be found in Table 3. The table summarises the most stringent limits on rare and forbidden decay channels of D + and D + s mesons that have been investigated to date. The BES III collaboration has published a search for numerous decay channels of D + and D 0 mesons into final states comprising two electrons in 2018. 99 These are decays of charged D + mesons to two electrons accompanied by a pair of a neutral and a charged pseudoscalar meson (e.g. D + → π 0 π + e + e − ), decays of neutral D 0 mesons into two electrons plus a neutral meson (e.g. D 0 → π 0 e + e − ) or two oppositely charged pseudoscalars (e.g. D 0 → π + π − e + e − ). The analysis uses e + e − collision data corresponding to an integrated luminosity of 2.93 fb −1 . The data has been recorded at a center-of-mass energy of √ s = 3.773 GeV, which is close to the D + D − or D 0 D 0 mass threshold. Upper limits in the range (3 − 41) × 10 −6 are determined using a double tagging approach 101, 102 designed to measure absolute branching fractions. Due to the limited luminosity and lower charm production cross section in e + e − annihilations, the achieved limits are less stringent compared to those set by LHCb on muonic modes. 46,103 The measured limits for all channels under study are summarised in Table 3 for decays of charged D + mesons, and limits on decays of neutral D 0 mesons are listed in Table 4. In addition, a search for heavy Majorana neutrino LNV decays with two electrons has been published by BES III in 2019 , 100 where best limits on D + → K 0 S π − e + e + and D + → K − π 0 e + e + have Table 4. Most stringent upper limits (UL) at a 90 % CL on the branching fractions of (left) rare and (right) forbidden decays of D 0 mesons.
been reported, which can also be found in Table 3.
The most stringent upper limits on forbidden decays of neutral D 0 mesons have been reported by the BaBar collaboration in two successive publications 104, 105 during 2020. A data set of e + e − annihilation corresponding to integrated luminosity 468 fb −1 recorded at or close to the Υ (4S) resonance has been analysed. The upper limits are set relative to decays of purely hadronic decays. Upper limits in the order of (1.0 − 30.6) × 10 −7 and (5.0 − 42.8) × 10 −7 are found for decay modes involving a pair of oppositely charged and neutral hadrons, respectively. The achieved results can be found in Table 4. With the exception of D 0 → e ± µ ∓ , all best limits on forbidden D 0 meson decays are set by BaBar to date.
Note that more than half of the limits summarised in Table 3 and Table 4 originate from measurements published during 2018-2020, proving the great experimental progress that the field has made in recent times.

Observation of resonance-dominated semi-leptonic decays
Instead of vetoing resonance-dominated regions of the decay phase space, a complementary approach has recently been applied for D 0 → P 1 P 2 + − decays, where P 1 P 2 stands for a pair of pseudoscalar mesons. The branching fraction measurements are done binned in regions of q 2 . While limited sensitivity to the short-distance component is still given in bins which are away from the resonances, the regions around the resonances are fully dominated by long-distance contributions. However, signal decays in those bins can be used to perform in-depth studies of SM null tests as will be explained later in Sect. 6.
The observation of the decay mode D 0 → K − π + [µ + µ − ] ρ 0 /ω has been reported by LHCb in 2016. 106 The measurement is limited to a region of dimuon mass of [675 − 875] MeV/c 2 around the ρ 0 /ω mass, where a significant signal is observed for the first time. The branching fraction is measured relative to D 0 → K − π + π − π + decays to be 106 where the uncertainties are statistical and systematic, respectively. Since the decay is a Cabibbo-favoured mode, no |∆c| = |∆u| = 1 FCNC contributions are involved. However, the decay provides an important reference channel for further measurements of D 0 meson decays to four-body final states and to test QCD methods.
Integrating over the full dimuon-mass region, the total branching fractions of D 0 → π + π − µ + µ − and D 0 → K + K − µ + µ − decays has been found to be consistent with SM expectations 61 and measured to be: 46 where the uncertainties are statistical, systematic and due to the limited knowledge of the normalisation branching fraction, respectively. Even though dominated by resonant contributions, these represent the rarest observed decays of charmed mesons to date.
As of today, a single decay mode with two electrons in the final state has been observed. The BaBar collaboration succeeded in making the first observation of   108 coinciding with the dimuon-mass region defined by LHCb when measuring the muonic mode. 46 The analysis uses a data set of e + e − annihilation data corresponding to 468 fb −1 . Fully hadronic decays of D 0 → K − π + π − π + decays are used as normalisation. The branching fraction has been determined to be 108 where the uncertainties are statistical, systematic and due to the limited knowledge of the normalisation branching fraction, respectively. Assuming a SM rate, the available statistics accumulated by BaBar, Belle and BES III is expected to be insufficient to make further observations of rare four-body dielectron modes. LHCb has not yet started to measure branching fractions of resonance-dominated decays involving two electrons in the final state, however, contributions are expected in the near future. 109

Studies of semi-leptonic baryonic decays
The most recent searches for rare baryonic decays are restricted to decays of Λ c baryons decaying into a pair of leptons accompanied by a (anti)proton. While the world's best limits on the dielectron decay Λ c → p e + e − , as well as LFV and LNV decays have been set in the order (2.7 − 19.0) × 10 −6 by BaBar already in 2011, 98 LHCb has investigated decays into two muons and a proton in the year 2018 using pp collision data corresponding to 3 fb −1 . 110 The new limit improves the previous result on the branching fraction in the non-resonant region of dimuon mass by more than two orders of magnitude. Ranges ± 40 MeV/c 2 around the known ω and φ meson masses have been vetoed. The result reveals the capability where the uncertainties are statistical, systematic and due to the limited knowledge of the normalisation branching fraction, respectively. The resonant decay mode proceeding through an intermediate φ meson Λ c → p [µ + µ − ] φ has been used as normalisation, using the known branching fractions of the decays Λ c → p φ and φ → µ + µ − . A single measurement investigating decays of Σ c baryons exists by the E653 collaboration, 88 dating back to the year 1995. We give a summary of the most stringent limits on rare charmed Λ c baryon decay modes that have been studied to date in Table 6.

Model-independent constraints on Wilson coefficients
Using the most recent available experimental results presented in the previous section, we report the most stringent constraints on Wilson coefficients in the charm system, separately assuming lepton-flavour conservation and allowing for LFV.

Lepton flavour conserving bounds
Currently, the best constraints on dimuon and dielectron Wilson coefficients are obtained from upper limits on B(D + → π + µ + µ − ) and B(D + → π + e + e − ). Using the experimental upper limits on B(D + → π + µ + µ − ) < 6.7 · 10 −8 at 90 % CL presented in Sect. 4, and neglecting SM contributions, we obtain the following bounds for muons, Notice that coefficients in Eq. (19) are slightly improved with respect to Ref. 41. Right-handed currents can be included replacing C → C + C . The decay D 0 → µ + µ − provides additional information on C (µ) ( ) S,P,10 via the experimental upper bound presented in Eq. (13). The short-distance contributions can be written as where f D is the D 0 meson decay constant, see Sect 2.2.1. Albeit the active helicity suppression in Eq. (20), it gives the best constraints on NP contributions to C where again the SM contributions have been neglected since the relevant Wilson coefficients are zero in the SM, as discussed in Sect. 2.1. The most sizeable contribution in the SM is estimated to originate from long-distance contributions in D 0 → γ * γ * → µ + µ − . 58, 91, 92 Using the current limit B(D 0 → γγ) < 8.5 · 10 −7 at 90 % CL, 111 we obtain 58 well below the current experimental limit (see Eq. (13)). Eq. (22) can be considered as an upper limit and any measurement significantly exceeding this bound would signal NP. Constraints on dielectron modes are weaker than dimuon ones, using the upper limits on B(D + → π + e + e − ) < 1.1 · 10 −6 and B(D 0 → e + e − ) < 7.9 · 10 −8 at 90 % CL presented in Sect. 4, one obtains 18 which are approximately a factor three weaker than the muonic ones, see Eqs. (19) and (21). Similar bounds on muon and electron Wilson coefficients are obtained with the study of the transverse momentum (p T ) spectrum of dilepton pairs produced in pp collisions, recorded by the ATLAS 112 and CMS 113 detectors. In particular, the high-p T range of the spectrum can be used to determine bounds on Wilson coefficients. Due to the limited phase space of charm decays, bounds on τ Wilson coefficients are only accessible using high-p T data 114 and the current bounds read, In summary, Eqs. (19), (21), (23) and (24) represent the strongest constraints on Wilson coefficients for lepton flavour conserving transitions to date.

Lepton flavour violating bounds
The phenomenology of these decays can be adopted from the lepton flavour conserving case, by introducing lepton flavour indices for operators and different Wilson coefficients. For clarity, we additionally indicate Wilson coefficients corresponding to LFV operators with K. We update the bounds on these LFV Wilson coefficients 41 where S,P can be better constrained from data on leptonic decays, using the measurement of B(D 0 → e ± µ ∓ ) presented in Eq. (13). One obtains K S − K S ± 0.04 K 9 − K 9 2 + K P − K P + 0.04 K 10 − K 10 2 0.01 . (26) As stated in Ref. 41, K eτ couplings can be constrained by measurements of B(D 0 → e ± τ ∓ ). However, no experimental measurement exists to date and only high-p T constraints are available for K µτ and K eτ , see for instance Ref. 115.

Charming opportunities to probe the Standard Model with null tests
Due to its unique properties, the charm system allows the definition of clean observables which are null in the SM. In the following section, we outline ways to search for NP by testing for lepton universality, investigating angular distributions, CP asymmetries and radiative decays.

Testing lepton universality
In the last years, lepton universality (LU) ratios have gained interest due to anomalies found in the beauty sector. 5 The corresponding ratios in the charm sector are defined as 18,28,41,49 Table 7. R D + π + in the SM and in NP scenarios for different q 2 regions. Ranges correspond to uncertainties from form factors and resonance parameters. 41 SM |C In the SM, the electroweak gauge bosons couple to leptons of different generations with equal strength, leading to a well-controlled SM prediction of 28, 41 Phase-space and electromagnetic corrections to R hc F are highly suppressed and cannot exceed the percent level. 116-118 NP extensions do not necessarily respect LU, so R hc F can significantly deviate from unity. These ratios greatly profit from the cancellation of the dominant theoretical and experimental uncertainties if the same kinematic cuts (q 2 max , q 2 min ) are applied to the dielectron and dimuon modes. 119 As an example, Table 7 indicates experimentally possible NP values of R D + π + and R D + s K + for the full, low and high q 2 -integrated intervals. 41 Notice that the largest BSM signals are expected in the low and high q 2 regions, where the SM contributions are smaller. However, due to the limited knowledge of the resonance contributions in the low q 2 region, an interpretation in terms of specific NP models is harder than in the high q 2 region, where resonance influences are minimal. BSM predictions of LU ratios in D → P 1 P 2 + − modes have been worked out in Ref. 28. Effects up to 15% are found on R D 0 π + π − in the q 2 -integrated range for D 0 → π + π − µ + µ − decays. Due to the lack of branching fraction measurements of dielectron decay modes, LU tests in c → u + − processes are largely unexplored and only limits exist. Naïve LU ratios in Ref. 28 find bounds on R D 0 π + π − 0.1 and R D 0 K + K − 0.01 using the available upper limits on the dielectron modes 99 and branching fraction measurements of the dimuon modes 46 discussed in Sect. 4.2 and 4.3, respectively. As resonant-dominated dielectron modes are expected within the reach of the LHCb future sensitivity, 109 we expect first model-independent LU tests in rare charm decays the near future.

Angular observables
In the SM, the absence of axial vector lepton currents in rare charm decays (corresponding to C SM 10 = 0, cf. Eq. (8)) leads to very specific angular distributions of the final state particles, and allows for clean null tests irrespective of form factors and their uncertainties.
The theory predictions of form factors have improved in recent years and a variety of |∆c| = |∆u| = 1 decay modes can be studied, see Table 1. For angular distributions the number of observables increases with number of particles in the final state, as well as with the spin of particles included. In that manner, studying angular observables of different decay modes helps to pin down combinations of NP Wilson coefficients in a complementary way. In the following, the angular distribution for key channels are studied. Further model-dependent studies are available in e.g. Refs. 54, 55, 56, 66, 75.
The differential distribution of D → P + − transitions reads 120 where θ denotes the angle between the − -momentum and the P -momentum in the dilepton rest frame. Its particular angular distribution provides two clean null tests of the SM: the lepton forward-backward asymmetry A FB , and the "flat" term F H , where Γ = Γ(q 2 min , q 2 max ) = with integration limits depending on the q 2 -bin. These observables are sensitive to operators with Lorentz structures not present in the SM. Neglecting lepton masses, it follows that Taking the bounds on Wilson coefficients from Eqs. (19), (21) and (23), one obtains 18 A FB is zero in the SM up to higher corrections, which are suppressed by powers of m D /M W and come from higher dimensional operators or by α e /(4π) (D → π γγ → π + − ). 21, 120 F H is also highly suppressed in the SM. For instance, , and both are even smaller at high q 2 . 41 In addition, F H (D + (s) → π + (K + ) e + e − ) is even further suppressed in the SM. Then, any non-zero measurement of these observables hints to NP generated by (pseudo)-scalar and (pseudo)-tensor operators, see Eq. (33).
The observables A FB and F H introduced for D → P + − in Sect. 6.2.1 can also be defined for baryonic decays Λ c → p + − . Ref. 39 performs an analysis of A FB and F H in Λ c → p µ + µ − , considering only contributions from C (µ) ( ) 7,9,10 . In contrast to F H (D → P + − ), the SM prediction of F H (Λ c → p µ + µ − ) is polluted by resonance contributions (mainly from C R 9 ), leading to values of F H (Λ c → p µ + µ − ) SM ∼ 0.6, such that it cannot be considered as clean null test.
On the other hand, A FB (Λ c → p µ + µ − ) SM vanishes because it is proportional to the product C (µ) 9 C (µ) 10 . Consequently, a small NP contribution to C (µ) 10 could be enhanced through the interference with the resonant contributions from C R 9 , producing a non-zero value of A FB (Λ c → p µ + µ − ), causing a clear sign of NP.
Here, θ is the angle between the + -momentum and the D 0 -momentum in the dilepton center-of-mass system (cms), while θ P denotes the angle between the P +momentum and the negative direction of flight of the D 0 -meson in the (P 1 P 2 )-cms, and φ is the angle between the normal of the (P 1 P 2 )-plane and the ( + − )-plane in the D 0 rest frame. The angular coefficients I i ≡ I i (q 2 , p 2 , cos θ P1 ) in terms of transversity amplitudes can be found in Ref. 28. In particular, since C SM( ) 10 = 0, I 5,6,7 are zero in the SM. Hence, they constitute a formidable set of null tests observables. 28,61 From the CP-odd angular coefficients I 5,6,8,9 , asymmetries can be defined as: with the CP-averaged decay rate Γ ave = Γ+Γ 2 and the angular coefficients I i (Ī i ) for D 0 (D 0 ) mesons. In the SM, the observables A 5,6,8,9 are negligible given the experimental sensitivities today and constitute (approximate) SM null tests. 28,61 Some coefficients I i can be obtained by symmetric or asymmetric integrations in the decay angles, as for example the coefficient I 6 , which is proportional to cos θ  and A FB . This allows for a simplified experimental determination with respect to a full angular analysis. One can define the p 2 and cos θ P integrated observables as and focus on the effects from axial-vector NP contributions. 28 Figure 4, adopted from Ref. 28, shows the integrated I 5,7 observables as a function of q 2 for different NP benchmarks. I 5 and I 6 have similar BSM-sensitivity. Notice that NP effects are enhanced around the ρ/ω and the φ resonances, referred to as resonance catalyzed . 122 This has the advantage to allow profiting from the full sample statistics, while NP searches in branching fraction measurements are often limited to restricted q 2 ranges which suffer from low statistics.

CP asymmetries
Since CP -violating effects from SM contributions are CKM-suppressed in rare charm decays, as seen in Sect. 2, studies of CP asymmetries are powerful possibilities to obtain information on NP Wilson coefficients. In general, the CP asymmetry in the dilepton mass distribution is defined as where Γ corresponds to the decay rate of the CP -conjugated process. As in Eq. (32), Γ andΓ denote the q 2 integrated decay rates. Similar to the angular observables in D 0 → P 1 P 2 + − shown in Figure 4, CP asymmetries generated by NP contributions show maximal signatures in the vicinity of the resonances. The strong phases associated with the resonance contributions are needed to induce non-vanishing  asymmetries in the interference with NP amplitudes. Fig. 5 (left) shows the CP asymmetry for the decay mode D + s → K + µ + µ − in the q 2 region around the φ resonance with exemplary values for the strong phase and the NP contributions. The CP asymmetry can reach the percent level. CP asymmetries of this kind can be defined for any of the decay modes listed in Table 1, and apparent from Fig. 5, q 2 binning might be necessary to measure a non-zero value. NP models can only contribute to CP -violating observables in rare charm decays if they extend the SM by additional complex-valued couplings. Correlations between CP asymmetries in rare semileptonic decays and purely hadronic modes, such as the recently observed difference of CP asymmetries in D 0 → K + K − and D 0 → π + π − decays, ∆A CP , 124 were studied in Ref. 123 in the context of a generic Z -extension with generationdependent charges. A SM prediction for ∆A CP is not well-established, since different theoretical approaches give predictions differing by a factor ∼ 10. [125][126][127][128][129][130][131][132][133][134] In Ref. 123 it was shown that NP contributions of order 10 −3 in ∆A CP can induce large imaginary parts of C (µ) 9,10 , and vice-versa. This leads to measurable CP asymmetries in semi-leptonic decays at the percent level, as shown in the left of Fig. 5 Additional studies on CP asymmetries, mainly within specific BSM models, can be seen in e.g. Refs. 54, 56, 59, 65, 66, 122.

Experimental investigations of angular and CP asymmetries
The first measurement of asymmetries in semi-leptonic rare charm decays has been carried out by LHCb measuring angular and CP asymmetries in D 0 → π + π − µ + µ − and D 0 → K + K − µ + µ − decays using a data set corresponding to 5 fb −1 recorded during the years 2011 to 2016. 135 Among the possible angular asymmetries that can  Fig. 6. Measured forward-backward asymmetry of the dimuon system A FB , triple-product asymmetry A 2φ and CP asymmetry A CP for D 0 → π + π − µ + µ − (top) and D 0 → K + K − µ + µ − (bottom) decays. Due to a lack of statistics, no measurement has been performed in the grey shaded areas of dimuon mass. The blue line and band show the average and its uncertainty, respectively. Figures are taken from Ref. 135. be constructed, the A FB of the dimuon system, defined as and the triple-product asymmetry, A 2φ , defined as have been investigated. 135 See Sect. 6.2 for the definitions of the angles g . Comparing Eq. (35) and (36), A FB and A 2φ are related to angular observables I 6 and A 9 , respectively, making them SM null tests. 28,61 In addition, the CP asymmetry as defined in Eq. (39) has been measured. Experimentally, the flavour of the D 0 mesons at the moment of their production is determined by selecting neutral charm mesons arising from the decay chain D * + → D 0 π + , where the charge of the accompanying low-momentum pion unambiguously indicates the flavour of the D 0 mesons. The measured asymmetry is corrected for nuisance charge asymmetries introduced by asymmetric detection efficiencies for positively and negatively charged pions, and for asymmetric production rates of D * + and D * − mesons in pp collisions. Furthermore, the asymmetries are corrected for phase-space dependent variations of the total reconstruction and selection efficiencies. The investigated asymmetries are found to be 135 for the decay mode D 0 → π + π − µ + µ − , and for D 0 → K + K − µ + µ − . The first uncertainty is statistic and the second one systematic. The asymmetries are consistent with zero and therefore compatible with the SM expectations. The asymmetries have also been investigated as a function of dimuon mass to enhance the sensitivity to NP contributions. The measured asymmetries in bins of dimuon mass are shown in Fig. 6. No dependence of the asymmetries on dimuon mass is found. The precision is limited due to low statistics and reaches a level where NP predictions start. Conceptually, it is the first measurement of asymmetries in semi-leptonic rare charm decays and we see a large potential in future measurements of a similar kind.

Rare radiative charm decays
Complementary information with respect to semi-leptonic decays on NP couplings can be obtained from radiative decays.  136,137 where leading power corrections ∼ 1/m c are computed. The approach is limited by large uncertainties on hadronic parameters, which could be constrained by a measurement of the branching fraction of D + → ρ + γ in the future. Second, a hybrid approach combining heavy quark effective theory and chiral Lagrangian, 43,44 which leads to predictions comparable to calculations of previous work, see Refs 138,139 for details. Table 8 shows the SM predictions for the branching fractions of D 0 → φγ, D 0 → K * 0 γ and D 0 → ρ 0 γ decays as obtained in the two main approaches of Ref. 77. The SM branching fractions are ∼ 10 −5 − 10 −4 and therefore approximately one to three orders of magnitude above those of resonant-dominated semi-leptonic decays discussed in Sect. 6.1 -6.4.
CP -asymmetries in radiative decays constitute SM null tests. NP models can induce large CP asymmetries up to ≤ 10%. 77,140,141 Furthermore, the angular distribution of the photon encodes information on its chirality 63, 142-144 and provides

Experimental investigations of radiative decays
The most recent experimental study of the previously discussed decay topologies where the uncertainties are statistical and systematic, respectively. The measured branching fractions of the decays D 0 → K * 0 γ and D 0 → φγ are consistent with previous measurements by BaBar 148 and the world average. The branching fraction of D 0 → ρ 0 γ is slightly larger than predicted in most theoretical calculations, 44,77,138 indicating the poor convergence of the 1/m c and α s expansion. See Table 8 for a comparison, where we also add the results of Ref. 148 for completeness. In addition, CP asymmetries in these decays have been measured by Belle, representing the first investigation of CP asymmetries in radiative charm decays to date. The flavour of the D 0 meson has been determined by selecting D 0 decays of charged D * + mesons and charge dependent detection and reconstruction efficiency effects are corrected for using data-driven approaches. The CP asymmetries are measured to be 147 where the first and second uncertainties are due to statistic and systematic nature, respectively. At the current level of statistical precision, the measured asymmetries are all compatible with SM expectations and leave large room for NP effects in future measurements. Angular distributions have not yet been studied experimentally.

Future prospects
We expect significant experimental improvements in the field in the near future, mainly driven by the LHCb and Belle II collaborations. To date, the LHCb detector has recorded a data set of pp collisions corresponding to 9 fb −1 during 2011-2018. At the moment, the LHCb detector is undergoing a major upgrade with data taking restart originally planned in 2021. It has been planned to record up to 50 (23) fb −1 by 2030 (2024). Due to the COVID-19 situation, currently, a restart in 2022 is planned. 150 In parallel, an upgrade phase II is in preparation, lead by the ambitious goal to collect up to 300 fb −1 by 2038. The Belle II detector has started data taking and plans to record a data set of e + e − collisions corresponding to 50 ab −1 by the end of 2030. 151,152 We expect updated measurements of many searches for rare and forbidden decay modes presented in Sect. 4 by the current flavour experiments BES III, Belle II and LHCb. As an example, we summarise future sensitivities of upgrade LHCb for the selected benchmark channels D 0 → µ + µ − , D + → π + µ + µ − and Λ + c → p µ + µ − as taken from Ref. 109 in Table 9, where conservatively possible improvements in the detector performance have been neglected in the calculations. Using the most recent experimental results, 45 we add projections for the decay channels D + → π + e + e − and D + → π + e + µ − to the table by scaling the observed limits to 23 fb −1 and 300 fb −1 of integrated luminosity h . In particular for semi-leptonic decays involving muons, the projected limits will come close the expected resonant contributions across the full decay phase space. Searches for forbidden decays will reach upper limits below the allowed parameter space of NP modes and therefore have the fantastic potential to find NP or significantly reduce its parameters. Preliminary results 153 of the BES III collaboration indicate the potential to make significant contributions to investigations of final states with two electrons. The future limits in Table 9 translate into the following model-independent limits on the Wilson coefficients: for 23 fb −1 (300 fb −1 ). Improved bounds on C (e) S, P require a more precise limit on the branching fraction of D 0 → e + e − . Projections for future D 0 → e + e − branching Table 9. Estimated upper limits (UL) of selected rare and forbidden decay modes at LHCb for future data sets, taken from Ref. 109. Limits for the decay channels D + → π + e + e − and D + → π + e + µ− have been obtained by scaling the observed limits taken from Ref. 45 to 23 fb −1 and 300 fb −1 of integrated luminosity, assuming the upper limit to scale with the square root of the integrated luminosity.
We summarise the expected uncertainties on CP asymmetries in rare radiative decays at Belle II in Table 10, which correspond to an improvement of the statistical precision of more than a factor of seven with respect to the current best measurements, which will help to set stringent limits on NP models. We also expect the Belle II collaboration to be capable to investigate angular distributions of rare radiative decays or to study additional topologies, such as D 0 → P 1 P 2 γ (P 1,2 = K, π). Both options have recently been suggested by theory, 63,77,143 however, experimental results are not yet available. Efforts including channels with photons and electrons are also expected to be intensified by LHCb in the future. 109 For processes involving LFV with electrons and muons, the future limits in Table 9 on D + → π + e + µ − are expected to improve the current bounds on Wilson coefficients (see Eq. (25)) by a factor of 0.5 (0.3) at 23 fb −1 ( 300 fb −1 ). Future im-provements on the limit of D 0 → e ± µ ∓ are also desired to improve bounds on K ( ) S,P . With the exception of the decay mode D 0 → τ ± e ∓ , information on Wilson coefficients involving τ leptons are kinematically not accessible in charm hadron decays. Analyses of high-p T data, 114 recorded by the ATLAS 112 and CMS 113 detectors, will provide additional and complementary constraints. Refs. 114 gives estimations for future bounds on τ Wilson coefficients extrapolating data collected by the ATLAS detector to a data set corresponding to 3 ab −1 (2038) 154 as: In the longer term, a future circular e + e − collider (FCC-ee) 155 is planned as the next particle collider generation at CERN. Given the benchmark number of ∼ 550× 10 9 produced cc pairs at a center-of-mass energy corresponding to the mass of the Z boson, 155 fantastic ways to probe the SM in rare charm will open. In particular, the clean environment of a e + e − collider is well suited for studies of missing energy decay modes, such as dineutrino final states. First studies of dineutrino decays are already possible today at current e + e − colliders, such as BES III and Belle II. Recently, Refs. 76, 156 proposed the possibility to test charged lepton flavour conservation (cLFC) and LU with dineutrino decays. Using the current bounds on Wilson coefficients from decays of charged leptons, the experimental measurement of a branching fraction of dineutrino decays h c → F νν above 10 −5 would imply a clear sign of cLFC violation and therefore NP. Limits at the order of 10 −5 or less could be possible at the current and future e + e − colliders. Further details can be found in Refs. 76, 156.

Conclusions
In the past, NP searches in rare decays have mainly focused on K and B systems, and less attention has been devoted to rare charm decays. In this review, we present promising opportunities to test the SM in |∆c| = |∆u| = 1 processes. We discuss the current theoretical status and the most recent experimental measurements. NP searches are currently still possible in branching fraction measurements in restricted regions of the decay phase space, a window that might close soon given the expected sensitivities of current and future flavour experiments. However, we stress and advertise the possibility to define clean null-test observables in resonance-dominated rare and radiative decays. These allow for very clean NP searches with minimal uncertainties from hadronic effects and permit to fully exploit the available statistics of the decays. Theoretical and experimental exploration of rare charm decays have only started, and the field is expected to play a key role in the future of flavour physics as a complementary testing ground for NP searches. Bright prospects for current and future experimental flavour facilities such as LHCb, Belle II, BES III and the FCC-ee are outlined and open the door to a large and exciting new program in flavour physics. The charm system offers the unique possibility to discover NP in the up-type sector. If the anomalies presently seen in B physics are caused by NP, investigation of rare charm decays will play a crucial and complementary role in the understanding of its origin and nature.