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

TL;DR

This paper reports the first observation of specific rare B meson decays into J/psi and proton-antiproton pairs, along with precise measurements of their masses, revealing unexpectedly high branching fractions for B_s decays.

Contribution

The study provides the first observation of $B^0_{(s)} o J/ heta p ar{p}$ decays and offers the most precise measurements of $B^0$ and $B_s$ meson masses to date.

Findings

First observation of $B^0_{(s)} o J/ heta p ar{p}$ decays.

Measured branching fractions are significantly higher than theoretical expectations for $B_s$ decays.

Achieved the most precise measurements of $B^0$ and $B_s$ meson masses.

Abstract

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 \mbox{LHCb}\xspace detector. These decays are suppressed due to limited available phase space, as well as due to Okubo-Zweig-Iizuka 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.39\; \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.

Tables

Table 1. Table 1: Signal yields and masses for B 0 superscript 𝐵 0 {B}^{0} and B s 0 subscript superscript 𝐵 0 𝑠 {B}^{0}_{s} mesons.

Mode	Yield	$B_{\text{scaleto}\left(6.6pts\text{scaleto}\right)6.6pt}^0$ mass [MeV]
$B^0\to J/\psi p\overline{p}$	$256\pm 22$	$5279.74\pm 0.30$
$B_s^0\to J/\psi p\overline{p}$	$609\pm 31$	$5366.85\pm 0.19$

Table 2. Table 2: Systematic uncertainties on the branching fraction measurements for Run 1 and Run 2. The total uncertainties on the branching fraction ratios ( BFR BFR {\rm BFR} ) are the sum of the systematic uncertainties, added in quadrature. The total uncertainties on the absolute branching fractions ( ℬ ℬ \mathcal{B} ) include the normalization and the uncertainties on the ratio f s / f d subscript 𝑓 𝑠 subscript 𝑓 𝑑 f_{s}/f_{d} from external measurments as well.

	$ℬ\left(B^0\to J/\psi p\overline{p}\right)$	$ℬ\left(B_s^0\to J/\psi p\overline{p}\right)$
	Run 1 (Run 2)	Run 1 (Run 2)
Fit model	1.0 (0.5)%	1.0 (0.9)%
Detector resolution	0.6 (0.5)%	0.4 (0.6)%
PID efficiency	5.0 (4.0)%	5.0 (4.0)%
Trigger	1.0 (1.0)%	1.0 (1.0)%
Tracking	5.0 (5.0)%	5.0 (5.0)%
Simulation weighting	0.4 (0.4)%	0.3 (0.3)%
Multiple candidates	0.1 (0.1)%	0.1 (0.1)%
Total on $\mathrm{BFR}$	7.2 (6.5)%	7.2 (6.6)%
Normalization	6.1 (6.1)%	6.1 (6.1)%
$f_s/f_d$	$-\left(4.3\right)\%$	5.8 (5.8)%
Total on $ℬ$	$\mathrm{\hspace{0.33em}\hspace{0.33em}9.4}\left(10.1\right)$%	11.1 (10.7)%

Table 3. Table 3: Systematic uncertainties of B 0 superscript 𝐵 0 {B}^{0} and B s 0 subscript superscript 𝐵 0 𝑠 {B}^{0}_{s} mass measurements.

	$B^0$	$B_s^0$
	[MeV]	[MeV]
Momentum scale	0.097	0.124
Mass fit model	0.020	0.020
Energy loss correction	0.030	0.030
Total	0.103	0.129

Equations

$$\mathcal{B}(B^0\to J/\psi p\overline{p})$$

$$\mathcal{B}(B_s^0\to J/\psi p\overline{p})$$

$$\frac{\mathcal{B}(B^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})\times f_s/f_d}$$

$$\frac{\mathcal{B}(B_s^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})}$$

$$\frac{\mathcal{B}(B^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})\times f_s/f_d}$$

$$\frac{\mathcal{B}(B_s^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})}$$

$$\mathcal{B}(B^0\to J/\psi p\overline{p})$$

