The Belle Collaboration

The charmless hadronic decay $B_{s}^{0}\rightarrow\eta^{\prime}\eta$ is suppressed in the Standard Model (SM) and proceeds only through transitions sensitive to Beyond-the-Standard-Model (BSM) physics Bevan:2014iga . BSM scenarios, such as a fourth generation of fermions, supersymmetry with broken R-parity, and a two-Higgs doublet model with flavor-changing neutral currents, could affect the branching fraction and CP asymmetry of this decay belleiiphysicsbook . The expected branching fraction for $B_{s}^{0}\rightarrow\eta^{\prime}\eta$ in the SM spans a range of $(2-4)\times 10^{-5}$ bf1 ; bf2 ; bf3 ; bf4 ; bf5 . Once branching fractions for two-body decays $B_{d,s}^{0}\to\eta\eta$, $\eta^{\prime}\eta$, and $\eta^{\prime}\eta^{\prime}$ are measured, it would be possible to extract CP-violating parameters using a formalism based on SU(3)/U(3) symmetry bf1 . To achieve this goal, at least four of these six branching fractions need to be measured. Only the branching fraction for $B_{s}^{0}\to\eta^{\prime}\eta^{\prime}$ has been measured so far  bsepep .

In this Letter, we report the results of the first search for the decay $B_{s}^{0}\rightarrow\eta^{\prime}\eta$ using the full Belle data sample of $121.4~\textrm{fb}^{-1}$ collected at the $\Upsilon(5S)$ resonance. The inclusion of the charge-conjugate decay mode is implied throughout. The Belle detector was a large-solid-angle magnetic spectrometer that operated at the KEKB asymmetric-energy $e^{+}e^{-}$ collider KEKB . The detector components relevant to our study include a tracking system comprising a silicon vertex detector (SVD) and a central drift chamber (CDC), a particle identification (PID) system that consists of a barrel-like arrangement of time-of-flight scintillation counters (TOF) and an array of aerogel threshold Cherenkov counters (ACC), and a CsI(Tl) crystal-based electromagnetic calorimeter (ECL). All these components were located inside a superconducting solenoid coil that provided a 1.5 T magnetic field. A detailed description of the Belle detector can be found elsewhere Belle .

