Search for charged lepton flavor violating decays of $\Upsilon(1S)$

We present a search for the charged lepton-flavor-violating decays $\Upsilon(1S)\rightarrow\ell\ell^\prime$ and radiative charged lepton-flavour-violating decays $\Upsilon(1S)\rightarrow\gamma\ell\ell^\prime$ [$\ell,\ell^\prime = e, \mu, \tau$] using the 158 million $\Upsilon(2S)$ sample collected by the Belle detector at the KEKB collider. This search uses $\Upsilon(1S)$ mesons produced in $\pi^+\pi^-\Upsilon(1S)$ transitions. We do not find any significant signal, so we provide upper limits on the branching fractions at the 90% confidence level.


Introduction
Observations of neutrino oscillations [1] imply that the accidental lepton family symmetry in the standard model (SM) Lagrangian is broken. The minimal extension of the SM that can explain neutrino oscillations requires the presence of a right-handed neutrino. In such a framework, the conservation of individual lepton flavor is violated, and charged lepton-flavor-violating (CLFV) transitions can occur, mediated by W ± bosons and massive neutrinos. However, the existence of such CLFV transitions would imply a minimal value of B (µ ± → e ± γ) ∼ 10 −54 [2,3]. Several new physics models inspired by grand unified theories, such as supersymmetry and those predicting leptoquarks, typically enhance decay rates of CLFV transitions [4,5].
The effective Lagrangian of new physics (NP) models can be expressed as the sum of a dipole term, four-fermionic interactions, and a gluonic interaction part. The Wilson coefficients of the NP operators can be determined via fits to measurements of phenomena those involve CLFV interactions [6]. Several classes of operator, such as vector, axialvector, and tensor operators involved in four-fermionic interactions, allow CLFV transitions. Precise measurement of two-body vector meson CLFV decays allows one to probe the vector and tensor operators effectively.
Radiative lepton-flavor-violating (RLFV) transitions allow one to probe the operators which are not easily accessible in the two-body decays [6]. Using three-body vector meson RLFV decays, one can put constraints on the corresponding Wilson coefficients of axialvector, scalar, and pseudoscalar operators. Thus, RLFV studies of Υ(nS) [n = 1, 2, 3] provide complementary access to the NP parameters. Currently, there are no existing results available for the Υ(nS) → γ ± ∓ decays. We perform the first search for RLFV in Υ(1S) → γ ± ∓ decays using the Υ(2S) data sample.

Belle experiment
The world's largest Υ(2S) sample, corresponding to 158 million Υ(2S) events, was collected with the Belle detector at the KEKB asymmetric-energy e + e − collider [9] operating at a center-of-mass of energy ( √ s) of 10.02 GeV. We study the e + e − → qq (q = u, d, s, c) background using the 80 fb −1 data sample collected at 10.52 GeV.
The Belle detector is a large-solid-angle spectrometer, which includes a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter (ECL) comprised of 8736 CsI(Tl) crystals located inside a superconducting solenoid coil that provides a 1.5T magnetic field. An iron flux return located outside the coil is instrumented to detect K 0 L mesons and identify muons (KLM). The detector is described in detail elsewhere [10].

