EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2022-273 LHCb-PAPER-2022-029 December 21, 2022

Search for rare decays of ${D}^{0}$ mesons into two muons

LHCb collaboration

A search for the very rare ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay
is performed using data collected by the LHCb experiment in proton-proton
collisions at $\sqrt{s}=7$, 8 and $13~{}\text{\,Te\kern-1.00006ptV}$,
corresponding to an integrated luminosity of $9\text{\,fb}^{-1}$. The search
is optimised for ${D}^{0}$ mesons from
${{D}^{*+}}\\!\rightarrow{{D}^{0}}{{\pi}^{+}}$ decays but is also sensitive to
${D}^{0}$ mesons from other sources. No evidence for an excess of events over
the expected background is observed. An upper limit on the branching fraction
of this decay is set at
$\mathcal{B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})<3.1\times 10^{-9}$ at
a 90% CL. This represents the world’s most stringent limit, constraining
models of physics beyond the Standard Model.

Submitted to Phys. Rev. Lett.

© 2024 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence.

Processes with a change in quark flavour without a change in electric charge
are forbidden at the lowest order in the Standard Model (SM) of particle
physics. Flavour Changing Neutral Currents (FCNC) are additionally suppressed
by the Glashow-Iliopoulos-Maiani (GIM) mechanism [1]. FCNC have been
extensively studied in strange- and beauty-quark hadrons. In the charm sector
the GIM suppression is stronger because the mass differences between down-type
quarks are smaller than the ones between up-type quarks. These processes can
be enhanced by several orders of magnitude in new physics (NP) scenarios when
compared to the SM.

The ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay is among the most
interesting charm-hadrons decays,111Charge conjugate processes are implied
throughout. being fully leptonic and additionally suppressed by helicity
reasons. Its SM short-distance contribution is extremely suppressed, yielding
a branching fraction on the order of $10^{-18}$ [2]. Long-distance
contributions dominate through an intermediate two-photon state, and can be
estimated to be ${\cal B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})\simeq
2.7\times 10^{-5}{\cal B}({{D}^{0}}\\!\rightarrow\gamma\gamma)$, leading to a
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ branching fraction of at least
$3\times 10^{-13}$ [2]. The best upper limit on the
${{D}^{0}}\\!\rightarrow\gamma\gamma$ decay rate was set by the Belle
collaboration to be $8.5\times 10^{-7}$ at a 90% C.L. [3]; using the same
relation this turns into an upper limit on the long-distance contribution to
the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ branching fraction of
$2.3\times 10^{-11}$. The ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay
rate can be enhanced in many NP models [4]. In particular, it is interesting
to note the importance of ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays
as testing ground for models with leptoquarks proposed to explain deviations
from the SM observed in $B$ physics measurements [5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. In some of these models rare $B$
decays receive contributions at the loop level, while new particles could
mediate the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay at tree level.
This phenomenon has been extensively discussed in the literature [24, 25, 26,
27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. In addition, while in the
Minimal Supersymmetric Standard Model no sizeable contribution would come to
the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ rate [4], some Supersymmetric
models with $R$-parity symmetry violation would allow tree level contributions
[2, 39]. Recent discussions of predictions for
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays in such models can be found
in Refs. [40, 41]. In other explanations of the deviations in the $B$ physics
observables, limits on the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay
are used as additional constraint for the generic flavour couplings. In those
where additional vector bosons ($Z^{\prime}$) are introduced, the
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay usually does not give strong
constraints [42, 43]. Similarly, in other explanatory models [44, 45, 46] the
bound from ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays can be avoided.
Instead in some models with vector-like fermions [47], the
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay gives the strongest
constraint. Model-independent bounds on the Wilson coefficients related to
charm physics and in particular to ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$
decays have been set in Ref. [48, 49]. Furthermore, the
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ rate is correlated to the rate of
${{D}^{0}}-{{\kern 1.79993pt\overline{\kern-1.79993ptD}}{}^{0}}$ mixing in
many NP models [4]. This is of uttermost importance given the recent first
observation of the mass difference between neutral charm-meson eigen-states
[50].

The current world’s best limit on this decay is
$\mathcal{B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})<7.6\times 10^{-9}$ at
95% CL, and was obtained by the LHCb experiment exploiting about 0.9
$\text{\,fb}^{-1}$ of 2011 data [51]. The data used in Ref. [51] are also used
in this analysis and those results are superseded by this Letter.

This Letter presents a search for the
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay based on data collected by
the LHCb experiment in $p$ $p$ collisions corresponding to $9\text{\,fb}^{-1}$
of integrated luminosity. The data have been collected in 2011, 2012 (Run 1)
and 2015–2018 (Run 2) at $\sqrt{s}=7$, 8 and $13\text{\,Te\kern-1.00006ptV}$,
respectively. The ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay is
searched for using ${{D}^{*+}}\\!\rightarrow{{D}^{0}}{{\pi}^{+}}$ decays, as
this improves the background rejection and allows the yield of the decay to be
obtained from a two-dimensional fit to the dimuon invariant mass,
$m({\mu^{+}\mu^{-}})$, and the difference between the ${{D}^{*+}}$ and
${{D}^{0}}$ candidate masses, $\Delta m$. The yield is converted to the decay
branching fraction by normalising to two hadronic decays,
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$ and
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$, selected concurrently to the
signal (collectively referred to as ${{D}^{0}}\\!\rightarrow h^{+}h^{-}$).

The LHCb detector is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, described in detail in Refs. [52, 53]. The
simulated events used in this analysis are produced with the software
described in Refs. [54, *Sjostrand:2006za, 56, 57, 58, 59, 60, 61].

Events are selected online by a trigger that consists of a hardware stage,
which is based on information from the calorimeter and muon systems, followed
by two software stages. At the hardware trigger stage, events are required to
have a muon candidate with high transverse momentum, $p_{\mathrm{T}}$, or a
hadron, photon or electron candidate with high transverse energy in the
calorimeters. A first stage of the software trigger selects events with a muon
candidate, or a high $p_{\mathrm{T}}$ charged particle, or a combination of
two tracks, each of these displaced from the primary $pp$ collision vertex
(PV). In the second stage of the software trigger, dedicated algorithms select
candidate ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$,
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$ and
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays, combining two
oppositely charged tracks with suitable loose particle identification (PID)
requirements that form a secondary vertex separated from any PV. The invariant
mass of the ${D}^{0}$ candidate must be in the range
$1800–1940~{}\text{\,Me\kern-1.00006ptV\\!/}c^{2}$. To keep the same trigger
selection for the signal and normalisation channels, aside from PID, a scale
factor of order 0.2–3.0% depending on the data taking period is applied to the
normalisation channels trigger selections to limit their rate, keeping
randomly a fraction of the events. In the offline selection, only candidates
associated with a muon hardware trigger, or those where the rest of the event
contained a high transverse energy hadron or electron are kept for the signal
channel. For the normalisation channels, only candidates associated with a
high transverse energy hadron are kept.

In the offline analysis, ${D}^{0}$ candidates satisfying the trigger
requirements are formed with similar but more stringent criteria than the
second software-trigger stage. These ${D}^{0}$ candidates are then combined
with a charged particle originating from the same PV and having
$p_{\mathrm{T}}>110~{}\text{\,Me\kern-1.00006ptV\\!/}c$ to form
${{D}^{*}}^{+}\rightarrow{{D}^{0}}\pi^{+}$ candidates. To improve the mass
resolution, the ${{D}^{*+}}$ meson decay vertex is constrained to coincide
with the PV [62]. The candidate $\Delta m$ is required to be in the range
139.6–151.6 $\text{\,Me\kern-1.00006ptV\\!/}c^{2}$. A multivariate selection
based on a boosted decision tree (BDT) algorithm [63, 64, 65] is used to
suppress background from random combinations of charged particles, using as
input: the $p_{\mathrm{T}}$ of the pion from the ${D}^{*+}$ decay, the
smallest $p_{\mathrm{T}}$ and impact parameter significance with respect to
the PV of the ${D}^{0}$ decay products, the angle between the ${D}^{0}$
momentum and the vector connecting the primary and secondary vertices, and the
quality of the ${D}^{0}$ vertex. The BDT is trained separately for each Run of
data taking, using simulated decays as signal and data candidates from the
dimuon sample with
$m({\mu^{+}\mu^{-}})>1894~{}\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ as
background. The $k$-folding technique, with $k=9$, is applied [66]. The BDT
output ranges from 0 to 1, from background-like candidates to more signal-like
candidates. The BDT output is used to define three search regions:
$\text{BDT}\in[0.15,0.33]$, $[0.33,0.66]$, $[0.66,1.]$. The output of the same
BDT algorithm is computed also for the
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ candidates for calibration
purposes, but is not required in the selection.

A second source of background is due to two- and three-body ${D}^{0}$ decays,
with one or two hadrons misidentified as muons (e.g. ${{D}^{0}}\\!\rightarrow
h^{+}h^{-}$ or semileptonic decays). The misidentification occurs mainly for
hadrons that decay into a muon before the muon sub-detector. Although this
process is relatively rare, the large branching fractions of these modes
produce a background peaking in the signal region of the $m({\mu^{+}\mu^{-}})$
distribution that is partially suppressed by a multivariate muon
identification discriminant combining information from the Cherenkov
detectors, the calorimeters, and the muon sub-detector [67]. In addition, the
muon candidates are required to have associated muon chamber hits that are not
shared with any other track in the event. Signal and background sources can
originate both from ${{D}^{*+}}\\!\rightarrow{{D}^{0}}{{\pi}^{+}}$ decays or
from other sources (the PV or $B$ decays) combined with an unrelated pion
(untagged); both are taken into account in the yield estimation. A requirement
on the output of the multivariate muon identification discriminant is
simultaneously chosen for the three BDT regions, by optimising the sensitivity
to the minimum visible cross-section, as defined by an extension of the figure
of merit defined in Ref. [68]. Roughly one percent of events contain more than
one signal candidate after all selection requirements.

The signal yield is converted to the decay branching fraction by normalising
to the hadronic decays ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ and
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$, with branching fractions of
$(1.490\pm 0.027)\times 10^{-3}$ and $(3.999\pm 0.045)\times 10^{-2}$,
respectively [69], as

$\mathcal{B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})=\frac{{N}_{{{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}}}{N_{{{D}^{0}}\\!\rightarrow
h^{+}h^{-}}}\cdot\frac{\varepsilon_{h^{+}h^{-}}}{\varepsilon_{\mu^{+}\mu^{-}}}\cdot
s\cdot\mathcal{B}({{D}^{0}}\\!\rightarrow h^{+}h^{-})\equiv\alpha
N_{{{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}}$ (1)

where $\varepsilon$ is the efficiency and $N$ is the yield of the given
channel, $s$ is the scale factor of the normalisation channel and $\alpha$ is
defined as the single event sensitivity.

The efficiencies in Eq. (1) are factorised into different steps for ease of
estimation and evaluated with respect to the previous steps: detector
acceptance, reconstruction and selection, PID, and trigger.

The reconstruction and selection efficiencies are obtained from simulated
samples. The simulated candidates are assigned weights with an iterative
procedure that improves the agreement with data using the following variables:
pseudorapidity of the ${D}^{0}$ meson, transverse momentum of the ${D}^{0}$
meson and number of tracks in the event. It is verified that after weighting,
all variables used in the selection agree well between data and simulation.
The weights obtained from the ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$
candidates are used to correct also the signal simulation.

Possible residual differences between data and simulation in the tracking
efficiencies are determined using control channels in data [70]. The PID
efficiencies are determined from data using samples of kinematically
identified charged particles from
${{B}^{+}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{{K}^{+}}$ and
${{D}^{*+}}\\!\rightarrow{{D}^{0}}(\\!\rightarrow{{K}^{-}}{{\pi}^{+}}){{\pi}^{+}}$
decays [67], weighted to match the kinematic properties of the signal and the
normalisation channels, respectively. The efficiencies are determined in bins
of the $p$ and $p_{\mathrm{T}}$ of the tracks. A total systematic uncertainty
of 1–3% is associated to the binning scheme and background determination in
the calibration samples.

The efficiency of the second level of the software trigger is unity with
respect to the offline-selected candidates by construction, as the selection
is tighter in every requirement. The hardware and first level software trigger
efficiencies are evaluated with the Tistos method [71] in data. For the signal
channel, the
${{B}^{+}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}(\\!\rightarrow{\mu^{+}}{\mu^{-}}){{K}^{+}}$
decay is used as the calibration channel, selected with the same requirements
as those used for the analysis of $B$ decays into two muons [72, 73]. The
calibration is performed in intervals of the
${J\mskip-3.0mu/\mskip-2.0mu\psi}$ $p_{\mathrm{T}}$ and pseudorapidity. For
each interval, a scaling factor between data and simulation is obtained and
applied to the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ simulation.
Compatible results are obtained repeating the calibration in intervals of
$p_{\mathrm{T}}$ and the maximum $\rm{IP}\chi^{2}$ of the muons, where
$\rm{IP}\chi^{2}$ is defined as the difference in the vertex-fit $\chi^{2}$ of
a given PV reconstructed with and without the track under consideration. The
typical scaling between data and simulation deviates from unity by 2–6%. The
normalisation channels, given their high yields, are self-calibrated. The
Tistos method is applied to the ${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$
and ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ channels and trigger
efficiencies are obtained. To minimise cross-correlation biases, only
candidates in events that satify a muon trigger independently of the candidate
are used as calibration sample. The calibration of the hadronic hardware
trigger is also validated with independent estimates based on control samples
in data, obtained with similar methods as in Ref. [74], from which a 15%
relative systematic uncertainty is assigned to the hadronic trigger efficiency
calibration.

The efficiency of the BDT requirement, and the signal fraction in the BDT
intervals, are calibrated in data by applying the same estimator to the
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decay, which is topologically
very similar to the signal. The distribution of the BDT output is obtained in
background subtracted ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays
in data and simulation, and found to be compatible, as shown in the
Supplemental Material [75]. A small correction is determined and applied to
the signal; its uncertainty is assigned as systematic uncertainty to the
signal efficiency.

The yields of the normalisation channels are obtained through a fit to the
$\Delta m$ distribution (Fig. 1), requiring the reconstructed ${D}^{0}$ mass
to be within $\pm 10\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ of the known
${D}^{0}$ mass. The signal probability distribution function is composed of a
sum of a Gaussian and a Crystal-Ball function [76] with power-law tails on
both sides. The background is described with a threshold function, as defined
in Ref. [51]. The parameters of the Crystal-Ball function are estimated with
simulation: the power of the tail is fixed, while the position where the power
tails start may vary freely in the fit. In addition, the signal width and all
background parameters are left free in the fit.

Figure 1: Distributions of $\Delta m$ for (left)
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$ and (right)
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ normalisation channels
candidates for (top) Run 1 and (bottom) Run 2 data. The distributions are
superimposed with the fit.

Using Eq. (1), values of $\alpha$ for both normalisation channels are
obtained, and found to be in good agreement with each other for each data
taking run and for the full sample. As an additional cross-check, the ratio of
the efficiency corrected yields of the two normalisation channels is obtained
and compared to the ratio of their branching fractions. The value is stable
across the data taking years and compatible with the world average [69]. The
average single event sensitivity is found to be $\alpha=(2.15\pm 0.34)\times
10^{-11}$, corresponding to at most one expected signal
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay under the SM hypothesis.

The signal yield is obtained through an unbinned maximum-likelihood fit to the
two-dimensional distribution of $m({\mu^{+}\mu^{-}})$ and $\Delta m$,
performed simultaneously in the three BDT intervals and in the two data taking
Runs. The distributions projected onto the two variables are shown in Fig. 2.
Each of the two projections is selected using only candidates in the signal
region of the other variable, where the signal regions are defined as
$m({{D}^{0}})\in[1840,1885]\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ and $\Delta
m\in[144.9,146.1]\text{\,Me\kern-1.00006ptV\\!/}c^{2}$, respectively. The full
distributions can be seen in the Supplemental Material [75]. The two variables
are treated as uncorrelated, as the correlations in all fit components are
verified to be negligible. After the full selection, only combinatorial and
misidentified hadronic ${D}^{0}$ decays are found to contribute to the
background. Other sources, such as semileptonic decays, are found to be
negligible and well modelled by the same shape as the combinatorial
background.

The shape of the signal and misidentified background
(${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$,
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$) distributions in the two fit
variables is obtained from simulation, reconstructed as
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays. The model parameters are
determined separately for Run 1 and Run 2; the resulting PDF describes the
distribution in each data taking year and BDT interval well. For signal, the
$m({\mu^{+}\mu^{-}})$ and $\Delta m$ distributions are both parametrised by
Crystal-Ball [76] functions with power tails on both sides. For
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays, a single Crystal-ball
function in $m({\mu^{+}\mu^{-}})$ and the sum of a Johnson [77] and Gaussian
function in $\Delta m$ are employed. For
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$ decays, a Johnson function is
used for the $m({\mu^{+}\mu^{-}})$ distribution, while the $\Delta m$
distribution is described by three Gaussian functions. The combinatorial
background is described by an exponential function in $m({\mu^{+}\mu^{-}})$
and a threshold function in $\Delta m$ [51]. Untagged signal and
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ components are included and
parametrised as their respective tagged component in $m({\mu^{+}\mu^{-}})$ and
with the same threshold function of the combinatorial background in $\Delta
m$. The fraction of this component is fixed to the value determined in each
BDT interval from a fit to ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$
data. The shape parameters obtained from the simulated samples are fixed in
the data fit, while the slope of the exponential of the combinatorial
background is left free to vary in each BDT interval.

A constraint on the expected number of misidentified
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays is determined from a
dedicated, high-statistics, simulation sample with the trigger and offline
selection applied. The most critical part of the simulated sample is the PID
efficiency due to the presence of a large fraction of $\pi\rightarrow\mu$
decays that mimic the signal and are not considered in the standard
calibration tools. The PID efficiency is obtained from simulation but it is
cross-checked using ${{D}^{+}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{+}}$
and ${{D}^{+}_{s}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{+}}$ control
samples in data where same-sign pions are weighted to match the kinematics of
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays. The agreement between
the PID efficiency determined with both methods is satisfactory over the full
range of the muon identification discriminant variable [75]. Therefore, no
systematic uncertainty is assigned on this estimate.

The yield of the misidentified ${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$
decays is constrained from an auxiliary fit to the $m({\mu^{+}\mu^{-}})$
sideband data, recomputed with the correct mass hypothesis. The fit is
performed using the $\Delta m$ distribution within a $\pm
10\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ region around the ${D}^{0}$ mass in
the ${{K}^{-}}{{\pi}^{+}}$ mass hypothesis. A correction is applied to take
into account this mass requirement. The correlation between this estimate and
the yield in the final fit is found not to influence the estimate of the
signal branching fraction.

The systematic uncertainties related to both the normalisation, through
$\alpha$, and the background shapes and yields, are included in the fit as
Gaussian constraints on the relevant parameters. The dominant systematic
uncertainty comes from the calibration of the hadronic trigger efficiency,
which is shared through auxiliary parameters among the normalisation channels,
and also with the misidentified
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ yields that depend on the same
estimate. The fit procedure is tested with pseudoexperiments. The values of
the floating shape parameters are obtained from the data fit. Unbiased
estimates of the branching fraction with correct coverage are obtained.

Figure 2: Distribution of (left) $m({\mu^{+}\mu^{-}})$ and (right) $\Delta m$
for the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ candidates in data from
(top) Run 1 and (bottom) Run 2, for the most sensitive BDT interval. The
distribution is superimposed with the fit to data. Each of the two
distributions is in the signal region of the other variable, see text for
details. Untagged and tagged decays are included in a single component for
signal and ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ background.

The $m({\mu^{+}\mu^{-}})$ and $\Delta m$ distributions in data are shown for
the most sensitive BDT interval in Fig. 2 and for all intervals in Ref. [75],
overlaid with the result of the fit. The data are consistent with the expected
background. The value obtained for the
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ branching fraction is
$\mathcal{B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})=(1.7\pm 1.0)\times
10^{-9}$, corresponding to $79\pm 45$ signal decays. The significance of this
signal is estimated comparing the test statistics in data with the
distribution of the test statistics in background-only pseudoexperiments, and
is found to have a $p$-value of $0.068$, corresponding to a significance of
$1.5\sigma$ (see also Ref. [75]). An upper limit on the branching fraction is
derived using the frequentist $\rm{CL_{s}}$ method [78] as implemented in the
GammaCombo framework [79, 80]. This yields

$\mathcal{B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})<3.1\,(3.5)\times
10^{-9}\text{ at 90 (95)\% CL }.$

The observed limit is larger than the one expected from background-only
pseudoexperiments,
$\mathcal{B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})<\mbox{$1.9\,(2.3)\times
10^{-9}$}\text{ at a 90 (95)\% CL}$, coherently with the central value for the
signal branching fraction.

The fit is repeated with different configurations: allowing the resolution of
the misidentified ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ background
to vary, using a double exponential function in place of a single one for the
combinatorial background, and reducing the range in the $\Delta m$ variable.
No significant change was found in the signal branching fraction with any
configuration.

In summary, a search for the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay
in data corresponding to $9\text{\,fb}^{-1}$ of $pp$ collision data collected
by the LHCb experiment is performed. No excess with respect to the background
expectation has been found and an upper limit of
$\mathcal{B}({{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}})<3.1\times 10^{-9}$ at
90% CL has been set. This result represents an improvement of more than a
factor two with respect to the previous LHCb result. This measurement
constitutes the most stringent limit on the relevant FCNC couplings in the
charm sector, allowing to set additional constraints on physics models beyond
the SM which predict the branching fractions of
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ and describe results from $B$
physics measurements.

## References

  * [1] S. L. Glashow, J. Iliopoulos, and L. Maiani, _Weak interactions with lepton-hadron symmetry_ , Phys. Rev. D2 (1970) 1285
  * [2] G. Burdman, E. Golowich, J. L. Hewett, and S. Pakvasa, _Rare charm decays in the Standard Model and beyond_ , Phys. Rev. D66 (2002) 014009, arXiv:hep-ph/0112235
  * [3] Belle collaboration, N. K. Nisar et al., _Search for the rare decay $D^{0}\rightarrow\gamma\gamma$ at Belle_, Phys. Rev. D93 (2016) 051102, arXiv:1512.02992
  * [4] E. Golowich, J. Hewett, S. Pakvasa, and A. A. Petrov, _Relating $D^{0}-\bar{D^{0}}$ mixing and $D^{0}\rightarrow l^{+}l^{-}$ with new physics_, Phys. Rev. D79 (2009) 114030, arXiv:0903.2830
  * [5] LHCb collaboration, R. Aaij et al., _Test of lepton universality with ${{B}^{0}}\\!\rightarrow{{K}^{*0}}\ell^{+}\ell^{-}$ decays_, JHEP 08 (2017) 055, arXiv:1705.05802
  * [6] LHCb collaboration, R. Aaij et al., _Test of lepton universality in beauty-quark decays_ , arXiv:2103.11769, to appear in Nature Physics
  * [7] LHCb collaboration, R. Aaij et al., _Test of lepton universality using ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow p{{K}^{-}}\ell^{+}\ell^{-}$ decays_, JHEP 05 (2020) 040, arXiv:1912.08139
  * [8] LHCb collaboration, R. Aaij et al., _Angular analysis of the $B^{+}\rightarrow K^{\ast+}\mu^{+}\mu^{-}$ decay_, Phys. Rev. Lett. 126 (2021) 161802, arXiv:2012.13241
  * [9] LHCb collaboration, R. Aaij et al., _Measurement of $C\\!P$-averaged observables in the ${{B}^{0}}\\!\rightarrow{{K}^{*0}}{\mu^{+}\mu^{-}}$ decay_, Phys. Rev. Lett. 125 (2020) 011802, arXiv:2003.04831
  * [10] LHCb collaboration, R. Aaij et al., _Angular analysis of the ${{B}^{0}}\\!\rightarrow{{K}^{*0}}{\mu^{+}\mu^{-}}$ decay using $3\text{\,fb}^{-1}$ of integrated luminosity_, JHEP 02 (2016) 104, arXiv:1512.04442
  * [11] LHCb collaboration, R. Aaij et al., _Angular analysis of the rare decay ${{B}^{0}_{s}}\rightarrow\phi{\mu^{+}}{\mu^{-}}$_, JHEP 11 (2021) 043, arXiv:2107.13428
  * [12] ATLAS collaboration, M. Aaboud et al., _Angular analysis of $B^{0}_{d}\rightarrow K^{*}\mu^{+}\mu^{-}$ decays in $pp$ collisions at $\sqrt{s}=8$ TeV with the ATLAS detector_, JHEP 10 (2018) 047, arXiv:1805.04000
  * [13] BaBar collaboration, B. Aubert et al., _Measurements of branching fractions, rate asymmetries, and angular distributions in the rare decays $B\rightarrow K\ell^{+}\ell^{-}$ and $B\rightarrow K^{*}\ell^{+}\ell^{-}$_, Phys. Rev. D73 (2006) 092001, arXiv:hep-ex/0604007
  * [14] BaBar collaboration, J. P. Lees et al., _Measurement of angular asymmetries in the decays $B\rightarrow K^{*}\ell^{+}\ell^{-}$_, Phys. Rev. D93 (2016) 052015, arXiv:1508.07960
  * [15] Belle collaboration, J.-T. Wei et al., _Measurement of the differential branching Fraction and forward-backward Asymmetry for $B\rightarrow K^{(*)}\ell^{+}\ell^{-}$_, Phys. Rev. Lett. 103 (2009) 171801, arXiv:0904.0770
  * [16] Belle collaboration, S. Wehle et al., _Lepton-flavor-dependent angular analysis of $B\rightarrow K^{\ast}\ell^{+}\ell^{-}$_, Phys. Rev. Lett. 118 (2017) 111801, arXiv:1612.05014
  * [17] CDF collaboration, T. Aaltonen et al., _Measurements of the angular distributions in the decays $B\rightarrow K^{(*)}\mu^{+}\mu^{-}$ at CDF_, Phys. Rev. Lett. 108 (2012) 081807, arXiv:1108.0695
  * [18] CMS collaboration, V. Khachatryan et al., _Angular analysis of the decay $B^{0}\rightarrow K^{*0}\mu^{+}\mu^{-}$ from pp collisions at $\sqrt{s}=8$ TeV_, Phys. Lett. B753 (2016) 424, arXiv:1507.08126
  * [19] CMS collaboration, A. M. Sirunyan et al., _Measurement of angular parameters from the decay $\mathrm{B}^{0}\rightarrow\mathrm{K}^{*0}\mu^{+}\mu^{-}$ in proton-proton collisions at $\sqrt{s}=$ 8 TeV_, Phys. Lett. B781 (2018) 517, arXiv:1710.02846
  * [20] LHCb collaboration, R. Aaij et al., _Measurements of the S-wave fraction in ${{B}^{0}}\\!\rightarrow{{K}^{+}}{{\pi}^{-}}{\mu^{+}\mu^{-}}$ decays and the ${{B}^{0}}\\!\rightarrow{{K}^{*}}(892)^{0}{\mu^{+}\mu^{-}}$ differential branching fraction_, JHEP 11 (2016) 047, Erratum ibid. 04 (2017) 142, arXiv:1606.04731
  * [21] LHCb collaboration, R. Aaij et al., _Branching fraction measurements of the rare $B^{0}_{s}\rightarrow\phi\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow f_{2}^{\prime}(1525)\mu^{+}\mu^{-}$ decays_, Phys. Rev. Lett. 127 (2021) 151801, arXiv:2105.14007
  * [22] LHCb collaboration, R. Aaij et al., _Differential branching fractions and isospin asymmetries of ${B}\\!\rightarrow K^{(*)}{\mu^{+}\mu^{-}}$ decays_, JHEP 06 (2014) 133, arXiv:1403.8044
  * [23] LHCb collaboration, R. Aaij et al., _Differential branching fraction and angular analysis of ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}{\mu^{+}\mu^{-}}$ decays_, JHEP 06 (2015) 115, Erratum ibid. 09 (2018) 145, arXiv:1503.07138
  * [24] M. Bauer and M. Neubert, _Minimal leptoquark explanation for the R ${}_{D^{(*)}}$ , RK , and $(g-2)_{g}$ anomalies_, Phys. Rev. Lett. 116 (2016) 141802, arXiv:1511.01900
  * [25] M. Freytsis, Z. Ligeti, and J. T. Ruderman, _Flavor models for $\bar{B}\rightarrow D^{(*)}\tau\bar{\nu}$_, Phys. Rev. D92 (2015) 054018, arXiv:1506.08896
  * [26] R. Benbrik and C.-H. Chen, _Leptoquark on $P\rightarrow\ell^{+}\nu$, FCNC and LFV_, Phys. Lett. B672 (2009) 172, arXiv:0807.2373
  * [27] S. Fajfer and N. Košnik, _Leptoquarks in flavor changing neutral current charm decays_ , Phys. Rev. D79 (2009) 017502
  * [28] S. de Boer and G. Hiller, _Flavor and new physics opportunities with rare charm decays into leptons_ , Phys. Rev. D93 (2016) 074001, arXiv:1510.00311
  * [29] D. Bečirević, N. Košnik, O. Sumensari, and R. Zukanovich Funchal, _Palatable leptoquark scenarios for lepton flavor violation in exclusive $b\rightarrow s\ell_{1}\ell_{2}$ modes_, JHEP 11 (2016) 035, arXiv:1608.07583
  * [30] Y. Cai, J. Gargalionis, M. A. Schmidt, and R. R. Volkas, _Reconsidering the One Leptoquark solution: flavor anomalies and neutrino mass_ , JHEP 10 (2017) 047, arXiv:1704.05849
  * [31] K. Kowalska, E. M. Sessolo, and Y. Yamamoto, _Constraints on charmphilic solutions to the muon g-2 with leptoquarks_ , Phys. Rev. D99 (2019) 055007, arXiv:1812.06851
  * [32] I. Bigaran, J. Gargalionis, and R. R. Volkas, _A near-minimal leptoquark model for reconciling flavour anomalies and generating radiative neutrino masses_ , JHEP 10 (2019) 106, arXiv:1906.01870
  * [33] I. Doršner et al., _Physics of leptoquarks in precision experiments and at particle colliders_ , Phys. Rept. 641 (2016) 1, arXiv:1603.04993
  * [34] I. Doršner, S. Fajfer, and M. Patra, _A comparative study of the $S_{1}$ and $U_{1}$ leptoquark effects in the light quark regime_, Eur. Phys. J. C80 (2020) 204, arXiv:1906.05660
  * [35] I. Doršner, S. Fajfer, and O. Sumensari, _Muon $g-2$ and scalar leptoquark mixing_, JHEP 06 (2020) 089, arXiv:1910.03877
  * [36] K. Kowalska, E. M. Sessolo, and Y. Yamamoto, _Flavor anomalies from asymptotically safe gravity_ , Eur. Phys. J. C81 (2021) 272, arXiv:2007.03567
  * [37] M. Bordone, O. Catà, T. Feldmann, and R. Mandal, _Constraining flavour patterns of scalar leptoquarks in the effective field theory_ , JHEP 03 (2021) 122, arXiv:2010.03297
  * [38] F. Abudinén, T. Blake, U. Egede, and T. Gershon, _Prospects for studies of ${D^{*0}}\rightarrow{\mu^{+}}{\mu^{-}}$ and ${B^{{*0}}_{(s)}}\rightarrow{\mu^{+}}{\mu^{-}}$ decays_, Eur. Phys. J. C 82 (2022) 459, arXiv:2202.03916
  * [39] G. Burdman and I. Shipsey, _$D^{0}$ \- $\bar{D}^{0}$ mixing and rare charm decays_, Ann. Rev. Nucl. Part. Sci. 53 (2003) 431, arXiv:hep-ph/0310076
  * [40] R.-M. Wang et al., _Decays $D^{+}_{(s)}\rightarrow\pi(K)^{+}\ell^{+}\ell^{-}$ and $D^{0}\rightarrow\ell^{+}\ell^{-}$ in the MSSM with and without R-parity_, Int. J. Mod. Phys. A30 (2015) 1550063, arXiv:1409.0181
  * [41] W. Altmannshofer, P. S. B. Dev, A. Soni, and Y. Sui, _Addressing R ${}_{D^{(*)}}$, R${}_{K^{(*)}}$, muon $g-2$ and ANITA anomalies in a minimal $R$-parity violating supersymmetric framework_, Phys. Rev. D102 (2020) 015031, arXiv:2002.12910
  * [42] Z. Kang and Y. Shigekami, _$(g-2)_{\mu}$ versus flavor changing neutral current induced by the light $(B-L)_{\mu\tau}$ boson_, JHEP 11 (2019) 049, arXiv:1905.11018
  * [43] A. K. Alok, A. Dighe, S. Gangal, and D. Kumar, _Predictions for $B_{s}\rightarrow{\bar{K}}^{*}\ell\,\ell$ in non-universal $Z^{\prime}$ models_, Eur. Phys. J. C80 (2020) 682, arXiv:1912.02052
  * [44] D. Guadagnoli, M. Reboud, and O. Sumensari, _A gauged horizontal $SU(2)$ symmetry and $R_{K^{(\ast)}}$_, JHEP 11 (2018) 163, arXiv:1807.03285
  * [45] R. Alonso et al., _A clockwork solution to the flavor puzzle_ , JHEP 10 (2018) 099, arXiv:1807.09792
  * [46] A. K. Alok et al., _New-physics signals of a model with a vector-singlet up-type quark_ , Phys. Rev. D92 (2015) 013002, arXiv:1504.00517
  * [47] K. Ishiwata, Z. Ligeti, and M. B. Wise, _New vector-like fermions and flavor physics_ , JHEP 10 (2015) 027, arXiv:1506.03484
  * [48] S. Fajfer and N. Kosnik, _Prospects of discovering new physics in rare charm decays_ , Eur. Phys. J. C75 (2015) 567, arXiv:1510.00965
  * [49] A. Bharucha, D. Boito, and C. Méaux, _Disentangling QCD and new physics in $D^{+}\rightarrow\pi^{+}\ell^{+}\ell^{-}$_, JHEP 04 (2021) 158, arXiv:2011.12856
  * [50] LHCb collaboration, R. Aaij et al., _Observation of the mass difference between neutral charm-meson eigenstates_ , Phys. Rev. Lett. 127 (2021) 111801, arXiv:2106.03744
  * [51] LHCb collaboration, R. Aaij et al., _Search for the rare decay ${{D}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}}$_, Phys. Lett. B725 (2013) 15, arXiv:1305.5059
  * [52] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [53] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [54] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [55] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [56] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [57] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [58] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [59] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [60] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [61] N. Davidson, T. Przedzinski, and Z. Was, _PHOTOS interface in C++: Technical and physics documentation_ , Comp. Phys. Comm. 199 (2016) 86, arXiv:1011.0937
  * [62] W. D. Hulsbergen, _Decay chain fitting with a Kalman filter_ , Nucl. Instrum. Meth. A552 (2005) 566, arXiv:physics/0503191
  * [63] L. Breiman, J. H. Friedman, R. A. Olshen, and C. J. Stone, Classification and regression trees, Wadsworth international group, Belmont, California, USA, 1984
  * [64] Y. Freund and R. E. Schapire, _A decision-theoretic generalization of on-line learning and an application to boosting_ , J. Comput. Syst. Sci. 55 (1997) 119
  * [65] H. Voss, A. Hoecker, J. Stelzer, and F. Tegenfeldt, _TMVA - Toolkit for Multivariate Data Analysis with ROOT_ , PoS ACAT (2007) 040
  * [66] A. Blum, A. Kalai, and J. Langford, _Beating the hold-out: bounds for k-fold and progressive cross-validation_ , in Proceedings of the twelfth annual conference on computational learning theory,  COLT ’99, (New York, NY, USA), 203, ACM, 1999
  * [67] R. Aaij et al., _Selection and processing of calibration samples to measure the particle identification performance of the LHCb experiment in Run 2_ , Eur. Phys. J. Tech. Instr. 6 (2019) 1, arXiv:1803.00824
  * [68] G. Punzi, _Sensitivity of searches for new signals and its optimization_ , eConf C030908 (2003) MODT002, arXiv:physics/0308063
  * [69] Heavy Flavor Averaging Group, Y. Amhis et al., _Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of 2018_, Eur. Phys. J. C81 (2021) 226, arXiv:1909.12524, updated results and plots available at https://hflav.web.cern.ch
  * [70] LHCb collaboration, R. Aaij et al., _Measurement of the track reconstruction efficiency at LHCb_ , JINST 10 (2015) P02007, arXiv:1408.1251
  * [71] S. Tolk, J. Albrecht, F. Dettori, and A. Pellegrino, _Data driven trigger efficiency determination at LHCb_ , LHCb-PUB-2014-039. CERN-LHCb-PUB-2014-039, CERN, Geneva, 2014
  * [72] LHCb collaboration, R. Aaij et al., _Analysis of neutral $B$-meson decays into two muons_, Phys. Rev. Lett. 128 (2022) 041801, arXiv:2108.09284
  * [73] LHCb collaboration, R. Aaij et al., _Measurement of the ${{B}^{0}_{s}}\rightarrow{\mu^{+}}{\mu^{-}}$ decay properties and search for the ${{B}^{0}}\rightarrow{\mu^{+}}{\mu^{-}}$ and ${{B}^{0}_{s}}\rightarrow{\mu^{+}}{\mu^{-}}\gamma$ decays_, Phys. Rev. D105 (2022) 012010, arXiv:2108.09283
  * [74] A. Martin Sanchez, P. Robbe, and M.-H. Schune, _Performances of the LHCb L0 Calorimeter Trigger_ , LHCb-PUB-2011-026, CERN-LHCb-PUB-2011-026, CERN, Geneva, 2012
  * [75] See Supplemental Material at the end of this Letter for additional figures.
  * [76] T. Skwarnicki, A study of the radiative cascade transitions between the Upsilon-prime and Upsilon resonances, PhD thesis, Institute of Nuclear Physics, Krakow, 1986, DESY-F31-86-02
  * [77] N. L. Johnson, _Systems of frequency curves generated by methods of translation_ , Biometrika 36 (1949) 149
  * [78] A. L. Read, _Presentation of search results: The CL S technique_, J. Phys. G28 (2002) 2693
  * [79] M. Kenzie et al., _GammaCombo: A statistical analysis framework for combining measurements, fitting datasets and producing confidence intervals_ ,  doi: 10.5281/zenodo.3371421
  * [80] LHCb collaboration, R. Aaij et al., _Measurement of the CKM angle $\gamma$ from a combination of LHCb results_, JHEP 12 (2016) 087, arXiv:1611.03076

## Supplemental Material

This section contains the additional figures mentioned in the main text.
Figures 3 and 4 present the same data shown in Figure 1 in the main body, but
for all the BDT intervals. Figures 5 and 6 represent the same distributions
and subdivisions of the previous plots, but without restricting each to the
signal region of the other variable, i.e. they contain the full data used for
the signal search. Figure 7 shows the data in the two-dimensional plane of the
$m({\mu^{+}\mu^{-}})$ and $\Delta m$ variables, as well as the signal regions
in each variable. Figure 8 displays the result of the calibration of the BDT
output described in the text. Figure 9 shows the test on the particle
identification variable mentioned in the main body. Finally, Figure 10 shows
the value of the $\rm{CL_{s}}$ estimator used to compute the upper limit on
the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ branching fraction, as a
function of the branching fraction itself.

Figure 3: Distributions of (left) $m({\mu^{+}\mu^{-}})$ and (right) $\Delta m$
for the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ candidates for Run 1 data
in, top to bottom, the three BDT intervals. The distributions are superimposed
with the fit to data. Each of the two distributions is in the signal region of
the other variable, see text for details. Untagged and tagged decays are
included in a single component for signal and
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ background.

Figure 4: Distributions of (left) $m({\mu^{+}\mu^{-}})$ and (right) $\Delta m$
for the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ candidates for Run 2 data
in, top to bottom, the three BDT bins. The distributions are superimposed with
the fit to data. Each of the two distributions is in the signal region of the
other variable, see text for details. Untagged and tagged decays are included
in a single component for signal and
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ background.

Figure 5: Distributions of (left) $m({\mu^{+}\mu^{-}})$ and (right) $\Delta m$
for the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ candidates for Run 1 data
in, top to bottom, the three BDT bins. The distributions are superimposed with
the fit to data. Untagged and tagged decays are included in a single component
for signal and ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ background.
Unlike the correspondent figures in the main body of the Letter, here all
events are shown.

Figure 6: Distributions of (left) $m({\mu^{+}\mu^{-}})$ and (right) $\Delta m$
for the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ candidates for Run 2 data
in, top to bottom, the three BDT bins. The distributions are superimposed with
the fit to data. Untagged and tagged decays are included in a single component
for signal and ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ background.
Unlike the correspondent figures in the main body of the Letter, here all
events are shown.

Figure 7: Two-dimensional distribution of $m({{D}^{0}})$ versus $\Delta m$ for
the most sensitive BDT interval in the (left) Run 1 and (right) Run 2 data.
Red lines represent the signal region of each variable as defined in the text.

Figure 8: Distribution of the BDT output. The
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ candidates in data and
simulation are shown and used to weight the simulated
${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays distribution used for the
search.  Figure 9: Efficiency for two pions to pass a requirement on the
$ProbNNmu$ variable (after additional requirements on its identification as
muon) in Run 2, as measured from
${{D}^{+}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{+}}$ and
${{D}^{+}_{s}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{+}}$ decays in data
and ${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays in simulation. Same
sign pions are used in data to avoid contamination for hadronic resonances
decaying to two real muons. Figure 10: Value of $\rm{CL_{s}}$ as a function of
the ${{D}^{0}}\\!\rightarrow{\mu^{+}}{\mu^{-}}$ branching fraction.