Search for the Flavor-Changing Neutral Current Decay $D^0 \to \mu^+\mu^-$ with the HERA-B Detector

We report on a search for the flavor-changing neutral current decay $D^0 \to \mu^+\mu^-$ using $50 \times 10^6$ events recorded with a dimuon trigger in interactions of 920 GeV protons with nuclei by the HERA-B experiment. We find no evidence for such decays and set a 90% confidence level upper limit on the branching fraction $Br(D^0 \to \mu^+\mu^-)<2.0 \times 10^{-6}$.


Introduction
The decay D 0 → µ + µ − 13 is sensitive to flavor-changing neutral currents (FCNC), which, due to the GIM mechanism [1], are forbidden at lowest order and strongly suppressed at second order in the Standard Model (SM). In the SM, the expected contribution of short distance processes to the branching fraction for this decay is of the order of 10 −19 [2,3], well below the sensitivity of current experiments. Long distance effects may enhance the branching fraction to roughly 10 −13 [3,4], still undetectable by foreseeable experiments. However, several extensions of the SM including the Minimal Supersymmetric Standard Model (MSSM), models with multiple Higgs doublets, with horizontal gauge bosons, and with extra fermions, predict an enhancement of the branching fraction by several orders of magnitude. A comprehensive study of this issue has been presented in Burdman et al. [4] (see also references therein). According to this report, an MSSM variant with R-parity violation predicts the highest estimated branching fraction, 3.5 × 10 −6 . Thus, the value of the branching fraction of the D 0 → µ + µ − decay mode is exquisitely sensitive to physics beyond the SM. Recently, the CDF Collaboration has published a new upper limit of 2.5 × 10 −6 at the 90% confidence level [5].
In this letter we report on a search for the D 0 → µ + µ − decay using 50 × 10 6 events recorded with a dimuon trigger in interactions of 920 GeV protons with nuclei in the experiment HERA-B [6]. of reconstructed J/ψ → µ + µ − events [7] since possible biases arising from the dimuon trigger and muon identification largely cancel. Thus, for proton interactions on a single target of atomic weight A, the branching fraction limit can be determined using the following inequality: where • n cl is the upper limit on the number of D 0 → µ + µ − decays; • N J/ψ is the number of observed J/ψ → µ + µ − events; • a D 0 and a J/ψ are the acceptances for D 0 → µ + µ − and J/ψ → µ + µ − after event quality and particle identification cuts are applied (i.e., cuts applied to both channels); • ǫ D 0 is the efficiency for D 0 → µ + µ − of cuts designed to select secondary vertices (i.e., cuts not applied to J/ψ → µ + µ − ); • σ pA D 0 and σ pA J/ψ are the production cross sections per target nucleus for D 0 and J/ψ; The terms a D 0 , a J/ψ and ǫ D 0 are evaluated with a complete Monte Carlo simulation. Equation (1) also requires knowledge of the relative production cross sections of D 0 and J/ψ. Large errors on the available measurements of the D 0 cross section dominate the systematic error.

