Observation of the $B_s^0 \rightarrow J/\psi \phi \phi$ decay

The $B_s^0 \rightarrow J/\psi \phi \phi$ decay is observed in $pp$ collision data corresponding to an integrated luminosity of 3 fb$^{-1}$ recorded by the LHCb detector at centre-of-mass energies of 7 TeV and 8 TeV. This is the first observation of this decay channel, with a statistical significance of 15 standard deviations. The mass of the $B_s^0$ meson is measured to be $5367.08\,\pm \,0.38\,\pm\, 0.15$ MeV/c$^2$. The branching fraction ratio $\mathcal{B}(B_s^0 \rightarrow J/\psi \phi \phi)/\mathcal{B}(B_s^0 \rightarrow J/\psi \phi)$ is measured to be $0.0115\,\pm\, 0.0012\, ^{+0.0005}_{-0.0009}$. In both cases, the first uncertainty is statistical and the second is systematic. No evidence for non-resonant $B_s^0 \rightarrow J/\psi \phi K^+ K^-$ or $B_s^0 \rightarrow J/\psi K^+ K^- K^+ K^-$ decays is found.


Introduction
The study of the B 0 s → J/ψ φ φ decay, previously unobserved, allows a precise measurement of the B 0 s meson mass and a search for possible resonances in the φ φ and J/ψ φ invariant mass spectra, similar to what has been reported for the B + → J/ψ φ K + decay mode [1][2][3] (the inclusion of charge conjugate processes is implied throughout). The most recent theoretical predictions for heavy hadron masses, based on lattice QCD calculations, can be found in Refs. [4][5][6]. The current experimental knowledge of the B 0 s mass, as summarized in Ref. [7], is dominated by results from the LHCb experiment [8], which were obtained with the B 0 s → J/ψ φ decay using a small fraction of the integrated luminosity collected in the 2010-2012 LHC run. The B 0 s mass measurement using this decay is limited by the precision of the momentum scale. The B 0 s → J/ψ φ φ decay mode is a good alternative to B 0 s → J/ψ φ since the kinetic energy available to the final-state particles (Q-value) is much lower, leading to a 65% reduction in the systematic uncertainty arising from the precision of the momentum scale.
The B 0 s → J/ψ φ φ decay is also of interest in searches for intermediate states in the B 0 s decay chain. In recent years, many new charmonium or charmonium-like states have been discovered, which are not easily accommodated in the quark model of hadrons [9,10]. In a study of B + → J/ψ φ K + decays, the CDF collaboration reported evidence for a state, in the J/ψ φ invariant mass spectrum, called Y (4140) with mass and width values of m = 4143.0 ± 2.9 (stat) ± 1.2 (syst) MeV/c 2 and Γ = 11.7 +8.3 −5.0 (stat) ± 3.7 (syst) MeV [3]. The Belle and BaBar collaborations searched for the Y (4140) using the same B + decay mode [1,11] and found no significant signal, although the upper limits on the production rate did not contradict the CDF measurement. Recently, the D0 collaboration reported a similar structure [12]. At the LHC, both the LHCb and CMS collaborations have searched for the state in question. The LHCb collaboration found no evidence with 0.37 fb −1 of pp collision data [2], in 2.4 σ disagreement with the CDF measurement. A CMS search for the same signature [13] supports the CDF observation. With two out of five experiments failing to observe the Y (4140) resonance the question of its existence still remains open. The search for resonances in the φ φ invariant mass spectrum is also of interest. Several experiments have reported a near-threshold enhancement in the φ φ invariant mass distribution from the J/ψ → γ φ φ decay [14][15][16]. A partial-wave analysis showed that the structure is dominated by a 0 −+ state called η(2225). This resonance is still controversial and its observation in a different decay mode would be conclusive.
Theoretical predictions of the B 0 s → J/ψ φ φ branching fraction are difficult due to the presence of three vector mesons in the final states. The B 0 s → J/ψ φ φ decay is the B 0 s counterpart of the measured B + → J/ψ φ K + and B 0 → J/ψ φ K 0 decays [17]. All these channels are strongly suppressed with respect to the similar decays without the additional φ meson in the final state. The suppression factors of the last two channels are 0.048 ± 0.004 and 0.057 ± 0.012 for the charged and neutral decays [7]. A qualitative comparison with these branching fractions can be done considering that the phase space of the decay B 0 s → J/ψ φ φ is smaller by a factor of seven, so the B 0 s → J/ψ φ φ branching fraction is expected to be ∼ 10 −5 . This paper presents the first observation of the decay B 0 s → J/ψ φ φ and the decay branching fraction measurement with respect to the reference decay B 0 s → J/ψ φ. A measurement of the B 0 s mass is also presented. The data sample corresponds to an integrated luminosity of 3.0 fb −1 pp collisions collected by the LHCb experiment. The data were recorded in the years 2011 and 2012 at centre-of-mass energies of 7 TeV and 8 TeV, respectively.