The $\Upsilon(5S)$ resonance decays into $B_{s}^{*0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{*0}$, $B_{s}^{*0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{0}$, and $B_{s}^{0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{0}$ pairs, where the relative fractions of the two former decays are $f_{B_{s}^{*0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{*0}}=(87.0\pm 1.7)\%$ and $f_{B_{s}^{*0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{0}}=(7.3\pm 1.4)\%$ frac , respectively. Signal $B_{s}^{0}$ mesons originate from the direct decays of $\Upsilon(5S)$ or from radiative decays of the excited vector state $B_{s}^{*0}$. The $\Upsilon(5S)$ production cross section is $340\pm 16$ pb frac . To present our nominal result for $\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)$ we use the world average value for the fraction of $B_{s}^{(*)0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{(*)0}$ in $b\bar{b}$ events $f_{s}=0.201\pm 0.031$ PDG , the data sample is therefore estimated to contain $(16.60\pm 2.68)\times 10^{6}$ $B_{s}^{0}$ mesons. We also report the results for $f_{s}\times\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)$.

To maximize discovery potential of the analysis and to validate the signal extraction procedure, we use a sample of background Monte Carlo (MC) simulated events equivalent to six times the data statistics. In addition, to estimate the overall reconstruction efficiency we use a high-statistics signal MC sample, where the other $B_{s}^{(\ast)0}$ meson decays according to known branching fractions PDG . Both samples are used to develop a model implemented in the unbinned extended maximum-likelihood (ML) fit to data. The MC-based model is validated with a control sample of $B^{0}\rightarrow\eta^{\prime}K_{S}^{0}$ decays reconstructed from 711 ${\rm fb^{-1}}$ of $\Upsilon(4S)$ data.

We reconstruct $\eta$ candidates using pairs of electromagnetic showers not matched to the projections of charged tracks to the ECL and therefore identified as photons. We require that the reconstructed energies of these showers exceed 50 (100) MeV in the barrel (endcap) region of the ECL. The larger energy threshold for the endcaps is due to the larger beam-related background in these regions. To reject hadronic showers mimicking photons, the ratio of the energies deposited by a photon candidate in the $(3\times 3)$ and $(5\times 5)$ ECL crystal arrays centered on the crystal with the largest deposited energy is required to exceed 0.75. The reconstructed invariant mass of the $\eta$ candidates is required to be $515\leq M(\gamma\gamma)\leq 580$ ${\rm MeV}/{\it c}^{2}$, which corresponds, approximately, to a $\pm 3\sigma$ Gaussian resolution window around the nominal $\eta$ mass PDG . To suppress combinatorial background arising due to low-energy photons, the magnitude of the cosine of the helicity angle ($\cos\theta_{\textrm{hel}}$) is required to be less than 0.97, where $\theta_{\textrm{hel}}$ is the angle in the $\eta$ candidate’s rest frame between the directions of its Lorentz boost from the laboratory frame and one of the photons.

The $\eta^{\prime}$ candidates are formed by combining pairs of oppositely charged pions with the $\eta$ candidates. We require the reconstructed $\eta^{\prime}$ invariant mass to be in the range $920\leq M(\pi^{+}\pi^{-}\eta)\leq 980$ ${\rm MeV}/{\it c}^{2}$, which corresponds, approximately, to the range $[-10,+6]\sigma$ of the Gaussian resolution, after performing a kinematic fit constraining the reconstructed mass of the $\eta$ candidate to the nominal $\eta$ mass PDG . To identify charged pion candidates, the ratios of PID likelihoods, $R_{i/\pi}={{\mathcal{L}}}_{i}/({\mathcal{L}}_{\pi}+{\mathcal{L}}_{i})$, are used, where $L_{\pi}$ is the likelihood for the track being a pion, while $L_{i}$ is the corresponding likelihood for the kaon ($i=K$) or electron ($i=e$) hypotheses. We require $R_{K/\pi}\leq 0.6$ and $R_{e/\pi}\leq 0.95$ for pion candidates. The likelihood for each particle species is obtained by combining information from CDC, TOF and ACC nakano_pid , and (for electrons only) ECL eid . According to MC studies, these requirements reject 28% of background, while the resulting efficiency loss is below 3%. Charged pion tracks are required to originate from near the interaction point (IP) by restricting their distance of closest approach to the $z$ axis to be less than 4.0 cm along the $z$ axis and 0.3 cm perpendicular to it, respectively. The $z$ axis is opposite to the direction of the $e^{+}$ beam. These selection criteria suppress beam-related backgrounds and reject poorly reconstructed tracks. To reduce systematic uncertainties associated with track reconstruction efficiency, the transverse momenta of charged pions are required to be greater than 100 ${\rm MeV}/{\it c}$.

To identify $B_{s}^{0}\rightarrow\eta^{\prime}\eta$ candidates we use (shown here in natural units) the beam-energy-constrained $B_{s}^{0}$ mass, $M_{\rm bc}=\sqrt{E_{\rm beam}^{2}-p_{B_{s}}^{2}}$, the energy difference, $\Delta E=E_{B_{s}}-E_{\rm beam}$, and the reconstructed invariant mass of the $\eta^{\prime}$, where $E_{\rm beam}$, $p_{B_{s}}$ and $E_{B_{s}}$ are the beam energy, the momentum and energy of the $B_{s}^{0}$ candidate, respectively. All these quantities are calculated in the $e^{+}e^{-}$ center-of-mass frame. To improve the $\Delta E$ resolution, the $\eta^{\prime}$ candidates are further constrained to the nominal mass of $\eta^{\prime}$, though most of the improvement comes from the $\eta$ mass constraint. Signal candidates are required to satisfy selection criteria $M_{\rm bc}>5.3$ ${\rm GeV}/{\it c}^{2}$ and $-0.4\leq\Delta E\leq 0.3$ GeV. In a Gaussian approximation, the $\Delta E$ resolution is approximately 40 MeV. Similarly, the $M_{\rm bc}$ resolution is 4 ${\rm MeV}/c^{2}$. To take advantage of all available information in case the data indicate signal presence, we include $M(\pi^{+}\pi^{-}\eta)$ in the three-dimensional (3D) ML fit used to statistically separate the signal from background. We define the signal region: $5.35<M_{\rm bc}<5.43$ ${\rm GeV}/{\it c}^{2}$, $-0.25<\Delta E<0.10$ GeV, and $0.94<M(\pi^{+}\pi^{-}\eta)<0.97$ ${\rm GeV}/{\it c}^{2}$. The area outside the signal region is considered as sideband. To optimize sensitivity we use a narrower signal region $5.39<M_{\rm bc}<5.43$ ${\rm GeV}/{\it c}^{2}$ which would contain the largest signal contribution.

Hadronic continuum events from $e^{+}e^{-}\to q\bar{q}$ ($q=u,d,c,s$) are the primary source of background. Because of large initial momenta of the light quarks, continuum events exhibit a “jetlike” event shape, while $B_{s}^{(*)0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{(*)0}$ events are distributed isotropically. We utilize modified Fox-Wolfram moments ksfw , used to describe the event topology, to discriminate between signal and continuum background. A likelihood ratio ($\mathcal{LR}$) is calculated using Fisher discriminant coefficients obtained in an optimization based on these moments. We suppress the background using a discovery-optimized selection on $\mathcal{LR}$ obtained by maximizing the value of Punzi’s figure of merit punzi :

The background events containing real $\eta^{\prime}$ mesons exhibit a peak in the $M(\pi^{+}\pi^{-}\eta)$ distribution, however, they are distributed smoothly in $M_{\rm bc}$ and $\Delta E$. The fraction of this peaking background is a free parameter in our ML fits.

About 14% of the reconstructed signal MC events contain multiple candidates primarily arising due to misreconstructed $\eta$ mesons. In such events we retain the candidate with the smallest value of $\sum{\chi^{2}_{\eta}}+\chi^{2}_{\pi^{+}\pi^{-}}$, where $\chi^{2}_{\eta}$ denotes the $\eta$ mass-constrained fit statistic, the summation is over the two $\eta$ candidates, and $\chi^{2}_{\pi^{+}\pi^{-}}$ quantifies the quality of the vertex fit for two pion tracks. Simulation shows that this procedure selects the correct $B_{s}^{0}$ candidate in 62% of such events. The overall reconstruction efficiency is 10%.

To extract the signal yield, we perform an unbinned extended ML fit to the 3D distribution of $M_{\rm bc}$, $\Delta E$, and $M(\pi^{+}\pi^{-}\eta)$. The likelihood function is

To validate our fitting model and adjust the PDF shape parameters used to describe the signal, we use the control sample of $B^{0}\rightarrow\eta^{\prime}K_{S}^{0}$ decays. We reconstruct $K_{S}^{0}$ candidates via secondary vertices associated with pairs of oppositely charged pions ks_reco using a neural network technique NN . The following information is used in the network: the momentum of the $K_{S}^{0}$ candidate in the laboratory frame; the distance along the $z$ axis between the two track helices at the point of their closest approach; the flight length in the $x-y$ plane; the angle between the $K_{S}^{0}$ momentum and the vector joining the $K_{S}^{0}$ decay vertex to the IP; the angle between the pion momentum and the laboratory-frame $K_{S}^{0}$ momentum in the $K_{S}^{0}$ rest frame; the distance-of-closest-approach in the $x-y$ plane between the IP and the two pion helices; and the pion hit information in the SVD and CDC. The selection efficiency is 87% over the momentum range of interest. We also require that the reconstructed $\pi^{+}\pi^{-}$ invariant mass is within 12 ${\rm MeV}/{\it c}^{2}$, which is about 3.5$\sigma$, of the nominal $K_{S}^{0}$ mass PDG . We require $5.20\leq M_{\rm bc}\leq 5.30$ ${\rm GeV}/{\it c}^{2}$ for $B^{0}$ candidates. The control-sample signal region is $5.27<M_{\rm bc}<5.29$ ${\rm GeV}/{\it c}^{2}$, $-0.20<\Delta E<0.10$ GeV, and $0.94<M(\pi^{+}\pi^{-}\eta)<0.97$ ${\rm GeV}/{\it c}^{2}$. All other selection criteria are the same as those used to select $B_{s}^{0}$ candidates. This control sample is used to validate the $\eta$ and $\eta^{\prime}$ reconstruction and its effect on the resolution functions and PDF shape parameters. The validation of $K_{S}^{0}$ reconstruction was performed previously in a similar $B_{s}^{0}$ analysis Pal:2015ghq .

The presence of four photons in the final state gives rise to a sizable misreconstruction probability for the signal events. We study these self-crossfeed (SCF) events using the signal MC sample. A large correlation between $M_{\rm bc}$ and $\Delta E$ for such signal events is taken into account by describing the correctly reconstructed signal and SCF components separately with two different PDF sets. The latter comprise approximately 14% of the reconstructed signal and are excluded from the estimate of its efficiency. The Pearson correlation coefficient for the region with largest correlations for SCF signal events is 27%.

A sum of a Gaussian and a Crystal Ball xbal function is used to model the correctly reconstructed signal in each of the three fit variables. For $M_{\rm bc}$ and $M(\pi^{+}\pi^{-}\eta)$ we use a sum of these two functions with the same mean but different widths, while for $\Delta E$ both the mean and width are different. A Bukin function bukin and an asymmetric Gaussian are used to model the SCF contribution in $M_{\rm bc}$ and $\Delta E$, respectively. For $M(\pi^{+}\pi^{-}\eta)$, we use a sum of a Gaussian and a first-order Chebyshev polynomial. In our nominal fit to data the fraction of correctly reconstructed signal is fixed to its MC value. The signal PDF shape parameters for $M_{\rm bc}$ and $\Delta E$ are validated using the $B^{0}\rightarrow\eta^{\prime}K_{S}^{0}$ control sample.

We use an ARGUS argus function to describe the background distribution in $M_{\rm bc}$ and a first-order Chebyshev polynomial for $\Delta E$. To model the peaking part in $M(\pi^{+}\pi^{-}\eta)$ we use the signal PDF, because the peak is due to real $\eta^{\prime}$ mesons, while an additional first-order Chebyshev polynomial is used for the non-peaking contribution. The projections of the fit to the $B^{0}\rightarrow\eta^{\prime}K_{S}^{0}$ control sample are shown in Fig. 1.

To further test and validate our fitting model, ensemble tests are carried out by generating MC pseudoexperiments. In these experiments we use PDFs obtained from full detector simulation and the $B^{0}\rightarrow\eta^{\prime}K_{S}^{0}$ data. We perform 1000 pseudoexperiments for each assumed number of signal events. An ML fit is executed for each sample prepared in these experiments. The signal yield distribution obtained from these fits exhibits good linearity. We use the results of pseudoexperiments to construct classical confidence intervals (without ordering) using a procedure due to Neyman frequentist_approach . For each ensemble of pseudoexperiments, the lower and upper ends of the respective confidence interval represent the values of fit signal yields for which 10% of the results lie below and above these values, respectively. These intervals are then combined to prepare a classical confidence belt belt_method ; belt_method_2 used to make a statistical interpretation of the results obtained from data. The confidence intervals prepared using this statistical method are known to slightly “overcover” for the number of signal events fc , therefore resulting in a conservative upper limit.

We apply the 3D model to the data and obtain $2.7\pm 2.5$ signal and $57.3\pm 7.8$ background events. The signal-region projections of the fit are shown in Fig. 2. We observe no significant signal and estimate a 90% confidence-level (CL) upper limit on the branching fraction for the decay $B_{s}^{0}\rightarrow\eta^{\prime}\eta$ using the following formula:

Sources of systematic uncertainties and their relative contributions are listed in Table 1. The relative uncertainties on $f_{s}$ and $\sigma(\Upsilon(5S))$ are 15.4% and 4.7%, respectively. The systematic uncertainty due to $\eta$ reconstruction is 2.1% per $\eta$ candidate eta_syst . Track reconstruction track_syst and PID systematic uncertainties are 0.35% and 2% per track, respectively. We estimate the systematic uncertainty due to the $\mathcal{LR}$ requirement to be 10%, which represents the relative change in efficiency when this requirement is varied by $\pm$0.02 about the nominal value of 0.95. This range of variation is defined by the statistics of the control sample which is used to validate the efficiency and its dependence on the $\mathcal{LR}$ requirement. Systematic uncertainty due to signal PDF shape is estimated by varying the fixed parameters within their statistical uncertainties determined with $B^{0}\rightarrow\eta^{\prime}K_{S}^{0}$ data. When varying these parameters, we observe an 11% change in the signal yield obtained from the data and use this number as an estimate of PDF parametrization systematics. Systematic uncertainty due to $f_{B_{s}^{(*)0}\accentset{\rule{4.91673pt}{0.6pt}}{B}_{s}^{(*)0}}$ is evaluated by varying relative fractions of possible contributions to signal PDF and is 1.3%. When varying the SCF contribution by $\pm 50$% of itself, we observe a 4% change in the results of the fit to data, which we use as an estimate of SCF PDF systematic uncertainty. The relative uncertainties on $\eta$ and $\eta^{\prime}$ branching fractions are 1% and 1.2%, respectively. The statistical uncertainty due to MC statistics is estimated to be 0.1%. The overall systematic uncertainties for $\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)$ and $f_{s}\times\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)$ are estimated by adding the individual contributions in quadrature and are 23.1% and 17.2%, respectively. These systematic uncertainties are included in the $N_{\textrm{UL}}^{90\%}$ estimates of approximately 7.0 and 6.9 events by smearing the fit yield distributions while constructing the confidence belt used to extract the results. We estimate the upper limits on the branching fraction $\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)<6.5\times 10^{-5}$ and on the product $f_{s}\times\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)<1.3\times 10^{-5}$ at 90% CL. Finally, using the number of signal events obtained from the fit we estimate $\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)=(2.5\pm 2.2\pm 0.6)\times 10^{-5}$ and $f_{s}\times\mathcal{B}(B_{s}^{0}\rightarrow\eta^{\prime}\eta)=(0.51\pm 0.44\pm 0.09)\times 10^{-5}$, where, for each of the two estimates, the first uncertainty is statistical and the second is systematic. We summarize the results in Table 2.