The Detector and Trigger
The HERA-B fixed-target spectrometer [6] operates at the 920 GeV proton beam of the HERA storage ring at DESY and features a vertex detector and extensive tracking and particle identification systems. It has a large geometrical coverage from 15 mrad to 220 mrad in the bending plane and 15 mrad to 160 mrad in the vertical plane. Figure 1 shows a plan view of the detector in the configuration of the 2002-2003 data run.
The target system [9] consists of two stations of four wires each. The wires are positioned above, below, and on either side of the beam and are made from various materials including carbon, titanium, and tungsten. The stations are separated by 40 mm along the beam direction. The wires can be individually moved into the halo of the HERA proton beam and the interaction rate for each inserted wire can be adjusted independently.
The Vertex Detector System (VDS) [10] is a forward microvertex detector integrated into the HERA proton ring. It provides a precise measurement of primary and secondary vertices. The VDS consists of 7 stations (4 stereo views) of double-sided silicon strip detectors (50 × 70 mm, 50 µm pitch) integrated into a Roman pot system inside a vacuum vessel and operated as near as 10 mm from the beam. An additional station is mounted immediately downstream of the 3 mm thick aluminum window of the vacuum vessel.
The first station of the main tracker is placed upstream of the 2.13 T-m spectrometer dipole magnet. The remaining 6 tracking stations extend from the downstream end of the magnet to the electromagnetic calorimeter (ECAL) at 13 m downstream of the target. Each tracking station is divided into inner and outer detectors. The large area Outer Tracker (OTR) [11] consists of 95,000 channels of honey-comb drift cells and covers the region starting from 200 mm from the beam center. The region starting from the beam pipe and extending up to the start of the OTR acceptance is covered by micro-strip gas chambers with GEM foils (Inner Tracker or ITR [12]).
Particle identification is performed by a Ring Imaging Cherenkov detector (RICH) [13], the ECAL [14] and a muon detector (MUON) [15]. The MUON detector is segmented into four superlayers. Iron and concrete shielding extends from just behind the ECAL to the penultimate MUON superlayer, except for gaps for the superlayers themselves. The first two superlayers consist of three layers of tube chambers with different stereo angles. The last two superlayers each consist of one layer of pad chambers.
For the sample considered here, triggers are initiated by coincidences of pad chamber hits ("pretriggers") from the third and last MUON superlayers [16]. Starting from the pretrigger coordinates, the First Level Trigger (FLT) [17] searches for tracks of charged particles in the MUON tracking layers and four of the main tracker stations (OTR only). The momentum of a candidate track is inferred by assuming the track originates at the target wires and calculating the bending angle within the magnetic field. Triggered events are required to have at least one FLT track and at least one additional pretrigger candidate.
The Second Level Trigger (SLT) [18] processor farm receives the track and pretrigger candidates from the FLT. Starting from regions of interest defined by the pretriggers, it searches for tracks in the main tracker, extrapolates those found through the magnet and attempts to follow them through the VDS. Events with at least two fully reconstructed tracks which form a good vertex are triggered. Events passing the SLT are assembled and transferred to a fourth level farm (no third level algorithm is applied), where a fraction of the events is fully reconstructed online to provide monitoring. No event selection is applied at the fourth level. The final archiving rate is about 100 Hz.

Monte Carlo Simulation
A Monte Carlo (MC) simulation is used to determine the D 0 → µ + µ − and J/ψ → µ + µ − acceptance and efficiency terms in equation (1) and to estimate the background due to D 0 → K − π + and D 0 → π + π − decays. The Monte Carlo events for pA → D 0 +X are generated in two steps. First, a cc pair is generated with PYTHIA 5.7 and hadronized with JETSET 7.4 [19] such that a D 0 is always produced. The fractional longitudinal momentum (x F ) distribution of the D 0 is forced into the form (1 − |x F |) n and its transverse momentum (p T ) distribution is forced into the form exp(−bp 2 T ) by an appropriate selection of events. For n and b, we use averages of measurements by E653 [20] and E743 [21]: n = 7.7 ± 1.4, b = 0.83 ± 0.07 (GeV/c) −2 . After the generation of the D 0 is completed, the remaining energy is given as input to FRITIOF 7.02 [22], which generates the underlying event taking into account further interactions inside the nucleus. The generated event conserves momentum.
A similar two-step procedure is used to generate the J/ψ → µ + µ − sample, differing only in the method for achieving agreement of J/ψ kinematic distributions with previous measurements. In this case, the generated events are reweighted according to the parameterizations of the J/ψ differential cross sections in proton-gold collisions at 800 GeV by E789 [23]: dσ/dp 2 The detector response is simulated with the GEANT 3.21 package [24]. Realistic detector efficiencies, readout noise and dead channels are taken into account. The simulated events are processed by the same trigger and reconstruction codes as the data. We have checked that the mis-identification probabilities of pions from K 0 S decay, kaons from φ(1020) decay and protons from Λ decay estimated by the Monte Carlo agree with those measured in the data.