Event selection
We use the EVTGEN package [11], with QED final-state radiation simulated by PHO-TOS [12], for the generation of Monte Carlo (MC) simulation events. We generate the signal events for two-body CLFV modes using both the vector to two leptons decay model, VLL and the phase space decay model, PHSP. The reconstruction efficiencies for the MC signal events generated with the VLL model are smaller (by approximately 8%) than for the PHSP model. We will use the MC signal events generated with the VLL decay model for twobody CLFV decays to quote the most conservative upper limits. We use the PHSP model to generate the signal events for RLFV modes. We are using TAUOLA [13] or PYTHIA [14] for generating signal events for the subsequent decays of τ leptons. A GEANT3-based [15] MC simulation is used to model the response of the detector. Thus, dedicated MC samples are generated for different signal modes to determine signal efficiencies and define selection criteria. Background studies and the optimization of those criteria are performed using an MC simulated sample of Υ(2S) events with a size corresponding to the integrated luminosity. Dominant backgrounds arise from Υ(2S) → π + π − Υ(1S) decays with Υ(1S) → ± ∓ [ = e, µ, τ ]. The MC samples for these decays, corresponding to about 20 times of data sample sizes, are used to study backgrounds. For two-body CLFV searches, Υ(1S) candidates are reconstructed in the e ± µ ∓ , µ ± τ ∓ , and e ± τ ∓ final states. We recon- Also, we reconstruct Υ(1S) → e + e − and Υ(1S) → µ + µ − decays, which are used to validate and calibrate the analysis. To validate the recoil Υ(1S) sample along with muon and electron identifications, we measure the branching fractions of Υ(1S) → e + e − , and Υ(1S) → µ + µ − decays. As taus are mostly reconstructed in the leptonic decays, validating lepton identification with the high momentum leptons from Υ(1S) is also important for Υ(1S) → ± τ ∓ and Υ(1S) → γ ± τ ∓ decays. This analysis follows a blind analysis procedure.
Charged tracks are required to originate from the vicinity of the interaction point (IP). The distance of the closest approach to the IP is required to be within 3.5 cm along the beam direction and within 1.5 cm in the transverse plane. The combined information from the CDC, TOF, and ACC is used to identify charged pions based on the pion likelihood ratio, L π = P π /(P π + P K ), where P π and P K are likelihood values for the pion and kaon hypotheses, respectively [16]. Pions are required to have L π > 0.6, which has an identification efficiency of 94%. Muon candidates are identified using a likelihood ratio, L µ , which is based on the difference between the range of the track calculated from the particle momentum and that measured in the KLM. This ratio includes the value of χ 2 formed from the KLM hit locations with respect to the extrapolated track. The muon identification efficiency for the applied selection, L µ > 0.95, is 89%, with a pion misidentification probability of 1.4% [17]. Identification of electrons uses an analogous likelihood ratio, L e , based on specific ionization from the CDC, the ratio of the energy deposited in the ECL to the momentum measured by the CDC and SVD combined, the shower shape in the ECL, hit information from the ACC, and matching between the position of the charged track and the ECL cluster. The electron identification efficiency for the applied selection, L e > 0.6, is 95%, with a pion misidentification probability 0.3% [18]. To recover the bremsstrahlung energy loss for electrons and positrons, we include the energy from the photon(s) within 50 mrad of each of the e ± tracks, which improves the efficiency for true signal events by 2.7%. For Υ(1S) → ± τ ∓ decay, most of the background comes from Υ(1S) → τ + τ − and Υ(1S) → ± ∓ decays. To suppress events coming from Υ(2S) decays, other than Υ(2S) → π + π − Υ(1S), we define the recoil mass of two pions as: Where E total , E ππ , and | − → p ππ | are the total energy of the colliding e + e − beams, the energy of the two pions from Υ(2S) and the magnitude of the 3-momentum of the pion pair, respectively in the center-of-mass (CM) frame. The distribution of M recoil ππ is shown in Fig. 1. The M recoil ππ distribution peaks at the Υ(1S) mass for signal events, while it is flat for the combinatorial background. We consider the events within 9.450 < M recoil ππ < 9.466 GeV/c 2 , corresponding to a ±3σ region around the nominal Υ(1S) mass [19]. To suppress the background from e + e − γ and µ + µ − γ, we remove the events with the cosine of the angle between the two pions in the Υ(2S) rest frame (cos θ ππ ) greater than 0.5. We define the visible tau momentum (p τ vis ) as the sum of the momentum carried by the daughter charged track(s) of τ in the lab frame. We select the τ candidates with p τ vis > 0.3 GeV/c. Furthermore, τ − → π − π + π − ν τ is reconstructed with the invariant mass of the three-pion lower than 1.8 GeV/c 2 and energy in the lab frame greater than 2.6 GeV. These τ selections are wide enough to account for the missing momentum from neutrinos. Also, we fit the three-pion vertex for the τ − → π − π + π − ν τ decay and events with fitted χ 2 < 15 are selected to reduce combinatorial backgrounds. We count the number of tracks identified as muons or electrons with energy in the lab frame greater than 1 GeV as prompt muons (N µ ) and prompt electrons (N e ), respectively. For Υ(1S) → µ ± τ ∓ decays, in order to reject the background coming from the Υ(1S) → µ + µ − decay, we select the events with N µ = 1 and Similarly, for Υ(1S) → e ± τ ∓ decays, we select the events with N e = 1 and N µ ≤ 1 (N µ = 0) for We suppress a large number of Υ(1S) → µ + µ − and Υ(1S) → e + e − backgrounds by the selections of prompt leptons.
For the Υ(1S) → e ± µ ∓ study, the distribution of lepton pair invariant mass (M eµ ) for π + π − recoil sample in Υ(2S) data is shown in Fig. 2. We select the events with M eµ within 9.09 to 9.65 GeV/c 2 by selecting a ±3σ region around the mean position and the Υ(1S) momentum in the lab frame (| − → p eµ |) less than 4.4 GeV/c to reduce the e + e − → qq events.
For RLFV decays, there is an extra photon in the final state. Therefore, we include all the selections which are used for the pion and lepton in the corresponding two-body CLFV decay previously discussed.
decay will be replaced by M γeµ and | − → p γeµ |, respectively. In addition to the above selections to the corresponding non-radiative mode, we select photons with energy in the lab frame greater than 200 MeV to remove soft photons and beam backgrounds. The photons used in the bremsstrahlung recovery are not considered in reconstructing the radiative candidates. Inside the signal search window, we find 3%, 8%, and 7% multiple Υ(2S) candidates for Υ(1S) → γe ± µ ∓ , Υ(1S) → γµ ± τ ∓ , and Υ(1S) → γe ± τ ∓ decays, respectively. Multiplicity due to misreconstructed charged particles is handled using a procedure similar to two-body CLFV decays. Multiplicity occurring from the multiple photon candidates is removed by selecting the event randomly. Best candidate selection efficiencies for Υ(1S) → γe ± µ ∓ , Υ(1S) → γµ ± τ ∓ , and Υ(1S) → γe ± τ ∓ are 93%, 84%, and 87%, respectively.   To study the calibration modes, we select events with lepton pair invariant mass (M ) within 9.09 to 9.65 GeV/c 2 and momentum of the reconstructed Υ(1S) in the lab frame (| − → p |) less than 4.4 GeV/c. To extract the signal for Υ(1S) → ± ∓ decays, we perform an unbinned maximum likelihood (UML) fits to ∆M = M ππ −M , where = e, µ. The signal probability density function (PDF) used is a sum of two Gaussians sharing a common mean. Backgrounds from all the sources are flat in the ∆M window and small compared to the signal yields. We fit the background with a first-order Chebyshev polynomial. To account for any resolution difference between data and MC, the mean and the width parameter of the primary Gaussian (σ 1 ) are floated in the fit, and the width of the secondary Gaussian is set to σ 2 = k × σ 1 , with the factor k fixed from MC.
Expected signal efficiencies for ee and µµ are estimated to be 28.3% and 35.6%, respectively. Fig. 3 shows the fits to Υ(2S) data. The signal yields obtained for the ee and µµ final states are 191353 ± 467 and 246255 ± 504 events, respectively. The data-MC differences of widths for the e ± e ∓ and µ ± µ ∓ final states are estimated to be 12% and 16%, respectively. One can calculate the branching fractions using the following relation: where, N Υ(2S) , N sig and are the number of Υ(2S) produced in e + e − collision, signal yield in data and the effective signal efficiency (after implementing all the systematic corrections) respectively. Using Equation (4.1), the calculated branching fractions including only statistical uncertainties for Υ(1S) → e + e − and Υ(1S) → µ + µ − are (2.40 ± 0.01) × 10 −2 and (2.46 ± 0.01) × 10 −2 , respectively. These are consistent with the world average values [19].
These results are discussed further in Section 8 after including the systematic uncertainty.  . ∆M fit to Υ(2S) data for Υ(1S) → e ± µ ∓ decay. The fitted signal PDF is represented by the filled red region, the dashed cyan line represents the flat background and the dotted magenta curve is the peaking background from lepton misidentification. The solid blue curve represents the overall fit to data. The long-dashed red curve represents the signal PDF corresponding to 5 hypothetical signal events.

