Observation of the decay B0s → η c ϕ and evidence for B 0s → η c π + π −

JHEP07(2017)021 Published for SISSA by Springer Received:February 28, 2017 Revised:May 3, 2017 Accepted:June 17, 2017 Published:July 5, 2017 Observation of the decay B0 s→ηcφand evidence for B0 s→ηcπ+π− The LHCb collaboration E-mail: [email protected] Abstract: A study of B0 s→ηcφand B0 s→ηcπ+π−decays is performed using pp collision data corresponding to an integrated luminosity of 3.0 fb−1, collected with the LHCb detector in Run 1 of the LHC. The observation of the decay B0 s→ηcφis reported, where the ηcmeson is reconstructed in the p¯p,K+K−π+π−,π+π−π+π−and K+K−K+K−decay modes and the φ(1020) in the K+K−decay mode. The decay B0 s→J/ψφ is used as a normalisation channel. Evidence is also reported for the decay B0 s→ηcπ+π−, where the ηcmeson is reconstructed in the p¯pdecay mode, using the decay B0 s→J/ψπ+π−as a normalisation channel. The measured branching fractions are B(B0 s→ηcφ) = (5.01 ±0.53 ±0.27 ±0.63) ×10−4, B(B0 s→ηcπ+π−) = (1.76 ±0.59 ±0.12 ±0.29) ×10−4, where in each case the first uncertainty is statistical, the second systematic and the third uncertainty is due to the limited knowledge of the external branching fractions. Keywords: B physics, Branching fraction, Flavor physics, Hadron-Hadron scattering (experiments) ArXiv ePrint: 1702.08048 Open Access, Copyright CERN, for the benefit of the LHCb Collaboration. Article funded by SCOAP3. https://doi.org/10.1007/JHEP07(2017)021 JHEP07(2017)021 Contents 1 Introduction 1 2 Detector and simulation 2 3 Analysis strategy 3 4 Event selection 4 5 Efficiency correction 6 6 Fit models 7 6.1 Model for B0 s→ηcπ+π−decays 7 6.2 Model for B0 s→ηcφdecays 8 7 Results 11 8 Systematic uncertainties 15 9 Conclusions 16 A Fit projections 17 B Correlation matrix 17 The LHCb collaboration 27 1 Introduction When a B0 smeson decays through the ¯ b→¯cc¯sprocess, interference between the direct decay amplitude, and the amplitude after B0 s−B0 soscillation, gives rise to a CP-violating phase, φs. This phase is well predicted within the Standard Model (SM) [1] and is sensitive to possible contributions from physics beyond the SM [2–5]. The φsphase is best measured using the “golden” channel1B0 s→J/ψφ [6–10] and the precision of this measurement is expected to be dominated by its statistical uncertainty until the end of LHC running. In addition to B0 s→J/ψφ, other modes have been used to constrain φs:B0 s→J/ψπ+π−[6], B0 s→D+ sD− s[11], and B0 s→ψ(2S)φ[12]. In this paper, the first study of B0 s→ηcφand B0 s→ηcπ+π−decays is presented.2 These decays also proceed dominantly through a ¯ b→¯cc¯stree diagram as shown in figure 1. Unlike in B0 s→J/ψφ decays, the ηcφfinal state is purely CP-even, so that no angular analysis is required to measure the mixing phase φs. However, the size of the data sample recorded by the LHCb experiment in LHC Run 1 is not sufficient to perform time-dependent 1The simplified notation φand ηcare used to refer to the φ(1020) and the ηc(1S) mesons throughout this article. 2The use of charge-conjugate modes is implied throughout this article. – 1 – JHEP07(2017)021 B0 s ηc(1S) φ(1020), f0(980) W ¯ b s s ¯c c ¯s Figure 1. Leading diagram corresponding to B0 s→ηcφand B0 s→ηcπ+π−decays, where the π+π−pair may arise from the decay of the f0(980) resonance. analyses of B0 s→ηcφand B0 s→ηcπ+π−decays. Instead, the first measurement of their branching fractions is performed. No prediction is available for either B(B0 s→ηcφ) or B(B0 s→ηcπ+π−). Assuming B(B0 s→ηcφ) B(B0 s→J/ψφ)=B(B0→ηcK0) B(B0→J/ψK0)=B(B0 s→ηcπ+π−) B(B0 s→J/ψπ+π−)(1.1) allows B(B0 s→ηcφ) and B(B0 s→ηcπ+π−) to be estimated. From the known values of B(B0→ηcK0), B(B0→J/ψK0), B(B0 s→J/ψφ) and B(B0 s→J/ψπ+π−) [13], one finds B(B0 s→ηcφ) = O(10−3),(1.2) B(B0 s→ηcπ+π−) = O(10−4).(1.3) The measurements presented in this paper are performed using a dataset corresponding to 3 fb−1of integrated luminosity collected by the LHCb experiment in pp collisions during 2011 and 2012 at centre-of-mass energies of 7 TeV and 8 TeV, respectively. The paper is organised as follows: section 2describes the LHCb detector and the procedure used to generate simulated events; an overview of the strategy for the measurements of B(B0 s→ ηcφ) and B(B0 s→ηcπ+π−) is given in section 3; the selection of candidate signal decays is described in section 4; the methods to determine the reconstruction and selection efficiencies are discussed in section 5. Section 6describes the fit models. The results and associated systematic uncertainties are discussed in sections 7and 8. Finally, conclusions are presented in section 9. 2 Detector and simulation The LHCb detector [14,15] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing bor cquarks. 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 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 (IP), is measured with – 2 – JHEP07(2017)021 a resolution of (15 + 29/pT)µm, where pTis 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 [16], 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. Samples of simulated events are used to determine the effects of the detector geometry, trigger, and selection criteria on the invariant-mass distributions of interest for this paper. In the simulation, pp collisions are generated using Pythia [17,18] with a specific LHCb configuration [19]. The decay of the B0 smeson is described by EvtGen [20], which generates final-state radiation using Photos [21]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [22,23] as described in ref. [24]. Data-driven corrections are applied to the simulation to account for the small level of mismodelling of the particle identification (PID) performance [25]. In the simulation the reconstructed momentum of every track is smeared by a small amount in order to better match the mass resolution of the data. 3 Analysis strategy In the analysis of B0 s→ηcφdecays, the φmeson is reconstructed in the K+K−final state and the ηcmeson is reconstructed in the p¯p,K+K−π+π−,π+π−π+π−and K+K−K+K− final states. For clarity, the three four-body final states are referred to as 4hthroughout the paper. In determining the branching fraction, the decay B0 s→J/ψφ is used as a normalisation channel, where the J/ψ meson is reconstructed in the same decay modes as the ηcmeson. A similar strategy is adopted for the measurement of the branching fraction of B0 s→ηcπ+π−decays. However, due to the higher expected level of combinatorial background compared to B0 s→ηcφdecays, the ηcand J/ψ mesons are reconstructed only in the p¯pfinal state in the measurement of B(B0 s→ηcπ+π−). In both analyses, a two-stage fit procedure is performed. In the first stage, unbinned extended maximum likelihood (UML) fits are performed to separate signal candidates from background contributions. For the B0 s→ηc(→p¯p)π+π−decay the fit is done to the p¯pπ+π−mass distribution, while for the decays B0 s→ηc(→p¯p)φ(→K+K−) and B0 s→ηc(→4h)φ(→K+K−) it is made to the two-dimensional p¯pK+K−versus K+K− or 4hK+K−versus K+K−mass distributions, respectively. The likelihood function is L(N,a) = e−PjNj n! n Y l=1  X j NjPj(m;a) ,(3.1) where jstands for the event species, Njis the corresponding yield and Nis the vector of yields Nj,ais the vector of fitted parameters other than yields, nis the total number – 3 – JHEP07(2017)021 of candidates in the sample, and Pj(m) is the probability density function (PDF) used to parametrise the set of invariant-mass distributions mconsidered. The RooFit package [26] is used to construct the negative log-likelihood function (NLL), which is minimised using Minuit [27]. Using information from these fits, signal weights for each candidate, ωl, are obtained using the sPlot technique [28]. In the second stage, for B0 s→p¯pπ+π−candidates a weighted UML fit is made to the p¯pinvariant-mass spectrum, and weighted UML fits of the p¯pand the 4hinvariant-mass spectra are done for B0 s→p¯pφ and B0 s→4hφ candidates, respectively, to disentangle ηc and J/ψ candidates from nonresonant (NR) and remaining background contributions, as described in section 6. For the weighted fits, the NLL function is given by −lnL(N,a) = ζX j Nj−ζX l ωlln  X j NjPj(m;a) + ln(n!),(3.2) where ζ=Plωl/Plω2 lensures proper uncertainty estimates from the weighted likelihood fit [29]. For the observed numbers of ηcand J/ψ candidates in final state f,Nηc,f and NJ/ψ ,f , the measured branching fraction is B(B0 s→ηcX) = Nηc,f NJ/ψ ,f ×B(B0 s→J/ψX)×B(J/ψ →f) B(ηc→f)×ε(J/ψ)f ε(ηc)f ,(3.3) where Xrefers to either the φmeson or the π+π−pair. The branching fractions B(B0 s→ J/ψφ), B(B0 s→J/ψπ+π−), B(J/ψ →f) and B(ηc→f) are taken from ref. [13], and the efficiency correction factors, ε, are obtained from simulation. In order to maximise the sensitivity to B(B0 s→ηcφ), a simultaneous fit to the p¯pand 4hinvariant-mass spectra is performed. 4 Event selection A common strategy for the event selection, comprising several stages, is adopted for all final states. First, online requirements are applied at the trigger level, followed by an initial offline selection in which relatively loose criteria are applied. Boosted decision trees (BDTs) [30], implemented using the TMVA software package [31], are then used to further suppress the combinatorial background arising from random combinations of tracks originating from any PV. Finally, the requirements on the output of the BDTs and on the PID variables are simultaneously optimised for each final state, to maximise the statistical significance of the signal yields. At the hardware trigger stage, events are required to have a muon with high pTor a hadron with high transverse energy in the calorimeters. The software trigger requires a two-, three- or four-track secondary vertex (SV) with a significant displacement from any PV. At least one charged particle must have a large transverse momentum and be inconsistent with originating from a PV. A multivariate algorithm [32] is used for the identification of secondary vertices consistent with the decay of a bhadron into charged hadrons. In addition, for the 4hfinal states, an algorithm is used to identify inclusive – 4 – JHEP07(2017)021 φ→K+K−production at a secondary vertex, without requiring a decay consistent with abhadron. In the initial stage of the offline selection, candidates for B0 s→p¯pπ+π−and B0 s→p¯pK+K−(B0 s→4hK+K−) decays are required to have four (six) good quality, high-pTtracks consistent with coming from a vertex that is displaced from any PV in the event. Loose PID criteria are applied, requiring the tracks to be consistent with the types of hadrons corresponding to the respective final states. In addition, the B0 scandidates, formed by the combination of the final-state candidates, are required to originate from a PV by requiring a small angle between the B0 scandidate momentum vector and the vector joining this PV and the B0 sdecay vertex, and a small χ2 IP, which is defined as the difference in the vertex-fit χ2of the considered PV reconstructed with and without the candidate. When forming the B0 scandidates for B0 s→p¯pπ+π−and B0 s→p¯pK+K−decays, the p¯p mass resolution is improved by performing a kinematic fit [33] in which the B0 scandidate is constrained to originate from its associated PV (that with the smallest value of χ2 IP for the B0 s), and its reconstructed invariant mass is constrained to be equal to the known value of the B0 smass [13]. No significant improvement of the 4hmass resolution is observed for B0 s→4hK+K−decays. In order to reduce the combinatorial background, a first BDT, based on kinematic and topological properties of the reconstructed tracks and candidates, is applied directly at the initial stage of the offline selection of candidate B0 s→4hK+K− decays. It is trained with events from dedicated simulation samples as signal and data from the reconstructed high-mass sidebands of the B0 scandidates as background. In the second step of the selection, the offline BDTs are applied. They are trained using the same strategy as that used for the training of the first BDT. The maximum distance of closest approach between final-state particles, the transverse momentum, and the χ2 IP of each reconstructed track, as well as the vertex-fit χ2per degree of freedom, the χ2 IP, and the pointing angle of the B0 scandidates are used as input to the BDT classifiers used to select candidate B0 s→p¯pπ+π−and B0 s→p¯pK+K−decays. For the p¯pK+K−final state, the direction angle, the flight distance significance and the χ2 IP of the reconstructed B0 s candidate are also used as input to the BDT, while the pTof the B0 scandidate is used for the p¯pπ+π−final state. The difference in the choice of input variables for the p¯pK+K− and the p¯pπ+π−final states is due to different PID requirements applied to pions and kaons in the first stage of the offline selection. The optimised requirements on the BDT output and PID variables for B0 s→p¯pπ+π−(B0 s→p¯pK+K−) decays retain ∼45% (40%) of the signal and reject more than 99% (99%) of the combinatorial background, inside the mass-fit ranges defined in section 6. Dedicated BDT classifiers are trained to select candidate B0 s→4hK+K−decays using the following set of input variables: the pTand the IP with respect to the SV of all reconstructed tracks; the vertex-fit χ2of the ηcand φcandidates; the vertex-fit χ2, the pT, the flight-distance significance with respect to the PV of the B0 scandidate, and the angle between the momentum and the vector joining the primary to the secondary vertex of the B0 scandidate. The optimised requirements on the BDT output and PID variables, for each of the 4hmodes, retain about 50% of the signal and reject more than 99% of the combinatorial background inside the mass-fit ranges defined in section 6. – 5 – JHEP07(2017)021 B0 s→2K2πφ B0 s→4πφ B0 s→4Kφ B0 s→p¯pφ B0 s→p¯pπ+π− ε(J/ψ ) ε(ηc)1.047 ±0.011 1.068 ±0.016 0.962 ±0.028 1.038 ±0.009 1.004 ±0.015 Table 1. Ratio of efficiencies between the normalisation and signal channels for each final state. From simulation, after all requirements for B0 s→4hK+K−decays, a significant contamination is expected from B0 s→D+ s3hdecays, where the D+ sdecays to φπ+and 3h is any combination of three charged kaons and pions. This background contribution has distributions similar to the signal in the 4hK+K−and K+K−invariant-mass spectra, while its distribution in the 4hinvariant-mass spectrum is not expected to exhibit any peaking structure. In order to reduce this background contamination, the absolute difference between the known value of the D+ smass [13] and the reconstructed invariant mass of the system formed by the combination of the φcandidate and any signal candidate track consistent with a pion hypothesis is required to be >17 MeV/c2. This requirement is optimised using the significance of B0 s→J/ψK+K−candidates with respect to background contributions. This significance is stable for cut values in the range [9,25] MeV/c2, with a maximum at 17 MeV/c2, which removes about 90% of B0 s→D+ s3hdecays, with no significant signal loss. 5 Efficiency correction The efficiency correction factors appearing in eq. (3.3) are obtained from fully simulated events. Since the signal and normalisation channels are selected based on the same requirements and have the same final-state particles with very similar kinematic distributions, the ratio between the efficiency correction factors for B0 s→ηcXand B0 s→J/ψX decays are expected to be close to unity. The efficiency correction factors include the geometrical acceptance of the LHCb detector, the reconstruction efficiency, the efficiency of the offline selection criteria, including the trigger and PID requirements. The efficiencies of the PID requirements are obtained as a function of particle momentum and number of charged tracks in the event using dedicated data-driven calibration samples of pions, kaons, and protons [34]. The overall efficiency is taken as the product of the geometrical acceptance of the LHCb detector, the reconstruction efficiency and the efficiency of the offline selection criteria. In addition, corrections are applied to account for different lifetime values used in simulation with respect to the known values for the decay channels considered. The effective lifetime for B0 sdecays to ηcφ(ηcπ+π−) final state, being purely CP-even (CP-odd), is obtained from the known value of the decay width of the light (heavy) B0 s state [35]. The effective lifetime of B0 s→J/ψφ (B0 s→J/ψπ+π−) decays is taken from ref. [35]. The lifetime correction is obtained after reweighting the signal and normalisation simulation samples. The final efficiency correction factors, given in table 1, are found to be compatible to unity as expected. – 6 – JHEP07(2017)021 6 Fit models In this section the fit models used for the measurement of the branching fractions are described, first the model used for B0 s→ηcπ+π−decays in section 6.1, then the model used for B0 s→ηcφdecays in section 6.2. 6.1 Model for B0 s→ηcπ+π−decays Candidates are fitted in two stages. First, an extended UML fit to the p¯pπ+π−invariantmass spectrum is performed in the range 5150–5540 MeV/c2, to discriminate B0 s→p¯pπ+π− events from combinatorial background, B0→p¯pπ+π−decays, and B0→p¯pKπ decays, where the kaon is misidentified as a pion. The p¯pπ+π−mass distribution of B0 s→p¯pπ+π− and B0→p¯pπ+π−candidates are described by Hypatia functions [36]. Both Hypatia functions share common core resolution and tail parameters. The latter are fixed to values obtained from simulation. The distribution of the misidentified B0→p¯pKπ background is described by a Crystal Ball function [37], with mode, power-law tail, and core resolution parameters fixed to values obtained from simulation. The combinatorial background is modelled using an exponential function. The mode and the common core resolution parameters of the Hypatia functions and the slope of the exponential functions, as well as all the yields, are allowed to vary in the fit to data. Using the information from the fit to the p¯pπ+π−spectrum, signal weights are then computed and the background components are subtracted using the sPlot technique [28]. Correlations between the p¯pand p¯pπ+π− invariant-mass spectra, for both signal and backgrounds, are found to be negligible. Second, a UML fit to the weighted p¯pinvariant-mass distribution is performed in the mass range 2900–3200 MeV/c2. In this region, three event categories are expected to populate the p¯pspectrum: the ηcand J/ψ resonances, as well as a possible contribution from nonresonant B0 s→(p¯p)NRπ+π−decays. The p¯pmass distribution of ηccandidates is described by the convolution of the square of the modulus of a complex relativistic Breit- Wigner function (RBW) with constant width and a function describing resolution effects. The expression of the RBW function is taken as Rres(m;mres,Γres)∝1 m2 res −m2−imresΓres ,(6.1) where mres and Γres are the pole mass and the natural width, respectively, of the resonance. From simulation, in the mass range considered, the p¯pinvariant-mass resolution is found to be a few MeV/c2, while Γηc= 31.8±0.8 MeV/c2[13]. Thus, the p¯pdistribution of ηccandidates is expected to be dominated by the RBW, with only small effects on the total ηclineshape from the resolution. On the other hand, due to the small natural width of the J/ψ resonance [13], the corresponding lineshape is assumed to be described to a very good approximation by the resolution function only. For the ηcand J/ψ lineshapes, Hypatia functions are used to parametrise the resolution, with tail parameters that are fixed to values obtained from simulation. A single core resolution parameter, σc¯c res, shared between these two functions, is free to vary in the fit to data. The ηcpole mass and the mode of the Hypatia function describing the J/ψ lineshape, which can be approximated – 7 – JHEP07(2017)021 by the pole mass of the resonance, are also free to vary, while the ηcnatural width is constrained to its known value [13]. The possible contribution from B0 s→(p¯p)NRπ+π− decays is parametrised by a constant. The angular distributions of P- and S-waves are characterised by a linear combination of odd- and even-order Legendre polynomials, respectively. In the case of a uniform acceptance, after integration over the helicity angles, the interference between the two waves vanishes. For a non-uniform acceptance, after integration, only residual effects from the interference between ηc(→p¯p)π+π−and J/ψ (→p¯p)π+π−amplitudes can arise in the p¯p invariant mass spectra. Due to the limited size of the current data sample, these effects are assumed to be negligible. Also, given the sample size and the small expected contribution of the NR p¯pcomponent, interference between the ηc(→p¯p)π+π−and (p¯p)NRπ+π− amplitudes is neglected. In order to fully exploit the correlation between the yields of ηcand J/ψ candidates, the former is parametrised in the fit, rearranging eq. (3.3), as Nηc=NJ/ψ ×B(B0 s→ηcπ+π−) B(B0 s→J/ψπ+π−)×B(ηc→p¯p) B(J/ψ →p¯p)×ε(ηc)p¯p ε(J/ψ)p¯p ,(6.2) where B(B0 s→ηcπ+π−) and NJ/ψ are free parameters. The yield of the NR p¯pcomponent is also free to vary. 6.2 Model for B0 s→ηcφdecays The procedure and the fit model used to measure B(B0 s→ηcφ) is based on that described in section 6.1. However, several additional features are needed to describe the data, as detailed below. The K+K−invariant mass is added as a second dimension in the first step fit, which here consists of a two-dimensional (2D) fit to the p¯pK+K−or 4hK+K−and K+K−invariant mass spectra. This allows the contributions from φ→K+K−decays and nonresonant K+K−pairs to be separated. Thus, the first step of the fitting procedure consists of four independent two-dimensional UML fits to the p¯pK+K−versus K+K−and 4hK+K−versus K+K−invariant-mass spectra in the ranges 5200–5500 MeV/c2and 990–1050 MeV/c2, respectively.3 Similar 2D fit models are used for each 4hmode. The 4hK+K−distributions of B0 s→ 4hφ signal and B0→4hφ background contributions, as well as those of B0 s→4hK+K− and B0→4hK+K−backgrounds, are described by Hypatia functions. The 4hK+K−distribution of the combinatorial background is parametrised using two exponential functions, one for when the K+K−pair arises from a random combination of two prompt kaons, and another for when the K+K−pair originates from the decay of a prompt φmeson. The K+K−distribution of each contribution including a φin the final state is described by the square of the modulus of a RBW with mass-dependent width convolved with a Gaussian function accounting for resolution effects. The K+K−distributions of the contributions 3In order to better constrain the combinatorial background shape, the upper limit of the p¯pK+K− invariant-mass range is extended to 5550 MeV/c2. – 8 – JHEP07(2017)021 2K2π4π4K p¯p FFJ/ψ 0.39 ±0.03 0.28 ±0.03 0.29 ±0.05 0.76 ±0.03 FFηc0.49 ±0.05 0.63 ±0.07 0.31 ±0.03 0.24 ±0.03 FFbkg 0.12 ±0.10 0.02 ±0.03 0.32 ±0.19 n/a FFNR 0.24 ±0.11 0.12 ±0.09 0.15 ±0.16 n/a PkFFk1.24 ±0.07 1.05 ±0.11 1.08 ±0.08 1.00 Table 5. Fit fractions obtained from the parameters of the simultaneous fit to the p¯pand 4hinvariant-mass spectra. Uncertainties are statistical only. Due to interference between B0 s→ηc(→4h)φand B0 s→(4h)NRφamplitudes, for the 4hfinal states the sum of fit fractions, PkFFk, may be different from unity. The abbreviation “n/a” stands for “not applicable”. 8 Systematic uncertainties As the expressions for B(B0 s→ηcπ+π−) and B(B0 s→ηcφ) are based on the ratios of observed quantities, only sources of systematic uncertainties inducing different biases to the number of observed ηcand J/ψ candidates are considered. The dominant source of systematic uncertainties is due to the knowledge of the external branching fractions. These are estimated by adding Gaussian constraints on the external branching fractions in the fits, with widths corresponding to their known uncertainties [13]. A summary of the systematic uncertainties can be found in table 6. To assign systematic uncertainties due to fixing of PDF parameters, the fits are repeated by varying all of them simultaneously. The resolution parameters, estimated from simulation, are varied according to normal distributions, taking into account the correlations between the parameters and with variances related to the size of the simulated samples. The external parameters are varied within a normal distribution of mean and width fixed to their known values and uncertainties [13]. This procedure is repeated 1000 times, and for each iteration a new value of the branching fraction is obtained. The systematic uncertainties on the branching fraction are taken from the variance of the corresponding distributions. The systematic uncertainty due to the fixing of the values of the efficiencies is estimated by adding Gaussian constraints to the likelihood functions, with widths that are taken from the uncertainties quoted in table 1. The presence of intrinsic biases in the fit models is studied using parametric simulation. For this study, 1000 pseudoexperiments are generated and fitted using the nominal PDFs, where the generated parameter values correspond to those obtained in the fits to data. The biases on the branching fractions are then calculated as the difference between the generated values and the mean of the distribution of the fitted branching fraction values. To assign a systematic uncertainty from the model used to describe the detector resolution, the fits are repeated for each step replacing the Hypatia functions by bifurcated Crystal Ball functions, the parameters of which are obtained from simulation. The difference from the nominal branching fraction result is assigned as a systematic uncertainty. – 15 – JHEP07(2017)021 Source Value [%] B(B0 s→ηcπ+π−)B(B0 s→ηcφ) Fixed PDF parameters 5.7 1.4 Efficiencies 3.4 0.8 Fit bias 1.7 1.4 Resolution model 0.6 4.4 φ(1020) barrier radius n/a 1.6 Acceptance (4h) n/a 1.6 Nonresonant p¯pn/a 1.0 Sum 6.8 5.4 External branching fractions 16.4 12.6 Table 6. Summary of systematic uncertainties. The “Sum” of systematic uncertainties is obtained from the quadratic sum of the individual sources, except the external branching fractions, which are quoted separately. All values are in % of the measured branching fractions. The abbreviation “n/a” stands for “not applicable”. The Blatt-Weisskopf parameter rof the φis arbitrarily set to 3 (GeV/c)−1. To assign a systematic uncertainty due to the fixed value of this rparameter, the fits are repeated for different values taken in the range 1.5–5.0 (GeV/c)−1. The maximum differences from the nominal branching fraction result are assigned as systematic uncertainties. To assign a systematic uncertainty due to the assumption of a uniform acceptance, the simultaneous fit is repeated after correcting the 4hinvariant-mass distributions for acceptance effects. A histogram describing the acceptance effects in each of the 4hinvariant- mass spectra is constructed from the ratio of the normalised 4hinvariant-mass distributions taken from simulated samples of B0 s→(4h)φphase space decays, obtained either directly from EvtGen, or after processing through the full simulation chain. The simultaneous fit is repeated after applying weights for each event from the central value of its bin in the 4hinvariant-mass distribution. The difference from the nominal branching fraction result is assigned as a systematic uncertainty. No significant dependence on the binning choice was observed. The systematic uncertainty due to neglecting the presence of a nonresonant p¯pcontribution in the p¯pspectrum for B0 s→p¯pφ candidates is estimated by repeating the simultaneous fit with an additional component described by an exponential function, where the slope and the yield are allowed to vary. The difference from the nominal branching fraction result is assigned as a systematic uncertainty. 9 Conclusions This paper reports the observation of B0 s→ηcφdecays and the first evidence for B0 s→ ηcπ+π−decays. The branching fractions are measured to be B(B0 s→ηcφ) = (5.01 ±0.53 ±0.27 ±0.63) ×10−4, B(B0 s→ηcπ+π−) = (1.76 ±0.59 ±0.12 ±0.29) ×10−4, – 16 – JHEP07(2017)021 where in each case the two first uncertainties are statistical and systematic, respectively, and the third uncertainties are due to the limited knowledge of the external branching fractions. The significance of the B0 s→ηcπ+π−decay mode, including systematic uncertainties, is 4.6σ. The results for B(B0 s→ηcπ+π−) and B(B0 s→ηcφ) are in agreement with expectations based on eqs. (1.1), (1.2) and (1.3). The data sample recorded by the LHCb experiment in Run 1 of the LHC is not sufficiently large to allow a measurement of the CP-violating phase φsfrom time-dependent analysis of B0 s→ηcφor B0 s→ηcπ+π−decays. However, in the future with significant improvement of the hadronic trigger efficiencies [40], these decay modes may become of interest to add sensitivity to the measurement of φs. Acknowledgments 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); FOM and NWO (The Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania); MinES and FASO (Russia); MinECo (Spain); SNSF and SER (Switzerland); NASU (Ukraine); STFC (United Kingdom); NSF (U.S.A.). We acknowledge the computing resources that are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN (Italy), SURF (The Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and Yandex LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), PL-GRID (Poland) and OSC (U.S.A.). 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 AvH Foundation (Germany), EPLANET, Marie Sk lodowska-Curie Actions and ERC (European Union), Conseil G´en´eral de Haute-Savoie, Labex ENIGMASS and OCEVU, R´egion Auvergne (France), RFBR and Yandex LLC (Russia), GVA, XuntaGal and GENCAT (Spain), Herchel Smith Fund, The Royal Society, Royal Commission for the Exhibition of 1851 and the Leverhulme Trust (United Kingdom). A Fit projections The p¯pπ+π−invariant mass distribution and the fit projection are shown in figure 5. The four p¯p(4h)K+K−and K+K−invariant-mass distributions and the corresponding twodimensional fit projections are shown in figures 6to 9. B Correlation matrix The statistical correlation matrix for the simultaneous fit to the p¯pand 4hinvariant-mass distributions for B0 s→p¯pφ and B0 s→4hφ candidates is given in table 7. – 17 – JHEP07(2017)021 ] 2 ) [MeV/c − π + π ppm( 5200 5300 5400 5500 ) 2 Candidates / (10.0 MeV/c 0 50 100 150 200 250 Data Total PDF − π + π pp → 0 s B − π + π pp → 0 B − π + π ppCombinatorial − π + K pp → 0 B LHCb Figure 5. Distribution of the p¯pπ+π−invariant mass. Points with error bars show the data. The solid curve is the projection of the total fit result. The short-dashed blue, the dashed-double-dotted green, the dashed-single-dotted yellow and medium-dashed red curves show the B0 s→p¯pπ+π−, B0→p¯pπ+π−,B0→p¯pK+π−and combinatorial background contributions, respectively. – 18 – JHEP07(2017)021 ] 2 ) [MeV/c − K + Kppm( 5200 5300 5400 5500 ) 2 Candidates / (6.0 MeV/c 0 10 20 30 40 50 60 70 80 90 Data Total PDF φ pp → s 0 B − K + Kp p → s 0 B φ pp → 0 B − K + Kp p → 0 B φ ppCombinatorial − K + Kp pCombinatorial − π + Kp p → 0 B LHCb ] 2 ) [MeV/c − K + Km( 1000 1020 1040 ) 2 Candidates / (1.2 MeV/c 0 20 40 60 80 100 120 LHCb Figure 6. Distribution of the p¯pK+K−(top) and K+K−(bottom) invariant masses. Points with error bars show the data. The solid black curve is the projection of the total fit result. The short-dashed and dotted blue curves show the B0 s→p¯pφ and B0 s→p¯pK+K−contributions, respectively. The long-dashed and medium-dashed red curves show the contributions of the combinatorial background with prompt φand NR K+K−, respectively. The dashed-double-dotted green, dashedtriple-dotted green and dashed-single-dotted-yellow curves show the B0→p¯pφ,B0→p¯pK+K− and B0→p¯pK+π−contributions, respectively. – 19 – JHEP07(2017)021 ] 2 ) [MeV/c − K + K π 2Km(2 5200 5300 5400 5500 ) 2 Candidates / (6.0 MeV/c 0 20 40 60 80 100 120 140 Data Total PDF φπ 2K 2→ s 0 B − K + K π 2K 2→ s 0 B φπ 2K 2→ 0 B − K + K π 2K 2→ 0 B φπ 2KCombinatorial 2 − K + K π 2KCombinatorial 2 LHCb ] 2 ) [MeV/c − K + Km( 1000 1020 1040 ) 2 Candidates / (1.2 MeV/c 0 20 40 60 80 100 120 140 160 180 LHCb Figure 7. Distribution of the 2K2πK+K−(top) and K+K−(bottom) invariant masses. Points with error bars show the data. The solid black curve is the projection of the total fit result. The short-dashed and dotted blue curves show the B0 s→2K2πφ and B0 s→2K2πK+K−contributions, respectively. The long-dashed and medium-dashed red curves show the contributions of the combinatorial background with prompt φand NR K+K−, respectively. The dashed-double-dotted and dashed-triple-dotted green curves show the B0→2K2πφ and B0→2K2πK+K−contributions, respectively. – 20 – JHEP07(2017)021 ] 2 ) [MeV/c − K + K π m(4 5200 5300 5400 5500 ) 2 Candidates / (6.0 MeV/c 0 20 40 60 80 100 Data Total PDF φπ 4→ s 0 B − K + K π 4→ s 0 B φπ 4→ 0 B − K + K π 4→ 0 B φπ Combinatorial 4 − K + K π Combinatorial 4 LHCb ] 2 ) [MeV/c − K + Km( 1000 1020 1040 ) 2 Candidates / (1.2 MeV/c 0 20 40 60 80 100 120 140 160 180 LHCb Figure 8. Distribution of the 4πK+K−(top) and K+K−(bottom) invariant masses. Points with error bars show the data. The solid black curve is the projection of the total fit result. The short-dashed and dotted blue curves show the B0 s→4πφ and B0 s→4πK+K−contributions, respectively. The long-dashed and medium-dashed red curves show the contributions of the combinatorial background with prompt φand NR K+K−, respectively. The dashed-double-dotted and dashed-triple-dotted green curves show the B0→4πφ and B0→4πK+K−contributions, respectively. – 21 – JHEP07(2017)021 ] 2 ) [MeV/c − K + KKm(4 5200 5300 5400 5500 ) 2 Candidates / (6.0 MeV/c 0 5 10 15 20 25 30 35 Data Total PDF φ K 4→ s 0 B − K + KK 4→ s 0 B φ K 4→ 0 B − K + KK 4→ 0 B φ KCombinatorial 4 − K + KKCombinatorial 4 LHCb ] 2 ) [MeV/c − K + Km( 1000 1020 1040 ) 2 Candidates / (1.2 MeV/c 0 10 20 30 40 50 LHCb Figure 9. Distribution of the 4KK+K−(top) and K+K−(bottom) invariant masses. Points with error bars show the data. The solid black curve is the projection of the total fit result. The short-dashed and dotted blue curves show the B0 s→4Kφ and B0 s→4KK+K−contributions, respectively. The long-dashed and medium-dashed red curves show the contributions of the combinatorial background with prompt φand NR K+K−, respectively. The dashed-double-dotted and dashed-triple-dotted green curves show the B0→4Kφ and B0→4KK+K−contributions, respectively. – 22 – JHEP07(2017)021 κ2K2π NR κ4K NR κ4π NR κ2K2π bkg κ4K bkg κ4π bkg ξ2K2π NR ξ4K NR ξ4π NR B(B0 s→ηcφ) +0.22 −0.00 −0.05 −0.19 +0.00 −0.00 +0.55 +0.04 +0.13 κ2K2π NR −0.00 −0.01 +0.65 −0.00 −0.00 +0.07 +0.01 +0.02 κ4K NR −0.01 +0.00 −0.06 −0.00 +0.00 +0.37 +0.02 κ4π NR −0.00 −0.01 +0.25 −0.10 +0.09 −0.72 κ2K2π bkg +0.00 −0.00 −0.47 −0.01 −0.01 κ4K bkg −0.00 +0.01 −0.21 +0.02 κ4π bkg −0.01 +0.01 −0.07 ξ2K2π NR −0.02 +0.17 ξ4K NR −0.12 ξ2K2π bkg ξ4K bkg ξ4π bkg mηcmJ/ψ δϕ2K2πδϕ4Kδϕ4π B(B0 s→ηcφ)−0.40 −0.04 −0.11 −0.02 +0.00 −0.23 +0.28 −0.44 κ2K2π NR −0.00 −0.01 −0.02 +0.01 +0.00 −0.02 +0.06 −0.09 κ4K NR −0.00 −0.26 −0.01 −0.03 −0.00 +0.01 −0.10 +0.00 κ4π NR +0.05 −0.05 +0.52 +0.45 −0.04 −0.15 +0.02 +0.09 κ2K2π bkg +0.55 +0.01 +0.01 −0.03 +0.00 +0.49 −0.05 +0.08 κ4K bkg −0.00 +0.04 −0.01 −0.03 −0.00 +0.01 −0.17 −0.00 κ4π bkg +0.00 −0.00 −0.33 +0.04 −0.00 −0.01 +0.00 +0.03 ξ2K2π NR −0.82 −0.00 −0.12 −0.18 +0.01 −0.23 +0.14 −0.24 ξ4K NR −0.00 −0.52 +0.07 +0.21 −0.01 −0.08 +0.07 −0.01 ξ4π NR −0.09 +0.06 −0.49 −0.61 +0.01 +0.18 −0.01 +0.05 ξ2K2π bkg +0.01 +0.07 +0.07 −0.02 +0.41 −0.11 +0.18 ξ4K bkg −0.03 −0.11 −0.02 +0.05 −0.26 +0.01 ξ4π bkg +0.34 −0.02 −0.09 −0.01 +0.28 mηc−0.01 −0.35 +0.07 +0.02 mJ/ψ +0.01 +0.01 −0.01 δϕ2K2π−0.09 +0.10 δϕ4K−0.12 Table 7. Statistical correlation matrix for the parameters from the simultaneous fit to the p¯pand 4hinvariant-mass spectra for B0 s→p¯pφ and B0 s→4hφ candidates. – 23 – JHEP07(2017)021 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. References [1] CKMfitter Group collaboration, J. Charles et al., CP violation and the CKM matrix: Assessing the impact of the asymmetric Bfactories,Eur. Phys. J. C 41 (2005) 1 [hep-ph/0406184] [INSPIRE]. [2] W. Altmannshofer, A.J. Buras, S. Gori, P. Paradisi and D.M. Straub, Anatomy and phenomenology of FCNC and CPV effects in SUSY theories,Nucl. Phys. B 830 (2010) 17 [arXiv:0909.1333] [INSPIRE]. [3] W. Altmannshofer, A.J. Buras and D. Guadagnoli, The MFV limit of the MSSM for low tan β: Meson mixings revisited,JHEP 11 (2007) 065 [hep-ph/0703200] [INSPIRE]. [4] A.J. Buras, Flavour theory: 2009,PoS(EPS-HEP 2009)024 [arXiv:0910.1032] [INSPIRE]. [5] C.-W. Chiang, A. Datta, M. Duraisamy, D. London, M. Nagashima and A. Szynkman, New physics in B0 s→J/ψφ: A general analysis,JHEP 04 (2010) 031 [arXiv:0910.2929] [INSPIRE]. [6] LHCb collaboration, Precision measurement of CP violation in B0 s→J/ψK+K−decays, Phys. Rev. Lett. 114 (2015) 041801 [arXiv:1411.3104] [INSPIRE]. [7] D0 collaboration, V.M. Abazov et al., Measurement of the CP-violating phase φJ/ψφ susing the flavor-tagged decay B0 s→J/ψφ in 8 fb−1of p¯pcollisions,Phys. Rev. D 85 (2012) 032006 [arXiv:1109.3166] [INSPIRE]. [8] CDF collaboration, T. Aaltonen et al., Measurement of the CP-violating phase βJ/ψφ sin B0 s→J/ψφ decays with the CDF II detector,Phys. Rev. D 85 (2012) 072002 [arXiv:1112.1726] [INSPIRE]. [9] CMS collaboration, Measurement of the CP-violating weak phase φsand the decay width difference ∆Γsusing the B0 s→J/ψφ(1020) decay channel in pp collisions at √s= 8 TeV, Phys. Lett. B 757 (2016) 97 [arXiv:1507.07527] [INSPIRE]. [10] ATLAS collaboration, Measurement of the CP-violating phase φsand the B0 smeson decay width difference with B0 s→J/ψφ decays in ATLAS,JHEP 08 (2016) 147 [arXiv:1601.03297] [INSPIRE]. [11] LHCb collaboration, Measurement of the CP-violating phase φsin ¯ B0 s→D+ sD− sdecays, Phys. Rev. Lett. 113 (2014) 211801 [arXiv:1409.4619] [INSPIRE]. [12] LHCb collaboration, First study of the CP-violating phase and decay-width difference in B0 s→ψ(2S)φdecays,Phys. Lett. B 762 (2016) 253 [arXiv:1608.04855] [INSPIRE]. [13] Particle Data Group collaboration, C. Patrignani et al., Review of particle physics,Chin. Phys. C 40 (2016) 100001 [INSPIRE]. [14] LHCb collaboration, The LHCb detector at the LHC,2008 JINST 3S08005 [INSPIRE]. [15] LHCb collaboration, LHCb detector performance,Int. J. Mod. Phys. A 30 (2015) 1530022 [arXiv:1412.6352] [INSPIRE]. [16] R. Aaij et al., The LHCb trigger and its performance in 2011,2013 JINST 8P04022 [arXiv:1211.3055] [INSPIRE]. – 24 – JHEP07(2017)021 37 Institute for High Energy Physics (IHEP), Protvino, Russia 38 ICCUB, Universitat de Barcelona, Barcelona, Spain 39 Universidad de Santiago de Compostela, Santiago de Compostela, Spain 40 European Organization for Nuclear Research (CERN), Geneva, Switzerland 41 Institute of Physics, Ecole Polytechnique F´ed´erale de Lausanne (EPFL), Lausanne, Switzerland 42 Physik-Institut, Universit¨at Z¨urich, Z¨urich, Switzerland 43 Nikhef National Institute for Subatomic Physics, Amsterdam, The Netherlands 44 Nikhef National Institute for Subatomic Physics and VU University Amsterdam, Amsterdam, The Netherlands 45 NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine 46 Institute for Nuclear Research of the National Academy of Sciences (KINR), Kyiv, Ukraine 47 University of Birmingham, Birmingham, United Kingdom 48 H.H. Wills Physics Laboratory, University of Bristol, Bristol, United Kingdom 49 Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom 50 Department of Physics, University of Warwick, Coventry, United Kingdom 51 STFC Rutherford Appleton Laboratory, Didcot, United Kingdom 52 School of Physics and Astronomy, University of Edinburgh, Edinburgh, United Kingdom 53 School of Physics and Astronomy, University of Glasgow, Glasgow, United Kingdom 54 Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom 55 Imperial College London, London, United Kingdom 56 School of Physics and Astronomy, University of Manchester, Manchester, United Kingdom 57 Department of Physics, University of Oxford, Oxford, United Kingdom 58 Massachusetts Institute of Technology, Cambridge, MA, United States 59 University of Cincinnati, Cincinnati, OH, United States 60 University of Maryland, College Park, MD, United States 61 Syracuse University, Syracuse, NY, United States 62 Pontif´ıcia Universidade Cat´olica do Rio de Janeiro (PUC-Rio), Rio de Janeiro, Brazil, associated to 2 63 University of Chinese Academy of Sciences, Beijing, China, associated to 3 64 School of Physics and Technology, Wuhan University, Wuhan, China, associated to 3 65 Institute of Particle Physics, Central China Normal University, Wuhan, Hubei, China, associated to 3 66 Departamento de Fisica , Universidad Nacional de Colombia, Bogota, Colombia, associated to 8 67 Institut f¨ur Physik, Universit¨at Rostock, Rostock, Germany, associated to 12 68 National Research Centre Kurchatov Institute, Moscow, Russia, associated to 32 69 Instituto de Fisica Corpuscular, Centro Mixto Universidad de Valencia - CSIC, Valencia, Spain, associated to 38 70 Van Swinderen Institute, University of Groningen, Groningen, The Netherlands, associated to 43 aUniversidade Federal do Triˆangulo Mineiro (UFTM), Uberaba-MG, Brazil bLaboratoire Leprince-Ringuet, Palaiseau, France cP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS), Moscow, Russia dUniversit`a di Bari, Bari, Italy eUniversit`a di Bologna, Bologna, Italy fUniversit`a di Cagliari, Cagliari, Italy gUniversit`a di Ferrara, Ferrara, Italy hUniversit`a di Genova, Genova, Italy iUniversit`a di Milano Bicocca, Milano, Italy jUniversit`a di Roma Tor Vergata, Roma, Italy kUniversit`a di Roma La Sapienza, Roma, Italy lAGH - University of Science and Technology, Faculty of Computer Science, Electronics and Telecommunications, Krak´ow, Poland – 31 – JHEP07(2017)021 mLIFAELS, La Salle, Universitat Ramon Llull, Barcelona, Spain nHanoi University of Science, Hanoi, Viet Nam oUniversit`a di Padova, Padova, Italy pUniversit`a di Pisa, Pisa, Italy qUniversit`a degli Studi di Milano, Milano, Italy rUniversit`a di Urbino, Urbino, Italy sUniversit`a della Basilicata, Potenza, Italy tScuola Normale Superiore, Pisa, Italy uUniversit`a di Modena e Reggio Emilia, Modena, Italy vIligan Institute of Technology (IIT), Iligan, Philippines wNovosibirsk State University, Novosibirsk, Russia †Deceased – 32 –