The LHCb detector
The LHCb detector [18,19] 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 silicon-strip vertex detector surrounding the pp interaction region, 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 placed downstream of the magnet. The polarity of the dipole magnet is reversed periodically throughout data-taking. The tracking system provides a measurement of momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary vertex (PV), the impact parameter, is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction.
Simulated events are used to determine trigger, reconstruction and selection efficiencies and reconstructed mass distributions. In addition, simulated samples are used to estimate possible peaking backgrounds from B meson decays that can mimic the B 0 s → J/ψ φ (φ) final states. In the simulation, pp collisions are generated using Pythia 6 [20] 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]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [24] as described in Ref. [25].

Event selection
The final states of the signal and reference channels differ only by the presence of an extra φ meson in the former mode. The selections of the B 0 s → J/ψ φ φ and B 0 s → J/ψ φ candidates are done in almost the same way, allowing a partial cancellation of systematic uncertainties in the evaluation of the efficiency ratio. The J/ψ meson is reconstructed in the J/ψ → µ + µ − decay while the φ meson is reconstructed in the φ → K + K − decay.
Events are selected by the hardware triggers requiring a single muon with transverse momentum p T > 1.48 GeV/c or a muon pair with product of transverse momenta greater than (1.3 GeV/c) 2 . At the first stage of the software trigger, events are selected that contain two muon tracks with p T > 0.5 GeV/c and invariant mass m(µ + µ − ) > 2.7 GeV/c 2 , or a single muon track with p T > 1 GeV/c and χ 2 IP > 16 with respect to any PV. The quantity χ 2 IP is the difference between the χ 2 values of a given PV reconstructed with and without the track considered. The second stage of the software trigger selects a muon pair with an invariant mass that is consistent with the known J/ψ mass [7]. The decay length significance of the reconstructed J/ψ candidate, S L , is required to be greater than 3, where S L is the distance between the J/ψ vertex and the PV, divided by its uncertainty.
The offline analysis uses a cut-based preselection, followed by a multivariate analysis. In the preselection all the tracks are required to have a good-quality track fit. In the φ → K + K − decay reconstruction, kaons are selected with p > 3 GeV/c and p T > 200 MeV/c, and the vertex is required to have a good-quality fit. Particle identification (PID) is performed using information from all the subdetectors. A loose requirement is applied to the PID discriminant of kaons with respect to the pion misidentification DLL Kπ > 0, where DLL xπ = ln L x − ln L π is the delta-log-likelihood for the x particle hypothesis with respect to the pion. For the J/ψ → µ + µ − decay, the two muons are required to have p > 5 GeV/c and to satisfy a loose PID selection, DLL µπ > −1. The invariant mass of the J/ψ candidate is required to be in the interval [3036, 3156] MeV/c 2 , corresponding to a ±4σ interval around the nominal mass of the J/ψ meson [7]. To select the final B 0 s → J/ψ φ (φ) decay, the φ and J/ψ meson candidates are required to pass the selection cuts p T (φ) > 300 MeV/c and p T (J/ψ ) > 400 MeV/c, and to form a good-quality displaced vertex. The collinearity angle, defined as the angle between the reconstructed B 0 s momentum and the flight direction determined from the secondary vertex, is required to be smaller than 1.8 • . In the B 0 s → J/ψ φ decay selection, to reduce the contamination from non-resonant B 0 s → J/ψ K + K − decays, the dikaon invariant mass is required to be in the range [980, 1080] MeV/c 2 . To improve the mass and decay-time resolutions, a kinematic fit [26] is applied to both B 0 s decays, constraining the mass of the J/ψ candidate to its known value [7] and the B 0 s momentum to point to the PV. Finally, the B 0 s candidate invariant mass is required to be in the interval [5250, 5490] MeV/c 2 .
Different multivariate selection algorithms, based on a boosted decision tree (BDT) [27,28] with the AdaBoost algorithm [29], are used to select the signal and the reference channel samples. The BDT is trained with simulated B 0 s samples for the signals, while for the background, a sample of 40 million simulated events containing inclusive B → J/ψ X decays is used. For the B 0 s → J/ψ φ φ decay channel, the simulated sample is generated according to phase space. The BDT input variables are the p T of the φ and J/ψ mesons and the vertex χ 2 , flight distance significance, S L , collinearity angle and the impact parameter of the B 0 s meson with respect to the PV. The BDT discriminant threshold is chosen to maximise the figure of merit, /(3/2 + √ b) [30], where is the signal efficiency determined using simulated events and b is the number of expected background candidates estimated using mass sideband events in the data. For the B 0 s → J/ψ φ channel the BDT discriminant is selected to maximize s/ √ s + b, where s and b are the expected signal and background yields, estimated from simulated events and sideband data, respectively.
In the B 0 s → J/ψ φ φ selection, no restriction is initially put on the K + K − system invariant mass, with both the resonant and the non- also passes the cuts, resulting in a duplicated candidate. So a genuine resonant B 0 s → J/ψ φ(K + K − ) φ(K + K − ) event will most of the time produce also a "fake" non-resonant candidate, given the low probability that the invariant mass of two wrongly-coupled kaons is around the φ mass. In order to remove these "fake" candidates, the K + K − system masses are required to satisfy |m(K + K − ) − m φ | < 15 MeV/c 2 . After this cut, 1.8% of events contain double candidates. For each of these events, one candidate is chosen at random. In the B 0 s → J/ψ φ decay selection, this ambiguity problem is not present, so a tight cut on the |m(K + K − ) − m φ | is not applied.

