Study of $B^0_s \rightarrow J/\psi\pi^+\pi^-K^+K^-$ decays

The decays $B^0_s \rightarrow J/\psi \pi^+\pi^- K^+ K^-$ are studied using a data set corresponding to an integrated luminosity of 9fb$^{-1}$, collected with the LHCb detector in proton-proton collisions at centre-of-mass energies of 7, 8 and 13TeV. The decays $B^0_s \rightarrow J/\psi K^{\ast0} \bar{K}^{\ast0}$ and $B^0_s \rightarrow \chi_{c1}(3872)K^+K^-$, where the $K^+K^-$ pair does not originate from a $\phi$ meson, are observed for the first time. Precise measurements of the ratios of branching fractions between intermediate $\chi_{c1}(3872)\phi$, $J/\psi K^{\ast0}\bar{K}^{\ast0}$, $\psi(2S)\phi$ and $\chi_{c1}(3872)K^+K^-$ states are reported. A structure, denoted as $X(4740)$, is observed in the $J/\psi\phi$ mass spectrum and, assuming a Breit-Wigner parameterisation, its mass and width are determined to be \begin{eqnarray*} m_{X(4740)}&=&4741 \pm 6 \pm 6\,{\mathrm{MeV}}/c^2 \,, \\ \Gamma_{X(4740)}&=&53 \pm 15 \pm 11\,{\mathrm{MeV}} \,, \end{eqnarray*} where the first uncertainty is statistical and the second is systematic. In addition, the most precise single measurement of the mass of the $B^0_s$ meson is performed and gives a value of $$ m_{B^0_s} = 5366.98 \pm 0.07 \pm 0.13\,{\mathrm{MeV}}/c^2\,. $$


Introduction
Decays of beauty hadrons to final states with charmonia provide a unique laboratory to study the properties of charmonia and charmonium-like states. A plethora of new states has been observed in such decays, including the χ c1 (3872) particle [1], pentaquark [2][3][4][5] and numerous tetraquark [5][6][7][8][9][10][11][12][13][14] candidates as well as conventional charmonium states, such as the tensor D-wave ψ 2 (3823) meson [15,16]. The nature of many exotic charmonium-like candidates remains unclear. A comparison of production rates with respect to those of conventional charmonium states in decays of beauty hadrons can shed light on their production mechanisms [17]. For example, the D * D rescattering mechanism [18,19] would give a large contribution to the χ c1 (3872) production and affect the pattern of decay rates of beauty hadrons. A modified pattern is also expected for a compact-tetraquark interpretation of the χ c1 (3872) state [20].
The decay chain B 0 s → (χ c1 (3872) → J/ψπ + π − ) (φ → K + K − ) is experimentally easiest to study in (quasi) two-body decays of a B 0 s meson with a χ c1 (3872) particle in the final state. This decay has recently been studied by the CMS collaboration, which found the ratio of branching fractions for the B 0 s → χ c1 (3872)φ and B 0 → χ c1 (3872)K 0 decays to be compatible with unity, and two times smaller than the ratio of branching fractions for B 0 s → χ c1 (3872)φ and B + → χ c1 (3872)K + decays [21]. The decay of the B 0 s meson into the J/ψπ + π − K + K − final state allows the mass spectrum of the J/ψφ system to be studied. Four tetraquark candidates have been observed by the LHCb collaboration using an amplitude analysis of B + → J/ψφK + decays [12]. These states are denoted by the PDG as χ c1 (4140), χ c1 (4274), χ c0 (4500) and χ c0 (4700) [22]. In B 0 s → J/ψπ + π − φ decays, the J/ψφ mass can be probed up to approximately 300 MeV/c 2 above the allowed kinematic limit in B + → J/ψφK + decays.
The B 0 s → J/ψφφ decay has been observed by the LHCb collaboration [23] using a data set collected in 2011-2012 LHC data taking. While the energy release in this decay is small, the measured branching fraction is large, possibly indicating a non-trivial decay dynamics that enhances the decay rate. A comparison of the rate of the B 0 s → J/ψφφ decay with that of B 0 s → J/ψη φ, B 0 s → J/ψη η and B 0 s → J/ψK * 0 K * 0 decay modes could clarify the dynamics: the first two modes have similar quark content while the third is formed by three vector mesons and results in the J/ψK + K − π + π − final state.
In this paper, a sample of B 0 s → J/ψπ + π − K + K − decays is analysed, with the J/ψ meson reconstructed in the µ + µ − final state. The study is based on proton-proton (pp) collision data, corresponding to integrated luminosities of 1, 2 and 6 fb −1 , collected with the LHCb detector at centre-of-mass energies of 7, 8 and 13 TeV, respectively. This data sample is used to measure the rates of the B 0 s → χ c1 (3872)φ, B 0 s → J/ψK * 0 K * 0 decays, where K * 0 denotes the K * (892) 0 resonance, and B 0 s → χ c1 (3872)K + K − decays, where the K + K − pair does not originate from a φ meson. The presence of B 0 s → (ψ(2S) → J/ψπ + π − ) (φ → K + K − ) decays in the same sample provides a convenient mode for normalising the observed rates of the different final states since the branching fraction of this decay is known [22]. This paper presents measurements of the following ratios of branching fractions (B), R J/ψK * 0 K * 0 ψ(2S)φ ≡ B B 0 s → J/ψK * 0 K * 0 × (B (K * 0 → K + π − )) 2 B (B 0 s → ψ(2S)φ) × B (ψ(2S) → J/ψπ + π − ) × B (φ → K + K − ) . (1c) The J/ψφ mass spectrum from B 0 s → J/ψπ + π − φ decays is investigated to search for resonant contributions. The large size of the analysed sample and the low level of background also allows for a precise determination of the mass of the B 0 s meson. The mass is measured using a subsample enriched in B 0 s → ψ(2S)φ decays, which have a small energy release.