Introduction
During data-taking, the interaction rate was maintained at approximately 5 MHz, resulting in an average of approximately 0.6 interactions per filled HERA bunch. A total of 50 × 10 6 triggers were recorded, with 38% from runs with one or two carbon targets, 55% from runs with carbon and tungsten targets operated simultaneously, 5% from runs with one titanium wire and 2% from runs with a single tungsten wire. After eliminating runs with problematic detector performance, poor beam conditions or non-standard trigger conditions, 47 × 10 6 events remain.
The events are reconstructed with the standard HERA-B analysis package [25,26]. The search for J/ψ and D 0 candidates is performed by analyzing the mass spectrum for unlike-sign dimuon candidates. D 0 candidates are further selected by requiring a detached vertex. The number of background events in the signal region is calculated from sidebands. The upper limit is obtained by comparing the observed number of events inside the signal region with the background estimate and normalizing to the number of observed J/ψ's, after appropriate corrections are applied.
The mass spectrum for unlike-sign dimuon candidates after minimal cuts is shown in Fig. 2. Peaks at the ω/ρ, φ, J/ψ and ψ(2S) masses are clearly visible. The fall-off towards low mass is caused by an implicit p T cut of about 0.7 GeV/c imposed by the SLT when defining regions of interest for subsequent tracking. A fit in the mass interval 2.4-3.5 GeV/c 2 to an exponential background plus a Gaussian with a radiative tail yields 147710 ± 520 J/ψ decays in a two standard deviation window centered on the mean value from the Gaussian fit: 3.095 GeV/c 2 . The fitted width is 44 MeV/c 2 . The mean value is 1.9 MeV/c 2 below the PDG value [8]. From these numbers we estimate the mass resolution in the D 0 region to be about 25 MeV/c 2 and that the reconstructed mass of a D 0 → µ + µ − signal would be within 1.1 MeV/c 2 of the PDG [8] value.

Event Selection
The main background for the D 0 → µ + µ − channel results from muon pairs from independent π ± or K ± decays which appear to form a secondary vertex displaced from the primary vertex. The cuts are designed to minimize this background while maintaining high efficiency for D 0 decays.
The data sample is first reduced by application of relatively loose cuts. In this first pass, a general vertex-finding algorithm is applied and events with at least one reconstructed primary vertex and at most one reconstructed primary vertex per target wire are selected. Muon candidate tracks are then selected by requiring that the muon probability, P µ , derived from the MUON hit information, be greater than 0.01, that the χ 2 per degree of freedom of the track fit (χ 2 tr /d.o.f.) be less than 20, and that the kaon identification likelihood probability from the RICH reconstruction be less than 0.4. The two muon tracks are then excluded from the primary vertex, and primary and dimuon vertices are fitted. The χ 2 probability of the primary vertex is required to be greater than 0.01 and that of the dimuon vertex greater than 0.2. The number of tracks per primary vertex is required to be less than 50. Approximately 0.93 million dimuons with mass in the D 0 region and 69000 J/ψ decays (estimated here and below using the fit described in Sec. 4.1) survive the initial cuts.
In the following, we call the muon pair with a fitted secondary vertex a "dimuon pseudoparticle". We consider two regions of µ + µ − invariant mass: 2.7−4.0 GeV/c 2 , which we refer to as the "J/ψ region", and 1.59−2.15 GeV/c 2 , which we refer to as the "D 0 region".
Further selection cuts are divided into two parts: • Common cuts applied both for D 0 and J/ψ regions. These cuts are mainly quality cuts, the efficiencies of which are nearly identical for muons from possible D 0 → µ + µ − and for muons from J/ψ → µ + µ − . Small differences in the efficiencies coming from the different momentum distributions of muons from D 0 decay and J/ψ decay are evaluated from the data and Monte Carlo. • Cuts applied only for the D 0 region (lifetime cuts). These cuts are intended to select a well-defined detached secondary vertex associated with the selected primary vertex.
The following common selection criteria are applied for D 0 and for J/ψ regions: • a cut on total track multiplicity (N tr ) to suppress multi-event pile-up which is enhanced by the two-muon requirement; • a cut on χ 2 tr /d.o.f. for each muon to suppress ghosts and π/K decays in flight; • a cut on P µ to reduce fake dimuon events; • a cut on the transverse momentum of each muon (p µ T ) to suppress muons from pion and kaon decays.
To optimize the common cuts, we employ a blind analysis technique: all dimuons from the D 0 signal region are masked and the cuts are chosen to maximize the quantity N J/ψ / √ B D 0 , where N J/ψ is the number of J/ψ candidates above background found in the µ + µ − invariant mass spectrum in a two standard deviation window around the J/ψ position and B D 0 is the expected background in the D 0 signal region (1.815 − 1.915 GeV/c 2 ), estimated from D 0 sidebands (1.59 − 1.79 GeV/c 2 and 1.94 − 2.14 GeV/c 2 ). The resulting cuts are: N tr < 45, χ 2 tr /d.o.f. < 7.5, P µ > 0.7, p µ T > 0.7 GeV/c.
After all common cuts are applied, about 238000 events in the D 0 region and 46000 events in the J/ψ peak remain. The dimuon mass distributions for D 0 and J/ψ regions with the above cuts are shown in Fig. 3.