$$\mathcal{B}(B_s^0\to J/\psi p\overline{p})$$

$$m_{B^0}$$

$$m_{B_s^0}$$

$$\epsilon (\vec{\phi })=\sum _{i,j,k,l}{c}_{i,j,k,l}P(\cos {\theta }_{\ell },i)P(\cos {\theta }_h,j)P({\chi }^{\prime },k)P(m_{p\overline{p}}^{\prime },l),$$

$$c_{i,j,k,l}=$$

$$\times P(\cos {\theta }_{\ell },i)P(\cos {\theta }_h,j)P({\chi }^{\prime },k)P(m_{p\overline{p}}^{\prime },l).$$

$$\frac{\mathcal{B}(B^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})\times f_s/f_d}$$

$$\frac{\mathcal{B}(B_s^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})}$$

$$\frac{\mathcal{B}(B^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})\times f_s/f_d}$$

$$\frac{\mathcal{B}(B_s^0\to J/\psi p\overline{p})}{\mathcal{B}(B_s^0\to J/\psi \varphi )\times \mathcal{B}(\varphi \to K^{+}{K}^{-})}$$

$$\mathcal{B}(B^0\to J/\psi p\overline{p})$$

$$\mathcal{B}(B_s^0\to J/\psi p\overline{p})$$

$$\mathcal{B}(B^0\to J/\psi p\overline{p})$$

$$\mathcal{B}(B_s^0\to J/\psi p\overline{p})$$

$$V=\left(\begin{array}{cc}{\sigma }_{\text{stat, Run 1}}^2+{\sigma }_{\text{syst, Run 1}}^2& {\sigma }_{\text{syst, Run 1}}\times {\sigma }_{\text{syst, Run 2}}\\ {\sigma }_{\text{syst, Run 1}}\times {\sigma }_{\text{syst, Run 2}}& {\sigma }_{\text{stat, Run 2}}^2+{\sigma }_{\text{syst, Run 2}}^2\end{array}\right),$$

$$V_{B^0}$$

$$V_{B_s^0}$$

$$\overline{x}$$

$${\sigma }_x^2$$

$$w_i=\frac{\sum _{k\in \{1,2\}}(V^{-1}{)}_{ik}}{\sum _{j,k\in \{1,2\}}(V^{-1}{)}_{jk}}.$$

Full text

EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

​​​

CERN-EP-2019-006

LHCb-PAPER-2018-046

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

LHCb collaboration†††Authors are listed at the end of this paper.

The first observation of the decays $B^{0}_{(s)}\rightarrow 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 Okubo-Zweig-Iizuka or Cabibbo suppression. The measured branching fractions are

[TABLE]

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})}\text{\,Me\kern-1.00006ptV}$, and the $B^{0}_{s}$ mass as ${5366.85\pm 0.19\;({\rm stat})\pm 0.13\;({\rm syst})}\text{\,Me\kern-1.00006ptV}$.

Published in Phys. Rev. Lett. 122, 191804 (2019)

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

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 ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}{{K}^{-}}$ decays, the LHCb collaboration has found states consistent with $|uudc\overline{c}\rangle$ pentaquarks decaying to ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}p$ [3, 4] (charge conjugation is implied throughout this Letter). The decays ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ are sensitive to pentaquark searches in the ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}p$ and ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ components and to glueball states [5, 6] in the $p\overline{p}$ system. Baryonic ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}$ 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 Cabibbo-Kobayashi-Maskawa element $V_{cd}$, while the ${{B}^{0}_{s}}$ mode is Okubo-Zweig-Iizuka suppressed [8, 2, 9]. The naïve theoretical expectation for the branching fraction $\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{})$ 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}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ decays to tensor glueball states via a possible resonant contribution of $f_{J}(2220)\rightarrow{p}{\overline{p}}$, which could enhance the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ decay branching fraction up to order $10^{-6}$. Hints towards such enhancements were noted in a previous LHCb measurement using 1$\text{\,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}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ decay [11].