Detector and simulation
The LHCb detector [24,25] 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 [26], 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 [27,28] placed downstream of the magnet.
The tracking system provides a measurement of the momentum of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The momentum scale is calibrated using samples of J/ψ → µ + µ − and B + → J/ψK + decays collected concurrently with the data sample used for this analysis [29,30]. The relative accuracy of this procedure is estimated to be 3 × 10 −4 using samples of other fully reconstructed b hadrons, Υ and K 0 S mesons. The minimum distance of a track to a primary pp-collision vertex (PV), the impact parameter (IP), 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 (RICH) [31]. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter [32]. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [33].
The online event selection is performed by a trigger [34], 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. The hardware trigger selects muon candidates with high transverse momentum or dimuon candidates with a high value of the product of the p T of the muons. In the software trigger two oppositely charged muons are required to form a good-quality vertex that is significantly displaced from every PV, with a dimuon mass exceeding 2.7 GeV/c 2 .
Simulated events are used to describe signal shapes and to compute the efficiencies needed to determine the branching fraction ratios. In the simulation, pp collisions are generated using Pythia [35] with a specific LHCb configuration [36]. Decays of unstable particles are described by the EvtGen package [37], in which final-state radiation is generated using Photos [38]. The χ c1 (3872) → J/ψπ + π − decays are simulated proceeding via an S-wave J/ψρ 0 intermediate state [39]. The model described in Refs. [40][41][42][43] is used to describe the ψ(2S) decays. The simulation is corrected to reproduce the transverse momentum and rapidity distributions of the B 0 s observed in data. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [44] as described in Ref. [45]. To account for imperfections in the simulation of charged-particle reconstruction, the track reconstruction efficiency determined from simulation is corrected using data-driven techniques [46].