D 0 Selection
To isolate possible D 0 mesons, cuts are applied to three quantities: the separation between primary and secondary vertices, the proper decay length, and the impact parameter of the dimuon pseudoparticles to the primary vertex. At HERA-B energies, nearly all D 0 mesons originate at the target wires, i.e., the fraction arising from B decays is negligible (< 0.1%).
Most D 0 mesons decay within a few millimeters of the primary vertex in the laboratory frame. This distance is comparable to the precision of the secondary vertex measurement in the longitudinal direction. To ensure that the secondary vertex is well separated from the associated primary vertex, we compare this distance with uncertainties of primary and secondary vertices and cut on the separation value ∆z = (z sec − z pr )/ σ 2 zsec + σ 2 zpr , where z pr and z sec are the z-coordinate (along the beam direction) of primary and secondary vertices, respectively, and σ zsec , σ zpr are their calculated errors. The average resolutions are σ zsec = 500 µm and σ zpr = 420 µm.
The proper decay length is given by cτ = mc · L/p, where m is the dimuon invariant mass, L is the decay length in the laboratory frame, and p is the reconstructed momentum of the dimuon pseudoparticle.
The impact parameter (I v ) is defined as the distance between the primary vertex and the point of intersection of the dimuon pseudoparticle flight direction with the xy plane at the z position of the primary vertex. The impact parameter resolution is typically in the range of 30 µm to 60 µm.
To optimize the cuts, we apply a blind analysis technique similar to that used to optimize the common cuts. Since these three quantities are correlated, a three-dimensional optimization is performed wherein the ratio is the number of reconstructed Monte Carlo events in the D 0 peak. The "experimental sensitivity",S, is the average 90% confidence level upper limit on the number of signal events obtained for an ensemble of experiments with the expected background estimated from the D 0 sidebands, assuming no signal from D 0 → µ + µ − is present (from Table XII of Ref. [27]).
This optimization procedure produces the following cuts: ∆z > 7.0, cτ > 250 µm and I v < 110 µm. The optimized experimental sensitivity is S opt = 5.53 events.