Results
The B 0 s → J/ψ φ φ decay branching fraction is measured with respect to the reference decay where N obs (B 0 s → J/ψ φ φ) and N obs (B 0 s → J/ψ φ) are the numbers of observed events and (B 0 s → J/ψ φ φ) and (B 0 s → J/ψ φ) are the selection efficiencies. Figure 1 shows the invariant mass of the reconstructed B 0 s → J/ψ φ φ decay, for all the candidates surviving the pre-selection, the BDT and the selection on the m(K + K − ) around the φ mass. In order to evaluate the number of signal decays, an unbinned extended maximum likelihood fit is performed assuming a Gaussian signal peak and an exponential combinatorial background. The observed signal yield is 128 ± 13 events, where the uncertainty is statistical. Using Wilks's theorem [31], the statistical significance is found to be 15 standard deviations. As expected, the mass resolution is good, σ = 3.05±0.41 MeV/c 2 , due to the low Q-value of the decay.
In the reference channel, in order to discriminate between the resonant B 0 s → J/ψ φ and the non-resonant B 0 s → J/ψ K + K − decays, a fit to the K + K − mass spectrum is made, where the combinatorial background in the B 0 s mass window is statistically removed using the sPlot technique [32]. Figure 2 (left) shows the mass distribution of the B 0 s → J/ψ K + K − candidates with the fit results superimposed. The B 0 s peak is described by a double Crystal Ball function [33], while the underlying combinatorial background is described by an exponential function plus a second-order polynomial.  and a second-order polynomial for the non-resonant component. The observed B 0 s → J/ψ φ yield is 82120 ± 330 events where the uncertainty is statistical.
All the efficiencies (detector acceptance, reconstruction, trigger and selection) are evaluated using the simulated samples together with a data-driven method [34] for tracking and PID. To check the reliability of the simulation, a comparison is made between data and simulation for all of the kinematic variables used in the selection; good agreement is found. Since the ratio of (B 0 s → J/ψ φ φ) over (B 0 s → J/ψ φ) is evaluated, many systematic effects, related to possible small deviation of simulation with respect to data, cancel or are significantly reduced.
The efficiency ratio (B 0 s → J/ψ φ φ)/ (B 0 s → J/ψ φ) is evaluated to be 0.2778 ± 0.0015, where the uncertainty is statistical, due to the limited simulated sample sizes. As expected, the efficiency of the B 0 s → J/ψ φ φ channel is lower than that of B 0 s → J/ψ φ one, due to the presence of the additional φ → K + K − decay and the fact that on average the decay products have a smaller transverse momentum.
From the event yields and the ratio of efficiencies, and using the known φ → K + K − branching fraction [7], the branching fraction ratio is measured to be The systematic uncertainty will be discussed in Section 5.
From the fit to the B 0 s invariant mass distribution in the B 0 s → J/ψ φ φ decay, the mass of the B 0 s meson is measured to be m(B 0 s ) = 5367.08 ± 0.38 (stat) ± 0.15 (syst) MeV/c 2 . The J/ψ φ and φ φ mass distributions are shown in Fig. 3 for both data and simulation. For the data, the sPlot techique is used to subtract the background from the signal. Since the B 0 s → J/ψ φ φ process is a decay of a pseudoscalar into three vector mesons, its accurate description is complex and affected by large theoretical uncertainty. Here, to simulate the B 0 s → J/ψ φ φ decay, a simple phase-space decay model is used, which turns out not to provide a satisfactory description of the data. The disagreement can be due to either intermediate resonances or the simplified description of the decay. More data are needed to resolve the issue. Presently, due to the low statistics and the unknown decay dynamics, it is difficult to draw any conclusions from the two mass distributions.