In summary, we have used the full data sample recorded by the Belle experiment at the $\Upsilon(5S)$ resonance to search for the decay $B_{s}^{0}\rightarrow\eta^{\prime}\eta$. We observe no statistically significant signal and set a 90% CL upper limit of $6.5\times 10^{-5}$ on its branching fraction. To date, our result is the only experimental information on $B_{s}^{0}\rightarrow\eta^{\prime}\eta$ and is twice as large as the most optimistic SM-based theoretical prediction. This decay can be probed further at the next-generation Belle II experiment belle2 at the SuperKEKB collider in Japan.

We thank the KEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group, and the Pacific Northwest National Laboratory (PNNL) Environmental Molecular Sciences Laboratory (EMSL) computing group for strong computing support; and the National Institute of Informatics, and Science Information NETwork 5 (SINET5) for valuable network support. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the Promotion of Science (JSPS), and the Tau-Lepton Physics Research Center of Nagoya University; the Australian Research Council including grants DP180102629, DP170102389, DP170102204, DP150103061, FT130100303; Austrian Federal Ministry of Education, Science and Research (FWF) and FWF Austrian Science Fund No. P 31361-N36; the National Natural Science Foundation of China under Contracts No. 11435013, No. 11475187, No. 11521505, No. 11575017, No. 11675166, No. 11705209; Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; the CAS Center for Excellence in Particle Physics (CCEPP); the Shanghai Pujiang Program under Grant No. 18PJ1401000; the Shanghai Science and Technology Committee (STCSM) under Grant No. 19ZR1403000; the Ministry of Education, Youth and Sports of the Czech Republic under Contract No. LTT17020; Horizon 2020 ERC Advanced Grant No. 884719 and ERC Starting Grant No. 947006 “InterLeptons” (European Union); the Carl Zeiss Foundation, the Deutsche Forschungsgemeinschaft, the Excellence Cluster Universe, and the VolkswagenStiftung; the Department of Atomic Energy (Project Identification No. RTI 4002) and the Department of Science and Technology of India; the Istituto Nazionale di Fisica Nucleare of Italy; National Research Foundation (NRF) of Korea Grant Nos. 2016R1D1A1B01010135, 2016R1D1A1B02012900, 2018R1A2B3003643, 2018R1A6A1A06024970, 2018R1D1A1B07047294, 2019K1A3A7A09033840, 2019R1I1A3A01058933; Radiation Science Research Institute, Foreign Large-size Research Facility Application Supporting project, the Global Science Experimental Data Hub Center of the Korea Institute of Science and Technology Information and KREONET/GLORIAD; the Polish Ministry of Science and Higher Education and the National Science Center; the Ministry of Science and Higher Education of the Russian Federation, Agreement 14.W03.31.0026, and the HSE University Basic Research Program, Moscow; University of Tabuk research grants S-1440-0321, S-0256-1438, and S-0280-1439 (Saudi Arabia); the Slovenian Research Agency Grant Nos. J1-9124 and P1-0135; Ikerbasque, Basque Foundation for Science, Spain; the Swiss National Science Foundation; the Ministry of Education and the Ministry of Science and Technology of Taiwan; and the United States Department of Energy and the National Science Foundation.