Event selection
Candidate B 0 s → J/ψπ + π − K + K − decays are reconstructed using similar selection criteria to those used in Refs. [47][48][49]. Muon and hadron candidates are identified using combined information from the RICH, calorimeter and muon detectors [50]. They are required to have a transverse momentum larger than 550, 200 and 400 MeV/c for muon, pion and kaon candidates, respectively. To ensure that the particles can be efficiently separated by the RICH detectors, kaons and pions are required to have a momentum between 3.2 and 150 GeV/c. To reduce the combinatorial background due to particles produced promptly in the pp interaction, only tracks that are inconsistent with originating from a primary vertex are used. Pairs of oppositely charged muons consistent with originating from a common vertex are combined to form J/ψ candidates. The mass of the dimuon candidate is required to be between 3.05 and 3.15 GeV/c 2 .
Selected J/ψ candidates are combined with two oppositely charged kaons as well as two oppositely charged pions to form B 0 s → J/ψπ + π − K + K − candidates. A requirement on the quality of the common six-prong vertex is imposed. To improve the mass resolution for the B 0 s candidates, the mass of the µ + µ − pair is constrained to the known mass of the J/ψ meson [22] and the B 0 s candidate is constrained to originate from its associated PV. 1 Finally, the decay time of the B 0 s candidates is required to be between 0.2 and 2.0 mm/c. The lower limit is used to reduce background from particles coming from the PV while the upper limit suppresses poorly reconstructed candidates.
A possible background from Λ 0 b → J/ψπ + π − pK − and B 0 → J/ψπ + π − K + π − decays, with the proton or a pion misidentified as a kaon, is suppressed using a veto. After assigning the proton or pion mass to one of the kaons, only candidates outside the mass intervals 5.606 < m J/ψπ + π − pK − < 5.632 GeV/c 2 and 5.266 < m J/ψπ + π − K + π − < 5.288 GeV/c 2 are retained in the analysis. The mass distribution of the selected B 0 s → J/ψπ + π − K + K − candidates is shown in Fig. 1. The data are fit with the sum of a modified Gaussian function with power-law tails on both sides of the distribution [51,52] and a linear combinatorial background component. The B 0 s signal yield is (26.5 ± 0.2) × 10 3 candidates.
Only candidates with 0.995 < m K + K − < 1.060 GeV/c 2 and 5.30 < m J/ψπ + π − K + K − < 5.48 GeV/c 2 are considered. To improve the resolution on the J/ψπ + π − mass and to eliminate a small correlation between the m J/ψK + K − π + π − LHCb Figure 1: Distribution of the J/ψπ + π − K + K − mass of selected B 0 s candidates shown as points with error bars. A fit, described in the text, is overlaid. and m J/ψπ + π − variables, the m J/ψπ + π − variable is computed using a kinematic fit [53] that constrains the mass of the B 0 s candidate to its known value [22]. In each region, the three-dimensional fit model is defined as a sum of eight components. Four of these components correspond to decays of B 0 s mesons: 1. a signal B 0 s → X cc φ component, described by the product of B 0 s , X cc and φ signal templates, discussed in detail in the next paragraph; 2. a component corresponding to B 0 s → X cc K + K − decays, where the K + K − pair does not originate from a φ meson, parameterised by the product of B 0 s and X cc signal templates and a slowly varying template describing the non-resonant K + K − distribution, referred to below as the non-resonant K + K − function; 3. a component corresponding to B 0 s → J/ψπ + π − φ decays, parameterised as a product of the B 0 s and φ signal templates and a slowly varying template describing the nonresonant J/ψπ + π − mass distribution, referred to as the non-resonant J/ψπ + π − function hereafter; 4. a component corresponding to the decay B 0 s → J/ψπ + π − K + K − with no narrow resonance in either the J/ψπ + π − or the K + K − systems, described by the product of the B 0 s signal template and a slowly varying function f bkg m J/ψπ + π − , m K + K − , described below.
Four additional components correspond to random X cc φ, X cc K + K − , J/ψπ + π − φ and J/ψπ + π − K + K − combinations. Their parameterisation uses a second-order polynomial function in m J/ψK + K − π + π − , denoted as F B 0 s . These four background components are: 1. a component corresponding to random combinations of X cc and φ signals, parameterised as a product of the F B 0 s function and the X cc and φ signal templates; 2. a component corresponding to random combinations of an X cc signal with a non-resonant K + K − pair, parameterised as a product of the F B 0 s function, the signal X cc template and the non-resonant K + K − function; 3. a component corresponding to random combinations of a φ signal with a non-resonant J/ψπ + π − combination, parameterised as a product of the F B 0 s function, the signal φ template and the non-resonant J/ψπ + π − function; 4. a component corresponding to random J/ψπ + π − K + K − combinations parameterised as a product of the F B 0 s function and the function f bkg m J/ψπ + π − , m K + K − . In the m J/ψπ + π − K + K − distribution, the B 0 s signal shape is modelled with a modified Gaussian function with power-law tails on both sides of the distribution [51,52]. The tail parameters are fixed from simulation, while the mass of the B 0 s meson is allowed to vary. The detector resolution taken from simulation is corrected by a scale factor, s B 0 s , that accounts for a small discrepancy between data and simulation [16] and is allowed to vary. The φ and X cc signal templates are modelled with relativistic P-wave and S-wave Breit-Wigner functions, respectively, convolved with the detector resolution functions described below. Due to the proximity of the χ c1 (3872) state to the D 0 D * 0 threshold, modelling this component with a Breit-Wigner function may not be adequate [18,[54][55][56][57]. However, the analyses in Refs. [16,58] demonstrate that a good description of data is obtained with a Breit-Wigner line shape when the detector resolution is included. The mass of the ψ(2S) state is allowed to vary, while the width is fixed to its known value [22]. The width of the χ c1 (3872) state and the mass difference m ψ(2S) − m χ c1 (3872) are constrained to their known values [16,58] using Gaussian constraints. The detector resolution is described by a symmetric modified Gaussian function with power-law tails on both sides of the distribution [51,52], with all parameters determined from simulation. The resolution functions for the X cc templates are corrected by a common scale factor, s X cc , to account for a small discrepancy in the detector resolution between data and simulation [16,58]. This factor is determined from data. The non-resonant K + K − and J/ψπ + π − distributions are modelled by the product of a linear function and two-body, Φ 2,5 (m K + K − ), and three-body, Φ 3,5 m J/ψπ + π − phase-space distributions for five-body B 0 s decays [59]. 2 The function f bkg is parameterised by 2 The phase-space mass distribution of an l-body combination of particles from a n-body decay is approximated by Φ l,n (x) ∝ x , and x min , x max denote the minimal and maximal values of x, respectively, where P bkg is a polynomial function that is linear in one variable for each fixed value of the other variable.
The fit is performed simultaneously to the two J/ψπ + π − mass regions, with the B 0 s and X cc masses and the resolution scale factors, s B 0 s and s X cc , as shared parameters. The J/ψπ + π − K + K − , J/ψπ + π − and K + K − mass distributions together with projections of the simultaneous fit are shown in Figs. 2 and 3. The fit procedure is tested using a large sample of pseudoexperiments, generated using the nominal model with parameters extracted from data. Biases of O(1%) on the yields of different fit components are observed and the results are corrected for these biases. The corrected yields of the B 0 s → χ c1 (3872)φ and B 0 s → ψ(2S)φ decays and the resolution scale factors are listed in Table 1. The statistical significance for the B 0 s → χ c1 (3872)φ signal is calculated to be in excess of 10 standard deviations using Wilks' theorem [60]. Apart from the signal B 0 s → ψ(2S)φ component, only the B 0 s → ψ(2S)K + K − , and the combinatorial J/ψπ + π − φ and J/ψπ + π − K + K − components are found to contribute in a non-negligible way to the ψ(2S) mass region. In the χ c1 (3872) region, the contribution from the B 0 s → χ c1 (3872)K + K − component is found to be small and the combinatorial components χ c1 (3872)φ and χ c1 (3872)K + K − are negligible. The resolution scale factors, s B 0 s and s X cc , are similar to those obtained in Refs. [16,58].
The results are cross-checked using a two-dimensional unbinned extended maximum-likelihood fit to the background-subtracted J/ψπ + π − and K + K − mass distributions, where the sPlot technique [61] is used with the J/ψπ + π − K + K − mass as a discriminating variable. The results of this fit are found to be in very good agreement with the results listed in Table 1.
The ratio of branching fractions defined in Eq. (1a) is calculated from where the signal yields, N B 0 s →χ c1 (3872)φ and N B 0 s →ψ(2S)φ , are taken from Table 1 and ε B 0 s →χ c1 (3872)φ and ε B 0 s →ψ(2S)φ are the efficiencies to reconstruct and select the B 0 s → χ c1 (3872)φ and B 0 s → ψ(2S)φ decays, respectively. The efficiencies are defined as the product of the detector geometric acceptance and the reconstruction, selection, hadron identification and trigger efficiencies. All of the efficiency contributions, except the hadron-identification efficiency, are determined using simulated samples. The hadron-identification efficiency is determined using large calibration samples of D * + → (D 0 → K − π + ) π + , K 0 S → π + π − and D + s → (φ → K + K − ) π + decays selected in data [31,62]. The efficiency ratio is found to be 0.66 ± 0.01, where the uncertainty is only  that due to the size of the simulated samples. The efficiency ratio differs from unity due to the harder p T spectrum of pions in the B 0 s → χ c1 (3872)φ decays. The resulting value of R where the uncertainty is statistical.
The decay B 0 s → χ c1 (3872)K + K − , where the K + K − pair does not originate from a φ meson, is studied using a sample of selected B 0 s → J/ψπ + π − K + K − candidates with the J/ψπ + π − and J/ψπ + π − K + K − masses in the ranges 3.85 < m J/ψπ + π − < 3.90 GeV/c 2 and 5.30 < m J/ψπ + π − K + K − < 5.48 GeV/c 2 . A two-dimensional unbinned extended maximum-likelihood fit is performed to the J/ψK + K − π + π − and J/ψπ + π − mass distributions. The fit function comprises the sum of four components:  2. a component corresponding to B 0 s → J/ψπ + π − K + K − decays, parameterised as a product of the B 0 s signal template and the non-resonant J/ψπ + π − function; 3. a component corresponding to random combinations of χ c1 (3872) particles with a K + K − pair, parameterised as a product of the χ c1 (3872) signal template and the F B 0 s function; 4. a component corresponding to random J/ψπ + π + K + K − combinations, parameterised as a product of the three-body phase-space function Φ 3,5 m J/ψπ + π − and a two-dimensional non-factorisable bilinear function.
The J/ψπ + π − K + K − and J/ψπ + π − mass distributions together with projections of the fit are shown in Fig. 4. The yield of B 0 s → χ c1 (3872)K + K − signal decays is  Figure 4: Distributions of the (left) J/ψπ + π − K + K − and (right) J/ψπ + π − mass of selected B 0 s → χ c1 (3872)K + K − candidates shown as points with error bars. A fit, described in the text, is overlaid.
which significantly exceeds the yield of N B 0 s →χ c1 (3782)φ shown in Table 1, pointing to a sizeable contribution from the B 0 s → χ c1 (3872)K + K − decays, where the K + K − pair does not originate from a φ meson.
The fraction of B 0 s → χ c1 (3872) (φ → K + K − ) decays is estimated using an unbinned maximum-likelihood fit to the background-subtracted K + K − mass distribution from signal B 0 s → χ c1 (3872)K + K − decays. The background-subtracted K + K − mass distribution is obtained by applying the sPlot technique [61] to the results of the two-dimensional fit to the B 0 s → χ c1 (3872)K + K − decays described above. The background-subtracted K + K − mass distribution is further corrected for the K + K − mass-dependent efficiency by applying a weight, to each candidate. The efficiencies ε B 0 s →χ c1 (3872)φ and ε B 0 s →χ c1 (3872)K + K − are calculated using simulated samples, where a phase-space decay model is used for the three-body B 0 s → χ c1 (3872)K + K − decays. The background-subtracted and efficiency-corrected K + K − mass distribution of the B 0 s → χ c1 (3872)K + K − candidates is shown in Fig. 5. In addition to a clear narrow structure, corresponding to B 0 s → χ c1 (3872) (φ → K + K − ) decays, a sizeable number of B 0 s → χ c1 (3872)K + K − decays, where the K + K − pair does not originate from the φ meson is visible. The K + K − mass distribution for m K + K − > 1.1 GeV/c 2 cannot be described by phase-space, and possibly contains contributions from the f 0 (980), f 2 (1270), f 0 (1370) and f 2 (1525) resonances decaying to a pair of kaons, as has been observed in B 0 s → J/ψK + K − decays [63,64]. An amplitude analysis of a larger data sample would be required to properly disentangle individual contributions. However, a narrow φ component can be separated from the non-φ components using an unbinned maximum-likelihood fit to the background-subtracted and efficiency-corrected K + K − mass distribution. The fit function comprises two components  the φ signal template (see Sec. 4) multiplied by the phase-space function Φ 2,3 (m K + K − ) for the three-body B 0 s → χ c1 (3872)K + K − decay; 2. a component that accounts for non-resonant B 0 s → χ c1 (3872)K + K − decays and decays via broad high-mass K + K − intermediate states, modelled by a product of a phase-space function Φ 2,3 (m K + K − ) for three-body B 0 s → χ c1 (3872)K + K − decays and a third-order polynomial function.
The shape of the second component is flexible enough to accommodate contributions from wide K + K − resonances. The projection of the fit is overlaid in Fig. 5. The fraction of the φ → K + K − signal component is found to be f φ = (38.9 ± 4.9) % .
This fraction is converted into the ratio of branching fractions R K + K − , defined in Eq. (1c), 5366.79 ± 0.06 where the uncertainty is statistical. This is the first observation of the decay B 0 s → χ c1 (3872)K + K − , where K + K − pair does not originate from a φ meson.
The fit model is similar to that used in Sec. 4 but with some modifications. First, the model is symmetric with respect to an interchange of K + π − and K − π + pairs. Second, for components that account for K * 0 K * 0 , K * 0 K − π + or K + π − K * 0 combinations, corrections are applied due to the limited phase space available in the decays. These shapes are derived from fits to simulated samples and comprise symmetric products of phase-space functions and linear polynomials. The K * 0 (K * 0 ) signal is parameterised by a relativistic P-wave Breit-Wigner function. The width of the K * 0 meson, 47.3 ± 0.5 MeV, is not small [22] and the fit ranges are wide, hence the correct determination of all components would require a full amplitude analysis that properly accounts for interference effects. Such an analysis is beyond the scope of this paper. However, fits to simulated samples of B 0 s → J/ψπ + π − K + K − decays with different compositions of intermediate states show that the simple model described here allows for a reliable determination of the B 0 s → J/ψK * 0 K * 0 component. The J/ψπ + π − K + K − , K + π − and K − π + mass distributions together with projections of the fit are shown in Fig. 6 and the parameters of interest are summarized in Table 2. A study of a large sample of pseudoexperiments generated and fitted with the nominal model, indicates a small bias of O(1%) on the signal yield. The quoted yield is corrected for this bias.

B 0 s mass measurement
The precision on the B 0 s mass value, reported in Table 1, is improved by imposing a constraint on the reconstructed mass of the ψ(2S) state [53]. Applying this constraint improves the B 0 s mass resolution and significantly decreases systematic uncertainties on the mass measurement, since the mass of the ψ(2S) meson is known with high precision [65]. The mass of the B 0 s meson is determined from an unbinned extended maximum-likelihood fit to the ψ(2S)K + K − mass distribution for a sample of B 0 s → J/ψK + K − π + π − decays with m K + K − < 1.06 GeV/c 2 and with the J/ψπ + π − mass within a narrow region around the known mass of the ψ(2S) meson, 3.679 < m J/ψπ + π − < 3.694 GeV/c 2 .
The ψ(2S)K + K − mass distribution is fitted with a two-component function comprising a signal component modelled with the B 0 s signal template and a background component modelled with a second-order polynomial function. The ψ(2S)K + K − mass distribution together with the fit results is shown in Fig. 7. The fit results are summarized in Table 3. Studies of simulated samples show that the selection requirements introduce a small bias in the measured mass of long-lived heavy-flavour hadrons [66][67][68]. The corrected value for  Figure 7: Distribution of the ψ(2S)K + K − mass for selected B 0 s → J/ψK + K − π + π − candidates, enriched in B 0 s → (ψ(2S) → J/ψπ + π − ) (φ → K + K − ) decays (points with error bars). A fit, described in the text, is overlaid.
where the uncertainty is statistical only.

J/ψφ mass spectrum
The J/ψφ mass spectrum in B 0 s → J/ψπ + π − φ decays is studied using a sample of selected B 0 s → J/ψπ + π − K + K − candidates with the K + K − mass in the range 0.995 < m K + K − < 1.06 GeV/c 2 and excluding the J/ψπ + π − mass regions around the narrow ψ(2S) and χ c1 (3872) states, i.e. 3.672 < m J/ψπ + π − < 3.700 GeV/c 2 and 3.864 < m J/ψπ + π − < 3.880 GeV/c 2 . A two-dimensional unbinned extended maximum-likelihood is performed to the J/ψπ + π − K + K − and K + K − mass distributions. The fit function comprises a sum of four components: 1. a component corresponding to B 0 s → J/ψπ + π − φ decays, parameterised by the product of the B 0 s and φ signal templates described in Sec. 4; 2. a component corresponding to B 0 s → J/ψπ + π − K + K − decays, parameterised by the product of the B 0 s signal template and the non-resonant K + K − function; 3. a component corresponding to random J/ψπ + π − φ combinations, parameterised by the product of the φ signal template and the F B 0 s function; 4. a component describing random J/ψπ + π − K + K − combinations, parameterised by the product of the phase-space function Φ 2,5 (m K + K − ) and the two-dimensional non-factorisable bilinear function described in Sec. 4.
The J/ψπ + π − K + K − and K + K − mass spectra together with the projections of the fit are shown in Fig. 8. The sPlot technique is applied to obtain a background-subtracted J/ψφ mass distribution of B 0 s → J/ψπ + π − φ decays. The resulting distribution is shown in Fig. 9 (left). It shows a prominent structure at a mass around 4.74 GeV/c 2 . No such structure is seen if the K + K − mass is restricted to the region of 1.06 < m K + K − < 1.15 GeV/c 2 . This structure cannot be explained by B 0 s → X cc φ decays via a narrow intermediate X cc resonance since contributions from B 0 s → ψ(2S)φ and B 0 s → χ c1 (3872)φ decays are explicitly vetoed. If no veto is applied, B 0 s → ψ(2S)φ decays would produce a broad structure in the J/ψφ mass spectrum that peaks around 4.76 GeV/c 2 and has a width that is approximately twice that of the observed structure. Studies with simulated samples indicate that after the veto is applied the remaining contributions from these decays are totally negligible. No sizeable contributions from decays via other narrow charmonium states are observed in the background-subtracted J/ψπ + π − mass spectrum. The background-subtracted π + π − mass distribution of candidates in the mass range 4.68 < m J/ψφ < 4.78 GeV/c 2 is found to have no structure. The background-subtracted φπ + π − mass spectrum of the B 0 s → J/ψπ + π − φ decays is shown in Fig. 9 (right). The spectrum exhibits significant deviations from the phase-space distribution, indicating possible presence of excited φ states, referred to as φ * states hereafter. The decays B 0 s → J/ψφ * via intermediate φ(1680), φ(1850) or φ(2170) states [22] are studied using simulated samples. It is found that the J/ψφ mass spectra from B 0 s → J/ψφ * decays exhibit no structure and for the J/ψφ mass exceeding 4.4 GeV/c 2 can be described by a monotonically decreasing function. If the intervals used to reject the B 0 s → ψ(2S)φ and B 0 s → χ c1 (3872)φ decays are significantly increased, in excess of 60 MeV/c 2 , it is possible to generate two wide regions with decreased yields around 4.65 GeV/c 2 and 4.82 GeV/c 2 in the J/ψφ mass spectrum. The positions and shapes of these dips depend on the assumed mass and width of the φ * state and for certain choices of the φ * states, two dips in the monotonically decreased spectrum could sculpt a bump. The complicated interference between several decay chains, including different intermediate φ * states, could result in a distorted J/ψφ mass spectrum. In order to ascertain if the structure at 4.74 GeV/c 2 , seen in  Figure 8: Distribution of the (left) J/ψπ + π − K + K − and (right) K + K − mass for selected B 0 s → J/ψπ + π − φ candidates. A fit, described in the text, is overlaid. : Background-subtracted (left) J/ψφ and (right) φπ + π − mass distributions from B 0 s → J/ψπ + π − φ decays (points with error bars). The expectation from simulated B 0 s → J/ψπ + π − φ decays is overlaid (green solid line). In the right figure, the backgroundsubtracted φπ + π − mass distribution in the region 4.68 < m J/ψφ < 4.78 GeV/c 2 is shown (red open circles with error bars) together with the corresponding expectation from simulated B 0 s → J/ψπ + π − φ decays (blue dashed line). Fig. 9 (left), is resonant and not due to interference an amplitude analysis, similar to that in Refs. [2,9,11,12] would be required. Such an analysis is beyond the scope of this paper. Under the assumption that this structure, referred to as X(4740) hereafter, has a resonant nature, its mass and width are determined through an unbinned extended maximum-likelihood fit to the background-subtracted J/ψφ mass distribution in the range 4.45 < m J/ψφ < 4.90 GeV/c 2 . The fit function comprises two components: GeV/c 2 Figure 10: Background-subtracted J/ψφ mass distribution from B 0 s → J/ψπ + π − φ decays (points with error bars). A fit, described in the text, is overlaid.
2. a component, corresponding to B 0 s → J/ψπ + π − φ decays, parameterised by the product of the Φ 2,4 m J/ψφ function and a third-order polynomial function.
The background-subtracted J/ψφ mass spectrum with superimposed results of the fit is shown in Fig. 10 and the results are listed in Table 4.
The statistical significance of the observed structure is estimated using Wilks' theorem [60] and found to be 5.5 standard deviations. The significance estimate is validated using a large number of pseudoexperiments comprising no X(4740) signal component.

Systematic uncertainties
Due to the similar decay topologies, systematic uncertainties largely cancel in the ratios R.
The remaining contributions to systematic uncertainties are summarized in Table 5 and discussed below. The largest source of systematic uncertainty on the ratios arise from imperfect knowledge of the shapes of signal and background components used in the fits. To estimate this uncertainty, several alternative models for the signal, non-resonant signal and background components are tested. For the B 0 s signal shape and the detector resolution functions in the X cc signal templates, the bifurcated Student's t-distribution is tested as an alternative model. For the Breit-Wigner functions describing the φ and K * 0 signal shapes, the meson radii in the Blatt-Weisskopf barrier factors [69] are varied between 1.5 and 5 GeV −1 . The mass and width of the K * 0 meson are varied within their uncertainties [22]. The degree of the polynomials used in the non-resonant J/ψπ + π − and K + K − functions, the F B 0 s and P bkg functions and all other polynomial functions used in the fits are increased by one. The largest systematic uncertainty for the ratio R K + K − is associated with the parameterisation of the fit component for the B 0 s → χ c1 (3872)K + K − decays, where the K + K − pair does not originate from a φ meson. The explicit inclusion of a B 0 s → χ c1 (3872)f 0 (980) component is considered, where the f 0 (980) state decays into a K + K − pair. The f 0 (980) line shape is modelled by a Flatté function [70], with parameters taken from Refs. [71,72]. The systematic uncertainty on the ratio R J/ψK * 0 K * 0 ψ(2S)φ due to the fit range for K ± π ∓ masses is studied by increasing this range to 0.63 < m K ± π ∓ < 1.25 GeV/c 2 . For each alternative model the ratio of event yields is remeasured, and the maximum deviation with respect to the nominal model, 1.8%, 2.6% and 7.3% for the R and R K + K − ratios, respectively, is assigned as a systematic uncertainty.
An additional systematic uncertainty on the ratios arises due to differences between data and simulation. In particular, there are differences in the reconstruction efficiency of charged-particle tracks that do not cancel completely in the ratio due to the different kinematic distributions of the final-state particles. The track-finding efficiencies obtained from simulated samples are corrected using data-driven techniques [46]. The uncertainties related to the efficiency correction factors, together with the uncertainty on the hadron-identification efficiency due to the finite size of the calibration samples [31,62], are propagated to the ratios of the total efficiencies using pseudoexperiments and account for 0.3%, 0.1% and 0.3% for the R χ c1 (3872)φ ψ(2S)φ , R J/ψK * 0 K * 0 ψ(2S)φ and R K + K − ratios, respectively. A systematic uncertainty on the ratios related to the knowledge of the trigger efficiencies is estimated by comparing the ratios of trigger efficiencies in data and simulation for   [73] and is taken to be 1.1% for all three ratios of branching fractions. Other data-simulation differences are investigated by varying selection criteria in data. The resulting variations in the efficiency ratios do not exceed 2%, which is taken as the corresponding systematic uncertainty.
The final systematic uncertainty considered on the ratios of branching fractions is due to the knowledge of the ratios of efficiencies in Eqs. (3), (6) and (9) due to the finite size of the simulated samples. It is determined to be 1.0%, 0.9% and 1.3% for the R and R K + K − ratios, respectively. No systematic uncertainty is included for the admixture of the CP -odd and CP -even B 0 s eigenstates in the decays, which is assumed to be the same for all four involved channels [74]. In the extreme case that one decay is only from the short-lifetime eigenstate and the other from the long-lifetime eigenstate, the corresponding ratio of branching fractions would change by 3.8%.
The systematic uncertainties on the B 0 s mass measurement are summarised in Table 6. The most important source of systematic uncertainty is related to the momentum scale calibration in data. This effect is evaluated by varying the scale within its known uncertainty [30]. The resulting change in the mass of 122 keV/c 2 is assigned as a systematic uncertainty. Other sources of uncertainty are related to energy loss corrections and the imprecise knowledge of the K ± and ψ(2S) meson masses [22]. The amount of material traversed in the tracking system by a particle is known to 10% accuracy, which leads to an uncertainty on the estimated energy loss of particles in the detector. This systematic uncertainty is calculated in Ref. [30] to be 15 keV/c 2 . The uncertainties on the known kaon and ψ(2S) masses [22,65] are propagated to the uncertainty in the B 0 s mass using simulated samples and are found to be 27 and 10 keV/c 2 , respectively. Using the ψ(2S) mass constraint significantly reduces the systematic uncertainties associated with the momentum scale and energy loss correction. The B 0 s → J/ψK * 0 K * 0 signal sample has a smaller statistical uncertainty, see Table 2, however, the systematic uncertainties due to momentum scaling and energy loss are twice as large, making this sample non-competitive for a precise measurement of the B 0 s mass. Systematic uncertainties on the mass and width of the X(4740) structure are summarised in Table 7. The uncertainty related to the imperfect knowledge of the signal and background shapes is estimated using alternative fit models. Relativistic P-and D-wave Breit-Wigner functions are used as alternative shapes for signal with the value of the Blatt-Weisskopf barrier factor meson radius varied between 1.5 and 5 GeV −1 .
A model comprising a product of an S-wave Breit-Wigner function with a phase-space function that accounts for the proximity of the upper edge for the J/ψφ mass spectrum from B 0 s → J/ψπ + π − φ decays is also used. For the background component, a convex monotonically decreasing third-order polynomial function and a product of the Φ 2,4 m J/ψφ function with a third-order polynomial function are tested as alternative models. The maximal deviations with respect to the baseline model of 2.8 MeV/c 2 and 8.4 MeV for the mass and width of the X(4740) state, respectively, are taken as corresponding systematic uncertainties. The contributions from B 0 s → ψ(2S)K + K − and B 0 s → χ c1 (3872)K + K − decays are explicitly suppressed in the analysis by excluding the mass regions 3.672 < m J/ψπ + π − < 3.700 GeV/c 2 and 3.864 < m J/ψπ + π − < 3.880 GeV/c 2 around the known masses of the ψ(2S) and χ c1 (3872) states [16,22,58,65]. Repeating the analysis using wider exclusion ranges, causes changes of 4.6 MeV/c 2 and 5.1 MeV in the mass and width of the X(4740) structure, respectively. These changes are taken as systematic uncertainties due to possible remaining contributions from B 0 s → ψ(2S)K + K − and B 0 s → χ c1 (3872)K + K − decays. Large interference effects between the signal and coherent part of the background can also distort the visible shape of the resonance. To probe the importance of this effect, the signal fit component F S is modelled with a coherent sum of an S-wave Breit-Wigner amplitude A m J/ψφ and a coherent background where the positive linear polynomial b(m J/ψφ ) stands for the magnitude of the coherent background amplitude and ϕ denotes the phase of the coherent background, chosen to be independent of the J/ψφ mass. The deviations of the mass and width of the X(4740) structure obtained from this fit are taken as systematic uncertainties related to neglecting possible interference effects between the signal and the coherent part of the background. The complicated interference pattern for the B 0 s → J/ψφ * decays via different φ * states also can distort the J/ψφ mass spectrum. However, to quantify this effect a full amplitude analysis, similar to Refs. [2,9,11,12] is needed, that is beyond the scope of this paper, and no systematic uncertainty is assigned. Other sources of systematic uncertainties on the mass and width of the X(4740) structure, namely the momentum scale and the background-subtraction procedure are found to be negligible with respect to the leading systematic uncertainties related to the fit model. For each choice of the fit Sum in quadrature 5.5 11.1 model, the statistical significance of the observed X(4740) structure is calculated from data using Wilks' theorem [60]. The smallest significance found is 5.3 standard deviations, taken as its significance including systematic uncertainties.

Summary
A study of B 0 s → J/ψπ + π − K + K − decays is made using pp collision data corresponding to an integrated luminosity of 1, 2 and 6 fb where the first uncertainty is statistical and the second systematic. The ratio R χ c1 (3872)φ ψ(2S)φ is consistent with but more precise than the value of (2.21 ± 0.29 ± 0.17) × 10 −2 recently reported by the CMS collaboration [21]. The decays B 0 s → J/ψK * 0 K * 0 and B 0 s → χ c1 (3872)K + K − , where the K + K − pair does not originate from a φ meson, are observed for the first time. A full amplitude analysis, similar to Refs. [63,64], is needed to resolve possible contributions from two-body decays via K + K − resonances, like B 0 s → χ c1 (3872)f 0 (980) and B 0 s → χ c1 (3872)f 2 (1525), that in turn could be useful for a better understanding of the nature of the χ c1 (3872) state.
A precise measurement of the B 0 s mass is performed using a sample of selected B 0 s → J/ψπ + π − K + K − candidates enriched in B 0 s → ψ(2S)φ decays. The mass of the B 0 s meson is determined to be m B 0 s = 5366.98 ± 0.07 ± 0.13 MeV/c 2 , which is the most precise single measurement of this observable. This result is combined with other precise measurements by the LHCb collaboration using B 0 s → J/ψφ [23], B 0 s → J/ψφφ [75], B 0 s → χ c2 K + K − [76] and B 0 s → J/ψpp [77] decays. The combined mass is calculated using the best linear unbiased estimator [78], accounting for correlations of systematic uncertainties between the measurements. The LHCb average for the mass of the B 0 s meson is found to be m LHCb B 0 s = 5366.94 ± 0.08 ± 0.09 MeV/c 2 .
A dedicated analysis using a larger data set is needed to resolve if this state is different from the χ c0 (4700) state, observed in the B + → J/ψφK + decays [11,12].