Observation of $B^0_{(s)} \to J/\psi p \overline{p}$ decays and precision measurements of the $B^0_{(s)}$ masses

The first observation of the decays $B^0_{(s)} \to J/\psi p \overline{p}$ is reported, using proton-proton collision data corresponding to an integrated luminosity of $5.2{\text{fb}}^{-1}$, collected with the LHCb detector. These decays are suppressed due to limited available phase space, as well as due to OZI or Cabibbo suppression. The measured branching fractions are \begin{align*} \mathcal{B}(B^0 \to J/\psi p \overline{p})&= (4.51\pm 0.40\; \text{(stat)} \pm 0.44\; \text{(syst)}) \times 10^{-7}, \\ \mathcal{B}(B^0_s \to J/\psi p \overline{p})&= (3.58\pm 0.19\; \text{(stat)} \pm 0.33\; \text{(syst)}) \times 10^{-6}. \end{align*} For the $B^0_s$ meson the result is much higher than the expected value of $ {\cal O} (10^{-9})$. The small available phase space in these decays also allows for the most precise single measurement of both the $B^0$ mass as ${5279.74 \pm 0.30\; ({\rm stat})\pm 0.10\; ({\rm syst})}$~MeV, and the $B^0_s$ mass as ${5366.85 \pm 0.19\; ({\rm stat})\pm 0.13\; ({\rm syst})}$~MeV.