Signal extraction for two-body CLFV decays
We extract the signal yield from a UML fit to the ∆M variable. ∆M should peak at the nominal mass difference between Υ(2S) and Υ(1S), approximately 560 MeV/c 2 [19]. A sum of two Gaussians sharing a common mean has been used as the signal PDF. To estimate the peaking background, the shape of the peaking background is considered to have the same shape as the signal PDF. The qq backgrounds are flat, and they are modeled with a first-order Chebyshev polynomial. The width of the signal PDF in the data is fixed at the MC width, corrected by the average of the data-MC difference for the µµ and ee samples.
A few Υ(1S) → µ + µ − (Υ(1S) → e + e − ) events mimic our signal when a µ ± (e ± ) is misidentified as an e ± (µ ± ). The amount of Υ(1S) → τ + τ − background is estimated to be negligible. The number of Υ(1S) → µ + µ − background events is estimated to be 3.5 ± 0.4 using a large MC sample. Such backgrounds are difficult to remove completely. To estimate the background from muon to electron misidentification in the data, we derive a correction factor for electron misidentification efficiency of muons using an e + e − → µ + µ − sample collected at √ s = 10.52 GeV with a tag-and-probe method. The data to MC correction factor for electron misidentification efficiency is estimated to be 2.5 ± 0.5, which leads to an estimation of this peaking background yield of 8.8 ± 2.0. The background from electron to muon misidentification is expected to be consistent with zero (0.1 ± 0.1 events). Fitted distribution of Υ(2S) data is shown in Fig. 4. To consider the peaking background, we include a fixed PDF of 8.8 events in the data fit (dotted magenta line), and uncertainty (2.0) will be added in systematic uncertainty. The estimated signal efficiency for the Υ(1S) → e ± µ ∓ mode is 32.5%. We finally obtain a yield of −1.3 ± 3.7 signal events for the Υ(1S) → e ± µ ∓ decay.