These decays also allow for high-precision mass measurements. The kinetic energies in the ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}$ rest systems of the decay products ($Q$-values) are approximately $306\text{\,Me\kern-1.00006ptV}$ for ${B}^{0}$ and $393\text{\,Me\kern-1.00006ptV}$ for ${B}^{0}_{s}$ decays. The small $Q$-values imply a very small contribution from momentum uncertainties to the ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}$ 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$\text{\,fb}^{-1}$ of $pp$ collision data collected by the LHCb experiment. As a normalization mode, the copious ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi(\rightarrow{{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<\eta<5$, designed for the study of particles containing $b$ or $c$ quarks. The detector includes a high-precision tracking system consisting of a silicon-strip 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.111Natural units with $\hbar=c=1$ are used throughout. 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\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ 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$$\text{\,fb}^{-1}) and 13 TeV (2.2$$\text{\,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, *Sjostrand:2007gs] with a specific LHCb configuration [21]. Decays of hadronic particles are described by EvtGen [22], in which final-state radiation is generated using Photos [23]. For the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ mode, simulation samples are generated according to a decay model based on results reported in Ref. [24], while the ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ 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 [25, *Agostinelli:2002hh] as described in Ref. [27].

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 $\chi^{2}_{\rm IP}$, defined as the difference in $\chi^{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\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ meson. Another pair of oppositely charged tracks, identified as protons and originating from a common vertex, is combined with the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ candidate to form a ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}$ candidate. The entire decay topology is submitted to a kinematic fit where the dimuon invariant mass is constrained to the known ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mass [28]. The ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ control mode candidates are reconstructed in a similar fashion, replacing the ${p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ combination with a pair of charged tracks identified as ${{K}^{+}}{{K}^{-}}$ candidates, required to have an invariant mass within $\pm 5\text{\,Me\kern-1.00006ptV}$ of the known $\phi$-meson mass [28]. All charged tracks are required to be of good quality and have $p_{\mathrm{T}}$ $>$ 300  MeV ($p_{\mathrm{T}}$ $>$ 550  MeV) for $p$ or $K$ ($\mu$). For the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ mode, the contamination from $B^{0}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}K^{+}\pi^{-}$ 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\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ meson candidate combined with two unrelated charged tracks.

At this stage, a multidimensional gradient-boosting algorithm [29] is used to weight the simulated ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ 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)}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}h^{+}h^{\prime-}$ decay topologies, $h^{+}$ and $h^{\prime-}$ 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}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{K}^{+}}{{\pi}^{-}}$, yielding similar results.

For further background suppression, two multivariate classifiers are applied, each employing a gradient-boosted decision tree (BDT) [31]. In the first stage, the BDTkin classifier, based on kinematical and topological variables of the ${B}^{0}_{s}$ candidate, is trained using the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ decays from simulation as signal proxy, and selected ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{K}^{+}}{{K}^{-}}$ candidates in the mass window $[5450,5700]$ MeV as background. For BDTkin, only kinematic variables whose distributions are similar between the signal and the control mode are employed. These include the $p$, $p_{\mathrm{T}}$, and $\chi^{2}_{\rm IP}$ values of the ${{B}^{0}_{s}}$ meson, the $\chi^{2}$ probability from a kinematic fit [32] to the decay topology, and the impact parameter (IP) of the muons with respect to the associated PV.

To choose the BDTkin selection cut, the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ signal figure of merit, $S/\sqrt{S+B}$, is required to exceed five in a 2$\sigma$ window around the ${B}^{0}_{s}$ mass peak. The background yield, $B$, is estimated from a fit to the ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ invariant mass distribution. To estimate the expected signal yield, $S$, the central value of the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ branching fraction quoted in Ref. [11] is used, along with the signal efficiency obtained from simulation.