Results
After applying all cuts, 31 events remain in the D 0 mass region as shown in Fig. 4. The sidebands are indicated by dashed lines and signal region is indicated by dotted lines. The signal region contains 3 events.
The number of background events in the D 0 mass region from D 0 → K − π + decays is estimated from Monte Carlo to be 1.8 ± 1.0. All other D 0 decay modes give a negligible contribution. The background shape is therefore not significantly influenced by charm decay and we estimate the background using the shape of the mass plot before D 0 selection cuts are applied (see Fig. 3a). After normalizing to the number of events in the sidebands, the expected background in the signal region is 6.0 ± 1.2. We have checked that the background shape does not significantly change relative to the uncut distribution when one of the three selection cuts is removed (three tests). The maximum difference between background estimates derived from the uncut distribution and the three partially cut distributions is 4%. We note that a simple linear interpolation between sidebands also predicts 6.0 ± 1.2 background events in the signal region.
The upper limit calculation for the D 0 → µ + µ − decay mode requires knowledge of the D 0 to J/ψ production cross section ratio in proton-nucleus interactions. Three experiments have published results on the D 0 production cross section at 800 GeV. The x F coverages of experiments E653 [20] and E743 [21] are similar to that of HERA-B: −0.2 < x F < 0.1. The x F coverage of experiment E789 [28] is relatively small and the E653 and E743 parameterizations have been used to extrapolate to the full x F range. For this reason we do not include it in the calculation of the D 0 cross section. From the remaining two measurements of 22 +9 −7 ± 5 µb/nucleon and 38 ± 3 ± 13 µb/nucleon we obtain σ pN D 0 = 27.3 ± 7.7 µb/nucleon.
The prompt J/ψ production cross section per nucleon, σ pN J/ψ , has been measured by two fixed target experiments [23,29] at 800 GeV. Both measurements were performed with nuclear targets and extrapolated to atomic weight (A) of 1, assuming σ pA J/ψ = σ pN J/ψ · A α J/ψ . After adjusting these measurements using the most recent measurement of α J/ψ = 0.955 ± 0.005 [30], and averaging them we obtain a prompt J/ψ production cross section of σ pN J/ψ = 333 ± 6 ± 26 nb/nucleon at 800 GeV. We assume that the ratio of D 0 and J/ψ cross sections does not change significantly between 800 GeV and 920 GeV.
The ratio of acceptances, not including cuts applied only for selecting D 0 → µ + µ − , is a D 0 /a J/ψ = 0.287 ± 0.028. The error is evaluated by studying the efficiencies of the applied cuts when detector characteristics (channel efficiencies, maskings, resolutions) are varied within ranges that reflect their uncertainties. The quadratic sum of contributions from all significant sources, including triggering, target geometry, muon and kaon identification and Monte Carlo statistics is 9.8%. The efficiency of the additional secondary vertex cuts applied in the D 0 region is ǫ D 0 = (6.83 ± 1.08) × 10 −2 .
Equation (1) can be rewritten as: is the detector sensitivity after summation over targets.
The values and associated errors for all terms in equation (3) are given in Table 1. Using these numbers, F sens = 1.57 × 10 6 . After combining statistical and systematic errors quadratically the total relative error on F sens from all contributing terms in equation (3) is 37%.  [27,32,33]. We choose an approach similar to that first advocated by Feldman and Cousins [27]. This method has been since extended by Conrad et al. [34] to incorporate uncertainties in detector sensitivity and the background estimate based on an approach described by Cousins and Highland [35]. A further refinement of the Conrad et al. method by Hill [36] results in more appropriate behavior of the upper limit when the observed number of events is less than the estimated background, as is the case for the present measurement. We have adopted this method but note that Table 2 contains all of the numbers needed to calculate an upper limit using any of the methods in the papers cited above. We assume that the probability density functions of F sens and background estimates are Gaussian-distributed.
Using equation (2) and taking systematic errors and background fluctuations into account with the Hill prescription [37], the upper limit at 90% confidence level on B(D 0 → µ + µ − ) is 2.0 × 10 −6 . Table 2 Summary of parameters entering into the upper limit calculation Sensitivity factor, F sens (1.57 ± 0.58) × 10 6 Number of events in the signal region 3 Expected background events 6.0 ± 1.2

Conclusions
We have investigated the dimuon mass spectrum in a search for the flavorchanging neutral current decay D 0 → µ + µ − . Using the values of D 0 and J/ψ production cross sections published in the literature we have set an upper limit on the branching fraction of B(D 0 → µ + µ − ) < 2.0 × 10 −6 at 90% confidence level. For comparison, the lowest previously published limit on B(D 0 → µ + µ − ) is 2.5 × 10 −6 at 90% confidence level [5].
According to the recent calculations by Burdman et al. [4], in the case of strong R-parity violation in the MSSM the D 0 → µ + µ − branching fraction is at the level of 3.5 × 10 −6 . Thus, our result disfavors this scenario and constrains the size of R-parity violating couplings which enter into the predicted branching fraction.