Υ(1S) → ± τ ∓ decay
For µτ and eτ decays of Υ(1S), we extract the signal from an UML fit to the recoil mass of ππ (M recoil ππ ), where = µ, e. M recoil ππ can be defined by replacing ππ with ππ in Eq. (3.1). As M recoil ππ is calculated from all the particles from the Υ(2S) except the τ , M recoil ππ should peak at the nominal τ mass (around 1.78 GeV/c 2 ) [19]. Signal events of µτ and eτ decays are modeled with a sum of one Gaussian and one bifurcated Gaussian, sharing a common mean. We obtain a difference of 7% (27%) for the resolution between data and MC for M recoil ππµ (M recoil ππe ) using the data-MC difference for the Υ(1S) → µ + µ − (Υ(1S) → e + e − ) mode. The width of the Υ(1S) → µ ± τ ∓ (Υ(1S) → e ± τ ∓ ) signal PDF in the data is fixed from the MC width corrected by the data-MC difference for the M recoil ππµ (M recoil ππe ) parameter. For Υ(1S) → µ ± τ ∓ decays, the main backgrounds come from Υ(1S) → τ + τ − and Υ(1S) → µ + µ − decays. For the τ τ background, a charged lepton or a pion from one of the tau decays is used as the signal muon. In the M recoil ππµ distribution, missing neutrino energy from the misidentified τ shifts such events away from the actual τ mass. Thus, the τ τ background increases exponentially, starting near the nominal τ mass value. We model the τ τ background using the following exponential threshold PDF starting near M th (in GeV/c 2 ), where A and B are the two slope parameters of the τ τ background PDF. We try other fitting models and find the current model describes the background the best. To get the proper shape of the τ τ background PDF, we use a large Υ(1S) → τ + τ − MC sample. In the data fit, M th and B are fixed from the MC background, and A is allowed to float. Fig. 5 shows the fitted distributions of data for Υ(1S) → ± τ ∓ decays. The background from Υ(1S) → µ + µ − is obtained from a large MC sample. The expected number of µµ backgrounds is less than the number of τ τ backgrounds, but it widely populates around the signal region: it peaks at the lower mass value (< 1 GeV/c 2 ) and has a broad tail. The PDF is presented by a sum of one bifurcated Gaussian and one threshold function starting from 0 GeV/c 2 using a large Υ(1S) → µ + µ − sample. To fit the data, we float the yield of the µµ background fixing the shape of the Gaussian of µµ background from MC corrected by the data-MC difference for the M recoil ππµ parameter. For Υ(1S) → e ± τ ∓ decays, potential backgrounds arise from Υ(1S) → τ + τ − and Υ(1S) → e + e − decays. These backgrounds are handled using a procedure similar to that used for the τ τ and µµ backgrounds to Υ(1S) → µ ± τ ∓ decays.
The expected peaking backgrounds in the Υ(2S) data for Υ(1S) → µ ± τ ∓ and Υ(1S) → e ± τ ∓ decays are estimated from MC to be 0.7 ± 4.1 and 4.6 ± 5.4, respectively. As no significant peaking background is found in the Υ(2S) decay MC sample, we do not include a peaking background component in the fit. Considering both the τ reconstruction modes, the effective signal efficiency for Υ(1S) → µ ± τ ∓ (Υ(1S) → e ± τ ∓ ) decay is 8.8% (7.1%). In Υ(2S) data, we find the yield of Υ(1S) → µ ± τ ∓ and Υ(1S) → e ± τ ∓ signals to be −1.5±4.3 and −3.5 ± 2.7, respectively. Hence, there is no evidence for Υ(1S) → ± τ ∓ transitions. 6 Signal extraction for RLFV decays Our RLFV signal extraction procedure is very similar to that used for the corresponding non-radiative transition. We perform an UML fit to the mass difference ∆M = M ππγeµ − M γeµ . The signal PDF used is a sum of two Gaussians sharing a common mean. To Events/ (2 MeV/c Figure 6. ∆M fit to Υ(2S) data for the Υ(1S) → γe ± µ ∓ decay. The fitted signal PDF is represented by the filled red region and the dashed cyan line represents the background. The solid blue curve represents the overall fit to data. The long-dashed red curve represents the signal PDF corresponding to 5 hypothetical signal events. estimate the peaking background from leptonic decays of the Υ(1S), we use the same shape as the signal PDF as the shape of the background PDF and a large Υ(1S) → ± ∓ MC sample is used to have a more precise estimation. Other backgrounds are flat on the ∆M window and modeled with a first-order Chebyshev polynomial. Fig. 6 shows the ∆M fit for Υ(2S) data. To fit the data, we fix the width of the signal PDF from the MC signal width corrected by the average of data-MC difference for Υ(1S) → e + e − and Υ(1S) → µ + µ − modes. The yield of peaking background is estimated to be 0.1 ± 0.1. The signal efficiency for γeµ decay is 24.6%. From the Υ(2S) data fit, the signal yield for the Υ(1S) → γe ± µ ∓ decay is estimated to be 0.8 ± 1.5.

Υ(1S) → γ ± τ ∓ decay
To extract the signal for γµτ and γeτ decays, we define the recoil mass of ππγ (M recoil ππ γ ) using Eq. (3.1). We perform a UML fit to M recoil ππ γ to extract the signal yield and estimate efficiency. For signal events, M recoil ππ γ should peak at the nominal τ mass. A Gaussian and a bifurcated Gaussian sum sharing a common mean is used to model the signal events for Υ(1S) → γ ± τ ∓ decays.
For both of the Υ(1S) → γ ± τ ∓ decays, the dominant background comes from τ τ decays and hadronic decays of the Υ(1S). The τ τ background is treated using an approach similar to that for the Υ(1S) → µ ± τ ∓ decay, with the background shape as described by the Eq. (5.1). Also, we find some background from the radiative hadronic decays of Υ(1S).

Systematic uncertainty and correction
We calculate the systematic uncertainty from various sources such as the number of Υ(2S), track reconstruction, photon reconstruction, identification of pions from Υ(2S), lepton identification, uncertainty in signal efficiency, secondary branching fraction, and the fitting model.
The uncertainty on the number of Υ(2S) events was determined from a study of hadronic decays to be 2.3% [20]. Reconstruction efficiency of charged particle tracks are studied using a partially reconstructed D * + → D 0 [K 0 S (π + π − )π + π − ]π + decay sample with p T > 200 MeV/c. Systematic uncertainty per track is estimated to be 0.35%. Due to correlation, uncertainties in charged track finding are added linearly. The efficiency of photon reconstruction is estimated with radiative Bhabha events, and the associated uncertainty is 2.0% [21].
Uncertainty from pion identification in Υ(2S) → π + π − Υ(1S) reconstruction may affect our results. In order to estimate it, we use the results of a dedicated study based on the D * + → D 0 (K − π + )π + decay. A correction for the difference in efficiency (between data and signal MC) is obtained from the same source. This correction is used to correct the efficiency, and its uncertainty is included as the systematic uncertainty due to pion identification. For all the decays, the efficiency correction factor and systematic uncertainty from pion pair reconstruction are estimated to be 1.00 and 1.9%, respectively. For the electron identification with L e > 0.6 and the muon identification with L µ > 0.95, systematic uncertainty are calculated from the comparison between data and MC for 2γ → ee/µµ decays. We calculate an efficiency correction factor and systematic uncertainty for all of the electrons and muons using the same approach. For τ − → π − π + π − ν τ decay, the systematic uncertainty due to pion identification is estimated using the D * + -based method described above. The efficiency correction factors associated with the leptons pair reconstructions for eµ, µτ , eτ , γeµ, γµτ , and γeτ decays are 0.99, 0.98, 0.97, 0.95, 0.94, and 0.97, respectively and corresponding systematic uncertainties are 1.9%, 2.1%, 2.3%, 2.6%, 2.8%, and 2.5%, respectively.
Due to the limited number of generated MC signal events, there is an uncertainty in the fitted number of signal events as well as in the signal efficiency ( ), and the corresponding uncertainty is included in the systematic uncertainty.
We fix some parameters (such as the mean, width, and fractions of the two Gaussians) of the signal PDF while fitting the data. The associated systematic uncertainties are estimated by varying each of the fixed parameters by ±1σ from their central values and repeating the fit. For Υ(1S) → e + e − and Υ(1S) → µ + µ − decays, these PDF systematic uncertainties are estimated to be 0.08% and 0.04%, respectively. In the absence of significant signal events for the CLFV modes, we take the average value of the control modes (0.06%) as the systematic uncertainty from the signal PDF for each of the CLFV modes. Similarly, the systematic uncertainty from the background PDF for Υ(1S) → µ ± τ ∓ , Υ(1S) → e ± τ ∓ , Υ(1S) → γµ ± τ ∓ , and Υ(1S) → γe ± τ ∓ are estimated to be 1.4%, 0.8%, 0.7%, and 1.4%, respectively. No parameters were fixed to estimate the background for Υ(1S) → e ± µ ∓ , Υ(1S) → γe ± µ ∓ , Υ(1S) → e + e − , and Υ(1S) → µ + µ − modes. Therefore, for the above decays, the systematic uncertainty from the background PDF is estimated as zero. Due to correlation, systematic uncertainties from signal PDF and background PDF are added linearly. Table 1 summarizes the systematic uncertainties from various sources for all the modes. Systematic uncertainties from the different sources are added in quadrature in order to get the total systematic uncertainty for a particular signal mode. The systematic uncertainty due to the uncertainty in the peaking background for Υ(1S) → e ± µ ∓ (8.8±2.0) is directly included in the estimated upper limit of the branching fraction. For other modes, the effect of possible peaking background lowers the upper limit, and we do not consider it to report conservative upper limits.

Results
Using equation (4.1), the branching fractions are calculated as B[Υ(1S) → e + e − ] = (2.40 ± 0.01(stat) ± 0.12(syst)) × 10 −2 and B[Υ(1S) → µ + µ − ] = (2.46 ± 0.01(stat) ± 0.11(syst)) × 10 −2 which agree within ±1σ with world average values [19]. All of the results for the branching fractions of CLFV modes are dominated by statistical uncertainty. In the absence of significant signal, we estimate the upper limits (UL) of the branching fractions with a frequentist approach [22]. One can calculate the UL of branching fractions using the following relation: +2.1 ± 5.9 10.0 6.1 × 10 −6 − Υ(1S) → γe ± τ ∓ 5.0 −9.5 ± 6.3 9.1 6.5 × 10 −6 − where N UL sig is the UL on the signal yield after including systematic uncertainty. We perform 5000 pseudo-experiments by generating the fixed background from the final PDF and varying the yield of the input signal within 1 to 20. We use the corresponding PDF that has been used to fit Υ(2S) data for generating the data sets for pseudo-experiments. The fraction of pseudo-experiments with a fitted yield greater than the estimated signal yield in data has been taken as the confidence level (CL). Systematic uncertainties of the CLFV modes are included by smearing the yield of the pseudo-experiments within the fluctuations. For Υ(1S) → e ± µ ∓ decay, the fitted signal yields of pseudo-experiments have been smeared within the corresponding uncertainty of peaking background to include the associated systematic uncertainty.