Multiquark hadronic states beyond the well-studied quark-antiquark (meson) and three-quark (baryon) combinations remain elusive even sixty years after their prediction in the quark model [1,2]. Employing an amplitude analysis of Λ 0 b → J/ψ pK − decays, the LHCb collaboration has found states consistent with |uudcc pentaquarks decaying to J/ψ p [3,4] (charge conjugation is implied throughout this Letter). The decays B 0 (s) → J/ψ pp are sensitive to pentaquark searches in the J/ψ p and J/ψ p components and to glueball states [5,6] in the pp system. Baryonic B 0 (s) decays are also interesting to study the dynamics of the final baryon-antibaryon system and its characteristic threshold enhancement, whose underlying origin has still to be completely understood [7].
In the leading Feynman diagrams shown in Fig. 1, the B 0 mode is Cabibbo suppressed due to the presence of the CKM element V cd , while the B 0 s mode is OZI suppressed [2,8,9]. The naïve theoretical expectation for the branching fraction B(B 0 s → J/ψ pp) is at the level of 10 −9 [10]. However, the presence of an intermediate pentaquark or glueball state can enhance the decay rate. The authors of Ref. [10] pointed out the potential sensitivity of B 0 s → J/ψ pp decays to tensor glueball states via a possible resonant contribution of f J (2220) → pp, which could enhance the B 0 s → J/ψ pp decay branching fraction up to order 10 −6 . Hints towards such enhancements were noted in a previous LHCb measurement using 1 fb −1 of pp collision data, where no observation for either mode was made, but a 2.8 standard deviation excess was seen for the B 0 s → J/ψ pp decay [11]. These decays also allow for high-precision mass measurements. The kinetic energies in the B 0 (s) rest systems of the decay products (Q-values) are approximately 306 MeV for B 0 and 393 MeV for B 0 s decays. The small Q-values imply a very small contribution from momentum uncertainties to the B 0 (s) mass measurements. In this Letter, the first observation of these modes along with their branching fraction and B 0 and B 0 s mass measurements are reported employing a data sample corresponding to 5.2 fb −1 of pp collision data collected by the LHCb experiment. As a normalization mode, the copious B 0 s → J/ψ φ(→ K + K − ) sample is used, which is similar in topology to the signal channels.
The LHCb detector [12,13] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region [14], a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes [15] placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. 1 Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [16]. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [17]. The online event selection is performed by a trigger [18], comprising a hardware stage based on information from the muon system, followed by a software stage that applies a full event reconstruction. The software trigger is a combination of event categories mostly relying on identifying J/ψ decays consistent with a B-meson decay topology with two muon tracks originating from a secondary decay vertex detached from the primary pp collision point.
The pp collision data used in this analysis were collected at center-of-mass energies of 7 and 8 TeV (3 fb −1 ) and 13 TeV (2.2 fb −1 ), during the Run 1 (2011 and 2012) and Run 2 (2015 and 2016) run periods, respectively. The data taking conditions differ enough between the two run periods, such that they are analyzed separately and the results combined at the end.
Samples of simulated events are used to study the properties of the signal and control channels. The pp collisions are generated using Pythia [19] with a specific LHCb configuration [20]. Decays of hadronic particles are described by EvtGen [21], in which final-state radiation is generated using Photos [22]. For the B 0 s → J/ψ φ mode, simulation samples are generated according to a decay model based on results reported in Ref. [23], while the B 0 (s) → J/ψ pp signal modes are generated uniformly in phase space. The interactions of the generated particles with the detector and its response are implemented using the Geant4 toolkit [24] as described in Ref. [25].
The event selection relies on the excellent vertexing and charged particle identification (PID) capabilities of the LHCb detector. For a given particle, the associated primary vertex (PV) corresponds to that with the smallest χ 2 IP , defined as the difference in χ 2 between the PV fit including and excluding the particle. Signal candidates are formed starting with a pair of charged tracks, consistent with muons originating from a common vertex significantly displaced from its associated PV and with an invariant mass consistent with the J/ψ meson. Another pair of oppositely charged tracks, identified as protons and originating from a common vertex, is combined with the J/ψ candidate to form a B 0 (s) candidate. The entire decay topology is submitted to a kinematic fit where the dimuon invariant mass is constrained to the known J/ψ mass [26]. The B 0 s → J/ψ φ control mode candidates are reconstructed in a similar fashion, replacing the pp combination with a pair of charged tracks identified as K + K − candidates, required to have an invariant mass within ±5 MeV of the known φ-meson mass [26]. All charged tracks are required to be of good quality and have p T > 300 MeV (p T > 550 MeV) for p or K (µ). For the B 0 s → J/ψ φ mode, the contamination from B 0 → J/ψ K + π − decays with a pion misidentified as a kaon is rejected by imposing a B 0 mass veto and using PID information. At this stage, the combinatorial background dominates, comprising a correctly reconstructed J/ψ meson candidate combined with two unrelated charged tracks.
For further background suppression, two multivariate classifiers are applied, each employing a gradient-boosted decision tree (BDT) [27]. In the first stage, the BDT kin classifier, based on kinematical and topological variables of the B 0 s candidate, is trained using the B 0 s → J/ψ φ decays from simulation as signal proxy, and selected J/ψ K + K − candidates in the mass window [5450, 5700] MeV as background. For BDT kin , only kinematic variables whose distributions are similar between the signal and the control mode are employed. These include the p, p T , and χ 2 IP values of the B 0 s meson, the χ 2 probability from a kinematic fit [28] to the decay topology, and the impact parameter (IP) of the muons with respect to the associated PV. Prior to the training, a multidimensional gradientboosting algorithm [29] is used to weight the simulated B 0 s → J/ψ φ events to match background-subtracted data distributions in all the training variables. These weights are denoted as GB-weights. The background-subtracted data distributions are obtained using the sPlot technique [30]. Under the assumption that the relative corrections between data and simulation are similar among different B 0 (s) → J/ψ h + h − decay topologies, h + and h − being charged hadrons, the GB-weights obtained from the control mode are applied to the signal mode. To validate this assumption, similar GB-weights are derived using another control mode, B 0 → J/ψ K + π − , yielding similar results.
To choose the BDT kin selection cut, the B 0 s → J/ψ pp signal figure of merit, S/ √ S + B, is required to exceed five in a 2σ window around the B 0 s mass peak. The background yield, B, is estimated from a fit to the J/ψ pp invariant mass distribution. To estimate the expected signal yield, S, the central value of the B 0 s → J/ψ pp branching fraction quoted in Ref. [11] is used, along with the signal efficiency obtained from simulation.
In the final selection stage, a second classifier, BDT PID , uses the hadron PID information from the Cherenkov detector system to distinguish between pions, kaons and protons. Aside from PID, the BDT PID training variables also include the p, p T and χ 2 IP values of the protons. The signal sample is taken as the B 0 s → J/ψ pp simulation incorporating the GB-weights for the kinematic variables, while the background sample is taken from events in data with m(J/ψ pp) ∈ [5450, 5500] MeV. The hadron PID variables in the simulation require further corrections to be representative of data. The PID variables are obtained from high-yield calibration samples of Λ + c → pK − π + and D * + → D 0 (→ K − π + )π + decays, which can be selected as a function of the p, p T and the number of tracks in the event using only kinematic information [31]. The optimal BDT PID selection criterion is chosen by maximizing the figure of merit S/ √ S + B, with the initial signal and background yields obtained from a fit to the m(J/ψ pp) distribution after the BDT kin selection.
For the B 0 s → J/ψ φ control mode, the selection is performed using a dedicated classifier, BDT CS , which includes the kinematic variables considered in BDT kin with the addition of the PID information.
After application of all selection requirements, the background is predominantly combinatorial. Approximately 1% of the selected events contain more than one candidate at this stage; a single candidate is selected randomly. The efficiency of the trigger, detector acceptance, reconstruction and selection procedure is approximately 1%, as estimated from simulation.
The B 0 and B 0 s signal and background yields are determined via an extended maximum likelihood fit to the J/ψ pp invariant mass distribution in the range [5220, 5420] MeV. Each signal shape is modeled as the sum of two Crystal Ball [32] functions sharing a common peak position, with tails on either sides of the peak to describe the radiative and misreconstruction effects. The background shape is modeled by a first-order polynomial  with parameters determined from the fit to data. The signal-model parameters are determined from simulation and only the B 0 and B 0 s central mass values are left as free parameters in the fit to data. The detector invariant-mass resolution is in agreement with simulations within a factor of 1.007 ± 0.004 as determined with the control mode. Residual discrepancies are accounted for in the systematic uncertainties. In order to validate the fit model, 1000 mass spectra are generated according to the model and fitted employing an alternative model comprising three Gaussian components for the signal and an exponential function for background. The difference between the input value of the yields and the mean of the fitted yields from the alternative model is assigned as a systematic uncertainty. The mass fit for the control mode uses a similar B 0 s signal lineshape, with the background modeled by an exponential function. The result of the fit to the combined Run 1 and Run 2 control mode yields a signal of 136,800 ± 400. The corresponding fit to the signal-mode candidates is shown in Fig. 2 with the results reported in Table 1, where clear signals of B 0 and B 0 s are observed.
The branching fractions measured with respect to the B 0 s → J/ψ φ control mode are where f s /f d is the ratio of the b-quark hadronization probabilities into B 0 s and B 0 mesons, and N corr denotes efficiency-corrected signal yields. For the signal modes, since the physics model is not known a priori, an event-by-event efficiency correction is applied to the data. It is derived from simulation as a function of the kinematic variables, which are given in detail in the Appendix.
Since the control mode has a topology very similar to that of the signal mode, most of the systematic uncertainties cancel in the branching-fraction ratio measurement. Residual systematic effects of the PID efficiency estimation are due to the correction procedure. An alternative PID correction is considered using proton calibration samples from decays of the long-lived Λ baryon to a proton and a pion, instead of prompt Λ + c decays. The difference between the two methods is assigned as a systematic uncertainty. The degree to which the simulation describes hadronic interactions with the detector material is less accurate for baryons than it is for mesons [21]. Following Ref. [33], a systematic uncertainty of 4% (1.1%) per proton (kaon) is assigned. Other systematic effects include the choice of the fit model, the weighting procedure, the trigger efficiency, and the presence of events with more than one candidate. The overall systematic uncertainties on the ratio of branching fractions are 7.2% (7.2%) and 6.5% (6.6%) for B 0 s (B 0 ) meson in Run 1 and Run 2, respectively, where the relevant contributions, listed in Table 2, are added in quadrature. Since the detector and the analysis methods remain the same between the two run periods, the systematic uncertainties are fully correlated while the statistical uncertainties are uncorrelated. The combination of the measurements is taken as a weighted mean to give the branching fraction ratios where the first uncertainty is statistical and the second is systematic.
For the absolute branching-fraction determination, the value [34] as the product of the two branching ratios, B(B 0 s → J/ψ φ) = (10.50 ± 0.13 ± 0.64) × 10 −4 and B(φ → K + K − ) = 0.489 ± 0.005, and the ratio of fragmentation probabilities f s /f d = 0.256 ± 0.020 [35]. For the B 0 s -meson normalization, the updated ratio f s /f d = 0.259 ± 0.015 [35] is used in Run 1, while for Run 2 it has been multiplied by an additional scale factor of 1.068 ± 0.046 [36] to take into account the dependence on the center of mass energy. The small S-wave K + K − fraction under the φ(1020) resonance, F S = 0.0070 ± 0.0005 [34], is accounted for as a correction. The absolute branching fractions are then combined to give   where the systematic uncertainty is the sum in quadrature of the overall systematic contribution on the ratio of branching fractions, the normalization mode uncertainty and the f s /f d uncertainty for the B 0 s signal. Table 2 summarizes the systematic uncertainties separately for the run periods. The dominant contributions are the normalization, the PID, and the tracking systematic uncertainties. For the B 0 meson, the external normalization measurement from Run 1, [34] is used, while for Run 2 the additional energy-dependent correction on f s /f d has an uncertainty of 4.3%. For the B 0 s meson, the measured B(B 0 , resulting in an uncertainty on f s /f d independent of the run condition. In addition, the small Q-values of the B 0 (s) → J/ψ pp decays also allow for precise measurements of the B 0 and B 0 s masses, with a resolution of 3.3 MeV (3.8 MeV) for the B 0 (B 0 s ) meson. The sources of systematic uncertainties include momentum scaling due to imperfections in the magnetic-field mapping, uncertainties on particle interactions with the detector material, and the choice of the signal model, as reported in Table 3. The with a correlation of 4 × 10 −4 in the statistical uncertainty. These represent the most precise single measurements for the B 0 and B 0 s masses. In summary, the first observation of the B 0 → J/ψ pp and B 0 s → J/ψ pp decays is reported. The measured branching fraction for the B 0 → J/ψ pp decay is consistent with theoretical expectations [10] while that for B 0 s → J/ψ pp is enhanced by two orders of magnitude with respect to predictions without resonant contributions [10]. More data are needed for glueball and pentaquark searches through a full Dalitz plot analysis. The world's best single measurements of the B 0 and B 0 s masses are also reported.
A. Efficiency parameterization for the signal mode The dihadron (in red) and dilepton (in blue) coordinate systems lie back-to-back with a common verticalŷ axis. The angle between the decay planes is χ ∈ (−π, π], while the two helicity angles, θ h and θ , are defined in the dihadron and dilepton rest frames, respectively. The 4-body phase-space of the decay B → J/ψ (→ µ + µ − )h + h − , where h ∈ {p, K}, is fully described by four independent kinematic variables. One of them is the dihadron invariant mass m h + h − . For a given m h + h − , the topology can be described by three angles, shown in Fig. 3: • θ and θ h : the helicity angles defined in the dimuon and dihadron rest frames, respectively; • χ: the azimuthal angle between the two decay planes of the dilepton and dihadron systems.
Since the final state is self-conjugate, the h − and the µ − particles are chosen to define the angles, for both B 0 (s) and B 0 (s) mesons. For the signal mode, the overall efficiency, including trigger, detector acceptance and selection procedure, is obtained from simulation as a function of the four kinematic variables, ϕ ≡ {m pp , cos θ , cos θ h , χ }. Here, m pp and χ are normalized such that all four variables in ϕ lie in the range (−1, 1]. The efficiency is parameterized as the product of Legendre polynomials ε( ϕ) = i,j,k,l c i,j,k,l P (cos θ , i)P (cos θ h , j)P (χ , k)P (m pp , l), where P (x, n) is a Legendre polynomial of order n in x ∈ (−1, 1]. Employing the maximum order of the polynomials as {3, 7, 7, 5} for {m pp , cos θ , cos θ h , χ }, respectively, was found to give a good parameterization. Simulation samples are employed, where B 0 (s) → J/ψ pp events are generated uniformly in phase space. The coefficients, c i,j,k,l , are determined from the simulation using a moments technique employing the orthogonality of Legendre polynomials c i,j,k,l = C Nrecon n=0 2i + 1 2 2j + 1 2 2k + 1 2 2l + 1 2 × P (cos θ , i)P (cos θ h , j)P (χ , k)P (m pp , l).
The sum is over the number of reconstructed decays, N recon , in the simulation sample after all selection criteria. The prefactor C ensures appropriate normalization. For a given data candidate, the corresponding kinematic variables, ϕ, are reconstructed and the efficiency, ε( ϕ), is computed according to the parameterization. The candidate is subsequently assigned a weight, 1/ε( ϕ), to account for the detector efficiency.

B. Combining the Run 1 and Run 2 branching fraction results
The ratio of branching fractions for Run 1 and Run 2 are provided here separately. For Run 1, To combine the results from the Run 1 and Run 2 data samples, the systematic uncertainties are taken as fully correlated between the run periods, since the data are collected employing the same detector and analysis techniques. The statistical uncertainties are considered uncorrelated since the datasets are disjoint. The covariance matrices are constructed as V = σ 2 stat, Run 1 + σ 2 syst, Run 1 σ syst, Run 1 × σ syst, Run 2 σ syst, Run 1 × σ syst, Run 2 σ 2 stat, Run 2 + σ 2 syst, Run 2 The weighted mean value and uncertainty are then calculated as respectively, where the variable x denotes the branching fraction and the weights are obtained from the inverse of the aforementioned covariance matrix as