In the final selection stage, a second classifier, BDTPID, uses the hadron PID information from the Cherenkov detector system to distinguish between pions, kaons and protons. Aside from PID, the BDTPID training variables also include the $p$, $p_{\mathrm{T}}$ and $\chi^{2}_{\rm IP}$ values of the protons. The signal sample is taken as the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ simulation incorporating the GB-weights for the kinematic variables, while the background sample is taken from events in data with $m({{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}p\bar{p})\in[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 ${{\mathchar 28931\relax}^{+}_{c}}\rightarrow{p}{{K}^{-}}{{\pi}^{+}}$ and $D^{\ast+}\rightarrow D^{0}(\rightarrow K^{-}{{\pi}^{+}})\pi^{+}$ decays, which can be selected as a function of the $p$, $p_{\mathrm{T}}$ and the number of tracks in the event using only kinematic information [33]. The optimal BDTPID selection criterion is chosen by maximizing the figure of merit $S/\sqrt{S+B}$, with the initial signal and background yields obtained from a fit to the $m({{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{})$ distribution after the BDTkin selection.

For the ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi$ control mode, the selection is performed using a dedicated classifier, BDTCS, which includes the kinematic variables considered in BDTkin 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\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ invariant mass distribution in the range $[5220,\;5420]$ MeV. Each signal shape is modeled as the sum of two Crystal Ball [34] 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\pm 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\pm 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}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi}$ control mode are

[TABLE]

where $f_{s}/f_{d}$ is the ratio of the $b$-quark hadronization probabilities into ${{B}^{0}_{s}}$ and ${{B}^{0}}$ mesons, and $N^{\rm 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 ${\mathchar 28931\relax}$ baryon to a proton and a pion, instead of prompt ${{\mathchar 28931\relax}^{+}_{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 [22]. Following Ref. [35], 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

[TABLE]

where the first uncertainty is statistical and the second is systematic. For the absolute branching-fraction determination, the value ${\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)\times\mathcal{B}(\phi\rightarrow{{K}^{+}}{{K}^{-}})\times f_{s}/f_{d}=(1.314\pm 0.016\pm 0.079)\times 10^{-4}}$ is obtained from Ref. [36] as the product of the two branching ratios, ${\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)=(10.50\pm 0.13\pm 0.64)\times 10^{-4}}$ and $\mathcal{B}(\phi\rightarrow{{K}^{+}}{{K}^{-}})=0.489\pm 0.005$, and the ratio of fragmentation probabilities ${f_{s}/f_{d}=0.256\pm 0.020}$ [37]. For the ${B}^{0}_{s}$-meson normalization, the updated ratio ${f_{s}/f_{d}=0.259\pm 0.015}$ [37] is used in Run 1, while for Run 2 it has been multiplied by an additional scale factor of $1.068\pm 0.046$ [38] to take into account the dependence on the center of mass energy. The small $S$-wave $K^{+}K^{-}$ fraction under the $\phi(1020)$ resonance, $F_{S}=0.0070\pm 0.0005$ [36], is accounted for as a correction. The absolute branching fractions are then combined to give

[TABLE]

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, $\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)\times\mathcal{B}(\phi\rightarrow{{K}^{+}}{{K}^{-}})\times f_{s}/f_{d}$ [36] 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 $\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)\times\mathcal{B}(\phi\rightarrow{{K}^{+}}{{K}^{-}})\times f_{s}/f_{d}$ is divided by $f_{s}/f_{d}$ to obtain the ${B}^{0}_{s}$ normalization, $\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)\times\mathcal{B}(\phi\rightarrow{{K}^{+}}{{K}^{-}})$, resulting in an uncertainty on $f_{s}/f_{d}$ independent of the run condition.

In addition, the small $Q$-values of the ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ 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 derived using well-known narrow resonances, uncertainties on particle interactions with the detector material, and the choice of the signal model, as reported in Table 3. The uncertainty on the proton mass is neglected. The final results are

[TABLE]