Systematic uncertainties
A summary of the systematic uncertainties on the measurement of the branching fraction ratio is given in Table 1. Since the various effects are uncorrelated, the total systematic uncertainty is evaluated by adding all terms in quadrature.
The average multiplicity of B 0 s → J/ψ φ φ candidates in the simulated sample is 1.006 compared to 1.018 of the data. The relative difference (1.2%) is assigned as a systematic uncertainty.
The use of a simplified decay model affects the determination of the detection efficiency and introduces some bias in the measurement. In order to evaluate the effect, the simulated sample is used to study the efficiency of the selection as a function of the two masses m(φ, φ) and m(J/ψ, φ). The efficiency is then evaluated in a simulated sample reweighted in such a way as to reproduce the mass distributions in the data. A relative difference ∆ / = 1.0% is found and is assigned as a systematic uncertainty due to the unknown decay model. Alternative functions for describing the signal component are tested: double Gaussian or double Crystal Ball function for the B 0 s → J/ψ φ and single Gaussian or single Crystal Ball function for B 0 s → J/ψ φ φ. In both cases a negligible change in yields is observed and therefore no systematic uncertainty is assigned. Conversely, different choices of the background parametrisation in the B 0 s → J/ψ φ data can lead to sizeable difference in the results. In order to estimate a systematic uncertainty, the fits are repeated using an exponential, a second-order polynomial and the sum of the two (the nominal fit). The largest difference in yield, 1.6%, between the nominal fit and the fit with the exponential, is taken as the systematic uncertainty. The same procedure applied to the signal channel results in a 0.8% change in yield. The two systematic uncertainties are added in quadrature to give an overall uncertainty of 1.8% in the modelling of the signal and backgrounds.
To evaluate the contamination from non-resonant B 0 s → J/ψ φ K + K − and B 0 s → J/ψ K + K − K + K − decays, a dedicated search is performed for these two channels in the whole allowed kinematic region (without any requirement on the K + K − mass). The yields are then extrapolated to the restricted kinematic region of the signal. For the B 0 s → J/ψ K + K − K + K − decay, the sPlot technique is first used to select the B 0 s decay and then the two K + K − mass spectra are fitted simultaneously to determine the yield of the fully resonant decay candidates and the non-resonant ones. The non-resonant component is the sum of true non-resonant decays plus the candidates obtained by exchanging the kaons pairings in the resonant decays. When the latter component is subtracted from the measured yield, the number of non-resonant candidates is found to be 22 ± 18. Extrapolating this number to the φ meson mass region and using the Feldman-Cousins method [35] gives an upper limit of 1.5 events in the signal region at 68.3% confidence level. A similar procedure is followed for the B 0 s → J/ψ φ K + K − decay. One K + K − pair is required to have the mass in the non-resonant range, m(K + K − ) > 1080 MeV/c 2 . In these events, no evidence of a mass peak is found in the mass spectrum m(J/ψ K + K − K + K − ) nor in the mass spectrum of the other kaon pair. Using the Feldman-Cousins method an estimated contamination of 6.2 events is found at 68.3% confidence level. The uncertainties on the two non-resonant modes are added linearly, resulting in an asymmetric relative uncertainty of −6%.
The data-driven method used to correct the tracking efficiency for the two additional kaons in the final state of B 0 s → J/ψ φ φ with respect to B 0 s → J/ψ φ decay has an uncertainty of 1.5% per track, resulting in an overall relative uncertainty of 3.0%. This term also takes into account the uncertainty of hadronic interactions in the detector material.
Due to the decay time requirement on the selected events, the lack of knowledge of the admixture of B 0 sH and B 0 sL eigenstates in the B 0 s → J/ψ φ φ decay is a further source of systematic uncertainty [36]. While for the B 0 s → J/ψ φ decay the simulation uses the measured fractions of B 0 sH and B 0 sL states [37], the B 0 s → J/ψ φ φ decay is simulated assuming a completely symmetric combination. In order to evaluate the systematic effect, the simulated sample is reweighted assuming the two extreme cases where the B s meson is a complete B 0 sH or a B 0 sL state. The observed difference in the efficiency is 2.1% and this number is assigned as the systematic uncertainty.
A detailed comparison between data and simulation is performed for all the variables used in the BDT selection. For both the B 0 s → J/ψ φ φ and B 0 s → J/ψ φ decay channels, all variables show good agreement and the relative branching fraction result is stable against changes in the threshold of the BDT response. The total systematic uncertainty on the ratio of branching fractions is found to be +4.4 −7.4 %. Table 2 gives a summary of the systematic uncertainties of the B 0 s mass measurement. For the B 0 s mass determination, the momentum scale calibration is the main source of systematic uncertainty. The momentum scale takes into account the limited knowledge of the detector alignment. By comparing measured mass values for several charmed mesons with precisely known values, an uncertainty of 0.03% on the momentum scale is estimated [38]. The corresponding uncertainty in the B 0 s mass value is ±0.12 MeV/c 2 . The uncertainty in the kaon mass [7] will affect the B 0 s mass determination, while the uncertainty on the J/ψ mass has a negligible effect. The effect is estimated by repeating the fit with the kaon mass shifted by ±σ, where σ is the uncertainty on the known kaon mass. The observed mass variation, ±0.06 MeV/c 2 , is assigned as a systematic uncertainty.
The fit model for the signal and background of the invariant mass distributions is another source of systematic uncertainty. The effect is estimated by comparing to the nominal case the fit results with those from alternative functions. The systematic uncertainty from this effect is ±0.02 MeV/c 2 .
The energy loss of the kaons in the detector is another possible source of bias in the mass measurement. A detailed study of this effect has been performed in Ref. [8] for the B 0 s → J/ψ φ decay. Following the same procedure in B 0 s → J/ψ φ φ decay, the effect is found to be ±0.06 MeV/c 2 , which is assigned as a systematic uncertainty. The bias for neglecting the QED radiative corrections in the final state is negligible due to the restricted phase space [8]. The uncertainty due to detector alignment is also negligible. Combining all of the above sources in quadrature, the total systematic uncertainty on the mass measurement is found to be ±0.15 MeV/c 2 .
As a cross check, the mass measurement is performed separately in the two data-taking periods and in two samples with opposite magnet polarity. All the measurements are consistent within the uncertainties.

Conclusions
This paper presents the first observation of the B 0 s → J/ψ φ φ decay channel, with a signal yield of 128 ± 13. Taking the B 0 s → J/ψ φ decay as the reference channel the relative branching fraction is measured to be B(B s → J/ψ φ φ) B(B s → J/ψ φ) = 0.0115 ± 0.0012 (stat) +0.0005 −0.0009 (syst) .
From a fit to the B 0 s invariant mass distribution in the B 0 s → J/ψ φ φ decay, the mass of the B 0 s meson is measured to be m(B 0 s ) = 5367.08 ± 0.38 (stat) ± 0.15 (syst) MeV/c 2 . This value is consistent with previous LHCb results [8] and with the world average [7]. The overall uncertainty is 20% larger than the current most precise measurement. As the systematic uncertainty is a factor of two smaller, further improvement can be expected when larger datasets become available.