with a correlation of $4\times 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}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ and ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ decays is reported. The measured branching fraction for the ${{B}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ decay is consistent with theoretical expectations [10] while that for ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ 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.

Acknowledgements

We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); MOST and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG (Germany); INFN (Italy); NWO (Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania); MSHE (Russia); MinECo (Spain); SNSF and SER (Switzerland); NASU (Ukraine); STFC (United Kingdom); NSF (USA). We acknowledge the computing resources that are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN (Italy), SURF (Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and Yandex LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), (France), KIT and DESY (Germany), INFN (Italy), SURF (Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and Yandex LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), PL-GRID (Poland) and OSC (USA). We are indebted to the communities behind the multiple open-source software packages on which we depend. Individual groups or members have received support from Fondazione Fratelli Confalonieri (Italy), AvH Foundation (Germany); EPLANET, Marie Skłodowska-Curie Actions and ERC (European Union); ANR, Labex P2IO and OCEVU, and Région Auvergne-Rhône-Alpes (France); Key Research Program of Frontier Sciences of CAS, CAS PIFI, and the Thousand Talents Program (China); RFBR, RSF and Yandex LLC (Russia); GVA, XuntaGal and GENCAT (Spain); the Royal Society and the Leverhulme Trust (United Kingdom); Laboratory Directed Research and Development program of LANL (USA).

Appendix

**A. Efficiency parameterization for the signal mode

**

The 4-body phase-space of the decay $B\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}(\rightarrow{\mu^{+}}{\mu^{-}})h^{+}h^{-}$, where $h\in\{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:

• •

$\theta_{\ell}$ and $\theta_{h}$: the helicity angles defined in the dimuon and dihadron rest frames, respectively;

• •

$\chi$: 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 $\mu^{-}$ particles are chosen to define the angles, for both ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}$ and ${\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}_{({s})}^{0}$ 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, $\vec{\varphi}\equiv\{m^{\prime}_{{p}\kern 0.70004pt\overline{\kern-0.70004ptp}{}},\cos\theta_{\ell},\cos\theta_{h},\chi^{\prime}\}$. Here, $m^{\prime}_{{p}\kern 0.70004pt\overline{\kern-0.70004ptp}{}}$ and $\chi^{\prime}$ are normalized such that all four variables in $\vec{\varphi}$ lie in the range $(-1,1]$. The efficiency is parameterized as the product of Legendre polynomials

[TABLE]

where $P(x,n)$ is a Legendre polynomial of order $n$ in $x\in(-1,1]$. Employing the maximum order of the polynomials as $\{3,7,7,5\}$ for $\{m^{\prime}_{{p}\kern 0.70004pt\overline{\kern-0.70004ptp}{}},\cos\theta_{\ell},\cos\theta_{h},\chi^{\prime}\}$, respectively, was found to give a good parameterization. Simulation samples are employed, where ${B}^{0}_{\scaleto{(}{6.6pt}s\scaleto{)}{6.6pt}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{p}\kern 1.00006pt\overline{\kern-1.00006ptp}{}$ 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

[TABLE]

The sum is over the number of reconstructed decays, $N_{\rm recon}$, in the simulation sample after all selection criteria. The prefactor $C$ ensures appropriate normalization. For a given data candidate, the corresponding kinematic variables, $\vec{\varphi}$, are reconstructed and the efficiency, $\varepsilon(\vec{\varphi})$, is computed according to the parameterization. The candidate is subsequently assigned a weight, $1/\varepsilon(\vec{\varphi})$, 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,

[TABLE]

while for Run 2,

[TABLE]

The absolute branching fractions obtained from our present knowledge of ${\mathcal{B}({{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\phi)\times\mathcal{B}(\phi\rightarrow{{K}^{+}}{{K}^{-}})\times f_{s}/f_{d}}=(1.314\pm 0.016\pm 0.079)\times 10^{-4}$ and ${f_{s}/f_{d}=0.259\pm 0.015}$ are, for Run 1,

[TABLE]

and for Run 2,

[TABLE]

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

[TABLE]

with the numerical values for the absolute branching fraction combination as

[TABLE]

The weighted mean value and uncertainty are then calculated as

[TABLE]

respectively, where the variable $x$ denotes the branching fraction and the weights are obtained from the inverse of the aforementioned covariance matrix as

[TABLE]